JTriggerEfficiency.cc File Reference

Auxiliary program to trigger Monte Carlo events. More...

#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 "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 "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 "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JRunByRun.hh"
#include "JSupport/JMeta.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Auxiliary program to trigger Monte Carlo events.

The options

• -a refers to the detector configuration that is used to convert Monte Carlo true information to raw data; and
• -b 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.
This allows for correlating the DAQ event to the Monte Carlo event.
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).
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.
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.
The two-fold (and higher) coincidence rates should still be provided on the command line.

Author

mdejong

Definition in file JTriggerEfficiency.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 100 of file JTriggerEfficiency.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTYPELIST<JAllTypes_t, JRootTypes_t>::typelist  typelist;

  JMultipleFileScanner<>   inputFile;
  JFileRecorder<typelist>  outputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 detectorFileA;
  string                 detectorFileB;
  int                    run;
  JTriggerParameters     parameters;
  bool                   triggeredEventsOnly;
  JPMTParametersMap      pmtParameters;
  JK40Rates              rates_Hz;
  JRunByRun              runbyrun;
  double                 sigma_ns;
  UInt_t                 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['R'] = make_field(run,                 "run number")  = -1;
    zap['r'] = make_field(runbyrun,            "option 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['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);
  }

  gRandom->SetSeed(seed);

  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);
  }

  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;

  if (runbyrun.is_valid()) {

    NOTICE("Using run-by-run:" << endl << runbyrun << endl);

    if (!runbyrun.hasNext()) {
      FATAL("Run-by-run simulation yields no input." << endl);
    }

    if (rates_Hz.getSinglesRate() != 0.0) {
      WARNING("Run-by-run simulation discards singles rate [Hz] " << rates_Hz.getSinglesRate() << endl);
    }

    try {    
      simbad.reset(new JK40RunByRunSimulator(summaryRouter, rates_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->getFilelist());

      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 = getLivetime(runbyrun->getFilelist());
    buffer.push(&JHead::DAQ);

    copy(buffer, header);

  } else {

    NOTICE("Using fixed rates [Hz]: " << rates_Hz << endl);

    try {    
      simbad.reset(new JK40DefaultSimulator(rates_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));

  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);

  const JPosition3D center = get<JPosition3D>(header);

  NOTICE("Apply detector offset from Monte Carlo run header (" << center << ")" << endl);

  detectorA -= center;
  detectorB -= center;

  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);

  JLimit_t     limit            = inputFile.getLimit();
  counter_type number_of_events = 0;
  int          trigger_counter  = 0;
  
  for (JMultipleFileScanner<>::const_iterator file = inputFile.begin(); file != inputFile.end(); ++file) {

    int mc_run_id = 0;
    
    try {

      const JHead head(getHeader(*file));

      mc_run_id = head.start_run.run_id;
    }
    catch(const JException& error) {
      FATAL(error);
    }
    
    JMultipleFileScanner<Evt> in(*file);

    limit.setLowerLimit(limit.getLowerLimit() - in.skip(limit.getLowerLimit()));
    
    for ( ; in.hasNext() && number_of_events != limit; ++number_of_events) {

      STATUS("event: " << setw(10) << number_of_events << '\r'); DEBUG(endl);
    
      Evt* event = in.next();

      event->mc_run_id = mc_run_id;
    
      DEBUG(*event << endl);

      bool trigger = false;

      if (!event->mc_hits.empty()) {
      
        int frame_index = (int) in.getCounter() + 1;                        // 1 event / slice

        if (runbyrun.is_valid() && runbyrun.hasNext()) {

          summaryRouter.update(runbyrun.next());

          summaryRouter.correct(dynamic_cast<const JPMTDefaultSimulatorInterface&>(simbad.getPMTSimulator()));
      
          frame_index = summaryRouter.getFrameIndex();
          run         = summaryRouter.getRunNumber();
        }

        // set the event time! 

        JTimeRange timeRange = getTimeRange(*event, period);

        if (!timeRange.is_valid()) {
          timeRange.setRange(0.0,0.0);
        }

        const double t0 = 0.5 * (timeRange.getLowerLimit() + timeRange.getUpperLimit());
        const double t1 = gRandom->Rndm() * getFrameTime();

        event->mc_t = getTimeOfFrame(frame_index)  +  t1 - t0;             // time since start of data taking run

        timeRange.add(event->mc_t);
        timeRange.add(period);

        JDAQUTCExtended utc(getTimeOfFrame(trigger_counter + 1));          // ensure incremental UTC times (e.g. for JDAQSplit.cc)

        if (event->mc_event_time != TTimeStamp(0)) {
          utc = getDAQUTCExtended(event->mc_event_time, t1);               // set UTC time according Monte Carlo event
        }

        const JDAQChronometer chronometer(detectorB.getID(), (run != -1 ? run : mc_run_id), 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));

            // L1
        
            timesliceL1.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
                                                  super_frame->getModuleIdentifier(),
                                                  module.getPosition()));

            buildL1(*timesliceL0.rbegin(), back_inserter(*timesliceL1.rbegin()));
        
            // L2
        
            timesliceL2.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
                                                  super_frame->getModuleIdentifier(),
                                                  module.getPosition()));
        
            buildL2(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceL2.rbegin()));

            // SN
        
            timesliceSN.push_back(JSuperFrame1D_t(super_frame->getDAQChronometer(),
                                                  super_frame->getModuleIdentifier(),
                                                  module.getPosition()));
        
            buildSN(buffer, *timesliceL1.rbegin(), back_inserter(*timesliceSN.rbegin()));

            DEBUG("L0 " << setw(8) << timesliceL0.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL0.rbegin()->size() << endl);
            DEBUG("L1 " << setw(8) << timesliceL1.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL1.rbegin()->size() << endl);
            DEBUG("L2 " << setw(8) << timesliceL2.rbegin()->getModuleID() << ' ' << setw(8) << timesliceL2.rbegin()->size() << endl);
            DEBUG("SN " << setw(8) << timesliceSN.rbegin()->getModuleID() << ' ' << setw(8) << timesliceSN.rbegin()->size() << 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));

          DEBUG("Event time: "
                << to->getFrameIndex() << ' ' 
                << eventTime           << ' '
                << timeRange           << endl);

          if (timeRange.overlap(eventTime)) {
        
            JTriggeredEvent tev(*to, 
                                timesliceRouter, 
                                moduleRouter, 
                                parameters.TMaxLocal_ns, 
                                getTimeRange(parameters));

            // Set the event counter to the index of the Monte Carlo event in the output.

            tev.setCounter(trigger_counter);

            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()) {
            outputFile.put(JSummaryslice(chronometer,simbad));
          }      
        }
      }

      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();
}