#include #include #include #include #include #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Hit.hh" #include "km3net-dataformat/offline/io_online.hh" #include "km3net-dataformat/tools/time_converter.hh" #include "JDAQ/JDAQEventIO.hh" #include "JAAnet/JAAnetToolkit.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JModuleRouter.hh" #include "JSupport/JTriggeredFileScanner.hh" #include "JSupport/JMonteCarloFileSupportkit.hh" #include "JSupport/JSupport.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { /** * Count number of coincident hits. * * \param data1 first data set * \param data2 second data set * \param Tmax_ns time window [ns] * \return number of hits */ inline int get_count(const std::set& data1, const std::set& data2, const double Tmax_ns) { int number_of_hits = 0; for (std::set::const_iterator p = data1.begin(), q = data2.begin(); ; ++p, ++q) { while (p != data1.end() && q != data2.end() && *p < *q - Tmax_ns) { ++p; } while (p != data1.end() && q != data2.end() && *q < *p - Tmax_ns) { ++q; } if (p == data1.end() || q == data2.end()) { break; } if (fabs(*p - *q) <= Tmax_ns) { ++number_of_hits; } } return number_of_hits; } } /** * \file * * Example program to test conversion between Monte Carlo and DAQ times. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; JTriggeredFileScanner<> inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string detectorFile; double Tmax_ns; int debug; try { JParser<> zap("Example program to test conversion between Monte Carlo and DAQ times."); zap['f'] = make_field(inputFile); zap['a'] = make_field(detectorFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['T'] = make_field(Tmax_ns) = 20.0; zap['d'] = make_field(debug) = 3; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } using namespace KM3NETDAQ; JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const JModuleRouter router(detector); //typedef JDAQSnapshotHit JHit_t; typedef JDAQTriggeredHit JHit_t; while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next(); const JDAQEvent* tev = ps; const Evt* event = ps; ASSERT(!tev->empty()); typedef map > map_type; map_type mc; map_type data; for (vector::const_iterator i = event->mc_hits.begin(); i != event->mc_hits.end(); ++i) { if (!is_noise(*i)) { DEBUG("Hit (MC): " << setw(5) << i->pmt_id << ' ' << FIXED(6,1) << getTime(*i) << endl); mc[i->pmt_id].insert(getTime(*i)); } } for (int test : {1, 2} ) { DEBUG("Test "<< test << endl); data.clear(); time_converter converter; if (test == 1) { converter = time_converter(*event, *tev); for (JDAQEvent::const_iterator i = tev->begin(); i != tev->end(); ++i) { const JModule& module = router.getModule(i->getModuleID()); const JPMT& pmt = module.getPMT(i->getPMT()); data[pmt.getID()].insert(converter.getTime(getTime(*i, pmt))); } } else if (test == 2) { Evt evt = *event; read(evt, *tev); converter = time_converter(evt); for (vector::const_iterator i = evt.hits.begin(); i != evt.hits.end(); ++i) { if (i->trig != 0) { const JModule& module = router.getModule(i->dom_id); const JPMT& pmt = module.getPMT(i->channel_id); data[pmt.getID()].insert(converter.getTime(getTime(i->tdc, pmt))); //data[pmt.getID()].insert(converter.getTime(getTime(i->tdc, pmt))); } } } for (map_type::const_iterator i = data.begin(); i != data.end(); ++i) { for (set::const_iterator hit = i->second.begin(); hit != i->second.end(); ++hit) { DEBUG("Hit (DAQ): " << setw(5) << i->first << ' ' << FIXED(6,1) << *hit << endl); } } int number_of_hits = 0; for (map_type::const_iterator p = mc.begin(), q = data.begin(); ; ++p, ++q) { while (p != mc.end() && q != data.end() && p->first < q->first) { ++p; } while (p != mc.end() && q != data.end() && q->first < p->first) { ++q; } if (p == mc.end() || q == data.end()) { break; } if (p->first == q->first) { number_of_hits += get_count(p->second, q->second, Tmax_ns); } } NOTICE("Number of hits " << setw(5) << number_of_hits << '/' << setw(5) << tev->size() << endl); ASSERT(number_of_hits != 0); } } ASSERT(inputFile.getCounter() != 0); return 0; }