53{
57
58
59
60
61
63 JLimit_t& numberOfEvents = inputFile.getLimit();
65 string detectorFile;
66 Double_t TMax_ns;
71 int filterLevel;
72 int muteChannel;
73 bool twoDim;
74 bool monitorOccupancy;
75 bool notJoin;
76
77 try {
78
79 JParser<> zap(
"Auxiliary program to estimate PMT and hit multiplicities.");
80
91 zap[
'F'] =
make_field(filterLevel,
"0 = raw data; 1 = filtered data; 2+ = only clean frames") = 1;
92 zap[
'D'] =
make_field(twoDim,
"2D mode for background subtraction");
93 zap[
'O'] =
make_field(monitorOccupancy,
"produces 2D PMT vs multiplicity plots");
94 zap[
'j'] =
make_field(notJoin,
"do not join consecutive hits on same PMT (diagnostic purpose)");
95
96 zap(argc, argv);
97 }
98 catch(const exception &error) {
99 FATAL(error.what() << endl);
100 }
101
102
103
104
105
106
108 bool hasTriggerParameters = true;
109
110 try {
112 NOTICE(
"Get trigger parameters from input." << endl);
113 }
115 ERROR(
"No trigger parameters from input." << endl);
116 hasTriggerParameters = false;
117
118 }
119
120
121 int prescale = 1;
122
123 if (hasTriggerParameters) {
124
125
126
127 if (parameters.writeL1.prescale > 0 && (selector == "JDAQTimeslice" || selector == "JDAQTimesliceL1") &&
128 parameters.TMaxLocal_ns < TMax_ns) {
129 FATAL(
"TMax_ns (-T) " << TMax_ns <<
" ns is larger than in the trigger " << parameters.TMaxLocal_ns <<
" ns." << endl);
130 }
131
132 if (selector == "JDAQTimeslice") {
133
134 prescale = (parameters.writeL1.prescale != 0) ? parameters.writeL1.prescale : parameters.writeL0.prescale;
135 }
136
137 if (selector == "JDAQTimesliceL0") {
138 prescale = parameters.writeL0.prescale;
139 }
140
141 if (selector == "JDAQTimesliceL1") {
142 prescale = parameters.writeL1.prescale;
143 }
144
145 if (selector == "JDAQTimesliceSN") {
146 prescale = parameters.writeSN.prescale;
147 }
148
149 }
150
151 if (prescale == 0) {
152 FATAL(
"[R] According to trigger parameters, the " << selector <<
" stream should be empty.");
153 } else {
154 NOTICE(
"[R] Prescale factor is " << prescale << endl);
155 }
156
158 FATAL(
"Invalid time over threshold " << totVeto_ns << endl);
159 }
160
161
162
163
164
166
167 try {
169 }
172 }
173
175 FATAL(
"Empty detector." << endl);
176
179
181
183
185
186
187
188
189
190
192 NOTICE(
"[R] Tmax [ns] " << TMax_ns << endl);
193 NOTICE(
"[R] Rate range [Hz] " << rateVeto_Hz << endl);
194 NOTICE(
"[R] ToT range [ns] " << totVeto_ns << endl);
195 NOTICE(
"[R] Timeslice stream " << selector << endl);
196 NOTICE(
"[R] Detector file has " << detSize <<
" DOMs." << endl);
197
198
199
200
201
202
203 TH1D* h_livetime = new TH1D("LiveTime", "Livetime", 1 + detSize, 0.0, 1.0 + detSize);
204
205 int ibin = 1;
206 h_livetime->GetXaxis()->SetBinLabel(ibin++, "Nominal");
207 for (JDetector::iterator module =
detector.begin(); module !=
detector.end(); ++module) {
208 h_livetime->GetXaxis()->SetBinLabel(ibin++, TString(
getLabel(*module)));
209 }
210
211 h_livetime->GetYaxis()->SetTitle("Livetime (s)");
212
213
214
216 const double xmin = -0.5;
217 const double xmax = nx - 0.5;
218
221 JManager_t
H(
new TH1D(
"%_H", NULL, nx, xmin, xmax));
H->Sumw2();
222 JManager_t P(new TH1D("%_P", NULL, nx, xmin, xmax)); P->Sumw2();
223
224 JManager2D_t HT(new TH2D("%_HT", NULL, nx, xmin, xmax, 100, -TMax_ns, +TMax_ns));
225
226 TH2D* CO_proto = new TH2D("%_CO", NULL, NUMBER_OF_PMTS, 0.5, 0.5 + NUMBER_OF_PMTS, NUMBER_OF_PMTS, -0.5, -0.5 + NUMBER_OF_PMTS);
227
228 if (monitorOccupancy) {
231 }
232
233 JManager2D_t CO(CO_proto);
234
235
236 const int sn_mul_min = 6;
237 TH1D* time_distr = new TH1D("T_distr", "Coincidence (M > 6) distribution in timeslice", 100, 0, 1e8);
238
239
240
241
242
243
245
246
247
248 int count_HRV = 0;
249 int count_FAF = 0;
250 int count_URV = 0;
251
252 int treeMismatchCount = 0;
253
254
255
256
257
259
261
263
267
269
270
271 bool hasMCHeader = true ;
272 bool isMupage = false;
273 double liveTimeMC = 0 ;
274 double liveTimeMCErr = 0 ;
275
277
278 try {
280 }
282 hasMCHeader = false;
283 }
284
285 NOTICE(
"[R] File source is " << (hasMCHeader ?
"MC" :
"DAQ") <<
"." << endl);
286 if (hasMCHeader) {
289 if (liveTimeMC > 0) {
290 isMupage = true;
291 NOTICE(
"[R] mupage live time is (" << liveTimeMC <<
" +/- " << liveTimeMCErr <<
") s." << endl);
292 } else {
293 NOTICE(
"[R] aanet header found but mupage live time is zero. Likely not a mupage file.");
294 }
295 }
296
297
298
299 int runNumber = 0;
300
301 if (summaryFile.hasNext()) {
302 runNumber = summaryFile.getEntry(0)->getRunNumber();
303 NOTICE(
"[R] Processing run #" << runNumber << endl);
304 summaryFile.rewind();
305 }
306
307
308
309
310
311 for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
312
313 STATUS(
"Entry: " << setw(10) << counter <<
'\r');
314
316 const Long64_t index = summaryFile.find(*timeslice);
317 JDAQSummaryslice* summary = (index != -1 ? summaryFile.getEntry(index) : NULL);
318
319 summaryRouter.
update(summary);
320
321
322 if (summary != NULL) {
324 ++treeMismatchCount;
325 ERROR(
"[R>T] Frame indices do not match [counter=" << counter <<
", timeslice=" << timeslice->
getFrameIndex() <<
", summaryslice=" << summary->
getFrameIndex() <<
"]" << endl);
326 if (treeMismatchCount > 1) {
327 FATAL(
"ROOT TTrees not aligned. Run #" << runNumber << endl);
328 }
329 continue;
330 }
331 }
332
333
334
335
336 for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
337
338
339
340
341
342 int moduleID = super_frame->getModuleID();
343
344 if (!router.hasModule(moduleID)) {
345
346 continue;
347 }
348
350
351
352
353
354
357
358 if (summary == NULL) {
359 WARNING(
"[R>T>F] Missing summary for timeslice." << endl);
360
363 }
364
365
367
368 } else {
369
370
371
373
374 if (!summaryRouter.
hasSummaryFrame(super_frame->getModuleIdentifier())) {
375 summaryRouter.print(cout);
376 FATAL(
"[R>D>F] Summary frame not found, module = " << super_frame->getModuleID() << endl);
377 } else {
378 summary_frame = summaryRouter.
getSummaryFrame(super_frame->getModuleID());
379 }
380
381
383 rate_Hz[i] = summary_frame.
getRate(i);
384 }
385
386
388
389 }
390
391
392
393
394
395
397 WARNING(
"[R>D>F] DAQ status not okay in timeslice #" << timeslice->
getFrameIndex() <<
", DOM=" << super_frame->getModuleID() << endl);
398 continue;
399 }
400
402 bool moduleLiveTime = 0;
403
404
405 for (vector<double>::const_iterator r = rate_Hz.begin(); !moduleLiveTime && r != rate_Hz.end(); ++r) {
406 if ((*r) > 0) {
407 moduleLiveTime = 1;
408 }
409 }
410
413 int frame_URV_count = 0;
414
416 if ((veto[i] = (veto[i] || !rateVeto_Hz(rate_Hz[i])))) {
417 frame_URV_count++;
418 }
419 }
420
421 count_HRV += frame_HRV_count;
422 count_FAF += frame_FAF_count;
423 count_URV += frame_URV_count;
424
425 int frame_VETO_count = frame_HRV_count + frame_FAF_count + frame_URV_count;
426
427
428 if (filterLevel >= 2 && frame_VETO_count > 0) {
429 moduleLiveTime = 0;
430 fill(veto.begin(), veto.end(), true);
431 } else if (filterLevel >= 1 && frame_HRV_count + frame_FAF_count > 0) {
434 }
435 }
436
437
438
439 if (muteChannel != -1) {
440
441
442 int mute_VETO_count =
445 !rateVeto_Hz(rate_Hz[muteChannel] );
446
447 if (mute_VETO_count == frame_VETO_count) {
448 moduleLiveTime = 1;
449 fill(veto.begin(), veto.end(), false);
450 }
451
452 veto[muteChannel] = true;
453
454 }
455
456
457 const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
459 const TString moduleLabel =
getLabel(module);
460
461
462 if (moduleLiveTime) {
463 h_livetime->Fill(moduleLabel,
getFrameTime() * 1e-9 * moduleLiveTime);
464 }
465
466
467
468
469
470 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, router.getModule(super_frame->getModuleID()));
471
472 if (!notJoin) {
473 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
474 i->join(match);
475 }
476 }
477
478 JSuperFrame1D_t&
data = JSuperFrame1D_t::multiplex(buffer);
479
480 filteredData.clear();
481
482 for (vector<JHitR0>::const_iterator h =
data.begin(); h !=
data.end(); ++h) {
483
484 const int pmt = h->getPMT();
485
486
487 if (!veto[pmt] &&
488 rateVeto_Hz(rate_Hz[pmt]) &&
489 totVeto_ns(h->getToT()) ) {
490 filteredData.push_back(*h);
491 }
492 }
493
494
495
496 if (filteredData.size() > 1) {
497
498 TH1D* h_hit =
H[moduleLabel];
499 TH1D* h_pmt = P[moduleLabel];
500
501 TH2D* h2_hit = HT[moduleLabel];
502 TH2D* h2_co = CO[moduleLabel];
503
504 sort(filteredData.begin(), filteredData.end());
505
506
507
508
509 for (vector<JHitR0>::const_iterator p = filteredData.begin(); p != filteredData.end(); ) {
510
511 vector<JHitR0>::const_iterator q = p;
512
513
514 while (++q != filteredData.end() && q->getT() - p->getT() <= TMax_ns ) {}
515
516
517 int hit_multiplicity =
distance(p,q);
519 for (vector<JHitR0>::const_iterator __p = p ; __p != q; ++__p ) { pmthit.insert(__p->getPMT()) ; }
520 int pmt_multiplicity = pmthit.size() ;
521
522
523 h_pmt->Fill(pmt_multiplicity);
524 h_hit->Fill(hit_multiplicity);
525
526 if (twoDim && hit_multiplicity > 1) {
527 const double W =
factorial(hit_multiplicity, 2);
528 for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
529 for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
531 h2_hit->Fill(hit_multiplicity, dt, 1.0/W);
532 }
533 }
534 }
535
536 if (monitorOccupancy) {
539 TString pmtLabel(pmtAddress.
toString());
540 h2_co->Fill(pmt_multiplicity, pmtLabel, 1.0/pmt_multiplicity);
541 }
542 }
543
544
545 if (hit_multiplicity >= sn_mul_min) {
546 time_distr->Fill(p->getT());
547 }
548
549 p = q;
550
551 }
552
553 }
554
555 }
556
557 }
558
559 NOTICE(
"[R] " << counter <<
" timeslices processed." << endl);
560
561 if (isMupage) {
562 h_livetime->Scale(liveTimeMC / (
getFrameTime() * 1e-9 * counter * prescale));
563 }
564
565
566 double nominalLiveTime = (isMupage ? liveTimeMC : counter *
getFrameTime() * 1e-9) / prescale;
567 h_livetime->Fill("Nominal", nominalLiveTime);
568
569 out.cd();
570
571 int nLiveDOMs =
H.size();
572
573 h_livetime->Write();
574
575
576
577
578 P.Write(out);
579
580 if (twoDim) {
581 HT.Write(out);
582 }
583
584 if (monitorOccupancy) {
585 CO.Write(out);
586 }
587
588 time_distr->Write();
589
590 double rate_HRV = count_HRV / (1.0 * counter);
591 double rate_FAF = count_FAF / (1.0 * counter);
592 double rate_URV = count_URV / (1.0 * counter);
593 double biolumIndex = (rate_HRV + rate_FAF) / (nLiveDOMs * NUMBER_OF_PMTS);
594 NOTICE(
"[R] Average HRV count per timeslice: " << rate_HRV << endl);
595 NOTICE(
"[R] Average FIFO-almost-full count per timeslice: " << rate_FAF << endl);
596 NOTICE(
"[R] Average user rate veto count per timeslice: " << rate_URV << endl);
597 NOTICE(
"[R] " <<
H.size() <<
" DOMs were active in the run." << endl);
598 NOTICE(
"[R] Bioluminescence proxy index: " << biolumIndex << endl);
599
601
602 out.Close();
603
605
606 return (counter ? 0 : 1);
607
608}
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
JAANET::livetime livetime
Lookup table for PMT addresses in detector.
const JModuleAddressMap & get(const int id) const
Get module address map.
virtual const JModuleAddressMap & getDefaultModuleAddressMap() const =0
Get default module address map.
Lookup table for PMT addresses in optical module.
const JPMTAddressTranslator & getAddressTranslator(const int tdc) const
Get PMT address translator.
Address of module in detector data structure.
Router for direct addressing of module data in detector data structure.
Data structure for a composite optical module.
Data structure for PMT physical address.
std::string toString() const
Convert PMT physical address to string.
Auxiliary class for multiplexing object iterators.
Utility class to parse command line options.
Auxiliary class to manage set of compatible ROOT objects (e.g. histograms) using unique keys.
General purpose class for object reading from a list of file names.
Router for fast addressing of summary data in KM3NETDAQ::JDAQSummaryslice data structure as a functio...
bool hasSummaryFrame(const JDAQModuleIdentifier &module) const
Has summary frame.
const JDAQSummaryFrame & getSummaryFrame(const JDAQModuleIdentifier &module) const
Get summary frame.
void update(const JDAQSummaryslice *ps)
Update router.
Template definition for direct access of elements in ROOT TChain.
1-dimensional frame with time calibrated data from one optical module.
2-dimensional frame with time calibrated data from one optical module.
int getFrameIndex() const
Get frame index.
bool testFIFOStatus() const
Test FIFO status.
int countHighRateVeto() const
Count high-rate veto status.
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
int countFIFOStatus() const
Count FIFO status.
bool testHighRateVeto() const
Test high-rate veto status.
bool testDAQStatus() const
Test DAQ status of packets.
Data storage class for rate measurements of all PMTs in one module.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
std::string getLabel(const JLocation &location)
Get module label for monitoring and other applications.
void load(const std::string &file_name, JDetector &detector)
Load detector from input file.
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
void setAxisLabels(TAxis *axis, const JModuleAddressMap &memo)
Set axis with PMT address labels.
long long int factorial(const long long int n)
Determine factorial.
static const double H
Planck constant [eV s].
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.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
static const JModuleCounter getNumberOfModules
Function object to count unique modules.
KM3NeT DAQ data structures and auxiliaries.
double getFrameTime()
Get frame time duration.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
double numberOfSeconds
Live time [s].
double errorOfSeconds
Uncertainty on live time [s].
Auxiliary class to select ROOT class based on class name.
Auxiliary class for defining the range of iterations of objects.
static counter_type max()
Get maximum counter value.