JKexing.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 "TF1.h"
 
#include "JSupernova.hh"
 
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JDAQ/JDAQEvaluator.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JROOT/JManager.hh"
 
#include "JROOT/JROOTClassSelector.hh"
#include "JTools/JQuantile.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 "JTrigger/JDAQHitToTSelector.hh"
 
using namespace std;
using namespace JPP;
using namespace KM3NETDAQ;
using namespace JSUPERNOVA;
 
/**
 * \author mlincett
 * \file
 * Example application to test supernova trigger on run files.
 */
int main(int argc, char **argv)
{
 
  typedef JRange<int>    JRange_t;
 
  JMultipleFileScanner<> inputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 outputFile;
  string                 detectorFile;
  double                 TMax_ns;
  double                 TVeto_ns;
  JROOTClassSelector     selector;
  JDAQHitToTSelector     totSelector_ns;
  int                    debug;
  int                    preTriggerThreshold;
  JRange_t               M;
  
 
  // Check input parameters
 
  try { 
 
    JParser<> zap("Example program test supernova triggers.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "kexing.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['M'] = make_field(M)                  = JRange<int>(6,10);
    zap['t'] = 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);
  }
 
  // Configure input streams
 
  JTreeScanner<>::debug = debug;
 
  JAutoTreeScanner<JDAQTimeslice, JDAQEvaluator> zmap = JType<JDAQTimesliceTypes_t>();
 
  JTreeScannerInterface<JDAQTimeslice>* pts = zmap[selector]; 
 
  pts->configure(inputFile);
  
  // Configure routers
 
  const JModuleRouter moduleRouter(detector);
  
  const JDAQHitRouter hitRouter(detector);
  
  // -----------------------------------
  // STEP 1: building vetoes from events
  // -----------------------------------
 
  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();
 
    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
  // -------------------------------
 
  // 0 -> no filter
  // 1 -> count once per track
  // 2 -> suppress track
  // 3 -> suppress based on trigger veto
  
  // Output data structures
  
  typedef JManager<int, TH1D>  JManager_t;
 
  JManager_t SNT(new TH1D("SNT_F%", NULL, 100, 0.0, 100));
 
  JManager_t MUL(new TH1D("MUL_F%", NULL, 1 + 31, -0.5, 31 + 1 - 0.5));
 
  const int nStages = 5;
  vector<vector<int> > trgHistory(nStages);
 
  // Loop
 
  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();
 
    JDataSN preTrigger(TMax_ns, preTriggerThreshold); 
 
    preTrigger(timeslice, moduleRouter, totSelector_ns);
 
    JTriggerSN trigger(TVeto_ns);
 
    trigger(preTrigger);
 
    // generate and configure event-based veto
 
    JVetoSet veto;
    
    if (triggeredEvents.count(fIndex)) { 
      veto = triggeredEvents.at(fIndex);
    }
 
    // trigger.setVeto(veto);
    
    // count trigger
 
    JSNFilterM trgF0(M, 0);
    JSNFilterM trgF1(M, 1);
 
    JSNFilterMV trgFV(M, veto);
 
    JRange<int> A = JRange<int>(2,31);
    JSNFilterMV trgAV(A, veto);
 
    int rawCount = count_if(preTrigger.begin(), preTrigger.end(), trgF0);
 
    int trgCountF0 = count_if(trigger.begin(), trigger.end(), trgF0);
    int trgCountF1 = count_if(trigger.begin(), trigger.end(), trgF1);
    int trgCountFV = count_if(trigger.begin(), trigger.end(), trgFV);
    int domCountF1 = trigger.getModules(trgF1).size();
 
    trgHistory[0].push_back(rawCount);
    trgHistory[1].push_back(trgCountF0);
    trgHistory[2].push_back(trgCountF1);
    trgHistory[3].push_back(trgCountFV);
    trgHistory[4].push_back(domCountF1);
 
    trigger.fill(MUL[1], trgF0);
    trigger.fill(MUL[2], trgF1);
    trigger.fill(MUL[3], trgFV);
    trigger.fill(MUL[4], trgAV);
 
  } 
 
  //-------------------------------------
  // STEP 3: generating trigger summaries
  // ------------------------------------
  
  for (int i = 0; i < nStages; i++) {
    for (unsigned j = 0; j < trgHistory[i].size(); j++) {
      SNT[i]->Fill(trgHistory[i][j]);
    }
  }
 
  SNT[0]->SetTitle("M[6,10] count before clustering");
  SNT[1]->SetTitle("M[6,10] count after clustering");
  SNT[2]->SetTitle("M[6,10] count after track self-veto");
  SNT[3]->SetTitle("M[6,10] count after track trigger-veto");
  SNT[4]->SetTitle("M[6,10] count after track self-veto, unique modules");
 
 
  if (outputFile != "") {
 
    TFile out(outputFile.c_str(), "RECREATE");
 
    h_vtr->Write();
    SNT.Write(out);
    MUL.Write(out);
 
    putObject(out, JMeta(argc,argv));
 
    out.Close();
  }
}