JRipple.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2F.h"
#include "TF1.h"
 
 
 
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JSummaryRouter.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JTimesliceRouter.hh"
#include "JTools/JCombinatorics.hh"
 
#include "JGizmo/JManager.hh"
#include "JLang/JObjectMultiplexer.hh"
#include "JMath/JMathToolkit.hh"
 
#include "JROOT/JROOTClassSelector.hh"
 
#include "JSupernova/JSupernova.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JSupportToolkit.hh"
 
 
#include "JRipple.hh"
 
 
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;
using namespace JSUPERNOVA;
 
/**
 * \author mlincett
 * \file
 * Example application to examine hit rates (and active channels) as a function of time.
 */
int main(int argc, char **argv)
{
 
  JMultipleFileScanner<> inputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 outputFile;
  string                 detectorFile;
  double                 TMax_ns;
  double                 TMaxSelf_ns;
  JROOTClassSelector     selector;
  int                    debug;
  int                    binWidth_ms;
  bool                   backVeto;
  bool                   globalOutputOnly;
  bool                   mode2D;
  JRange<int>            multiplicityRange;
  JRange<double>         totRange_ns;
  // Check input parameters
 
  try { 
 
    JParser<> zap("Example program to examine rates as a function of time on ms-level timescales.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "ripple.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['C'] = make_field(selector)            = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['T'] = make_field(TMax_ns,       "Time window for local coincidences (if 0 run in L0 mode)")  =    0.0;
    zap['S'] = make_field(TMaxSelf_ns,   "Time window for merge of double pulses")                            =   10.0;
    zap['B'] = make_field(binWidth_ms,   "Bin width (experimental)")                                     =      1;
    zap['V'] = make_field(backVeto,      "Apply retroactive veto.");
    zap['D'] = make_field(mode2D,        "L1 mode: create 2D histogram with time differences of coincidences (heavy memory usage, ignored if TMax_ns = 0)");
    zap['O'] = make_field(globalOutputOnly,  "Write only aggregate histograms");
    zap['M'] = make_field(multiplicityRange, "L1 mode: multiplicity range (ignored if TMax_ns = 0)")   = JRange<int>(2,31);
    zap['t'] = make_field(totRange_ns, "ToT acceptance range") = JRange<double>(0.0, 255.0);
    zap['d'] = make_field(debug)               = 1;
 
            
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  cout.tie(&cerr);
 
  if (((int)(getFrameTime() / 1e6)) % binWidth_ms) {
    FATAL("Frame time must be an integer multiple of bin width");
  }
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  // Configure input streams
 
  JTreeScanner<>::debug = debug;
 
  JAutoTreeScanner<JDAQTimeslice, JDAQEvaluator> zmap = JType<JDAQTimesliceTypes_t>();
 
  JTreeScannerInterface<JDAQTimeslice>* pts = zmap[selector]; 
 
  pts->configure(inputFile);
 
  int fEnd   = pts->rbegin()->getFrameIndex();
  int fStart = pts->begin( )->getFrameIndex();
 
  if (fEnd > inputFile.getUpperLimit()) {
    fEnd = fStart + inputFile.getUpperLimit();
  }
 
  double tEnd_ms   = fEnd   * getFrameTime() / 1.0e6;
  double tStart_ms = (fStart - 1) * getFrameTime() / 1.0e6;
 
  int     runNumber   = pts->begin()->getRunNumber();
 
  TString runTag      = Form("%d" , runNumber);
 
  double tRun_ms = tEnd_ms - tStart_ms; 
 
  NOTICE("START/END/DURATION [s]: " << tStart_ms / 1000 << " " << tEnd_ms / 1000 << " " << tRun_ms / 1000 << endl);
 
  // Configure data structure and routers
 
  typedef JSuperFrame2D<JHitR0>   JSuperFrame2D_t;
 
  typedef JSuperFrame1D<JHitR0>   JSuperFrame1D_t;
 
  const JMatchL0<JHitR0>      match(TMaxSelf_ns);     // time window self-coincidences [ns]
 
  const JModuleRouter moduleRouter(detector);
  
  const JDAQHitRouter hitRouter(detector);
    
  //JTimesliceRouter timesliceRouter;
  JSummaryRouter summaryRouter;
 
  // Configure output histograms
 
  typedef JManager<string, TH1F>  JManager_t;
 
  int nx = (tEnd_ms - tStart_ms) / binWidth_ms;
 
  int x_inf = tStart_ms;
 
  JManager_t RT_DOM(new TH1F("RT_%", NULL, nx      , x_inf, tEnd_ms));  
  JManager_t NC_DOM(new TH1F("NC_%", NULL, nx / 100, x_inf, tEnd_ms));  
  
  TString rt_tag = runTag + TString(".RT_DET_%");
  TString nc_tag = runTag + TString(".NC_DET_%"); 
 
  JManager_t RT_DET(new TH1F(rt_tag, NULL, nx,       x_inf, tEnd_ms));  
  JManager_t NC_DET(new TH1F(nc_tag, NULL, nx / 100, x_inf, tEnd_ms));  
 
  h2d_t* hT = NULL;
   
 
  if (mode2D) {     
 
    const int ny = 50;
 
    const int max_size = floor(h2d_limit / (sizeof(h2d_bintype) * ny));
 
    if (nx >= max_size) {
      FATAL("2D histogram size not supported by ROOT file output; limit input size (-n) below " << floor(max_size / 100.0) << endl);
    } 
 
    TString rt2d_tag = runTag + TString(".RT2D_DET"); 
 
    hT = new h2d_t(rt2d_tag, NULL, nx, x_inf, tEnd_ms, ny, -TMax_ns, +TMax_ns);
    hT->SetDirectory(0);
 
  }
 
  //---------------------
  // Timeslice processing
  //---------------------
 
  counter_type counter = 0;
  
  // preload first timeslice
 
  JDAQTimeslice curr;
 
  if (pts->hasNext()) {
    curr = *(pts->next());
    ++counter;
  } else {
    FATAL("Input file is too short.");
  }
 
  // loop over timeslices
 
  for ( ; pts->hasNext() && counter != inputFile.getLimit(); ++counter) {
 
    STATUS("timeslice: " << setw(10) << counter << '\r'); DEBUG(endl);
 
    JDAQTimeslice next = *(pts->next());
 
    const int ic = curr.getFrameIndex();
 
    const int in = next.getFrameIndex();
    
    const double tTimeslice_ms = getTimeOfRTS(ic) / 1.0e6;
 
    // look up next summaryslice
    
    JDAQSummaryslice* nextSummary = NULL;
 
    if (backVeto && ((in - ic) == 1)) { 
      nextSummary = new JDAQSummaryslice(next);
      summaryRouter.update(nextSummary);
    }
 
    // loop over frames
 
    for (JDAQTimeslice::const_iterator frame = curr.begin(); frame != curr.end(); ++frame) {
 
      const int      DOMID       = frame->getModuleID();
      const JModule& module      = moduleRouter.getModule(DOMID);
      const string   moduleLabel = getModuleLabel(module); 
     
      TH1F* RD = RT_DOM[moduleLabel];
 
      vector<bool> veto(NUMBER_OF_PMTS, false);
 
      // current veto
 
      if (frame->testHighRateVeto() || frame->testFIFOStatus()) {
        for (int pmt = 0; pmt < NUMBER_OF_PMTS; pmt++) {
          veto[pmt] = ( frame->testHighRateVeto(pmt) || frame->testFIFOStatus(pmt) );
        }        
      }  
 
     
      // retroactive veto
 
      if (nextSummary != NULL) {
 
        JDAQSummaryFrame nextFrame = summaryRouter.getSummaryFrame(DOMID);
 
        if (nextFrame.testHighRateVeto() || nextFrame.testFIFOStatus()) {
          for (int pmt = 0; pmt < NUMBER_OF_PMTS; pmt++) {
            veto[pmt] = veto[pmt] || ( nextFrame.testHighRateVeto(pmt) || nextFrame.testFIFOStatus(pmt) );
          }      
        }  
     
      }
 
      
      // determine number of active channels
 
      int nActiveChannels = count(veto.begin(), veto.end(), false);
      
      NC_DOM[moduleLabel]->Fill(tTimeslice_ms, nActiveChannels);
 
      
      // merge double pulses
 
      JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, module);
 
      for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) {
        if (veto[i->getPMTAddress()]) {
          i->reset();
        } else {
          i->join(match);
        } 
      }
 
      JSuperFrame1D_t& buffer1D = JSuperFrame1D_t::multiplex(buffer);
      buffer1D.pop_back();
      
      // process hit data
 
      vector<JHitR0> data;
 
      for (vector<JHitR0>::const_iterator h = buffer1D.begin(); h != buffer1D.end(); h++) {
        if (totRange_ns(h->getToT())) {
          data.push_back(*h);
        }
      }
      
      if (data.size() > 1) {
 
        // data need to be sorted for L1 processing
 
        if (TMax_ns > 0) {
          sort(data.begin(), data.end());
        }
 
 
        // loop over hits in frame
 
        for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); p++) {
          
          const double tHit = p->getT();
          const double tHit_ms = tTimeslice_ms + (tHit / 1.0e6);
 
          if (TMax_ns == 0.0) {
 
            // L0 mode
 
            RD->Fill(tHit_ms);
 
          } else {
 
            // L1 mode
            
            vector<JHitR0>::const_iterator q = p;
            
            while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
 
            int M = distance(p, q);
 
            if (multiplicityRange(M)) {
 
              RD->Fill(tHit_ms);
            
              if (mode2D) {
              
                // L1 2D mode
 
                const double W = factorial(M, 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());
 
                    hT->Fill(tHit_ms, dt, 1.0 / W);
        
                  }
                }
              }
            
            }
 
            p = (q - 1);
 
          } 
 
        }
      }
    }
 
    // update pointers
 
    curr = next;
 
    if (nextSummary != NULL) { delete nextSummary; }
 
  }
 
 
  //--------------------------------------
  // Histogram processing
  //--------------------------------------
 
  NOTICE("Processing histograms." << endl);
 
  for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
    
    string moduleLabel = getModuleLabel(*module);
 
    if (RT_DOM.count(moduleLabel) && NC_DOM.count(moduleLabel)) {
 
      TH1F* rt = RT_DOM.at(moduleLabel);
      TH1F* nc = NC_DOM.at(moduleLabel);
 
      for (int b = 1; b <= rt->GetXaxis()->GetNbins(); b++) {
 
        double r = rt->GetBinContent(b);
        double t = rt->GetBinCenter( b);
 
        RT_DET["SUM"]->Fill(t, r);
 
      }
 
 
      for (int b = 1; b <= nc->GetXaxis()->GetNbins(); b++) {
 
        double n = nc->GetBinContent(b);
        double t = nc->GetBinCenter(b );
 
        NC_DET["SUM"]->Fill(t, n);
 
      }
 
    } else {
 
      DEBUG(moduleLabel << " not active." << endl);
 
    }
 
  }
 
  NOTICE("Writing output file" << endl);
 
  if (outputFile != "") {
 
    TFile out(outputFile.c_str(), "RECREATE");
 
    NOTICE("Writing 1D histograms" << endl);
 
    RT_DET.Write(out);
    NC_DET.Write(out);
 
    NOTICE("Writing 2D histogram" << endl);
 
    if (hT != NULL) {
      hT->Write();
    }
 
    if (!globalOutputOnly) {
 
      TString dir_tag = runTag + TString(".Modules");
 
      NOTICE("Writing individual modules histograms" << endl);
 
      TDirectory* dir = out.mkdir(dir_tag);
      RT_DOM.Write(*dir);
      NC_DOM.Write(*dir);
    }
 
    NOTICE("Closing file" << endl);
 
    out.Close();
  }
 
}