54 using namespace KM3NETDAQ;
61 JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
66 JRange<double> rateVeto_Hz;
67 JRange<double> totVeto_ns;
68 JROOTClassSelector selector;
73 bool monitorOccupancy;
78 JParser<> zap(
"Auxiliary program to estimate PMT and hit multiplicities.");
83 zap[
'n'] =
make_field(numberOfEvents) = JLimit::max();
85 zap[
'V'] =
make_field(rateVeto_Hz) = JRange<double>(0, 1e5);
86 zap[
't'] =
make_field(totVeto_ns) = JRange<double>(0, 1e4);
87 zap[
'C'] =
make_field(selector) = getROOTClassSelection<JDAQTimesliceTypes_t>();
90 zap[
'F'] =
make_field(filterLevel,
"0 = raw data; 1 = filtered data; 2+ = only clean frames") = 1;
91 zap[
'D'] =
make_field(twoDim,
"2D mode for background subtraction");
92 zap[
'O'] =
make_field(monitorOccupancy,
"produces 2D PMT vs multiplicity plots");
93 zap[
'j'] =
make_field(notJoin,
"do not join consecutive hits on same PMT (diagnostic purpose)");
97 catch(
const exception &error) {
98 FATAL(error.what() << endl);
108 bool hasTriggerParameters =
true;
112 NOTICE(
"Get trigger parameters from input." << endl);
114 catch(
const JException& error) {
115 ERROR(
"No trigger parameters from input." << endl);
116 hasTriggerParameters =
false;
123 if (hasTriggerParameters) {
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);
132 if (selector ==
"JDAQTimeslice") {
134 prescale = (parameters.writeL1.prescale != 0) ? parameters.writeL1.prescale : parameters.writeL0.prescale;
137 if (selector ==
"JDAQTimesliceL0") {
138 prescale = parameters.writeL0.prescale;
141 if (selector ==
"JDAQTimesliceL1") {
142 prescale = parameters.writeL1.prescale;
145 if (selector ==
"JDAQTimesliceSN") {
146 prescale = parameters.writeSN.prescale;
152 FATAL(
"[R] According to trigger parameters, the " << selector <<
" stream should be empty.");
154 NOTICE(
"[R] Prescale factor is " << prescale << endl);
157 if (!totVeto_ns.is_valid()) {
158 FATAL(
"Invalid time over threshold " << totVeto_ns << endl);
168 load(detectorFile, detector);
170 catch(
const JException& error) {
174 if (detector.empty())
175 FATAL(
"Empty detector." << endl);
177 const JModuleRouter router(detector);
178 JSummaryRouter summaryRouter;
182 int detID = detector.getID();
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);
203 TH1D* h_livetime =
new TH1D(
"LiveTime",
"Livetime", 1 + detSize, 0.0, 1.0 + detSize);
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(
getModuleLabel(*module)));
211 h_livetime->GetYaxis()->SetTitle(
"Livetime (s)");
216 const double xmin = -0.5;
217 const double xmax = nx - 0.5;
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();
224 JManager2D_t HT(
new TH2D(
"%_HT", NULL, nx, xmin, xmax, 100, -TMax_ns, +TMax_ns));
228 if (monitorOccupancy) {
229 JModuleAddressMap memo = getDetectorAddressMap<JKM3NeT_t>().getDefaultModuleAddressMap();
233 JManager2D_t CO(CO_proto);
236 const int sn_mul_min = 6;
237 TH1D* time_distr =
new TH1D(
"T_distr",
"Coincidence (M > 6) distribution in timeslice", 100, 0, 1e8);
252 int treeMismatchCount = 0;
258 JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
260 JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
266 const JMatch_t match(TMax_ns);
271 bool hasMCHeader = true ;
272 bool isMupage =
false;
273 double liveTimeMC = 0 ;
274 double liveTimeMCErr = 0 ;
285 NOTICE(
"[R] File source is " << (hasMCHeader ?
"MC" :
"DAQ") <<
"." << endl);
287 liveTimeMC = mc_header.livetime.numberOfSeconds;
288 liveTimeMCErr = mc_header.livetime.errorOfSeconds;
289 if (liveTimeMC > 0) {
291 NOTICE(
"[R] mupage live time is (" << liveTimeMC <<
" +/- " << liveTimeMCErr <<
") s." << endl);
293 NOTICE(
"[R] aanet header found but mupage live time is zero. Likely not a mupage file.");
301 if (summaryFile.hasNext()) {
302 runNumber = summaryFile.getEntry(0)->getRunNumber();
303 NOTICE(
"[R] Processing run #" << runNumber << endl);
304 summaryFile.rewind();
311 for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
313 STATUS(
"Entry: " << setw(10) << counter <<
'\r');
316 const Long64_t index = summaryFile.find(*timeslice);
317 JDAQSummaryslice* summary = (index != -1 ? summaryFile.getEntry(index) : NULL);
319 summaryRouter.update(summary);
322 if (summary != NULL) {
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);
336 for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
342 int moduleID = super_frame->getModuleID();
344 if (!router.hasModule(moduleID)) {
349 JModuleAddressMap moduleAddressMap = detectorAddressMap.get(moduleID);
358 if (summary == NULL) {
359 WARNING(
"[R>T>F] Missing summary for timeslice." << endl);
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);
378 summary_frame = summaryRouter.getSummaryFrame(super_frame->getModuleID());
383 rate_Hz[i] = summary_frame.
getRate(i);
397 WARNING(
"[R>D>F] DAQ status not okay in timeslice #" << timeslice->
getFrameIndex() <<
", DOM=" << super_frame->getModuleID() << endl);
402 bool moduleLiveTime = 0;
413 int frame_URV_count = 0;
416 if (veto[i] = (veto[i] || !rateVeto_Hz(rate_Hz[i]))) {
421 count_HRV += frame_HRV_count;
422 count_FAF += frame_FAF_count;
423 count_URV += frame_URV_count;
425 int frame_VETO_count = frame_HRV_count + frame_FAF_count + frame_URV_count;
428 if (filterLevel >= 2 && frame_VETO_count > 0) {
430 fill(veto.begin(), veto.end(),
true);
431 }
else if (filterLevel >= 1 && frame_HRV_count + frame_FAF_count > 0) {
439 if (muteChannel != -1) {
442 int mute_VETO_count =
445 !rateVeto_Hz(rate_Hz[muteChannel] );
447 if (mute_VETO_count == frame_VETO_count) {
449 fill(veto.begin(), veto.end(),
false);
452 veto[muteChannel] =
true;
457 const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
458 const JModule& module = detector.getModule(address);
462 if (moduleLiveTime) {
463 h_livetime->Fill(moduleLabel,
getFrameTime() * 1e-9 * moduleLiveTime);
470 JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, router.getModule(super_frame->getModuleID()));
482 filteredData.clear();
486 const int pmt = h->getPMT();
490 rateVeto_Hz(rate_Hz[pmt]) &&
491 totVeto_ns(h->getToT()) ) {
492 filteredData.push_back(*h);
498 if (filteredData.size() > 1) {
500 TH1D* h_hit =
H[moduleLabel];
501 TH1D* h_pmt = P[moduleLabel];
503 TH2D* h2_hit = HT[moduleLabel];
504 TH2D* h2_co = CO[moduleLabel];
506 sort(filteredData.begin(), filteredData.end());
516 while (++q != filteredData.end() && q->getT() - p->getT() <= TMax_ns ) {}
519 int hit_multiplicity = distance(p,q);
522 int pmt_multiplicity = pmthit.size() ;
525 h_pmt->Fill(pmt_multiplicity);
526 h_hit->Fill(hit_multiplicity);
528 if (twoDim && hit_multiplicity > 1) {
529 const double W =
factorial(hit_multiplicity, 2);
532 double dt = JCombinatorics::getSign(__p->getPMT(), __q->getPMT()) * (__q->getT() - __p->getT());
533 h2_hit->Fill(hit_multiplicity, dt, 1.0/W);
538 if (monitorOccupancy) {
540 JPMTPhysicalAddress pmtAddress = moduleAddressMap.getAddressTranslator(*r);
541 TString pmtLabel(pmtAddress.toString());
542 h2_co->Fill(pmt_multiplicity, pmtLabel, 1.0/pmt_multiplicity);
547 if (hit_multiplicity >= sn_mul_min) {
548 time_distr->Fill(p->getT());
561 NOTICE(
"[R] " << counter <<
" timeslices processed." << endl);
564 h_livetime->Scale(liveTimeMC / (
getFrameTime() * 1e-9 * counter * prescale));
568 double nominalLiveTime = (isMupage ? liveTimeMC : counter *
getFrameTime() * 1e-9) / prescale;
569 h_livetime->Fill(
"Nominal", nominalLiveTime);
573 int nLiveDOMs =
H.size();
586 if (monitorOccupancy) {
592 double rate_HRV = count_HRV / (1.0 * counter);
593 double rate_FAF = count_FAF / (1.0 * counter);
594 double rate_URV = count_URV / (1.0 * counter);
595 double biolumIndex = (rate_HRV + rate_FAF) / (nLiveDOMs *
NUMBER_OF_PMTS);
596 NOTICE(
"[R] Average HRV count per timeslice: " << rate_HRV << endl);
597 NOTICE(
"[R] Average FIFO-almost-full count per timeslice: " << rate_FAF << endl);
598 NOTICE(
"[R] Average user rate veto count per timeslice: " << rate_URV << endl);
599 NOTICE(
"[R] " <<
H.size() <<
" DOMs were active in the run." << endl);
600 NOTICE(
"[R] Bioluminescence proxy index: " << biolumIndex << endl);
608 return (counter ? 0 : 1);
Utility class to parse command line options.
Data structure for all trigger parameters.
double getRate(const int tdc, const double factor=1.0) const
Get count rate.
int countFIFOStatus() const
Count FIFO status.
Auxiliary class to manage set of histograms.
Long64_t counter_type
Type definition for counter.
JSuperFrame2D< hit_type > JSuperFrame2D_t
long long int factorial(const long long int n)
Determine factorial.
int getFrameIndex() const
Get frame index.
Head getHeader(const JMultipleFileScanner_t &file_list)
Get Monte Carlo header.
JLimit JLimit_t
Type definition of limit.
Exception for null pointer operation.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
double getFrameTime()
Get frame time duration.
Data storage class for rate measurements of all PMTs in one module.
void setAxisLabels(TAxis *axis, const JModuleAddressMap &memo)
initialize axis with PMT address labels
void load(const JString &file_name, JDetector &detector)
Load detector from input file.
JDetectorAddressMap & getDetectorAddressMap()
Get detector address map.
bool testHighRateVeto() const
Test high-rate veto status.
std::string getModuleLabel(const JModuleLocation &location)
Get module label (DU-floor) for JMonitor applications.
int countHighRateVeto() const
Count high-rate veto status.
std::vector< frame_type >::iterator iterator
2-dimensional frame with time calibrated data from one optical module.
const JLimit & getLimit() const
Get limit.
static const int NUMBER_OF_PMTS
Total number of PMTs in module.
bool testDAQStatus() const
Test DAQ status of packets.
const JDAQFrameStatus & getDAQFrameStatus() const
Get DAQ frame status.
int getNumberOfModules(const JDetector &detector)
Get number of modules.
JTriggerParameters getTriggerParameters(const JMultipleFileScanner_t &file_list)
Get trigger parameters.
JSuperFrame1D< hit_type > JSuperFrame1D_t
bool testFIFOStatus() const
Test FIFO status.
bool putObject(TDirectory *dir, const T &object)
Write object to ROOT directory.