JKexing2D.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 "TH2D.h"
#include "TF1.h"
#include "TParameter.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "JDAQ/JDAQEvaluator.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JROOT/JManager.hh"
 
#include "JROOT/JROOTClassSelector.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "JLang/JObjectMultiplexer.hh"
#include "JMath/JMathToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JTools/JCombinatorics.hh"
#include "JTools/JQuantile.hh"
 
#include "km3net-dataformat/online/JDAQHit.hh"
 
#include "JTrigger/JHit.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JDAQHitToTSelector.hh"
 
#include "JSupernova/JSupernova.hh"
#include "JSupernova/JKexing2D.hh"
 
 
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;
using namespace JSUPERNOVA;
using namespace JROOT;
 
/**
 * \author mlincett
 * \file
 * Example application to analyse intra-run supernova background
 */
int main(int argc, char **argv)
{
 
  typedef JRange<int> JRange_t;
  typedef JRange<JDAQHit::JTOT_t>    JToTRange_t;
 
  JMultipleFileScanner<> inputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 outputFile;
  string                 detectorFile;
  double                 TMax_ns;
  double                 TVeto_ns;
  JToTRange_t            totRange_ns;
  JDAQHitToTSelector     totSelector_ns;
  JROOTClassSelector     selector;
  bool                   timeDifferences;
  int                    debug;
  int                    preTriggerThreshold;
 
  // Check input parameters
 
  try { 
 
    JParser<> zap("Example application to study supernova detection background");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "kexing2D.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['T'] = make_field(TMax_ns)             = 10.0;
    zap['V'] = make_field(TVeto_ns)            = 1000.0;  
    zap['D'] = make_field(timeDifferences);
    zap['C'] = make_field(selector)            = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['d'] = make_field(debug)               = 1;
    zap['P'] = make_field(preTriggerThreshold) = 4;
    zap['t'] = make_field(totSelector_ns) = JDAQHitToTSelector(4, 255);
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  // -----------------------------------
  // STEP 0: prepare
  // -----------------------------------
 
  // configure input streams
 
  JTreeScanner<>::debug = debug;
 
  JAutoTreeScanner<JDAQTimeslice, JDAQEvaluator> zmap = JType<JDAQTimesliceTypes_t>();
 
  JTreeScannerInterface<JDAQTimeslice>* pts = zmap[selector]; 
 
  pts->configure(inputFile);
 
  // detect input stream size
 
  
  int    fEnd   = pts->rbegin()->getFrameIndex();
  int    fStart = pts->begin( )->getFrameIndex(); 
 
 
  // crop histogram if -n is specified
  if (fEnd > inputFile.getUpperLimit()) {
    fEnd = fStart + inputFile.getUpperLimit();
  }
 
  int fLength = 1 + fEnd - fStart;
 
  NOTICE("begin | end | length = " << fStart << " | " << fEnd << " | " << fLength << endl);
 
 
  if (fStart > fEnd) {
    FATAL("FATAL: inconsistent TTree indexing" << endl);
  }  
 
  // detector routers
 
  const JModuleRouter moduleRouter(detector);
  
  const JDAQHitRouter hitRouter(detector);
 
  // data structures
 
  const int    nx   = fLength;
  const double xmin = fStart;
  const double xmax = fEnd + 1;
 
  const int    ny   = NUMBER_OF_PMTS + 1; 
  const double ymin = -0.5;
  const double ymax = ymin + ny;
 
  typedef JManager <string, TH2D>  JManager2D_t;
  typedef JManager <string, TH1D>  JManager1D_t;
 
  JManager2D_t MT(new TH2D(mul_p   , NULL, nx, xmin, xmax, ny, ymin, ymax));
  JManager1D_t ST(new TH1D(status_p, NULL, nx, xmin, xmax));
 
  // -----------------------------------
  // STEP 1: building vetoes from events
  // -----------------------------------
 
  int runNumber = 0;
 
  TH1D* h_vtr = new TH1D("VetoTimeRange","VetoTimeRange", 10000, 0, 10000);
 
  JTreeScanner<JDAQEvent, JDAQEvaluator> evIn(inputFile);
 
  map<int, JVetoSet> triggeredEvents; 
 
  for (; evIn.hasNext(); ) {
 
    STATUS("event: " << setw(10) << evIn.getCounter() << '\r'); DEBUG(endl);
    
    JDAQEvent* event = evIn.next();
 
    if (!runNumber) { runNumber = event->getRunNumber(); }
 
    JVeto vp(*event, hitRouter);
 
    triggeredEvents[event->getFrameIndex()].push_back(vp);
 
    h_vtr->Fill(vp.getLength());
 
  } 
 
  STATUS(triggeredEvents.size() << " JDAQEvents loaded in veto buffer." << endl);
 
  TParameter<int> RUNNR("RUNNR", runNumber);
  // TParameter<int> TSTART("TSTART", tStart)
 
  //-----------------------------
  // STEP 2: timeslice processing
  // ----------------------------
 
  counter_type counter = 0;
 
  for ( ; pts->hasNext() && counter != inputFile.getLimit(); ++counter) {
 
    STATUS("timeslice: " << setw(10) << counter << '\r'); DEBUG(endl);
 
    const JDAQTimeslice* timeslice = pts->next();
 
    int fIndex = timeslice->getFrameIndex();
 
    if (counter == 0) {
      NOTICE("Start frame index  = " << fIndex << " @ T " << timeslice->getTimesliceStart() << endl);
    }
        
    JVetoSet veto;
    
    if (triggeredEvents.count(fIndex)) { 
      veto = triggeredEvents.at(fIndex);
    }
 
    // L0-L1 processing
 
    typedef JHitR0 hit_type;
      
    typedef JSuperFrame2D<hit_type> JSuperFrame2D_t;
 
    typedef JSuperFrame1D<hit_type> JSuperFrame1D_t;
 
    const JMatchL0<hit_type> match(TMax_ns);
 
    double fractionActivePMTs = 0;
    
    int nDOMs = timeslice->size();
 
    for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
 
        int moduleID = frame->getModuleID();
 
        fractionActivePMTs += ((double) frame->countActiveChannels());
 
        JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame,
                                                               moduleRouter.getModule(moduleID));
        //                                                             totSelector_ns);
 
        buffer.preprocess(JPreprocessor::join_t, match);
 
        JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
 
        // sort(data.begin(), data.end()); // JSuperFrame1D is sorted
 
        // histogram for time differences, may not be used
 
        TH1D* h2dt = new TH1D("H2DT", "H2DT", 100, -TMax_ns, +TMax_ns); 
        
        for (vector<JHitR0>::const_iterator p = data.begin(); p != data.end(); ) {
 
          vector<JHitR0>::const_iterator q = p;
        
          while (++q != data.end() && q->getT() - p->getT() <= TMax_ns) {}
            
          int m = distance(p, q);
 
          // raw multiplicity count //
 
          MT["RAW"]->Fill(fIndex, m);
 
          // time difference distribution
            
          if (selector != "JDAQTimesliceSN" && timeDifferences && m > 1) { 
            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());
                h2dt->Fill(dt, 1.0/W);
              }
            }
          }
          
          // slide iterator
          
          p = q;
            
        } // end hit processing
        
        if (h2dt->GetEntries() > 0) {
          
          TF1 f("f", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
 
          f.SetParameter(0, h2dt->GetMaximum());
          f.SetParameter(1, h2dt->GetMean());
          f.SetParameter(2, h2dt->GetRMS() * 0.25);
          f.SetParameter(3, h2dt->GetMinimum());
          
          h2dt->Fit(&f, "Q", "same");
 
          double nb = h2dt->GetNbinsX();
          double bg_v = f.GetParameter(3) * nb;
          double sg = h2dt->GetSumOfWeights() - bg_v;
          
          // inclusive 2-fold after subtraction of fitted background //
          
          MT["FIT"]->Fill(fIndex, 2, sg);
          
        }
 
        delete h2dt;
 
    }
 
    // end timeslice processing
 
    fractionActivePMTs /= (NUMBER_OF_PMTS * nDOMs);
 
    ST["PMT"]->Fill(fIndex, fractionActivePMTs);
 
    // STANDARD PROCESSING
 
    JDataSN preTrigger(TMax_ns, preTriggerThreshold);// totSelector_ns); 
 
    preTrigger(timeslice, moduleRouter, totSelector_ns);
 
    JTriggerSN trigger(TVeto_ns);
 
    trigger(preTrigger);
 
    // generate and configure event-based veto
    
    // count trigger
 
    JRange_t A = JRange_t(4,31);
 
    // select single-module clusters with max multiplicity in range A
    JSNFilterM  trgA1(A, 1);
 
    // select non-vetoed clusters with max multiplicity in range A
    JSNFilterMV trgAV(A, veto);
    
    vector<double> m_a1 = trigger.getMultiplicities(trgA1);
 
    for (vector<double>::const_iterator p = m_a1.begin(); p != m_a1.end(); p++) {
      MT["TA1"]->Fill(fIndex, *p);
    }
 
    vector<double> m_av = trigger.getMultiplicities(trgAV);
 
    for (vector<double>::const_iterator p = m_av.begin(); p != m_av.end(); p++) {
      MT["TAV"]->Fill(fIndex, *p);
    }
    
  }
 
  //-------------------------------------
  // STEP 3: generating trigger summaries
  // ------------------------------------
 
  if (outputFile != "") {
 
    TFile out(outputFile.c_str(), "RECREATE");
 
    putObject(out, JMeta(argc,argv));
 
    RUNNR.Write();
    // TSTART.Write();
 
    h_vtr->Write();
        
    MT.Write(out);
    ST.Write(out);
 
    out.Close();
 
  }
 
}