JMonitorMultiplicity.cc

Go to the documentation of this file.

#include <iostream>
#include <string>
 
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"
#include "TROOT.h"
#include "TVectorD.h"
 
#include "JAAnet/JHead.hh"
#include "JCalibrate/JCalibrateK40.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetectorAddressMap.hh"
#include "JDetector/JDetectorSupportkit.hh"
#include "JDetector/JModuleAddressMap.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "JMath/JMathToolkit.hh"
#include "Jeep/JTimer.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JLang/JString.hh"
#include "JROOT/JManager.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JROOT/JROOTClassSelector.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JTrigger/JHit.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JSummaryRouter.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTools/JRange.hh"
 
/**
 * \file
 *
 * Auxiliary program to analyse coincidence multiplicity;
 * \author mlincetto, mkarel
 */
int main(int argc, char **argv)
{
  using namespace std;  
  using namespace KM3NETDAQ;
  using namespace JPP;  
  
  //-----------------------------------------------------------
  // parameter interface
  //-----------------------------------------------------------
  
  JMultipleFileScanner<JDAQTimesliceTypes_t> inputFile;
  JLimit_t& numberOfEvents =          inputFile.getLimit();
  string             outputFile;
  string             detectorFile;
  Double_t           TMax_ns;
  JRange<double>     rateVeto_Hz;
  JRange<double>     totVeto_ns;
  JROOTClassSelector selector;
  int                debug;
  int                filterLevel;
  int                muteChannel;
  bool               twoDim;
  bool               monitorOccupancy;
  bool               notJoin;
 
  try { 
 
    JParser<> zap("Auxiliary program to estimate PMT and hit multiplicities.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "monitormultiplicity.root";
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)       = JLimit::max();
    zap['T'] = make_field(TMax_ns)              = 10.0;                        // [ns]
    zap['V'] = make_field(rateVeto_Hz)          = JRange<double>(0, 1e5);      // [Hz]
    zap['t'] = make_field(totVeto_ns)           = JRange<double>(0, 1e4);      // min and max value of ToT of hits to be monitored [ns]
    zap['C'] = make_field(selector)             = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['d'] = make_field(debug)                = 1;
    zap['M'] = make_field(muteChannel) = -1;
    zap['F'] = make_field(filterLevel, "0 = raw data; 1 = filtered data; 2+ = only clean frames") = 1;
    zap['D'] = make_field(twoDim, "2D mode for background subtraction");
    zap['O'] = make_field(monitorOccupancy, "produces 2D PMT vs multiplicity plots");
    zap['j'] = make_field(notJoin, "do not join consecutive hits on same PMT (diagnostic purpose)");
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  
  //-----------------------------------------------------------
  // check the input parameters
  //-----------------------------------------------------------
 
  JTriggerParameters parameters;
  bool hasTriggerParameters = true;
 
  try {
    parameters = getTriggerParameters(inputFile);
    NOTICE("Get trigger parameters from input." << endl);
  }
  catch(const JException& error) {
    ERROR("No trigger parameters from input." << endl);
    hasTriggerParameters = false;
 
  }
 
  // keep track of prescale factor for proper MC normalization.
  int prescale = 1;
 
  if (hasTriggerParameters) {
 
    // for L1 (most common for sea analysis) check that time window is not larger than in trigger
    
    if (parameters.writeL1.prescale > 0  && (selector == "JDAQTimeslice" || selector == "JDAQTimesliceL1") &&
        parameters.TMaxLocal_ns < TMax_ns) {
      FATAL("TMax_ns (-T) " << TMax_ns << " ns is larger than in the trigger " << parameters.TMaxLocal_ns << " ns." << endl);
    }
    
    if (selector == "JDAQTimeslice") {
      // Jpp v8 legacy compatibility
      prescale = (parameters.writeL1.prescale != 0) ? parameters.writeL1.prescale : parameters.writeL0.prescale;
    }
 
    if (selector == "JDAQTimesliceL0") {
      prescale = parameters.writeL0.prescale;
    }
 
    if (selector == "JDAQTimesliceL1") {
      prescale = parameters.writeL1.prescale;
    }
    
    if (selector == "JDAQTimesliceSN") {
      prescale = parameters.writeSN.prescale;
    }
 
  }
 
  if (prescale == 0) { 
    FATAL("[R] According to trigger parameters, the " << selector << " stream should be empty.");
  } else {
    NOTICE("[R] Prescale factor is " << prescale << endl);
  }
 
  if (!totVeto_ns.is_valid()) {
    FATAL("Invalid time over threshold " << totVeto_ns   << endl); 
  }
    
  //-----------------------------------------------------------
  // load the detector file
  //-----------------------------------------------------------
  
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  if (detector.empty())
    FATAL("Empty detector." << endl);
 
  const JModuleRouter router(detector);
  JSummaryRouter summaryRouter;
  
  int detSize = getNumberOfModules(detector);
 
  int detID = detector.getID();
 
  JDetectorAddressMap& detectorAddressMap = getDetectorAddressMap(detID);
 
 
  //-----------------------------------------------------------
  // summary of input parameters
  //-----------------------------------------------------------
  
  NOTICE("[R] Frame time [ms] "       << getFrameTime() * 1e-6 << endl);
  NOTICE("[R] Tmax       [ns] "       << TMax_ns               << endl);
  NOTICE("[R] Rate range [Hz] "       << rateVeto_Hz           << endl);
  NOTICE("[R] ToT range  [ns] "       << totVeto_ns            << endl);
  NOTICE("[R] Timeslice stream      " << selector              << endl);
  NOTICE("[R] Detector file has " << detSize << " DOMs." << endl);
  
  //-----------------------------------------------------------
  // initialise histograms
  //-----------------------------------------------------------
 
  // nominal livetime + effective livetime for each DOM
  TH1D* h_livetime = new TH1D("LiveTime", "Livetime", 1 + detSize, 0.0, 1.0 + detSize); 
  
  int ibin = 1;
  h_livetime->GetXaxis()->SetBinLabel(ibin++, "Nominal");
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
    h_livetime->GetXaxis()->SetBinLabel(ibin++, TString(getLabel(*module)));
  }
 
  h_livetime->GetYaxis()->SetTitle("Livetime (s)");
 
 
  // multiplicity histograms
  const int    nx   = 1 + NUMBER_OF_PMTS; 
  const double xmin = -0.5;
  const double xmax = nx - 0.5;
    
  typedef JManager<TString, TH1D>  JManager_t;
  typedef JManager<TString, TH2D>  JManager2D_t;
  JManager_t H(new TH1D("%_H", NULL, nx, xmin, xmax)); H->Sumw2();
  JManager_t P(new TH1D("%_P", NULL, nx, xmin, xmax)); P->Sumw2();
  
  JManager2D_t HT(new TH2D("%_HT", NULL, nx, xmin, xmax, 100, -TMax_ns, +TMax_ns));
 
  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);
 
  if (monitorOccupancy) {
    JModuleAddressMap memo = getDetectorAddressMap<JKM3NeT_t>().getDefaultModuleAddressMap();
    setAxisLabels(CO_proto->GetYaxis(), memo);
  } 
 
  JManager2D_t CO(CO_proto);
        
  // benchmark histogram for time distribution of high-multiplicities
  const int sn_mul_min = 6;
  TH1D* time_distr = new TH1D("T_distr", "Coincidence (M > 6) distribution in timeslice", 100, 0, 1e8); 
 
  //-----------------------------------------------------------
  // output file and auxiliary data structures
  //-----------------------------------------------------------
  
  // output file
  TFile out(outputFile.c_str(), "recreate");
  
  // auxiliary counters
  
  int count_HRV = 0; // high-rate-veto counter
  int count_FAF = 0; // fifo-almost-full counter
  int count_URV = 0; // user-rate-veto counter
 
  int treeMismatchCount = 0;
  
  //-----------------------------------------------------------
  // file pre-processing
  //-----------------------------------------------------------
  
  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> summaryFile(inputFile);
 
  JObjectMultiplexer<JDAQTimesliceTypes_t, JDAQTimeslice> in(inputFile, selector);
 
  vector<JHitR0> filteredData;
  
  typedef JSuperFrame2D<JHitR0> JSuperFrame2D_t;
  typedef JSuperFrame1D<JHitR0> JSuperFrame1D_t;
  const JMatchL0<JHitR0> match(TMax_ns);
 
  counter_type counter = 0;
  
  // process Monte Carlo aanet header
  bool   hasMCHeader   = true ;
  bool   isMupage      = false;
  double liveTimeMC    = 0    ;
  double liveTimeMCErr = 0    ;
  
  JHead mc_header;
  
  try {
    mc_header = (JHead)getHeader(inputFile);
  }
  catch (const JException& error) {
    hasMCHeader = false;
  }
  
  NOTICE("[R] File source is " << (hasMCHeader ? "MC" : "DAQ") << "." << endl);
  if (hasMCHeader) {
    liveTimeMC = mc_header.livetime.numberOfSeconds;
    liveTimeMCErr = mc_header.livetime.errorOfSeconds;
    if (liveTimeMC > 0) {
      isMupage = true;
      NOTICE("[R] mupage live time is (" << liveTimeMC << " +/- " << liveTimeMCErr << ") s." << endl);
    } else {
      NOTICE("[R] aanet header found but mupage live time is zero. Likely not a mupage file.");
    } 
  }
 
  // retrieving and exporting run number to tree metadata;
 
  int runNumber = 0;
 
  if (summaryFile.hasNext()) { 
    runNumber = summaryFile.getEntry(0)->getRunNumber();
    NOTICE("[R] Processing run #" << runNumber << endl);
    summaryFile.rewind();
  }
  
  //-----------------------------------------------------------
  // file processing: loop over timeslices
  //-----------------------------------------------------------
  
  for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) {
 
    STATUS("Entry: " << setw(10) << counter  << '\r');
 
    const JDAQTimeslice*    timeslice = in.next();
    const Long64_t          index     = summaryFile.find(*timeslice);
    JDAQSummaryslice*       summary   = (index != -1 ? summaryFile.getEntry(index) : NULL);
    
    summaryRouter.update(summary);
 
    // check correspondence between summary data and timeslice data
    if (summary != NULL) {
      if (timeslice->getFrameIndex() != summary->getFrameIndex()) {
        ++treeMismatchCount;
        ERROR("[R>T] Frame indices do not match [counter=" << counter << ", timeslice=" << timeslice->getFrameIndex() << ", summaryslice=" << summary->getFrameIndex() << "]" << endl);
        if (treeMismatchCount > 1) {
          FATAL("ROOT TTrees not aligned. Run #" << runNumber << endl);
        }
        continue;
      }
    }
 
    //-----------------------------------------------------------
    // timeslice processing: loop over superframes
    //-----------------------------------------------------------
    for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) {
      
      // empty superframes are not skipped because SN timeslices could have empty superframes even if the module was active;
      // if DAQ error or FIFO full -> frame counts as dead time;
      // if frame is missing but no errors -> frame counts as live time (SN stream has many empty frames);
  
      int moduleID = super_frame->getModuleID();
 
      if (!router.hasModule(moduleID)) {
        // module is not mapped in detector file, nothing can be done
        continue;
      }
        
      JModuleAddressMap moduleAddressMap = detectorAddressMap.get(moduleID);
 
      //----------------------------------------------------------
      // retrieval of summary data from superframe or summaryframe
      //----------------------------------------------------------
      
      vector<double> rate_Hz(NUMBER_OF_PMTS, 0.0);
      JDAQFrameStatus status;
      
      if (summary == NULL) { 
        WARNING("[R>T>F] Missing summary for timeslice." << endl);
        // PMT rates estimated from timeslice data
        for (JDAQSuperFrame::const_iterator i = super_frame->begin(); i != super_frame->end(); ++i) {
          rate_Hz[i->getPMT()] += 1e9 / getFrameTime();
        }
        
        // frame status retrieved from superframe
        status = super_frame->getDAQFrameStatus();
        
      } else {
 
        // summary frame lookup
 
        JDAQSummaryFrame summary_frame;
 
        if (!summaryRouter.hasSummaryFrame(super_frame->getModuleIdentifier())) {
          summaryRouter.print(cout); 
          FATAL("[R>D>F] Summary frame not found, module = " << super_frame->getModuleID() << endl);
        } else {
          summary_frame = summaryRouter.getSummaryFrame(super_frame->getModuleID());
        }
        
        // PMT rates recovered from summaryslice data
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
          rate_Hz[i] = summary_frame.getRate(i);
        }
        
        // frame status retrieved from summaryframe
        status = summary_frame.getDAQFrameStatus();
        
      }
      
      //---------------------------
      // processing of summary data
      //---------------------------      
      
      // DAQ error (e.g. missing UDP trailer) -> discard frame
      if (!status.testDAQStatus()) {
          WARNING("[R>D>F] DAQ status not okay in timeslice #" << timeslice->getFrameIndex() << ", DOM=" << super_frame->getModuleID() << endl);
          continue;
      }
     
      vector<bool> veto(NUMBER_OF_PMTS, false);
      bool moduleLiveTime = 0;
      
      // default status of module is alive if at least one PMT is active
      for (vector<double>::const_iterator r = rate_Hz.begin(); !moduleLiveTime && r != rate_Hz.end(); ++r) {
        if ((*r) > 0) {
          moduleLiveTime = 1;
        } 
      }
     
      int frame_HRV_count = status.countHighRateVeto();
      int frame_FAF_count = status.countFIFOStatus();
      int frame_URV_count = 0; // user rate veto
 
      for (int i = 0; i != NUMBER_OF_PMTS; ++i) {
        if ((veto[i] = (veto[i] || !rateVeto_Hz(rate_Hz[i])))) {
          frame_URV_count++;
        }
      }
 
      count_HRV += frame_HRV_count;
      count_FAF += frame_FAF_count;
      count_URV += frame_URV_count;
 
      int frame_VETO_count = frame_HRV_count + frame_FAF_count + frame_URV_count;
 
      
      if (filterLevel >= 2 && frame_VETO_count > 0) {
        moduleLiveTime = 0;
        fill(veto.begin(), veto.end(), true);
      } else if (filterLevel >= 1 && frame_HRV_count + frame_FAF_count > 0) {
        for (int i = 0; i != NUMBER_OF_PMTS; ++i) { 
          veto[i] = (veto[i] || (status.testFIFOStatus(i) || status.testHighRateVeto(i)));
        }
      }
 
      // recover DOM in case the dead channel corresponds to the muted channel
      
      if (muteChannel != -1) {
 
        // veto conditions are not exclusive
        int mute_VETO_count = 
          status.testHighRateVeto(muteChannel) +
          status.testFIFOStatus(muteChannel  ) +
          !rateVeto_Hz(rate_Hz[muteChannel]  );
        
        if (mute_VETO_count == frame_VETO_count) {
          moduleLiveTime = 1;
          fill(veto.begin(), veto.end(), false);
        }
        
        veto[muteChannel] = true;
        
      }   
      
      // get address and module for current superframe DOM;
      const JModuleAddress& address = router.getAddress(super_frame->getModuleID());
      const JModule&        module  = detector.getModule(address);
      const TString         moduleLabel   = getLabel(module);
 
      // acknowledge DOM is live in this timeslice for analysis purpose;
      if (moduleLiveTime) {
        h_livetime->Fill(moduleLabel, getFrameTime() * 1e-9 * moduleLiveTime);
      }  
      
      //--------------------------------------------------------------------
      // superframe processing: hit joining, calibration, filter, processing
      //--------------------------------------------------------------------
      
      JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, router.getModule(super_frame->getModuleID())); 
 
      if (!notJoin) {
        for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
          i->join(match); 
        }
      }
      
      JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
 
      filteredData.clear();
      
      for (vector<JHitR0>::const_iterator h = data.begin(); h != data.end(); ++h) {
        
        const int pmt = h->getPMT();
 
        // veto and calibration;
        if (!veto[pmt]                && 
            rateVeto_Hz(rate_Hz[pmt]) &&
            totVeto_ns(h->getToT()) )     {
          filteredData.push_back(*h);
        } 
      }
      
      // processing signal;
      
      if (filteredData.size() > 1) {
 
        TH1D* h_hit =  H[moduleLabel];
        TH1D* h_pmt =  P[moduleLabel];
        
        TH2D* h2_hit = HT[moduleLabel];
        TH2D* h2_co  = CO[moduleLabel];
 
        sort(filteredData.begin(), filteredData.end());
      
        //-----------------------------------------------------------
        // superframe processing: loop on calibrated hits
        //-----------------------------------------------------------
        for (vector<JHitR0>::const_iterator p = filteredData.begin(); p != filteredData.end(); ) {
 
          vector<JHitR0>::const_iterator q = p;
        
          // coincidence building using a fixed time window
          while (++q != filteredData.end() && q->getT() - p->getT() <= TMax_ns ) {}
        
          // multiplicity estimation
          int hit_multiplicity = distance(p,q);
          set<int> pmthit;
          for (vector<JHitR0>::const_iterator __p = p ; __p != q; ++__p ) { pmthit.insert(__p->getPMT()) ; }
          int pmt_multiplicity = pmthit.size() ;
        
          // histogram filling
          h_pmt->Fill(pmt_multiplicity);
          h_hit->Fill(hit_multiplicity);
 
          if (twoDim && hit_multiplicity > 1) { 
            const double W = factorial(hit_multiplicity, 2);
            for (vector<JHitR0>::const_iterator __p = p; __p != q; ++__p) {
              for (vector<JHitR0>::const_iterator __q = __p; ++__q != q; ) {
                double dt = JCombinatorics::getSign(__p->getPMT(), __q->getPMT()) * (__q->getT() - __p->getT());
                h2_hit->Fill(hit_multiplicity, dt, 1.0/W);
              }
            }
          }
 
          if (monitorOccupancy) {
            for (set<int>::const_iterator r = pmthit.begin(); r != pmthit.end(); ++r) {
              JPMTPhysicalAddress pmtAddress = moduleAddressMap.getAddressTranslator(*r);
              TString pmtLabel(pmtAddress.toString());
              h2_co->Fill(pmt_multiplicity, pmtLabel, 1.0/pmt_multiplicity);
            }
          }
          
          // benchmarking
          if (hit_multiplicity >= sn_mul_min) {
            time_distr->Fill(p->getT());
          }
 
          p = q; 
 
        }
 
      }      
 
    }
    
  }
    
  NOTICE("[R] " << counter << " timeslices processed." << endl);
  
  if (isMupage) {
    h_livetime->Scale(liveTimeMC / (getFrameTime() * 1e-9 * counter * prescale));
  }
 
  // nominal live time;
  double nominalLiveTime = (isMupage ? liveTimeMC : counter * getFrameTime() * 1e-9) / prescale;
  h_livetime->Fill("Nominal", nominalLiveTime);
  
  out.cd();
  
  int nLiveDOMs = H.size();
  
  h_livetime->Write();
  
  // writing hit multiplicities is redundant given new options
  // H.Write(out);
 
  P.Write(out);
  
  if (twoDim) {
    HT.Write(out);
  }
 
  if (monitorOccupancy) {
    CO.Write(out);
  }
 
  time_distr->Write();
 
  double rate_HRV = count_HRV / (1.0 * counter);
  double rate_FAF = count_FAF / (1.0 * counter);
  double rate_URV = count_URV / (1.0 * counter);
  double biolumIndex = (rate_HRV + rate_FAF) / (nLiveDOMs * NUMBER_OF_PMTS);
  NOTICE("[R] Average HRV count per timeslice: " << rate_HRV << endl);
  NOTICE("[R] Average FIFO-almost-full count per timeslice: " << rate_FAF << endl);
  NOTICE("[R] Average user rate veto count per timeslice: " << rate_URV << endl);
  NOTICE("[R] " << H.size() << " DOMs were active in the run." << endl);
  NOTICE("[R] Bioluminescence proxy index: " << biolumIndex << endl);
 
  putObject(out, JMeta(argc,argv));
  
  out.Close();
 
  NOTICE(endl);
  
  return (counter ? 0 : 1);
 
}