#include #include #include #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 inputFile; JFileRecorder ::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* 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(inputFile); JTreeScanner::debug = debug; } else { // unsorted access; scan = new JTreeScanner(inputFile); JTreeScanner::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); }