#include #include #include #include #include #include #include "TROOT.h" #include "TFile.h" #include "TH1D.h" #include "antares-dataformat/TimeSlice.hh" #include "km3net-dataformat/online/JDAQSummaryslice.hh" #include "JDAQ/JDAQSummarysliceIO.hh" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JModuleRouter.hh" #include "JGeometry3D/JPosition3D.hh" #include "JGeometry3D/JRotation3D.hh" #include "JGeometry3D/JOmega3D.hh" #include "JTools/JQuantile.hh" #include "JSupport/JMultipleFileScanner.hh" #include "JSupport/JSupport.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" namespace { using JTOOLS::JQuantile; using JDETECTOR::JModuleRouter; using KM3NETDAQ::JDAQSummaryslice; /** * Auxiliary data structure to store rate as a function of index of module in detector. */ struct JRates : public std::vector { /** * Constructor. * * \param router module router * \param summary summary data */ JRates(const JModuleRouter& router, const JDAQSummaryslice& summary) : std::vector(router.getReference().size()) { const double factor = 1.0e-3; // [kHz] for (JDAQSummaryslice::const_iterator frame = summary.begin(); frame != summary.end(); ++frame) { if (router.hasModule(frame->getModuleID())) { const int index = router.getIndex(frame->getModuleID()); // index of module in detector for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) { (*this)[index].put(frame->getRate(pmt, factor)); // rate [kHz] } } } } /** * Constructor. * * \param router module router * \param summary summary data */ JRates(const JModuleRouter& router, const ExtendedSummary_TimeSlice& summary) : std::vector(router.getReference().size()) { ANTARES::setClock(); // Antares DAQ clock const double factor = 2.0e6 / getFrameTime(); // 2 ARS / PMT [kHz] for (ExtendedSummary_TimeSlice::const_iterator frame = summary.begin(); frame != summary.end(); ++frame) { if (router.hasModule(frame->lcm_id())) { const int index = router.getIndex(frame->lcm_id()); // index of module in detector const int count = frame->numberOfItemsOrg() << 4; // decompress data (*this)[index].put(count * factor); // rate [kHz] } } } }; /** * Main method. * * \param argc number of command line arguments * \param argv list of command line arguments * \return status */ template bool do_main(int argc, char **argv) { using namespace std; using namespace JPP; using namespace KM3NETDAQ; JMultipleFileScanner inputFile; JLimit_t& numberOfEvents = inputFile.getLimit(); string outputFile; string detectorFile; double roadWidth_m; double gridAngle_deg; int numberOfPMTs; int debug; try { JParser<> zap("Example program to analyse summary data."); zap['f'] = make_field(inputFile); zap['o'] = make_field(outputFile) = "monopole.root"; zap['a'] = make_field(detectorFile); zap['n'] = make_field(numberOfEvents) = JLimit::max(); zap['R'] = make_field(roadWidth_m); zap['G'] = make_field(gridAngle_deg); zap['N'] = make_field(numberOfPMTs); zap['d'] = make_field(debug) = 2; zap(argc, argv); } catch(const exception& error) { FATAL(error.what() << endl); } if (roadWidth_m <= 0.0) { FATAL("Invalid road width [m] " << roadWidth_m << endl); } if (gridAngle_deg <= 0.0) { FATAL("Invalid grid angle [deg] " << gridAngle_deg << endl); } if (gridAngle_deg >= 90.0) { FATAL("Invalid grid angle [deg] " << gridAngle_deg << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } DEBUG("Number of modules in detector " << detector.size() << endl); const JModuleRouter router(detector); // look-up table of module in detector const JOmega3D omega(gridAngle_deg * PI/180.0); // set of assumed monopole directions typedef vector< set > buffer1D; // index of module in detector -> indices of modules in detector typedef vector< buffer1D > buffer2D; // index of direction -> index of module in detector -> indices of modules in detector buffer2D M2(omega.size(), buffer1D(detector.size())); // look-up table of modules within road width around given module for each direction for (size_t i = 0; i != omega.size(); ++i) { const JRotation3D R(omega[i]); // rotates z-axis to assumed monopole direction for (size_t m1 = 0; m1 != detector.size(); ++m1) { // index of first module in detector for (size_t m2 = 0; m2 != m1; ++m2) { // index of second module in detector JPosition3D p1 = detector[m1].getPosition(); JPosition3D p2 = detector[m2].getPosition(); p1.rotate(R); p2.rotate(R); const double x = p2.getX() - p1.getX(); const double y = p2.getY() - p1.getY(); if (sqrt(x*x + y*y) <= roadWidth_m) { // application of road width M2[i][m1].insert(m2); } } } } TFile out(outputFile.c_str(), "recreate"); TH1D h1("h1", NULL, 1000, 0.0, 1.0e3); // average rate [kHz] while (inputFile.hasNext()) { STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); const T* summary = inputFile.next(); const JRates buffer(router, *summary); // rate as a function of index of module in detector for (size_t i = 0; i != omega.size(); ++i) { // scan directions JQuantile Q; // rate per direction for (size_t m1 = 0; m1 != buffer.size(); ++m1) { // index of module in detector Q += buffer[m1]; const set& M1 = M2[i][m1]; // indices of modules in detector for (set::const_iterator m2 = M1.begin(); m2 != M1.end(); ++m2) { Q += buffer[*m2]; } } if (Q.getCount() >= numberOfPMTs) { h1.Fill(Q.getMean()); // rate per PMT } } } STATUS(endl); out.Write(); out.Close(); return (inputFile.getCounter() != 0 ? true : false); } } /** * \file * * Example program to analyse summary data. * \author mdejong */ int main(int argc, char **argv) { if (do_main(argc, argv)) { return 0; } if (do_main (argc, argv)) { return 0; } }