JEventTimesliceWriter.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <fstream>

#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 "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtEvaluator.hh"

#include "JDetector/JCLBDefaultSimulator.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSimulator.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JK40DefaultSimulator.hh"
#include "JDetector/JPMTDefaultSimulator.hh"

#include "JDAQ/JDAQSummarysliceIO.hh"
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JTimeslice/JEventTimeslice.hh"
#include "JTimeslice/JRandomTimeslice.hh"

#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JTreeScannerInterface.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 * Auxiliary program to convert multiple Monte Carlo events to DAQ timeslices.
 * Events are timed according to a given rate (if > 0), otherwise they are timed according to Evt::mc_t. 
 * \author mdejong, mlincett
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;


  JMultipleFileScanner<JAAnetTypes_t>  inputFile;
  JFileRecorder       <JTYPELIST<JAAnetTypes_t, JDAQTimesliceL0, JMeta, JRootTypes_t>::typelist>  outputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 detectorFile;
  int                    run;
  JPMTParametersMap      pmtParameters;
  JK40Rates              rates_Hz;
  double                 eventRate_Hz;
  UInt_t                 seed;
  int                    debug;

  try { 

    JParser<> zap("Auxiliary program to convert multiple Monte Carlo events to time slices.");
    
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "timeslice.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['R'] = make_field(run)                 = 1;
    zap['P'] = make_field(pmtParameters)       = JPARSER::initialised();
    zap['B'] = make_field(rates_Hz)            = JPARSER::initialised();
    zap['r'] = make_field(eventRate_Hz);
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 0;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  gRandom->SetSeed(seed);


  setDAQLongprint(debug >= JEEP::debug_t);
  
  if (pmtParameters.getQE() != 1.0) {

    WARNING("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);
    
    rates_Hz.correct(pmtParameters.getQE());
  }

  // parameters

  DEBUG("PMT paramaters: " << endl << pmtParameters << endl);
  DEBUG("K40 rates:      " << endl << rates_Hz      << endl);
 
  if (eventRate_Hz < 0.0) {
    FATAL("Invalid event rate " << eventRate_Hz << "; consider using JRandomTimesliceWriter." << endl);
  }
  
  // detector
  
  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  
  // header
  
  Head header;
  
  try {
    header = getHeader(inputFile);
  }
  catch (const exception& error) {
    FATAL("Monte Carlo header is invalid");
  }
  
  double liveTimeMC = JHead(header).livetime.numberOfSeconds;

  if (liveTimeMC < 0) {
    FATAL("Monte Carlo live time is negative; input file may be corrupted.");
  }

  // background simulator

  JPMTParametersMap::Throw(false);

  JDetectorSimulator simbad(detector);

  simbad.reset(new JPMTDefaultSimulator(pmtParameters, detector));
  simbad.reset(new JK40DefaultSimulator(rates_Hz));
  simbad.reset(new JCLBDefaultSimulator());

  // output file
  
  outputFile.open();

  if (!outputFile.is_open()) {
    FATAL("Error opening file " << outputFile << endl);
  }
  
  outputFile.put(*gRandom);
  outputFile.put(JMeta(argc, argv));
  outputFile.put(header);

  // input stream;
  
  bool absTime = false;
  
  JTreeScannerInterface<Evt>* scan;

  if (eventRate_Hz == 0.0 && liveTimeMC == 0.0) {
    
    NOTICE("Event will be timed according to absolute MC time." << endl);
    
    absTime = true;
    
    // sorted access;
    scan = new JTreeScanner<Evt, JEvtEvaluator>(inputFile);
    JTreeScanner<Evt, JEvtEvaluator>::debug = debug;
    
  } else {
    
    // unsorted access;
    scan = new JTreeScanner<Evt>(inputFile);
    JTreeScanner<Evt>::debug = debug;
    
    if (eventRate_Hz == 0.0 && liveTimeMC > 0.0) {
      int nEntries = scan->getEntries();
      eventRate_Hz = nEntries / liveTimeMC;
      DEBUG(nEntries << " events to be written." << endl);
      NOTICE("MC live time is " << liveTimeMC << " s" << endl);
      NOTICE("Event rate set to " <<  eventRate_Hz << " Hz from MC header." << endl);
    } 
  }
  
  // begin timeslice generation;
  
  // state variables;
  bool pendingEvt       = false;
  bool pendingTimeslice = false;
  
  Evt* event = NULL; // pending event;
  
  JTimeRange timeRange;
  double tDAQ = 0; // DAQ event time;
  double tOff = 0; // hit production time offset;
  
  int frame_index = 0;
  JDAQChronometer chronometer;
  JRandomTimeslice timeslice; 
  
  counter_type evtCount = 0;

  while (pendingEvt || scan->hasNext()) {

    if (!pendingTimeslice) {
      // update background timeslice;
      frame_index++;
      chronometer = JDAQChronometer(detector.getID(), run, frame_index, JDAQUTCExtended(getTimeOfFrame(frame_index)));
      timeslice = JRandomTimeslice(chronometer, simbad);
      DEBUG("evt count: " << setw(10) << evtCount << endl); 
      STATUS("frame index: " << setw(10) << frame_index << " | evt count: " << setw(10) << evtCount << "\r"); DEBUG(endl);
      pendingTimeslice = true;
    }
    
    if (!pendingEvt) {
      // get new event;
      event = scan->next();
      timeRange = getTimeRange(*event);
      tOff = timeRange.is_valid() ? timeRange.getLowerLimit() : 0;
      // update DAQ time according to absolute or random;
      if (absTime)
        tDAQ  = getEvtValue(*event) + tOff;
      else
        tDAQ += gRandom->Exp(1.0e9 / eventRate_Hz);
      
      DEBUG("event time [s] " << setprecision(5) << tDAQ * 1.0e-9 << endl);
      pendingEvt = true;
    }
    
    // DEBUG("tDAQ : " << tDAQ << " - frame index: " << frame_index << " - DAQ frame index " << getFrameIndex(tDAQ) << endl);
    if (getFrameIndex(tDAQ) == frame_index) {
      
      // process event;
      if (timeRange.is_valid()) {
        DEBUG(*event << endl);
        // real start time of the hits will be always tDAQ;
        event->mc_t = tDAQ - tOff;
        timeslice.add(JEventTimeslice(chronometer, simbad, *event));
        evtCount++;
      }
      // event processed;
      outputFile.put(*event);
      pendingEvt = false;
    } else {
      // event time is past end of timeslice;
      DEBUG(timeslice << endl);
      outputFile.put(timeslice);
      pendingTimeslice = false;
    }
    
  }
  
  if (pendingTimeslice) {
    DEBUG(timeslice << endl);
    outputFile.put(timeslice);
    pendingTimeslice = false;
  }
  
  STATUS(endl);

  outputFile.put(*gRandom);
  outputFile.close();
  
  NOTICE(evtCount << " events written over " << frame_index << " timeslices. " << endl);
}