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.