#include <string> #include <iostream> #include <iomanip> #include <vector> #include <set> #include <TH1D.h> #include "antcc/Header.hh" #include "antcc/Event.hh" #include "antcc/Hit.hh" #include "JDAQ/JDAQTimeslice.hh" #include "JDAQ/JDAQEvent.hh" #include "JDAQ/JDAQSummaryslice.hh" #include "JTimeslice/JEventTimeslice.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JDetectorSimulator.hh" #include "JDetector/JModuleRouter.hh" #include "JDetector/JTimeRange.hh" #include "JDetector/JK40DefaultSimulator.hh" #include "JDetector/JPMTDefaultSimulator.hh" #include "JDetector/JCLBDefaultSimulator.hh" #include "JTrigger/JHit.hh" #include "JTrigger/JHitToolkit.hh" #include "JTrigger/JFrame.hh" #include "JTrigger/JTimeslice.hh" #include "JTrigger/JHitL0.hh" #include "JTrigger/JHitL1.hh" #include "JTrigger/JBuildL1.hh" #include "JTrigger/JBuildL2.hh" #include "JTrigger/JTrigger3DShower.hh" #include "JTrigger/JTrigger3DMuon.hh" #include "JTrigger/JEventOverlap.hh" #include "JTrigger/JTimesliceRouter.hh" #include "JTrigger/JTriggeredEvent.hh" #include "JTrigger/JTriggerParameters.hh" #include "JTrigger/JEventToolkit.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JFileRecorder.hh" #include "JSupport/JMonteCarloToolkit.hh" #include "JSupport/JSupport.hh" #include "JSirene/JSirene.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * Auxiliary program to trigger Monte Carlo events. */ int main(int argc, char **argv) { using namespace std; using namespace JSUPPORT; using namespace JTRIGGER; JMultipleFileScanner<JMonteCarloTypes_t> inputFile; JFileRecorder <JLANG::JAppend<JAllTypes_t, TRandom>::typelist> outputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string detectorFile; int run; JTriggerParameters parameters; set<int> triggerBits; bool triggeredEventsOnly; bool writeRawHits; bool generateBackground; UInt_t seed; int debug; try { JParser<> zap; zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = "trigger_efficieny.root"; zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['a'] = make_field(detectorFile); zap['R'] = make_field(run) = 1; zap['@'] = make_field(parameters); zap['T'] = make_field(triggerBits); zap['O'] = make_field(triggeredEventsOnly); zap['s'] = make_field(writeRawHits); zap['B'] = make_field(generateBackground); zap['S'] = make_field(seed) = 0; zap['d'] = make_field(debug) = 0; zap['@'] = JPARSER::initialised(); zap['T'] = JPARSER::initialised(); if (zap.read(argc, argv) != 0) return 1; } catch(const exception &error) { FATAL(error.what() << endl); } gRandom->SetSeed(seed); using namespace JSIRENE; using namespace KM3NETDAQ; using namespace JTRIGGER; using JTOOLS::make_range; cout.tie(&cerr); setLongprint(cout, debug >= 3); JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } parameters.set(getMaximalDistance(detector)); const JTrigger3DMuon trigger3DMuon (parameters); const JTrigger3DShower trigger3DShower(parameters); JModuleRouter moduleRouter(detector); JDetectorSimulator simbad(detector); if (generateBackground) simbad.reset(new JK40DefaultSimulator()); simbad.reset(new JPMTDefaultSimulator()); simbad.reset(new JCLBDefaultSimulator()); const double Tmax = getMaximalTime(detector); const JTimeRange period(-Tmax, +Tmax); typedef double hit_t; typedef vector <hit_t> JFrameL1_t; typedef JSuperFrame<hit_t> JSuperFrame_t; typedef JTimeslice <hit_t> JTimeslice_t; typedef JBuildL1 <hit_t> JBuildL1_t; typedef JBuildL2 <hit_t> JBuildL2_t; const JBuildL1_t buildL1(JHitToolkit<hit_t>::getT(parameters.TMaxLocal_ns)); Header header; try { header = inputFile.getHeader(); } catch(const JException& error) { FATAL(error); } if (header.coord_origin) { detector -= JPoint3D(header.coord_origin.x, header.coord_origin.y, header.coord_origin.z); } outputFile.open(); if (!outputFile.is_open()) FATAL("Error opening file " << outputFile << endl); outputFile.put(*gRandom); TH1D h1("Trigger bits", NULL, NUMBER_OF_TRIGGER_BITS, -0.5, NUMBER_OF_TRIGGER_BITS - 0.5); outputFile.put(header); outputFile.put(parameters); JTriggerCounter_t counter = 0; for (JMultipleFileScanner<Event>& in = inputFile; in.hasNext(); ) { STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl); Event* event = in.next(); const vector<Hit>& hitlist = event->HitList(); for (vector<Hit>::const_iterator i = hitlist.begin(); i != hitlist.end(); ++i) DEBUG(setw(5) << i->pm_id() << ' ' << fixed << setw(7) << setprecision(2) << i->pure_dt() << endl); // set the event time! const int frame_index = (int) in.getCounter() + 1; // 1 event / slice JTimeRange timeRange = getTimeRange(*event); if (!timeRange.is_valid()) timeRange.setRange(0.0,0.0); const double t0 = 0.5 * (timeRange.getLowerLimit() + timeRange.getUpperLimit()); const double t1 = getTimeOfFrame(frame_index) + 0.5 * getFrameTime(); event->eventTime(t1 - t0); // time since start of data taking run timeRange.add(event->eventTime()); const JDAQChronometer chronometer(run, frame_index, JDAQUTCExtended(getTimeOfFrame(frame_index))); const JEventTimeslice timeslice(chronometer, simbad, *event, period); DEBUG(timeslice << endl); const JTimesliceRouter timesliceRouter(timeslice, parameters.numberOfBins); JFrameL1_t frameL1; JTimeslice_t timesliceL1(timeslice.getDAQChronometer()); JSuperFrame_t buffer; for (JDAQTimeslice::const_iterator super_frame = timeslice.begin(); super_frame != timeslice.end(); ++super_frame) { // calibration const JModule& module = moduleRouter.getModule(super_frame->getModuleID()); buffer(*super_frame, module); // L1 frameL1.clear(); buildL1(buffer, back_inserter(frameL1)); // L2 const JBuildL2_t buildL2(buffer, 2, parameters.TMaxLocal_ns, parameters.ctMin); timesliceL1.push_back(JFrame<hit_t>(super_frame->getDAQChronometer(), JDAQPMTIdentifier(super_frame->getModuleID(), -1), module.getPosition())); buildL2(frameL1, back_inserter(timesliceL1.rbegin()->getData())); } // Trigger JTriggerInput trigger_input(timesliceL1); JTriggerOutput trigger_output; if (triggerBits.count(trigger3DMuon .getTriggerBit())) trigger3DMuon (trigger_input, back_inserter(trigger_output)); if (triggerBits.count(trigger3DShower.getTriggerBit())) trigger3DShower(trigger_input, back_inserter(trigger_output)); trigger_output.merge(JEventOverlap(parameters.TMaxEvent_ns)); bool trigger = false; for (JTriggerOutput::const_iterator to = trigger_output.begin(); to != trigger_output.end(); ++to) { for (int i = 0; i != h1.GetNbinsX(); ++i) { if (checkTriggerMask(*to,i)) h1.Fill((double) i); } JTimeRange eventTime = getTimeRange(*to).add(getTimeOfRTS(*to)); DEBUG("Event time: " << to->getFrameIndex() << ' ' << eventTime.getLowerLimit() << ' ' << eventTime.getUpperLimit() << ' ' << timeRange.getLowerLimit() << ' ' << timeRange.getUpperLimit() << endl); if (timeRange.overlap(eventTime)) { JTriggeredEvent tev(*to, timesliceRouter, moduleRouter, parameters.TMaxLocal_ns, parameters.getTimeRange()); // Set the event counter to the index of the Monte Carlo event in the output. tev.setCounter(counter); outputFile.put(tev); trigger = true; } } if (!triggeredEventsOnly || trigger) { if (writeRawHits) { const double t0 = event->eventTime() - getTimeOfRTS(timeslice.getFrameIndex()); vector<RawHit>& rawHits = const_cast<vector<RawHit>&>(event->RawHitList()); rawHits.clear(); for (JDAQTimeslice::const_iterator super_frame = timeslice.begin(); super_frame != timeslice.end(); ++super_frame) { const JModule& module = moduleRouter.getModule(super_frame->getModuleID()); for (JDAQSuperFrame::const_iterator hit = super_frame->begin(); hit != super_frame->end(); ++hit) { const JPMT& pmt = module.getPMT(hit->getPMT()); const double tdc = hit->getT() - t0; const double tot = hit->getToT(); const unsigned int id = rawHits.size() + 1; rawHits.push_back(RawHit(id, pmt.getID(), tdc, tot)); } } } outputFile.put(*event); ++counter; } } STATUS(endl); outputFile.put(*gRandom); outputFile.close(); }