JTriggerEfficiency.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
 
#include "TH1D.h"
#include "TRandom3.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Hit.hh"
#include "km3net-dataformat/definitions/pmt_status.hh"
#include "km3net-dataformat/definitions/weightlist.hh" 
 
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQToolkit.hh"
#include "JTimeslice/JEventTimeslice.hh"
#include "JSummaryslice/JSummaryslice.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JPhysics/JK40Rates.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSimulator.hh"
#include "JDetector/JModuleMapper.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JTimeRange.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JK40DefaultSimulator.hh"
#include "JDetector/JPMTDefaultSimulator.hh"
#include "JDetector/JCLBDefaultSimulator.hh"
 
#include "JGeometry3D/JCylinder3D.hh"
 
#include "JTrigger/JHit.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JBuildL1.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JTrigger3DShower.hh"
#include "JTrigger/JTriggerMXShower.hh"
#include "JTrigger/JTrigger3DMuon.hh"
#include "JTrigger/JTriggerBits.hh"
#include "JTrigger/JEventOverlap.hh"
#include "JTrigger/JTimesliceRouter.hh"
#include "JTrigger/JTriggeredEvent.hh"
#include "JTrigger/JTimesliceL1.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JEventToolkit.hh"
#include "JTrigger/JSummaryRouter.hh"
#include "JTrigger/JTriggerToolkit.hh"
#include "JTrigger/JK40RunByRunSimulator.hh"
#include "JTrigger/JPMTRunByRunSimulator.hh"
#include "JTrigger/JCLBRunByRunSimulator.hh"
 
#include "JSirene/JVisibleEnergyToolkit.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "JROOT/JRandom.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
namespace {
 
  using JDETECTOR::JDetector;
 
 
  /**
   * Get time offset of detector.
   *
   * \param  detector   detector
   * \return            time offset [ns]
   */
  inline double getT0(const JDetector& detector)
  {
    using namespace std;
    using namespace JPP;
 
    double T0_ns = numeric_limits<double>::max();
 
    for (const JModule& module : detector) {
      for (const JPMT& pmt : module) {
        if (!pmt.has(PMT_DISABLE)) {
          if (pmt.getT0() < T0_ns) {
            T0_ns = pmt.getT0();
          }
        }
      }
    }
 
    return T0_ns;
  }
}
 
/**
 * \file
 * Auxiliary program to trigger Monte Carlo events.
 *
 * The options 
 *   - <tt>-a</tt> refers to the detector configuration that is used to convert Monte Carlo true information to raw data; and
 *   - <tt>-b</tt> refers to the detector configuration that is used to convert raw data to calibrated data.
 *
 * The event counter of the triggered events, as returned by KM3NETDAQ::JDAQEvent::getCounter(),
 * is set to the corresponding index of the Monte Carlo event in the output.\n
 * This allows for correlating the DAQ event to the Monte Carlo event.\n
 * The time of the DAQ hits can be correlated to the time of the Monte Carlo hits using the frame index 
 * and the time of the Monte Carlo event (Evt::mc_t), respectively (see also time_converter).\n
 * The universal time of the DAQ event is set to the universal time of the Monte Carlo event 
 * with a correction for the time of the event (read hits) within the time slice.\n
 * In run-by-run mode, the trigger parameters as well as the singles rates and status of the PMTs
 * are read from the given raw data file.\n
 * The two-fold (and higher) coincidence rates should still be provided on the command line.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JTYPELIST<JAllTypes_t, JRootTypes_t>::typelist  typelist;
 
  JMultipleFileScanner<Evt> inputFile;
  JFileRecorder<typelist>   outputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 detectorFileA;
  string                 detectorFileB;
  double                 margin_m;
  int                    run;
  JTriggerParameters     parameters;
  bool                   triggeredEventsOnly;
  JPMTParametersMap      pmtParameters;
  JK40Rates              rates_Hz;
  JMixedK40Rates         mixed_Hz;
  string                 runbyrun;
  double                 sigma_ns;
  JRandom                seed;
  int                    debug;
 
  try { 
 
    JParser<> zap("Auxiliary program to trigger Monte Carlo events.");
    
    zap['f'] = make_field(inputFile,           "input file (output of detector simulation)");
    zap['o'] = make_field(outputFile,          "output file")                                                      = "trigger_efficieny.root";
    zap['n'] = make_field(numberOfEvents)                     = JLimit::max();
    zap['a'] = make_field(detectorFileA,       "detector used for conversion from Monte Carlo truth to raw data.");
    zap['b'] = make_field(detectorFileB,       "detector used for conversion of raw data to calibrated data.")     = "";
    zap['m'] = make_field(margin_m,            "margin to the can for the computation of visible energy [m]")      = 20.0;
    zap['R'] = make_field(run,                 "run number")  = -1;
    zap['r'] = make_field(runbyrun,            "DAQ file for run-by-run mode")                                     = JPARSER::initialised();
    zap['@'] = make_field(parameters,          "Trigger parameters (or corresponding file name)")                  = JPARSER::initialised();
    zap['O'] = make_field(triggeredEventsOnly, "optionally write only triggered events.");
    zap['P'] = make_field(pmtParameters,       "PMT simulation data (or corresponding file name)")                 = JPARSER::initialised();
    zap['B'] = make_field(rates_Hz,            "background rates [Hz]")                                            = JPARSER::initialised();
    zap['X'] = make_field(mixed_Hz,            "mixed-L1/L0 background rates [Hz]")                                = JPARSER::initialised();
    zap['s'] = make_field(sigma_ns,            "intrinsic time smearing of K40 coincidences [ns]")                 = JK40DefaultSimulatorInterface::getSigma();
    zap['S'] = make_field(seed,                "seed")        = 0;
    zap['d'] = make_field(debug,               "debug")       = 0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  seed.set(gRandom);
 
  JK40DefaultSimulatorInterface::setSigma(sigma_ns);
 
  setDAQLongprint(debug >= JEEP::debug_t);
 
  if (detectorFileB == "") {
    detectorFileB = detectorFileA;
  }
 
 
  JDetector detectorA;
  JDetector detectorB;
 
  try {
    load(detectorFileA, detectorA);
    load(detectorFileB, detectorB);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  const double TA_ns = getT0(detectorA);
  const double TB_ns = getT0(detectorB);
 
  JPMTParametersMap::Throw(true);
 
  if (!pmtParameters.is_valid()) {
    FATAL("Invalid PMT parameters " << pmtParameters << endl);
  }
 
  if (pmtParameters.getQE() != 1.0) {
 
    WARNING("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
    
    rates_Hz.correct(pmtParameters.getQE());
  }
 
  Head header;
 
  try {
    header = getHeader(inputFile);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  
  const JModuleRouter moduleRouter(detectorB);
  JDetectorSimulator  simbad      (detectorA);
  JSummaryRouter      summaryRouter;
 
  JTreeScanner<JDAQSummaryslice, JDAQEvaluator> scanner;                   // run-by-run
 
  set<int> fs;                                                             // frame indices of summary data already written
 
  JDAQUTCExtended UTC = JDAQUTCExtended::getInstance();                    // start of run
 
  if (runbyrun != "") {
 
    NOTICE("Using run-by-run: " << runbyrun << endl);
 
    scanner.configure(runbyrun);
 
    if (!scanner.empty()) {
      UTC = scanner.begin()->getTimesliceStart();
    } else {
      FATAL("Run-by-run simulation misses summary data." << endl);
    }
 
    try {    
      simbad.reset(new JK40RunByRunSimulator(summaryRouter, rates_Hz, mixed_Hz));
      simbad.reset(new JPMTRunByRunSimulator(summaryRouter, pmtParameters, detectorA));
      simbad.reset(new JCLBRunByRunSimulator(summaryRouter));
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
 
    try {
 
      parameters = getTriggerParameters(runbyrun);
 
      NOTICE("Set trigger parameters from run-by-run input." << endl);
    }
    catch(const JException& error) {
      WARNING("No trigger parameters from run-by-run input;\nrun with default/user input." << endl);
    }
 
    // set live time
 
    JHead buffer(header);
 
    buffer.DAQ.livetime_s = getTimeDuration(getUTCTimeRange(scanner));
    buffer.push(&JHead::DAQ);
 
    copy(buffer, header);
 
  } else {
 
    NOTICE("Using fixed rates [Hz]: " << rates_Hz << endl);
 
    try {    
      simbad.reset(new JK40DefaultSimulator(rates_Hz, mixed_Hz));
      simbad.reset(new JPMTDefaultSimulator(pmtParameters, detectorA));
      simbad.reset(new JCLBDefaultSimulator());
    }
    catch(const JException& error) {
      FATAL(error.what() << endl);
    }
  }
 
  // detector
 
  if (parameters.disableHighRateVeto) {
 
    NOTICE("Disabling high-rate veto of all PMTs." << endl);
 
    detectorB.setPMTStatus(HIGH_RATE_VETO_DISABLE);
  }
 
  parameters.set(getMaximalDistance(detectorB));
 
  DEBUG("Trigger:"        << endl << parameters   << endl);
  DEBUG("PMT parameters:" << endl << pmtParameters << endl); 
  
  const double     Tmax = max(getMaximalTime(detectorA),
                              getMaximalTime(detectorB));
 
  const JTimeRange period(-(Tmax + parameters.TMaxLocal_ns), 
                          +(Tmax + parameters.TMaxLocal_ns));
 
  const JTimeRange frame_time(0.0, getFrameTime());
 
  // hit type should be identical to that in JDataFilter.cc
 
  typedef double        hit_type;
 
  typedef JSuperFrame1D<hit_type>   JSuperFrame1D_t;
  typedef JSuperFrame2D<hit_type>   JSuperFrame2D_t;
  typedef JTimeslice   <hit_type>   JTimeslice_t;
  typedef JBuildL1     <hit_type>   JBuildL1_t;
  typedef JBuildL2     <hit_type>   JBuildL2_t;
 
  const JBuildL1_t buildL1(parameters);
  const JBuildL2_t buildL2(parameters.L2);
  const JBuildL2_t buildSN(parameters.SN);
 
  JTimesliceRouter timesliceRouter(parameters.numberOfBins);
 
  const JTrigger3DMuon   trigger3DMuon  (parameters);
  const JTrigger3DShower trigger3DShower(parameters);
  const JTriggerMXShower triggerMXShower(parameters, detectorB);
 
 
  TH1D h1("Trigger bits", NULL, NUMBER_OF_TRIGGER_BITS, -0.5, NUMBER_OF_TRIGGER_BITS - 0.5);
 
 
  outputFile.open();
 
  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }
 
  outputFile.put(*gRandom);
  outputFile.put(JMeta(argc, argv));
  outputFile.put(header);
  outputFile.put(parameters);
 
  JHead head;
  int   trigger_counter  =  0;
 
  JCylinder3D can(detectorA.begin(), detectorA.end());
  can.addMargin(margin_m);
  
  for (string file_name; inputFile.hasNext(); ) { 
 
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
    
    Evt* event = inputFile.next();
 
    DEBUG(*event << endl);
 
    if (has_primary(*event)) {
      
      if (event->w.size() <= WEIGHTLIST_VISIBLE_ENERGY_HADRON) {
        event->w.resize(WEIGHTLIST_VISIBLE_ENERGY_HADRON + 1, 0.0);
      }
      
      const double Evis = getVisibleEnergy             (*event, can);
      const double Elep = (has_leading_lepton(*event) ? getVisibleEnergyLeadingLepton(*event, can) : 0.0);
      const double Ehad = Evis - Elep;
      
      event->w[WEIGHTLIST_VISIBLE_ENERGY_LEPTON] = Elep;
      event->w[WEIGHTLIST_VISIBLE_ENERGY_HADRON] = Ehad;
    }
 
    if (file_name != inputFile.getFilename()) { 
 
      file_name  =  inputFile.getFilename();
      head       =  getHeader(file_name);
    }
 
    event->mc_run_id = head.start_run.run_id;
 
    if (run == -1) {
      run = event->mc_run_id;
    }
 
    bool   trigger     = false;
    int    frame_index = inputFile.getCounter() + 1;                     // event / slice
    double t1          = gRandom->Rndm() * getFrameTime();               // time of event within time slice
 
    JDAQUTCExtended utc(UTC.getTimeNanoSecond() + getTimeOfFrame(frame_index));
 
    if (runbyrun != "") {
 
      if (event->mc_event_time == TTimeStamp(0)) {
        FATAL("Monte Carlo event time undefined." << endl);
      }
 
      utc = getDAQUTCExtended(event->mc_event_time);      
      t1  = getTimeDifference(summaryRouter.getTimesliceStart(), utc) * 1.0e9;
 
      if (!frame_time(t1)) {
 
        const Long64_t    index = scanner.find(utc);
        JDAQSummaryslice* p     = scanner.getEntry(index);
 
        if (getTimeDifference(p->getTimesliceStart(), utc) < 0.0 && index != 0) {
          p = scanner.getEntry(index - 1);
        }
 
        summaryRouter.update(p);
        summaryRouter.correct(dynamic_cast<const JPMTDefaultSimulatorInterface&>(simbad.getPMTSimulator()));
      }
 
      t1  = getTimeDifference(summaryRouter.getTimesliceStart(), utc) * 1.0e9;
 
      run         = summaryRouter.getRunNumber();
      frame_index = summaryRouter.getFrameIndex();
      utc         = summaryRouter.getTimesliceStart();
 
      DEBUG("event: "
            << FIXED(15,3) << event->mc_event_time.AsDouble() << " [s] "
            << utc                                            << ' '
            << FIXED(12,0) << t1                              << " [ns] "
            << frame_time(t1)                                 << endl);
    }
 
    if (!event->mc_hits.empty() && frame_time(t1)) {
 
      JTimeRange timeRange = getTimeRange(*event, period);
 
      double t0 = 0.5 * (timeRange.getLowerLimit() + timeRange.getUpperLimit());
 
      DEBUG("Start time: " << FIXED(12,2) << t0 << ' ' << FIXED(12,2) << t1 << ' ' << FIXED(9,2) << TA_ns << endl);
 
      t1 += TA_ns;                                                       // add time offset of detector
 
      event->mc_t = getTimeOfFrame(frame_index)  +  t1 - t0;             // set to time since start of data taking run
 
      timeRange.add(event->mc_t);
      timeRange.add(period);
 
      const JDAQChronometer chronometer(detectorB.getID(), run, frame_index, utc);
 
      const JEventTimeslice timeslice(chronometer, simbad, *event, period);
 
      DEBUG(timeslice << endl);
 
    
      timesliceRouter.configure(timeslice);
 
      JTimeslice_t    timesliceL0(timeslice.getDAQChronometer());
      JTimeslice_t    timesliceL1(timeslice.getDAQChronometer());
      JTimeslice_t    timesliceL2(timeslice.getDAQChronometer());
      JTimeslice_t    timesliceSN(timeslice.getDAQChronometer());
 
      for (JDAQTimeslice::const_iterator super_frame = timeslice.begin(); super_frame != timeslice.end(); ++super_frame) {
 
        if (moduleRouter.hasModule(super_frame->getModuleID())) {
 
          // calibration
 
          const JModule&         module = moduleRouter.getModule(super_frame->getModuleID()); 
          const JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*super_frame, module);
 
          // L0
        
          timesliceL0.push_back(JSuperFrame1D_t(buffer));
 
          DEBUG("L0 " << setw(8) << timesliceL0.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL0.rbegin()->size() << LAMBDA([ps = timesliceL0.rbegin()](ostream& out) { for (const auto& i : *ps) { out << " " << i; } }) << endl); 
 
          // L1
        
          timesliceL1.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
                                                super_frame->getModuleIdentifier(),
                                                module.getPosition()));
 
          buildL1(*timesliceL0.rbegin(), back_inserter(*timesliceL1.rbegin()));
        
          DEBUG("L1 " << setw(8) << timesliceL1.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL1.rbegin()->size() << LAMBDA([ps = timesliceL1.rbegin()](ostream& out) { for (const auto& i : *ps) { out << " " << i; } }) << endl); 
 
          // L2
        
          timesliceL2.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
                                                super_frame->getModuleIdentifier(),
                                                module.getPosition()));
        
          buildL2(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceL2.rbegin()));
 
          DEBUG("L2 " << setw(8) << timesliceL2.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL2.rbegin()->size() << LAMBDA([ps = timesliceL2.rbegin()](ostream& out) { for (const auto& i : *ps) { out << " " << i; } }) << endl); 
 
          // SN
          {
            JTimeslice_t::value_type tv(super_frame->getDAQChronometer(),
                                        super_frame->getModuleIdentifier(),
                                        module.getPosition());
 
            buildSN(buffer, *timesliceL1.rbegin(), back_inserter(tv));
 
            if (!tv.empty()) {
 
              timesliceSN.push_back(tv);
 
              DEBUG("SN " << setw(8) << timesliceSN.rbegin()->getModuleID() << ' ' << setw(8) << timesliceSN.rbegin()->size() << endl);
              DEBUG("L2 " << setw(8) << timesliceL2.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL2.rbegin()->size() << LAMBDA([ps = timesliceL2.rbegin()](ostream& out) { for (const auto& i : *ps) { out << " " << i; } }) << endl); 
            }
          }
        }
      }
 
 
      // Trigger
 
      JTriggerInput  trigger_input(timesliceL2);
      JTriggerOutput trigger_output;
 
      trigger3DMuon  (trigger_input, back_inserter(trigger_output));
      trigger3DShower(trigger_input, back_inserter(trigger_output));
      triggerMXShower(trigger_input, timesliceL0, back_inserter(trigger_output));
 
      trigger_output.merge(JEventOverlap(parameters.TMaxEvent_ns));
 
      for (JTriggerOutput::const_iterator to = trigger_output.begin(); to != trigger_output.end(); ++to) {
 
        for (int i = 0; i != h1.GetNbinsX(); ++i) {
          if (to->hasTriggerBit(i)) {
            h1.Fill((double) i);
          }
        }
 
        JTimeRange eventTime = getTimeRange(*to).add(getTimeOfRTS(*to));
 
        eventTime.add(TA_ns);
        eventTime.sub(TB_ns);
 
        if (debug >= status_t) {
          cout << "Event time: "
               << to->getFrameIndex() << ' ' 
               << eventTime           << ' '
               << timeRange           << ' '
               << (timeRange.overlap(eventTime) ? "Y" : "N") << endl;
        }
 
        if (timeRange.overlap(eventTime)) {
        
          JTriggeredEvent tev(*to, 
                              timesliceRouter, 
                              moduleRouter, 
                              parameters.TMaxLocal_ns, 
                              getTimeRange(parameters));
 
          tev.setCounter(trigger_counter);                               // set the event counter to the index of the Monte Carlo event in the output.
 
          outputFile.put(tev);
        
          trigger = true;
        }
      }
   
 
      if (!triggeredEventsOnly || trigger) {
      
        if (parameters.writeL0()) {
          outputFile.put(timeslice);
        }
      
        if (parameters.writeL1()) {
          outputFile.put(JTimesliceL1<JDAQTimesliceL1>(timesliceL1, timesliceRouter, moduleRouter, parameters.TMaxLocal_ns));
        }
      
        if (parameters.writeL2()) {
          outputFile.put(JTimesliceL1<JDAQTimesliceL2>(timesliceL2, timesliceRouter, moduleRouter, parameters.L2.TMaxLocal_ns));
        }
      
        if (parameters.writeSN()) {
          outputFile.put(JTimesliceL1<JDAQTimesliceSN>(timesliceSN, timesliceRouter, moduleRouter, parameters.SN.TMaxLocal_ns));
        }
 
        if (parameters.writeSummary()) {
 
          if (fs.count(chronometer.getFrameIndex()) == 0) {
 
            outputFile.put(JSummaryslice(chronometer,simbad));
 
            fs.insert(chronometer.getFrameIndex());
          }
        }
      }
    }
 
    if (!triggeredEventsOnly || trigger) {
      
      outputFile.put(*event);
 
      ++trigger_counter;
    }
  }  
  STATUS(endl);
 
  {
    JMultipleFileScanner<JTYPELIST<JMetaTypes_t, MultiHead>::typelist> io(inputFile);
 
    io >> outputFile;
  }
 
  outputFile.put(h1);
  outputFile.put(*gRandom);
  outputFile.close();
}