Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
JPulsar.cc
Go to the documentation of this file.
1#include <string>
2#include <iostream>
3#include <iomanip>
4#include <map>
5
6#include "TROOT.h"
7#include "TFile.h"
8#include "TH1D.h"
9#include "TF1.h"
10#include "TFitResult.h"
11
12#include "JROOT/JMinimizer.hh"
13
15
20
21#include "JTrigger/JHit.hh"
22#include "JTrigger/JFrame.hh"
24#include "JTrigger/JHitR0.hh"
25#include "JTrigger/JBuildL0.hh"
26
27#include "JCalibrate/JTDC_t.hh"
28
30#include "JSupport/JSupport.hh"
31#include "JSupport/JMeta.hh"
32
33#include "JTools/JRange.hh"
36
37#include "Jeep/JProperties.hh"
38#include "Jeep/JParser.hh"
39#include "Jeep/JMessage.hh"
40
41
42namespace {
43
44 /**
45 * Get time within given period.
46 *
47 * The returned time is constraint to the interval <tt>[-T/2,+T/2]</tt>.
48 *
49 * \param t1 time
50 * \param T period
51 * \return time
52 */
53 inline double get_time(double t1, const double T)
54 {
55 for (int buffer[] = { 1000, 100, 10, 1, 0 }, *i = buffer; *i != 0; ++i) {
56
57 const double xmin = -0.5 * (*i) * T;
58 const double xmax = +0.5 * (*i) * T;
59
60 while (t1 > xmax) { t1 -= (*i) * T; }
61 while (t1 < xmin) { t1 += (*i) * T; }
62 }
63
64 return t1;
65 }
66}
67
68
69/**
70 * \file
71 *
72 * Application for dark room time calibration using external laser.
73 *
74 * The time calibration is made for the specified reference PMTs (option <tt>-! "<module identifier> <PMT channel>"</tt>.\n
75 * In this, the wild card value <tt>-1</tt> can be used.
76 * A ROOT histogram corresponding to each reference PMT is created, filled, fitted and stored to the specified output file.\n
77 * The detector can be updated using option <tt>-A</tt>\n
78 * The laser should be synchronised with the same clock system as the detector and
79 * produce single pulses at a fixed frequency (option <tt>-l <laserFrequency_Hz></tt>.
80 *
81 * In addition, the transit-time distribution of the reference PMTs can be monitored (see JLegolas.cc).\n
82 * For this, the expectation values should be less than one and the time-over-threshold distribution nominal.
83 *
84 * \author mbouwhuis
85 */
86int main(int argc, char **argv)
87{
88 using namespace std;
89 using namespace JPP;
90 using namespace KM3NETDAQ;
91
93
95 JLimit_t& numberOfEvents = inputFile.getLimit();
96 string outputFile;
97 string detectorFile;
98 JTDC_t TDC;
99 double laserFrequency_Hz;
100 bool overwriteDetector;
101 JROOTClassSelector selector;
102 string option;
103 double Wmin = 1e3; // Minimal signal intensity
104 JRange_t T_ns;
105 JRange_t X;
106 int debug;
107
108 try {
109
110 JProperties properties;
111
112 properties.insert(gmake_property(Wmin));
113
114 JParser<> zap("Application for dark room time calibration.");
115
116 zap['f'] = make_field(inputFile, "input file (time slice data from laser calibration).");
117 zap['o'] = make_field(outputFile, "output file.") = "pulsar.root";
118 zap['a'] = make_field(detectorFile, "detector file.");
119 zap['!'] = make_field(TDC,
120 "Set reference PMTs, e.g."
121 "\n-! \"808969848 0 808982077 23\" sets PMT 0 of module 808969848 and PMT 23 of module 808982077 as references.");
122 zap['n'] = make_field(numberOfEvents) = JLimit::max();
123 zap['l'] = make_field(laserFrequency_Hz, "laser frequency [Hz]") = 10000;
124 zap['A'] = make_field(overwriteDetector, "overwrite detector file provided through '-a' with correct time offsets.");
125 zap['O'] = make_field(option, "ROOT fit option, see TH1::Fit.") = "LS";
126 zap['@'] = make_field(properties) = JPARSER::initialised();
127 zap['C'] = make_field(selector, "timeslice selector, e.g. JDAQTimesliceL1.") = getROOTClassSelection<JDAQTimesliceTypes_t>();
128 zap['T'] = make_field(T_ns, "time window for time-over-threshold monitor") = JRange_t(-10.0, +10.0);
129 zap['x'] = make_field(X, "time window for histogram") = JRange_t();
130 zap['d'] = make_field(debug, "debug.") = 1;
131
132 zap(argc, argv);
133 }
134 catch(const exception& error) {
135 FATAL(error.what() << endl);
136 }
137
138
139 TDC.is_valid(true);
140
141 if (laserFrequency_Hz <= 0.0) {
142 FATAL("Invalid laser frequency " << laserFrequency_Hz << endl);
143 }
144
145 const double laserPeriod_ns = 1.0e9 / laserFrequency_Hz;
146
147 if (option.find('R') == string::npos) { option += 'R'; }
148 if (option.find('S') == string::npos) { option += 'S'; }
149 if (option.find('Q') == string::npos && debug < JEEP::debug_t) { option += 'Q'; }
150
151
153
154 try {
155 load(detectorFile, detector);
156 }
157 catch(const JException& error) {
158 FATAL(error);
159 }
160
161 const JModuleRouter moduleRouter(detector);
162
163
164 // ROOT I/O
165
166 TFile out(outputFile.c_str(), "recreate");
167
168 putObject(out, JMeta(argc, argv));
169
171
172 map_type zmap;
173
174 const double xmin = (X == JRange_t() ? -0.5 * laserPeriod_ns : X.getLowerLimit());
175 const double xmax = (X == JRange_t() ? +0.5 * laserPeriod_ns : X.getUpperLimit());
176 const int nx = 2 * (int) (xmax - xmin);
177
178 TH1D h0("h0", NULL, nx, xmin, xmax);
179 TH1D h1("h1", NULL, 256, -0.5, +255.5);
180
182
183 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
184
185 const JTDC_t::range_type range = TDC.equal_range(module->getID());
186
187 for (JTDC_t::const_iterator i = range.first; i != range.second; ++i) {
188
189 NOTICE("Reference PMT at"
190 << " string " << setw(3) << module->getString()
191 << " floor " << setw(2) << module->getFloor()
192 << " module " << setw(8) << module->getID()
193 << " channel " << setw(2) << i->second << endl);
194
195 const JPMTIdentifier id(module->getID(),i->second);
196
197 ostringstream os;
198
199 os << getLabel(module->getLocation()) << ' ' << getLabel(id);
200
201 zmap.insert(make_pair(id, new TH1D(os.str().c_str(), NULL, nx, xmin, xmax)));
202 }
203 }
204
205
206 typedef vector<JHitR0> JDataL0_t;
207
208 const JBuildL0<JHitR0> buildL0;
209
210
212
213 counter_type counter = 0;
214
215 for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
216
217 STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl);
218
219 const JDAQTimeslice* timeslice = in.next();
220
221 for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
222
223 const JTDC_t::range_type range = TDC.equal_range(frame->getModuleID());
224
225 if (range.first != range.second) {
226
227 const double t0 = get_time(getTimeOfFrame(frame->getFrameIndex()), laserPeriod_ns);
228
229 JDataL0_t dataL0;
230
231 buildL0(*frame, moduleRouter.getModule(frame->getModuleID()), back_inserter(dataL0));
232
233 for (JDataL0_t::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
234
235 const JPMTIdentifier id(frame->getModuleID(), hit->getPMT());
236
237 map_type::const_iterator p = zmap.find(id);
238
239 if (p != zmap.end()) {
240
241 const double t1 = get_time(t0 + hit->getT(), laserPeriod_ns);
242
243 p->second->Fill(t1);
244
245 h0.Fill(t1);
246
247 if (T_ns(t1)) {
248
249 h1.Fill(hit->getToT());
250
251 counts[id] += 1;
252 }
253 }
254 }
255 }
256 }
257 }
258 STATUS(endl);
259
260
261 map<JPMTIdentifier, double> t0; // fitted time offsets
262
263
264 TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2.0*TMath::Pi())*[2]) + [3]");
265
266 for (map_type::iterator it = zmap.begin(); it != zmap.end(); ++it) {
267
268 JPMTIdentifier pmt = it->first;
269 TH1D* h1 = it->second;
270
271 if (h1->GetEntries() == 0) {
272 WARNING("Histogram " << h1->GetName() << " empty" << endl);
273 continue;
274 }
275
276 STATUS("--- PMT = " << pmt << "; histogram " << h1->GetName() << endl);
277
278 // start values
279
280 Double_t ymax = 0.0;
281 Double_t x0 = 0.0; // [ns]
282 Double_t sigma = 2.0; // [ns]
283 Double_t ymin = 0.0;
284
285 for (int i = 1; i != h1->GetNbinsX(); ++i) {
286
287 const Double_t x = h1->GetBinCenter (i);
288 const Double_t y = h1->GetBinContent(i);
289
290 if (y > ymax) {
291 ymax = y;
292 x0 = x;
293 }
294 }
295
296 f1.SetParameter(0, ymax/sigma);
297 f1.SetParameter(1, x0);
298 f1.SetParameter(2, sigma);
299 f1.SetParameter(3, ymin);
300
301 for (Int_t i = 0; i != f1.GetNpar(); ++i) {
302 f1.SetParError(i, 0.0);
303 }
304
305
306 // fit
307
308 TFitResultPtr result = h1->Fit(&f1, option.c_str(), "same", x0 - 3 * sigma, x0 + 3 * sigma);
309
310 if (result.Get() == NULL) {
311 FATAL("Invalid TFitResultPtr " << h1->GetName() << endl);
312 }
313
314 if (debug >= notice_t || !result->IsValid() || f1.GetParameter(0) < Wmin) {
315 cout << "Histogram " << h1->GetName() << " fit " << (result->IsValid() ? "ok" : "failed") << endl;
316 cout << "\tw = " << FIXED(12,3) << f1.GetParameter(0) << endl;
317 cout << "\tx0 = " << FIXED(12,3) << x0 << endl;
318 cout << "\tt0 = " << FIXED(12,3) << f1.GetParameter(1) << endl;
319 }
320
321
322 // check for multiple peaks
323
324 int number_of_peaks = 0;
325
326 Double_t dx = 2.0 * f1.GetParameter(2);
327 Double_t W = 0.0;
328 Double_t Y = f1.GetParameter(3);
329
330 if (dx < (xmax - xmin) / nx) {
331 dx = (xmax - xmin) / nx; // minimal width
332 }
333
334 for (Int_t il = 1, ir = il; ir <= nx; ) {
335
336 for ( ; ir <= nx && h1->GetBinCenter(ir) <= h1->GetBinCenter(il) + dx; ++ir) {
337 W += h1->GetBinContent(ir) - Y;
338 }
339
340 if (W >= Wmin) {
341
342 number_of_peaks += 1;
343
344 W = 0.0; // reset
345 il = ir;
346 ir = il;
347
348 } else {
349
350 W -= h1->GetBinContent(il) - Y; // shift
351 il += 1;
352 }
353 }
354
355 if (number_of_peaks != 1) {
356 WARNING("Number of peaks " << h1->GetName() << ' ' << number_of_peaks << " != 1" << endl);
357 }
358
359 if (result->IsValid() && f1.GetParameter(0) >= Wmin) {
360 t0[pmt] = f1.GetParameter(1);
361 }
362 }
363
364 out.Write();
365 out.Close();
366
367
368 if (counter != 0) {
369
370 const double W = laserFrequency_Hz * counter * getFrameTime() * 1.0e-9;
371
372 NOTICE("Expection values [npe]" << endl);
373
374 for (map<JPMTIdentifier, int>::const_iterator i = counts.begin(); i != counts.end(); ++i) {
375 NOTICE(i->first << ' ' << FIXED(7,3) << i->second / W << endl);
376 }
377 }
378
379
380 if (overwriteDetector) {
381
382 int errors = 0;
383
384 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
385
386 const JTDC_t::range_type range = TDC.equal_range(module->getID());
387
388 if (range.first != range.second) {
389
390 const JPMTIdentifier id(module->getID(), range.first->second); // select t0 of the first reference PMT in this optical module
391
393
394 if (p0 != t0.end()) {
395
396 const double t1 = p0->second;
397
398 for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
399
400 map<JPMTIdentifier, double>::const_iterator p1 = t0.find(JPMTIdentifier(module->getID(),pmt));
401
402 if (p1 != t0.end())
403 module->getPMT(pmt).subT0(p1->second); // offset with t0 of this reference PMT
404 else
405 module->getPMT(pmt).subT0(t1); // offset with t0 of first reference PMT in same optical module
406 }
407
408 } else {
409
410 if (!module->getPMT(id.getPMTAddress()).has(PMT_DISABLE)) {
411 ++errors;
412 }
413
414 ERROR("Module/PMT "
415 << setw(10) << module->getID() << "/" << FILL(2,'0') << id.getPMTAddress() << FILL() << ' '
416 << "missing or insufficient signal." << endl);
417 }
418 }
419 }
420
421 if (errors == 0) {
422
423 NOTICE("Store calibration data on file " << detectorFile << endl);
424
425 detector.comment.add(JMeta(argc, argv));
426
427 store(detectorFile, detector);
428
429 } else {
430
431 FATAL("Number of errors " << errors << " != 0" << endl);
432 }
433 }
434}
string outputFile
Data structure for detector geometry and calibration.
TPaveText * p1
Basic data structure for L0 hit.
General purpose messaging.
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
#define STATUS(A)
Definition JMessage.hh:63
#define NOTICE(A)
Definition JMessage.hh:64
#define FATAL(A)
Definition JMessage.hh:67
int debug
debug level
Definition JSirene.cc:72
#define WARNING(A)
Definition JMessage.hh:65
ROOT I/O of application specific meta data.
Direct access to module in detector data structure.
Scanning of objects from multiple files according a format that follows from the extension of each fi...
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
Utility class to parse parameter values.
#define gmake_property(A)
macros to convert (template) parameter to JPropertiesElement object
int main(int argc, char **argv)
Definition JPulsar.cc:86
Auxiliary class to define a range between two values.
ROOT TTree parameter settings of various packages.
Basic data structure for time and time over threshold information of hit.
Detector data structure.
Definition JDetector.hh:96
Router for direct addressing of module data in detector data structure.
const JModule & getModule(const JObjectID &id) const
Get module parameters.
Utility class to parse parameter values.
General exception.
Definition JException.hh:24
Auxiliary class for multiplexing object iterators.
virtual bool hasNext() override
Check availability of next element.
virtual const pointer_type & next() override
Get next element.
Utility class to parse command line options.
Definition JParser.hh:1698
General purpose class for object reading from a list of file names.
Range of values.
Definition JRange.hh:42
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
Template L0 hit builder.
Definition JBuildL0.hh:38
const double xmax
const double xmin
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
Definition JLocation.hh:247
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
void store(const std::string &file_name, const JDetector &detector)
Store detector to output file.
@ debug_t
debug
Definition JMessage.hh:29
@ notice_t
notice
Definition JMessage.hh:32
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
bool putObject(TDirectory &dir, const TObject &object)
Write object to ROOT directory.
std::set< JROOTClassSelector > getROOTClassSelection(const bool option=false)
Get ROOT class selection.
Long64_t counter_type
Type definition for counter.
return result
Definition JPolint.hh:862
KM3NeT DAQ data structures and auxiliaries.
Definition DataQueue.cc:39
double getFrameTime()
Get frame time duration.
Definition JDAQClock.hh:162
double getTimeOfFrame(const int frame_index)
Get start time of frame in ns since start of run for a given frame index.
Definition JDAQClock.hh:185
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
Definition JDAQ.hh:26
std::map< int, range_type > map_type
static const int PMT_DISABLE
KM3NeT Data Definitions v3.6.0 https://git.km3net.de/common/km3net-dataformat.
Definition pmt_status.hh:12
Auxiliary data structure for sequence of same character.
Definition JManip.hh:330
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
Type definition of range.
Definition JHead.hh:43
Detector file.
Definition JHead.hh:227
Auxiliary class for TDC constraints.
Definition JTDC_t.hh:39
range_type equal_range(const int id)
Get range of constraints for given module.
Definition JTDC_t.hh:101
bool is_valid(const bool option=false) const
Check validity of TDC constrants.
Definition JTDC_t.hh:171
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
Auxiliary class to select ROOT class based on class name.
Auxiliary class for defining the range of iterations of objects.
Definition JLimit.hh:45
static counter_type max()
Get maximum counter value.
Definition JLimit.hh:128
Auxiliary class for ROOT I/O of application specific meta data.
Definition JMeta.hh:72