JMonopole.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
#include <vector>
#include <set>

#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<JQuantile>
  {
    /**
     * Constructor.
     *
     * \param  router        module router
     * \param  summary       summary data
     */
    JRates(const JModuleRouter& router, const JDAQSummaryslice& summary) :
      std::vector<JQuantile>(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<JQuantile>(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<class T>
  bool do_main(int argc, char **argv)
  {
    using namespace std;
    using namespace JPP;
    using namespace KM3NETDAQ;

    JMultipleFileScanner<T> 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);
    }


    cout.tie(&cerr);

    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<int> >  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<int>& M1 = M2[i][m1];                            // indices of modules in detector

          for (set<int>::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<KM3NETDAQ::JDAQSummaryslice>(argc, argv)) { return 0; }
  if (do_main<ExtendedSummary_TimeSlice>  (argc, argv)) { return 0; }
}