#include #include #include #include #include #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" /** * \file * 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.\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::typelist typelist; JMultipleFileScanner<> inputFile; JFileRecorder 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 JSuperFrame1D_t; typedef JSuperFrame2D JSuperFrame2D_t; typedef JTimeslice JTimeslice_t; typedef JBuildL1 JBuildL1_t; typedef JBuildL2 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(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 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(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(timesliceL1, timesliceRouter, moduleRouter, parameters.TMaxLocal_ns)); } if (parameters.writeL2()) { outputFile.put(JTimesliceL1(timesliceL2, timesliceRouter, moduleRouter, parameters.L2.TMaxLocal_ns)); } if (parameters.writeSN()) { outputFile.put(JTimesliceL1(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::typelist> io(inputFile); io >> outputFile; outputFile.put(h1); outputFile.put(*gRandom); outputFile.close(); }