52 int main(
int argc,
char **argv)
62 JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
67 JRange<double> rateVeto_Hz;
68 JRange<double> totVeto_ns;
69 JROOTClassSelector selector;
74 bool monitorOccupancy;
79 JParser<> zap(
"Auxiliary program to estimate PMT and hit multiplicities.");
84 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
86 zap[
'V'] =
make_field(rateVeto_Hz) = JRange<double>(0, 1e5);
87 zap[
't'] =
make_field(totVeto_ns) = JRange<double>(0, 1e4);
88 zap[
'C'] =
make_field(selector) = getROOTClassSelection<JDAQTimesliceTypes_t>();
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)");
98 catch(
const exception &error) {
99 FATAL(error.what() << endl);
109 bool hasTriggerParameters =
true;
113 NOTICE(
"Get trigger parameters from input." << endl);
115 catch(
const JException& error) {
116 ERROR(
"No trigger parameters from input." << endl);
117 hasTriggerParameters =
false;
124 if (hasTriggerParameters) {
128 if (parameters.
writeL1.
prescale > 0 && (selector ==
"JDAQTimeslice" || selector ==
"JDAQTimesliceL1") &&
130 FATAL(
"TMax_ns (-T) " << TMax_ns <<
" ns is larger than in the trigger " << parameters.
TMaxLocal_ns <<
" ns." << endl);
133 if (selector ==
"JDAQTimeslice") {
138 if (selector ==
"JDAQTimesliceL0") {
142 if (selector ==
"JDAQTimesliceL1") {
146 if (selector ==
"JDAQTimesliceSN") {
153 FATAL(
"[R] According to trigger parameters, the " << selector <<
" stream should be empty.");
155 NOTICE(
"[R] Prescale factor is " << prescale << endl);
158 if (!totVeto_ns.is_valid()) {
159 FATAL(
"Invalid time over threshold " << totVeto_ns << endl);
169 load(detectorFile, detector);
171 catch(
const JException& error) {
175 if (detector.empty())
176 FATAL(
"Empty detector." << endl);
178 const JModuleRouter router(detector);
179 JSummaryRouter summaryRouter;
183 int detID = detector.getID();
193 NOTICE(
"[R] Tmax [ns] " << TMax_ns << endl);
194 NOTICE(
"[R] Rate range [Hz] " << rateVeto_Hz << endl);
195 NOTICE(
"[R] ToT range [ns] " << totVeto_ns << endl);
196 NOTICE(
"[R] Timeslice stream " << selector << endl);
197 NOTICE(
"[R] Detector file has " << detSize <<
" DOMs." << endl);
204 TH1D* h_livetime =
new TH1D(
"LiveTime",
"Livetime", 1 + detSize, 0.0, 1.0 + detSize);
207 h_livetime->GetXaxis()->SetBinLabel(ibin++,
"Nominal");
208 for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
209 h_livetime->GetXaxis()->SetBinLabel(ibin++, TString(
getModuleLabel(*module)));
212 h_livetime->GetYaxis()->SetTitle(
"Livetime (s)");
217 const double xmin = -0.5;
218 const double xmax = nx - 0.5;
222 JManager_t
H(
new TH1D(
"%_H", NULL, nx, xmin, xmax));
H->Sumw2();
223 JManager_t P(
new TH1D(
"%_P", NULL, nx, xmin, xmax)); P->Sumw2();
225 JManager2D_t HT(
new TH2D(
"%_HT", NULL, nx, xmin, xmax, 100, -TMax_ns, +TMax_ns));
229 if (monitorOccupancy) {
230 JModuleAddressMap memo = getDetectorAddressMap<JKM3NeT_t>().getDefaultModuleAddressMap();
234 JManager2D_t CO(CO_proto);
237 const int sn_mul_min = 6;
238 TH1D* time_distr =
new TH1D(
"T_distr",
"Coincidence (M > 6) distribution in timeslice", 100, 0, 1e8);
253 int treeMismatchCount = 0;
259 JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
261 JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
265 typedef JSuperFrame2D<JHitR0> JSuperFrame2D_t;
266 typedef JSuperFrame1D<JHitR0> JSuperFrame1D_t;
267 const JMatchL0<JHitR0> match(TMax_ns);
272 bool hasMCHeader = true ;
273 bool isMupage =
false;
274 double liveTimeMC = 0 ;
275 double liveTimeMCErr = 0 ;
286 NOTICE(
"[R] File source is " << (hasMCHeader ?
"MC" :
"DAQ") <<
"." << endl);
290 if (liveTimeMC > 0) {
292 NOTICE(
"[R] mupage live time is (" << liveTimeMC <<
" +/- " << liveTimeMCErr <<
") s." << endl);
294 NOTICE(
"[R] aanet header found but mupage live time is zero. Likely not a mupage file.");
302 if (summaryFile.hasNext()) {
303 runNumber = summaryFile.getEntry(0)->getRunNumber();
304 NOTICE(
"[R] Processing run #" << runNumber << endl);
305 summaryFile.rewind();
312 for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
314 STATUS(
"Entry: " << setw(10) << counter <<
'\r');
317 const Long64_t index = summaryFile.find(*timeslice);
318 JDAQSummaryslice* summary = (index != -1 ? summaryFile.getEntry(index) : NULL);
320 summaryRouter.update(summary);
323 if (summary != NULL) {
326 ERROR(
"[R>T] Frame indices do not match [counter=" << counter <<
", timeslice=" << timeslice->
getFrameIndex() <<
", summaryslice=" << summary->
getFrameIndex() <<
"]" << endl);
327 if (treeMismatchCount > 1) {
328 FATAL(
"ROOT TTrees not aligned. Run #" << runNumber << endl);
337 for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
343 int moduleID = super_frame->getModuleID();
345 if (!router.hasModule(moduleID)) {
350 JModuleAddressMap moduleAddressMap = detectorAddressMap.get(moduleID);
359 if (summary == NULL) {
360 WARNING(
"[R>T>F] Missing summary for timeslice." << endl);
375 if (!summaryRouter.hasSummaryFrame(super_frame->getModuleIdentifier())) {
376 summaryRouter.print(cout);
377 FATAL(
"[R>D>F] Summary frame not found, module = " << super_frame->getModuleID() << endl);
379 summary_frame = summaryRouter.getSummaryFrame(super_frame->getModuleID());
384 rate_Hz[i] = summary_frame.
getRate(i);
398 WARNING(
"[R>D>F] DAQ status not okay in timeslice #" << timeslice->
getFrameIndex() <<
", DOM=" << super_frame->getModuleID() << endl);
403 bool moduleLiveTime = 0;
414 int frame_URV_count = 0;
417 if (veto[i] = (veto[i] || !rateVeto_Hz(rate_Hz[i]))) {
422 count_HRV += frame_HRV_count;
423 count_FAF += frame_FAF_count;
424 count_URV += frame_URV_count;
426 int frame_VETO_count = frame_HRV_count + frame_FAF_count + frame_URV_count;
429 if (filterLevel >= 2 && frame_VETO_count > 0) {
431 fill(veto.begin(), veto.end(),
true);
432 }
else if (filterLevel >= 1 && frame_HRV_count + frame_FAF_count > 0) {
440 if (muteChannel != -1) {
443 int mute_VETO_count =
446 !rateVeto_Hz(rate_Hz[muteChannel] );
448 if (mute_VETO_count == frame_VETO_count) {
450 fill(veto.begin(), veto.end(),
false);
453 veto[muteChannel] =
true;
458 const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
459 const JModule& module = detector.getModule(address);
463 if (moduleLiveTime) {
464 h_livetime->Fill(moduleLabel,
getFrameTime() * 1e-9 * moduleLiveTime);
471 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, router.getModule(super_frame->getModuleID()));
474 for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
479 JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
483 filteredData.clear();
487 const int pmt = h->getPMT();
491 rateVeto_Hz(rate_Hz[pmt]) &&
492 totVeto_ns(h->getToT()) ) {
493 filteredData.push_back(*h);
499 if (filteredData.size() > 1) {
501 TH1D* h_hit =
H[moduleLabel];
502 TH1D* h_pmt = P[moduleLabel];
504 TH2D* h2_hit = HT[moduleLabel];
505 TH2D* h2_co = CO[moduleLabel];
507 sort(filteredData.begin(), filteredData.end());
517 while (++q != filteredData.end() && q->getT() - p->getT() <= TMax_ns ) {}
520 int hit_multiplicity =
distance(p,q);
523 int pmt_multiplicity = pmthit.size() ;
526 h_pmt->Fill(pmt_multiplicity);
527 h_hit->Fill(hit_multiplicity);
529 if (twoDim && hit_multiplicity > 1) {
530 const double W =
factorial(hit_multiplicity, 2);
533 double dt = JCombinatorics::getSign(__p->getPMT(), __q->getPMT()) * (__q->getT() - __p->getT());
534 h2_hit->Fill(hit_multiplicity, dt, 1.0/W);
539 if (monitorOccupancy) {
541 JPMTPhysicalAddress pmtAddress = moduleAddressMap.getAddressTranslator(*
r);
542 TString pmtLabel(pmtAddress.toString());
543 h2_co->Fill(pmt_multiplicity, pmtLabel, 1.0/pmt_multiplicity);
548 if (hit_multiplicity >= sn_mul_min) {
549 time_distr->Fill(p->getT());
562 NOTICE(
"[R] " << counter <<
" timeslices processed." << endl);
565 h_livetime->Scale(liveTimeMC / (
getFrameTime() * 1e-9 * counter * prescale));
569 double nominalLiveTime = (isMupage ? liveTimeMC : counter *
getFrameTime() * 1e-9) / prescale;
570 h_livetime->Fill(
"Nominal", nominalLiveTime);
574 int nLiveDOMs =
H.size();
587 if (monitorOccupancy) {
593 double rate_HRV = count_HRV / (1.0 * counter);
594 double rate_FAF = count_FAF / (1.0 * counter);
595 double rate_URV = count_URV / (1.0 * counter);
596 double biolumIndex = (rate_HRV + rate_FAF) / (nLiveDOMs *
NUMBER_OF_PMTS);
597 NOTICE(
"[R] Average HRV count per timeslice: " << rate_HRV << endl);
598 NOTICE(
"[R] Average FIFO-almost-full count per timeslice: " << rate_FAF << endl);
599 NOTICE(
"[R] Average user rate veto count per timeslice: " << rate_URV << endl);
600 NOTICE(
"[R] " <<
H.size() <<
" DOMs were active in the run." << endl);
601 NOTICE(
"[R] Bioluminescence proxy index: " << biolumIndex << endl);
609 return (counter ? 0 : 1);