JKaboom.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 "TTree.h"
 
#include "JDAQ/JDAQEvaluator.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JROOT/JROOTClassSelector.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JAutoTreeScanner.hh"
#include "JSupport/JSupport.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JTrigger/JDAQHitToTSelector.hh"
 
#include "JSupernova.hh"
 
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;
using namespace JSUPERNOVA;
using namespace JROOT;
 
/**
 * \authors selhedri,iagoos
 * \file
 *  Application to analyse buffered L1 timeslices for the online CCSN analysis
 */
 
/**
 * Initialize the output ROOT file with CCSN candidate information
 * \param obs               Single-DOM observables
 * \param frame_number      Timeslice frame number
 * \param tr                Output ROOT tree
 */
inline int initialize_root_tree(JCoincidenceSN &obs, int *frame_number, TTree *tr){
    tr->Branch("dom_id", &obs.moduleID, "dom_id/I");
    tr->Branch("vetoIndex", &obs.vetoIndex, "vetoIndex/I");
    tr->Branch("event_time", &obs.time, "event_time/D");
    tr->Branch("timeslice_time_s", &obs.timeslice_time, "timeslice_time_s/D");
    tr->Branch("frame_number", frame_number, "frame_number/I");
    tr->Branch("multiplicity", &(obs.multiplicity), "multiplicity/I");
    return 0;
}
 
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;
  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)          = "kaboom.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['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);
  }
 
  cout.tie(&cerr);
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  TTree *tr = new TTree("singledom", "SN coincidence info");
  JCoincidenceSN SNdata(0,0,0);
  int frame_number;
  initialize_root_tree(SNdata, &frame_number, tr);
 
  // -----------------------------------
  // 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);
 
  // detector routers
 
  const JModuleRouter moduleRouter(detector);
 
  const JDAQHitRouter hitRouter(detector);
 
 
  // -----------------------------------
  // 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);
 
  //-----------------------------
  // STEP 2: timeslice processing
  // ----------------------------
 
  counter_type counter = 0;
 
  TFile *out = new TFile(outputFile.c_str(), "RECREATE");
  tr->SetDirectory(out);
 
  for ( ; pts->hasNext() && counter != inputFile.getLimit(); ++counter) {
 
    STATUS("timeslice: " << setw(10) << counter << '\r'); DEBUG(endl);
 
    const JDAQTimeslice* timeslice = pts->next();
 
    frame_number = timeslice->getFrameIndex();
 
    if (counter == 0) {
      NOTICE("Start frame index  = " << frame_number << " @ T " << timeslice->getTimesliceStart() << endl);
    }
 
    JVetoSet veto;
 
    if (triggeredEvents.count(frame_number)) {
      veto = triggeredEvents.at(frame_number);
    }
 
    // 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(preTriggerThreshold, NUMBER_OF_PMTS);
 
    // 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);
 
    JDAQUTCExtended tstart = timeslice->getTimesliceStart();
    double event_time = tstart.getTimeNanoSecond() * 1e-9;
 
    for(auto &p: trigger){
          SNdata = p.getPeak();
          SNdata.vetoIndex = trgA1(p) ? (trgAV(p) ? 0 : 2) : (trgAV(p) ? 1 : 3); // 0 if no veto, 1 if only DOM correlation veto, 2 if only muon veto, 3 if both
          SNdata.timeslice_time = event_time;
          tr->Fill();
    }
  }
 
  //-------------------------------------
  // STEP 3: Writing data
  // ------------------------------------
 
    tr->Write();
    out->Close();
}