#include #include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TF1.h" #include "JDAQ/JDAQEventIO.hh" #include "JDAQ/JDAQTimesliceIO.hh" #include "JDAQ/JDAQEvaluator.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JDAQHitRouter.hh" #include "JDetector/JPMTParametersMap.hh" #include "JTrigger/JBuildL0.hh" #include "JTrigger/JBuildL1.hh" #include "JLang/JSinglePointer.hh" #include "JROOT/JManager.hh" #include "JROOT/JROOTClassSelector.hh" #include "JSupport/JSingleFileScanner.hh" #include "JSupport/JTreeScannerInterface.hh" #include "JSupport/JTreeScanner.hh" #include "JSupport/JAutoTreeScanner.hh" #include "JSupport/JSupport.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to search for out of sync shifts around integral timeslices evolving over time. * \author shallmann */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; JSingleFileScanner<> inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; string detectorFile; double TMax_ns; int numberOfTimeslices; double binWidth_ns; double binWidth_min_timeEvolution; JROOTClassSelector selector; int debug; try { JParser<> zap("Example program to search for out of sync shifts around integral timeslices evolving over time."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = ""; zap['a'] = make_field(detectorFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['T'] = make_field(TMax_ns) = 20.0; zap['N'] = make_field(numberOfTimeslices) = 40; zap['W'] = make_field(binWidth_ns) = 10e3; zap['M'] = make_field(binWidth_min_timeEvolution) = 2; zap['C'] = make_field(selector) = "", getROOTClassSelection(); zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } map MASK; // mask PMTs with permament high-rate veto. const int WR = 0x80000000; // White Rabbit MASK[808969848] = 0x00000020; // floor 03, pmt 05 MASK[809544061] = 0x00000080; // floor 05, pmt 07 MASK[808432835] = 0x00004000; // floor 18, pmt 14 JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const JDAQHitRouter router(detector); typedef double hit_type; JBuildL0 buildL0; JBuildL1 buildL1(JBuildL1Parameters(TMax_ns, true)); vector data; JTreeScanner<>::debug = debug; JSinglePointer< JTreeScannerInterface > ps; JAutoTreeScanner zmap = JType(); if (selector == "") { if ((ps = new JTreeScanner< JAssertConversion >(inputFile))->getEntries() != 0 || (ps = new JTreeScanner< JAssertConversion >(inputFile))->getEntries() != 0 || (ps = new JTreeScanner< JAssertConversion >(inputFile))->getEntries() != 0 || (ps = new JTreeScanner< JAssertConversion >(inputFile))->getEntries() != 0 || (ps = new JTreeScanner< JAssertConversion >(inputFile))->getEntries() != 0) { } else { FATAL("No timeslice data." << endl); } NOTICE("Selected class " << ps->getClassName() << endl); } else { ps = zmap[selector]; ps->configure(inputFile); } typedef JManager JManager_t; // shift in integer time slices const Double_t ymin = -(numberOfTimeslices + 0.5); const Double_t ymax = +(numberOfTimeslices + 0.5); const Int_t ny = (Int_t) (ymax - ymin); // time in M minute binning const Double_t xmin = inputFile.getLimit().min()*getFrameTime()*1e-9;// time of first event - offset const Double_t xmax = min(inputFile.getLimit().max(), ps->getEntries())*getFrameTime()*1e-9; // time of last event const Int_t nx = (Int_t) ((xmax - xmin) / (60*binWidth_min_timeEvolution)) + 1; JManager_t manager(new TH2D("M2D_%", ";time [s];shift [timeslices]", nx, xmin, xmax, ny, ymin, ymax)); typedef map > map_type; // frame index -> event times map_type buffer; for (JTreeScanner in(inputFile.getFilename(), inputFile.getLimit()); in.hasNext(); ) { STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl); JDAQEvent* event = in.next(); // Average time of triggered hits double t0 = 0.0; for (JDAQEvent::const_iterator hit = event->begin(); hit != event->end(); ++hit) { t0 += getTime(*hit, router.getPMT(*hit)); } t0 /= event->size(); t0 += event->getFrameIndex() * getFrameTime(); buffer[event->getFrameIndex()].push_back(t0); } STATUS(endl); while (ps->hasNext()) { STATUS("event: " << setw(10) << ps->getCounter() << '\r'); DEBUG(endl); JDAQTimeslice* timeslice = ps->next(); map_type::const_iterator p = buffer.lower_bound(timeslice->getFrameIndex() - numberOfTimeslices); map_type::const_iterator q = buffer.upper_bound(timeslice->getFrameIndex() + numberOfTimeslices); int number_of_events = 0; for (map_type::const_iterator i = p; i != q; ++i) { number_of_events += i->second.size(); } if (number_of_events == 0) { continue; } for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) { // Test that none of the PMTs has high-rate veto. if ((frame->getStatus() & ~MASK[frame->getModuleID()] & ~WR) == 0) { data.clear(); buildL1(*frame, router.getModule(frame->getModuleID()), back_inserter(data)); TH2D* h2 = manager[frame->getModuleID()]; for (vector::const_iterator hit = data.begin(); hit != data.end(); ++hit) { const double t1 = *hit + frame->getFrameIndex() * getFrameTime(); for (map_type::const_iterator i = p; i != q; ++i) { for (map_type::mapped_type::const_iterator j = i->second.begin(); j != i->second.end(); ++j) { const double t0 = *j; // fill only shift is near integer time slice const double timeFromInt = (t1 - t0)/getFrameTime() - round((t1 - t0)/getFrameTime()); if ( abs(timeFromInt)*getFrameTime() < binWidth_ns/2. ){ h2->Fill(timeslice->getFrameIndex() * getFrameTime() * 1e-9, round(t1 - t0)/getFrameTime()); } } } } } } } STATUS(endl); if (outputFile != "") { manager.Write(outputFile.c_str()); } }