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 "JDAQ/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"
#include "JDAQ/JDAQSummaryslice.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/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 JMatch_t      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();
  }

}