#include "km3net-dataformat/online/JDAQ.hh" #include "km3net-dataformat/online/JDAQTimeslice.hh" #include "km3net-dataformat/online/JDAQSummaryslice.hh" #include "JDAQ/JDAQTimesliceIO.hh" #include "JDAQ/JDAQSummarysliceIO.hh" #include "JDAQ/JDAQEvaluator.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JModuleRouter.hh" #include "JTrigger/JHitL1.hh" #include "JTrigger/JBuildL2.hh" #include "JTrigger/JPMTSelector.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JTreeScanner.hh" #include "JSupport/JSupport.hh" #include "JLang/JObjectMultiplexer.hh" #include "JROOT/JROOTClassSelector.hh" #include "Jeep/JParser.hh" #include "JROOT/JManager.hh" #include "JTools/JCombinatorics.hh" #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TMath.h" #include #include using namespace JTOOLS; namespace MONITORL1DT { /** * Data structure for hit time and DOM identifier. */ class JElement { public: JElement() : id(0), t (0.0) {} JElement(const int __id, const double __t) : id(__id), t (__t) {} int id; double t; }; /** * Sort JElement in time. * * \param first first element * \param second second element * \return true if second later than first; else false */ inline bool operator<(const JElement& first, const JElement& second) { return first.t < second.t; } /** * Auxiliary data structure for histogram management. */ struct JHistogram { /** * Default constructor. */ JHistogram() : h2s(NULL), h1l(NULL) {} /** * Constructor. * * \param __h2s 2D histogram for signal * \param __h1l 1D histogram for background */ JHistogram(TH2D* __h2s, TH1D* __h1l) : h2s(__h2s), h1l(__h1l) { h2s->Sumw2(); h1l->Sumw2(); } TH2D* h2s; TH1D* h1l; }; // Taken from C. Pieterse nanobeacon code typedef JCombinatorics::pair_type pair_type; inline bool comparepair(const pair_type& A, const pair_type& B) { return( abs(A.first - A.second) > abs(B.first - B.second) ) ; } } /** * \file * * Auxiliary program to plot L1 time difference correlations between DOMs. * * \author kmelis, lnauta */ int main(int argc, char **argv) { using namespace std; using namespace MONITORL1DT; using namespace JPP; JMultipleFileScanner inputFile; string outputFile; string detectorFile; double Timewindow_ns; int binwidth; // since the data is taken on a ns level, having bins of non-integer value is non-sensical double TmaxL1_ns; JROOTClassSelector selector; unsigned int multiplicity; bool correct_time; double livetime_s; int debug; try { JParser<> zap("Program to create L1 hit time difference histograms from raw data."); zap['f'] = make_field(inputFile, "input file"); zap['o'] = make_field(outputFile, "output file") = "monitor.root"; zap['a'] = make_field(detectorFile, "detector file"); zap['t'] = make_field(TmaxL1_ns, "max time between L1 hits [ns]") = 1000.0; zap['T'] = make_field(Timewindow_ns, "time window around t=0 [ns]") = 2400.0; zap['w'] = make_field(binwidth, "binwidth [ns]") = 1; zap['C'] = make_field(selector, "datastream selector") = getROOTClassSelection(); zap['m'] = make_field(multiplicity, "minimal multiplicity of the L1 hits") = 2; zap['c'] = make_field(correct_time, "subtract expected arrival time from delta-t"); zap['L'] = make_field(livetime_s, "livetime of the data, set to positive value") = -1.0; // set to positive value if you want to set the livetime manually (i.e. MC) zap['d'] = make_field(debug) = 1; if (zap.read(argc, argv) != 0) return 1; } catch(const exception &error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } if (detector.empty()) { FATAL("Empty detector." << endl); } const JModuleRouter router(detector); typedef vector JFrameL1_t; const double ctMin = -1; // full angle over the DOM const JBuildL2 buildL2(JL2Parameters(multiplicity, TmaxL1_ns, ctMin)); // Manager histograms for single-DOM output vector zmap(detector.size()); const int nx = detector.size(); const double xmin = -0.5; const double xmax = nx - 0.5; const double ymin = -floor(Timewindow_ns) + 0.5; const double ymax = +floor(Timewindow_ns) + 0.5; // changed - to + to get an even number for Nbins, rebin needs numbers that have divisors. const int ny = (int) ((ymax - ymin) / binwidth); // Manage histograms for single-DU output Int_t num_of_floors = getNumberOfFloors(detector); JCombinatorics c(num_of_floors); Int_t npairs = c.getNumberOfPairs(); c.sort(MONITORL1DT::comparepair); JManager* manager; manager = new JManager (new TH2D("h%", "", npairs, 0.5, npairs+0.5, ny, ymin, ymax)); NOTICE("Running JMonitorL1dt: Monitoring L1 time differences and creating histograms." << endl); for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) { const JModuleAddress& address = router.getAddress(module->getID()); STATUS("Booking histograms for module " << module->getID() << endl); const JString title(module->getID()); JString titleString1D; JString titleString2D; titleString1D = title + ".1L"; titleString2D = title + ".2S"; zmap[address.first] = JHistogram(new TH2D((titleString2D).c_str(), NULL, nx, xmin, xmax, ny, ymin, ymax), new TH1D((titleString1D).c_str(), NULL, nx, xmin, xmax)); for (JDetector::iterator mod = detector.begin(); mod != detector.end(); ++mod) { zmap[address.first].h2s->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID())); zmap[address.first].h1l->GetXaxis()->SetBinLabel(distance(detector.begin(), mod)+1, Form("%i", mod->getID())); } } JTreeScanner summaryFile(inputFile); JLANG::JObjectMultiplexer in(inputFile, selector); int counter = 0; for ( ; in.hasNext() && counter != inputFile.getLimit(); ++counter) { STATUS("event: " << setw(10) << counter << '\r'); DEBUG(endl); const JDAQTimeslice* timeslice = in.next(); JFrameL1_t frameL1; vector buffer; vector DOM_OK(detector.size(), true); for (JDAQTimeslice::const_iterator super_frame = timeslice->begin(); super_frame != timeslice->end(); ++super_frame) { if (router.hasModule(super_frame->getModuleID()) && !super_frame->empty()) { const JModuleAddress& address = router.getAddress(super_frame->getModuleID()); const JModule& module = detector.getModule(address); frameL1.clear(); buildL2(*super_frame, module, back_inserter(frameL1)); for (JFrameL1_t::iterator L1hit = frameL1.begin(); L1hit != frameL1.end(); ++L1hit) { buffer.push_back(JElement(address.first, L1hit->begin()->getT())); } } } for (vector::iterator h1 = zmap.begin(); h1 != zmap.end(); ++h1) { if (!DOM_OK[distance(zmap.begin(), h1)]) { continue; } for (unsigned int i = 0; i < detector.size(); ++i) { if (DOM_OK[i]) { h1->h1l->Fill(i, getFrameTime() * 1e-9); } // fill the th1 with the frametime (how long did it measure) in seconds, all weight will give total measured time } } // fill histograms with correlations sort(buffer.begin(), buffer.end()); for (vector::const_iterator p = buffer.begin(); p != buffer.end(); ) { vector::const_iterator q = p; const JModule& module_p = detector.getModule((JModuleAddress)p->id); while (++q != buffer.end() && q->t - p->t <= Timewindow_ns ) { // if (p->id == q->id) { continue; } const JModule& module_q = detector.getModule((JModuleAddress)q->id); double dom_distance = getDistance(module_p.getPosition(), module_q.getPosition()); double time_correction = (correct_time ? (dom_distance / getSpeedOfLight()) : 0); zmap[p->id].h2s->Fill(q->id, q->t - p->t - time_correction); zmap[q->id].h2s->Fill(p->id, p->t - q->t + time_correction); if (module_p.getString() == module_q.getString()) { // indices are pairs running from 0-17 between 0-152, while floors are 1-18 and bins are 1-153 int xbin = c.getIndex(module_p.getFloor() - 1, module_q.getFloor() - 1) + 1; (*manager)[module_p.getString()]->Fill(xbin, q->t - p->t - time_correction); } } p++; // move to the next L1 hit } buffer.clear(); } STATUS(endl); // livetime in case set manually if (livetime_s > 0.0) { for (vector::iterator i = zmap.begin(); i != zmap.end(); ++i) { TH1D* hl = i->h1l; for (int ibin = 1; ibin <= hl->GetNbinsX(); ++ibin) { hl->SetBinContent(ibin, livetime_s); hl->SetBinError(ibin, 0.0000001); } } } TFile out(outputFile.c_str(), "recreate"); for (vector::iterator i = zmap.begin(); i != zmap.end(); ++i) { i->h2s->Write(); i->h1l->Write(); } manager->Write(out); out.Write(); out.Close(); }