#include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "JTools/JQuantile.hh" #include "JTools/JRange.hh" #include "JROOT/JManager.hh" #include "JROOT/JRootToolkit.hh" #include "JDAQ/JDAQEventIO.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JDAQHitRouter.hh" #include "JPhysics/JConstants.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JSupport.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to histogram string and floor time difference. * \author mdejong */ int main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; typedef JRange JRange_t; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; string detectorFile; JRange_t T_ns; double Wmin; double Qmin; int qaqc; int debug; try { JParser<> zap("Example program to histogram string and floor time difference."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = "squid.root"; zap['a'] = make_field(detectorFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['T'] = make_field(T_ns, "Time window for coincidences [ns]") = JRange_t(-50.0, +200.0); zap['W'] = make_field(Wmin, "Minimal number of entries") = 10.0; zap['q'] = make_field(Qmin, "Minimal fraction of coincidences") = 0.5; zap['Q'] = make_field(qaqc) = 0; zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const JDAQHitRouter router(detector); JManager H1(new TH1D("%", NULL, 200, -1.0e3, +1.0e3)); while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); JDAQEvent* event = inputFile.next(); map buffer; for (JDAQEvent::const_iterator hit = event->begin(); hit != event->end(); ++hit) { buffer[router.getModule(*hit)].put(getTime(*hit, router.getPMT(*hit)) + router.getModule(*hit).getZ() * getInverseSpeedOfLight()); } if (buffer.size() >= 2u) { for (map::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++p, ++q) { if (p->first.getString() == q->first.getString() && p->first.getFloor() + 1 == q->first.getFloor()) { H1[p->first]->Fill(q->second.getMean() - p->second.getMean()); } } } } STATUS(endl); TFile out(outputFile.c_str(), "recreate"); out << H1; out.Write(); out.Close(); int nin = 0; int nout = 0; for (const auto& i : H1) { Double_t W = 0.0; TH1D* h1 = i.second; if (h1->GetSumOfWeights() > Wmin) { for (Int_t ix = 1; ix <= h1->GetXaxis()->GetNbins(); ++ix) { const Double_t x = h1->GetBinCenter (ix); const Double_t y = h1->GetBinContent(ix); if (T_ns(x)) { W += y; } } DEBUG(i.first << ' ' << FIXED(6,0) << W << '/' << FIXED(6,0) << h1->GetSumOfWeights() << endl); if (W / h1->GetSumOfWeights() >= Qmin) nin += 1; else nout += 1; } } NOTICE("Number of modules out/in micro-sync " << nout << '/' << nin << endl); QAQC(nin << ' ' << nout << endl); return 0; }