JRunAnalyzer.hh

Go to the documentation of this file.

#ifndef __JRUNANALYZER__
#define __JRUNANALYZER__
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
 
#include "JSupport/JSupport.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JROOT/JManager.hh"
#include "JTools/JQuantile.hh"
#include "JTools/JRange.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"
#include "JMath/JConstants.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "TFile.h"
#include "TH1D.h"
#include "JRunHistograms.hh"
 
/**
 * \author rgruiz, adomi
 */
 
using JSUPPORT::JSingleFileScanner;
 
/**
 * Class dedicated to the analysis of %KM3NeT runs. 
 *
 */
class JRunAnalyzer {
 
  JSingleFileScanner<>  inputFile;
  JModuleRouter         router;
  TFile                 *outputFile;
  JRA_Histograms        histograms;
  JLimit_t              numberOfTimeslices;
  JLimit_t              numberOfSummaryslices;
  JLimit_t              numberOfEvents;
  int                   analysis_level;
  
public :
 
  /**
   * Constructor.
   *
   * \param file            file name
   * \param detector        detector
   * \param out             pointer to output file
   * \param nTimeslices     number of timeslices to be read
   * \param nSummaryslices  number of summary slices to be read
   * \param nEvents         number of events to be read
   * \param analysislevel   option for analysis of PMT or reco data
   */
  JRunAnalyzer (const JSingleFileScanner<>& file, const JDetector& detector, TFile *out, 
                JLimit_t nTimeslices, JLimit_t nSummaryslices, JLimit_t nEvents, int analysislevel):
    inputFile                       (file),
    router                      (detector),
    outputFile                       (out),
    numberOfTimeslices       (nTimeslices),
    numberOfSummaryslices (nSummaryslices),
    numberOfEvents               (nEvents),
    analysis_level         (analysislevel)
  {
    using namespace JPP;
    using namespace std;
    
    histograms     = JRA_Histograms(detector);
  } 
 
  /**
   * Destructor.
   */
  ~JRunAnalyzer(){}
  
  /*
   * Function template used to read events from a tree.
   *
   * \param scanner A JTreeScanner
   * \param 
   */
  void iterateEventTree (JTreeScanner < JDAQEvent > & scanner) {
 
    while (scanner.hasNext()) {
 
      JDAQEvent event = *(scanner.next());
 
      histograms.h_trigger.h_Snapshot_hits      -> Fill((Double_t) event.size<JDAQSnapshotHit> ());
      histograms.h_trigger.h_Triggered_hits     -> Fill((Double_t) event.size<JDAQTriggeredHit> ());
      histograms.h_trigger.h_Triggered_over_Snapshot_hits -> Fill((Double_t) event.size<JDAQTriggeredHit> () / event.size<JDAQSnapshotHit > ());
      histograms.h_trigger.h_Number_of_overlays -> Fill(event.getOverlays());
 
      for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
 
        if (event.hasTriggerBit(i)) {
 
          histograms.h_trigger.h_Trigger_bit_event -> Fill((Double_t) i);
        }
      }
 
      int counter_3dmuon = 0;
      typedef JRange<double>  JRange_t;
      JRange_t range(JRange_t::DEFAULT_RANGE());
      
      for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event.begin<JDAQTriggeredHit>(); hit != event.end<JDAQTriggeredHit>(); ++hit) {
 
        const JModule& module = router.getModule(hit->getModuleID());
        const double   t1     = getTime(hit->getT(), module.getPMT(hit->getPMT()));
        
        range.include(t1);
        
        JDAQTriggerMask trigger_mask(event.getTriggerMask(*hit));
 
        if(trigger_mask.hasTriggerBit(getTriggerBit<JTrigger3DMuon>())) counter_3dmuon++; 
         
        if (router.hasModule(hit->getModuleID())) {
 
          const JModule& module = router.getModule (hit->getModuleID());
        
          histograms.h_trigger.h_pmt_distribution_triggered_hits->Fill(hit->getPMT());
          histograms.h_trigger.h_tot_distribution_triggered_hits->Fill(hit->getToT());
 
          int String = module.getString();
          int Floor  = module.getFloor();
 
          histograms.h_trigger.h_Triggered_hits_per_module        -> Fill(String , Floor);
 
          if(trigger_mask.hasTriggerBit(getTriggerBit<JTrigger3DMuon>())){
            histograms.h_trigger.h_Triggered_hits_3dmuon_per_module -> Fill(String , Floor);
          }
          
          for (unsigned int i = 0; i != NUMBER_OF_TRIGGER_BITS; ++i) {
          
            if (hit -> hasTriggerBit(i)) {
 
              histograms.h_trigger.h_Trigger_bit_hit->Fill((Double_t) i);
            }
          }
        }else{
          FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
        }
      }
 
      histograms.h_trigger.h_event_duration->Fill(range.getLength());
 
      histograms.h_trigger.h_pmt_distribution_triggered_hits->Scale(1. / (Double_t) event.size<JDAQTriggeredHit> ());
     
      if(counter_3dmuon != 0) histograms.h_trigger.h_Triggered_hits_3dmuon->Fill(counter_3dmuon);
          
      for (JDAQEvent::const_iterator<JDAQSnapshotHit> hit = event.begin<JDAQSnapshotHit>() ; hit != event.end<JDAQSnapshotHit>() ; ++hit){
 
        if (router.hasModule(hit->getModuleID())) {
 
          const JModule& module = router.getModule (hit->getModuleID());
        
          int String = module.getString();
          int Floor  = module.getFloor();
          int pmt    = hit->  getPMT();
 
          if(analysis_level == 1){
            
            (*histograms.h_trigger.m_Snapshot_hits_per_pmt)[MAKE_STRING("Detector/DU" + to_string(String))] -> Fill(pmt, Floor);
 
          }
 
          histograms.h_trigger.h_pmt_distribution_snapshot_hits -> Fill(hit->getPMT()); 
          histograms.h_trigger.h_tot_distribution_snapshot_hits -> Fill(hit->getToT());
          histograms.h_trigger.h_Snapshot_hits_per_module       -> Fill(String, Floor);
          
          }else{
            FATAL("JModuleRouter trying to access non existing identifier: "<< hit->getModuleID());
        }
      }
    }
  }
 
  /*
   * Function template used to read reconstructed events from a tree.
   *
   * \param scanner A ParallelFileScanner for JDAQEvent and JEvt
   * \param 
   */
  void iterateRecoEventTree (JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > & scanner, TFile & f) {
 
    TH1D* h_tres = new TH1D("h_tres", ";Time residuals [ns]; Entries", 100, -50, 150);
    TH1D* h_likelihood = new TH1D ("h_likelihood", " ; Likelihood; Reconstructed Events ", 100, -1000, 1000);
    TH1D* h_beta0 = new TH1D("h_beta0", "; beta0; Reconstructed Events", 20, 0, 1);
    TH1D* h_energy = new TH1D ("h_energy", " ; Energy [GeV]; Reconstructed Events ", 65, 0, 9);
    setLogarithmicX(h_energy);
    TH1D* h_zenith = new TH1D("h_zenith", "; cosZenith; Reconstructed Events", 20, -1, 1);
    TH1D* h_azimuth = new TH1D("h_azimuth", "; cosAzimuth; Reconstructed Events", 20, -1, 1);
    TH1D* h_radial_position = new TH1D ("h_radial_position", "; Radial Position [m]; Reconstructed Events", 60, 0, 4.2);
    setLogarithmicX(h_radial_position);
    TH1D* h_z_position = new TH1D ("h_z_position", "; Z Position [m] ; Reconstructed Events", 50, 0, 3.2);
    setLogarithmicX(h_z_position);
 
    while (scanner.hasNext()) {
 
      JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> >::multi_pointer_type  ps  = scanner.next();
 
      JEvt* evt = ps;
      JDAQEvent* event = ps;
      
      if (!evt->empty()) {
 
        JEvt::iterator best = evt->begin();  
        h_energy -> Fill(best->getE());
        h_likelihood -> Fill(best->getQ());
        h_z_position -> Fill(best->getZ());
        h_radial_position -> Fill(sqrt(best->getX()*best->getX() + best->getY()*best->getY()));
        
        JTrack3D track = getTrack(*best); 
        h_zenith -> Fill(cos(track.getDirection().getTheta()));
        h_azimuth -> Fill(cos(track.getDirection().getPhi()));
        h_beta0 -> Fill(best->getW(JGANDALF_BETA0_RAD));
 
        std::vector<JHitL0> dataL0;
        const JBuildL0<JHitL0> buildL0;
        buildL0(JDAQTimeslice(*event, false),  router, back_inserter(dataL0));
 
        for (std::vector<JHitL0>::const_iterator hit = dataL0.begin(); hit != dataL0.end(); ++hit) {
 
          h_tres -> Fill(hit->getT() - track.getT(hit->getPosition()));
        }
      }
    }
 
    f.mkdir("Reco");                                                                                                          
    f.cd("Reco");                                                                                                            
 
    h_energy->Write();
    h_likelihood -> Write();
    h_z_position -> Write();
    h_radial_position -> Write();
    h_zenith -> Write();
    h_azimuth -> Write();
    h_beta0 -> Write();
    h_tres -> Write();
  }
 
  struct dom_type :
    public std::vector< JQuantile>
  {
    dom_type() :
      std::vector<JQuantile>(NUMBER_OF_PMTS)
    {}
  };
 
  /*
   * Function template to read and process time ordered summary slices.
   * 
   * \param scanner A JTreeScanner from which time ordered summary slices are accessed. 
   */
  void iterateSummarysliceTree (JTreeScanner < JDAQSummaryslice > & scanner) {
 
    using namespace std;
    using namespace JPP;
 
    map<int, map<int, dom_type> >       PMT_rate_quantiles;
    map<int, map<int, JQuantile> >      DOM_rate_quantiles;
 
    while (scanner.hasNext()){
 
      JDAQSummaryslice slice = *(scanner.next());
 
      for (JDAQSummaryslice::const_iterator frame = slice.begin() ; frame != slice.end() ; ++ frame) {
 
        if (router.hasModule(frame->getModuleID())) {
 
          const JModule& module = router.getModule (frame->getModuleID());
 
          int string = module.getString();
          int floor  = module.getFloor ();
 
          JDAQFrameStatus status = frame -> getDAQFrameStatus();
          int nFIFOcount         = status.countFIFOStatus();
          int nHRVcount          = status.countHighRateVeto();
 
          if (nHRVcount > 0) {
            histograms.h_summary.h_hrv_per_dom->Fill(string , floor, nHRVcount);
          }
        
          if (status.testDAQStatus() == false) {
            histograms.h_summary.h_daq_status_per_dom->Fill(string , floor, (1./distance(scanner.begin(), scanner.end()))*100);   
          }
        
          histograms.h_summary.h_fifo_per_dom->Fill(string , floor , nFIFOcount); 
 
          if(analysis_level == 1){
          
            TH2D*  h2   = (*histograms.h_summary.m_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))];
 
            double rate = 0;    
            const double factor = 1.0e-3;   // [kHz]
              
            for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
 
              rate += frame->getRate(i, factor);
 
              h2->Fill(i , frame->getRate(i, factor));
 
              PMT_rate_quantiles[string][floor][i].put(frame->getRate(i, factor)); 
              
              histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i, factor), frame->getWeight(i, factor));
 
            }
 
            histograms.h_summary.h_dom_rate_distribution->Fill(rate);
            
            DOM_rate_quantiles[string][floor].put(rate);
            
          } else {
 
            double rate = 0;    
            const double factor = 1.0e-3;   // [kHz]
                  
            for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
 
              rate += frame->getRate(i, factor);
 
              histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i, factor), frame->getWeight(i, factor));
            }
 
            histograms.h_summary.h_dom_rate_distribution->Fill(rate);
 
            DOM_rate_quantiles[string][floor].put(rate);
            
          }
          
        }else{
          FATAL("JModuleRouter trying to access non existing identifier: "<< frame->getModuleID());
        }      
      }      
    }
 
    for (const auto& i1 : DOM_rate_quantiles) {
      for (const auto& i2 : i1.second) {
        if (i2.second.getCount() > 0) {
          histograms.h_summary.h_rate_summary->Fill(i1.first, i2.first, i2.second.getMean());
        }
      }
    }
 
    if (analysis_level == 1){
      for (const auto& i1 : PMT_rate_quantiles){
 
        TH2D* h2 = (*histograms.h_summary.m_mean_summary_rate)             [MAKE_STRING("Detector/DU" + to_string(i1.first))];
        TH1D* h1 = (*histograms.h_summary.m_mean_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(i1.first))];
 
        for (const auto& i2 : i1.second) {
          for (int i3 = 0 ; i3 != NUMBER_OF_PMTS ; ++i3) {
            if (i2.second[i3].getCount() > 0){
              h2->Fill(i3, i2.first, i2.second[i3].getMean());
              h1->Fill(i2.second[i3].getMean());
            }
          }
        }
      }
    }
  }
  
  /*
   * Function template to read and process time ordered timeslices.
   * 
   * \param scanner A JTreeScanner from which time ordered timeslices are accessed. 
   */
  template <class T>
  void iterateTimesliceTree (JTreeScanner < T > & scanner) {
    
    int ts_type = JIndexOf<JDAQTimesliceTypes_t , T>::value;
 
    std::map < int , std::map <int , JQuantile> > DOM_rate_quantiles;
 
    double inverseFrameTimeSec = 1. / (1.0e-9 * getFrameTime());
 
    while (scanner.hasNext()){
      
      T slice = *(scanner.next());
 
      for(auto & frame : slice) {
        if (router.hasModule(frame.getModuleID())) {
 
          const JModule& module = router.getModule (frame.getModuleID());
          double  rate   = frame.numberOfHits * inverseFrameTimeSec;
          int     string = module.getString();
          int     floor  = module.getFloor ();
 
          DOM_rate_quantiles[string][floor].put(rate);
 
          vector <int>  pmt_hit_count (NUMBER_OF_PMTS , 0) ;
 
          if(analysis_level == 1){
          
            TH2D* h2 = (*histograms.h_timeslice.m_pmt_tot_distributions[ts_type])[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor) + "/" + to_string(frame.getModuleID()))];
            
            for (JDAQSuperFrame::const_iterator hit = frame.begin() ; hit != frame.end() ; ++hit){
          
            h2 -> Fill(hit->getPMT() , hit->getToT()) ;
            
            pmt_hit_count[hit->getPMT()]++;
            
            }
            
            for (int pmt = 0 ; pmt != NUMBER_OF_PMTS ; ++pmt) {
            
              (*histograms.h_timeslice.m_pmt_rate_distributions[ts_type])[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))] -> Fill(pmt , 1e-3 * pmt_hit_count[pmt] * inverseFrameTimeSec);
                  
            }
 
          }
            
          }else{
            FATAL("JModuleRouter trying to access non existing identifier: "<< frame.getModuleID());
          }      
        }      
    }
          
    for (std::map< int , std::map<int,JQuantile>>::const_iterator i = DOM_rate_quantiles.begin() ; i != DOM_rate_quantiles.end() ; i++){
 
      for (std::map<int , JQuantile>::const_iterator j = i->second.begin() ; j != i->second.end() ; j++){
 
        if (j->second.getCount() > 0) histograms.h_timeslice.h_dom_mean_rates[ts_type] -> Fill (i->first , j->first , j->second.getMean() ) ;
      }
    }
  }
 
  /*
   * Checks through the use of a JTreeScanner, if summary slices are present in the run file.
   * In case they are, the corresponding tree is iterated. 
   *
   */
  void readSummaryData() {
 
    JTreeScanner <JDAQSummaryslice> scanner(inputFile, numberOfSummaryslices);
    if (scanner.hasNext()) {
      histograms.initialize_summary_histograms();
      scanner.rewind();
      iterateSummarysliceTree(scanner);
    }
  }
 
  /*
   * Function template that checks through the use of a JTreeScanner, if objects of different timeslice classes are present in the run file.
   * In case they are, the corresponding tree is iterated. 
   *
   */
  template <class T>
  void readTimesliceData() {
 
    JTreeScanner <T> scanner(inputFile, numberOfTimeslices);
    if(scanner.hasNext()) {
      
      histograms.initialize_timeslice_histograms <T>();
      scanner.rewind();
      iterateTimesliceTree(scanner);
    }
  }
 
  /*
   * Function template that checks through the use of a JTreeScanner, if JDAQEvent objects are present in the run file.
   * In case they are, the corresponding tree is read. 
   *
   */
  void readEvents() {
 
    JTreeScanner <JDAQEvent> scanner(inputFile, numberOfEvents);
        
    if(scanner.hasNext()) {
 
      histograms.initialize_trigger_histograms();
      scanner.rewind();
      iterateEventTree(scanner);
    }
  }
 
  /*
   * Function template that checks through the use of a JTreeScanner, if JEvt objects are present in the run file.
   * In case they are, the corresponding tree is read. 
   */
  void readRecoEvents() {
 
    JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> >  scanner(inputFile, numberOfEvents);
        
    if(scanner.hasNext()) {
 
      scanner.rewind();
      iterateRecoEventTree(scanner, *outputFile);
    }
  }
 
   /*                                                                   
   * Writes the histograms to a root file.                                 
   * \param f The root file.                                    
   */
  void writeToFile(TFile *f,int analysis_level){
 
    f->mkdir("Detector");
    f->cd("Detector");
 
    if (histograms.h_summary.h_fifo_per_dom)          histograms.h_summary.h_fifo_per_dom          -> Write();
    if (histograms.h_summary.h_daq_status_per_dom)    histograms.h_summary.h_daq_status_per_dom    -> Write();
    if (histograms.h_summary.h_hrv_per_dom)           histograms.h_summary.h_hrv_per_dom           -> Write();
    if (histograms.h_summary.h_rate_summary)          histograms.h_summary.h_rate_summary          -> Write();
    if (histograms.h_summary.h_pmt_rate_distribution) histograms.h_summary.h_pmt_rate_distribution -> Write();
    if (histograms.h_summary.h_dom_rate_distribution) histograms.h_summary.h_dom_rate_distribution -> Write();
 
    histograms.Write_histogram_table_to_file(*f , MAKE_STRING("Detector"), histograms.h_timeslice.h_dom_mean_rates);
 
    for (typename vector < JManager < string , TH2D >* >::const_iterator i = histograms.h_timeslice.m_pmt_tot_distributions.begin();
         i != histograms.h_timeslice.m_pmt_tot_distributions.end() ; ++i){
 
      if ((*i)){
 
        for (typename JManager < string , TH2D >::const_iterator j = (*i) -> begin() ; j != (*i) -> end() ; ++j){
 
          j->second->Scale(1., "width");
        }
      }
    }
 
    histograms.Write_manager_table_in_key_dir (*f , histograms.h_timeslice.m_pmt_tot_distributions);
    histograms.Write_manager_table_in_key_dir (*f , histograms.h_timeslice.m_pmt_rate_distributions);
 
    if (histograms.h_summary.m_summary_rate_distribution) histograms.Write_manager_in_key_dir (*f , histograms.h_summary.m_summary_rate_distribution);
 
    if (histograms.h_summary.m_mean_summary_rate)              histograms.Write_manager_in_key_dir(*f , histograms.h_summary.m_mean_summary_rate);
    if (histograms.h_summary.m_mean_summary_rate_distribution) histograms.Write_manager_in_key_dir(*f , histograms.h_summary.m_mean_summary_rate_distribution);
 
    histograms.h_timeslice.Fill_mean_ToT_histograms();
 
    histograms.Write_manager_table_in_key_dir (*f , histograms.h_timeslice.m_mean_ToT);
    histograms.Write_manager_table_in_key_dir (*f , histograms.h_timeslice.m_mean_ToT_distribution);
 
    f->mkdir ( MAKE_STRING ("JDAQEvent").c_str());
    f->cd    ("JDAQEvent");
 
    if (histograms.h_trigger.h_Trigger_bit_event) histograms.h_trigger.h_Trigger_bit_event -> Write();
    if (histograms.h_trigger.h_Trigger_bit_hit)   histograms.h_trigger.h_Trigger_bit_hit   -> Write();
    if (histograms.h_trigger.h_Triggered_hits)        histograms.h_trigger.h_Triggered_hits        -> Write();
    if (histograms.h_trigger.h_Triggered_hits_3dmuon) histograms.h_trigger.h_Triggered_hits_3dmuon -> Write();
    if (histograms.h_trigger.h_Triggered_hits_3dmuon_per_module) histograms.h_trigger.h_Triggered_hits_3dmuon_per_module -> Write();
    if (histograms.h_trigger.h_Snapshot_hits)         histograms.h_trigger.h_Snapshot_hits         -> Write();
    if (histograms.h_trigger.h_Triggered_over_Snapshot_hits)       histograms.h_trigger.h_Triggered_over_Snapshot_hits  -> Write();
    if (histograms.h_trigger.h_event_duration)  histograms.h_trigger.h_event_duration -> Write();
    if (histograms.h_trigger.h_Number_of_overlays)  histograms.h_trigger.h_Number_of_overlays -> Write();
    if (histograms.h_trigger.h_pmt_distribution_triggered_hits) {histograms.h_trigger.h_pmt_distribution_triggered_hits -> Write();}
    if (histograms.h_trigger.h_pmt_distribution_snapshot_hits)  {histograms.h_trigger.h_pmt_distribution_snapshot_hits  -> Write();}
    if (histograms.h_trigger.h_tot_distribution_triggered_hits) { histograms.h_trigger.h_tot_distribution_triggered_hits->Scale(1, "width") ; histograms.h_trigger.h_tot_distribution_triggered_hits-> Write();}
    if (histograms.h_trigger.h_tot_distribution_snapshot_hits)  {histograms.h_trigger.h_tot_distribution_snapshot_hits->Scale(1, "width") ; histograms.h_trigger.h_tot_distribution_snapshot_hits -> Write();}
    if (histograms.h_trigger.h_Triggered_hits_per_module) { histograms.h_trigger.h_Triggered_hits_per_module-> Write();}
    if (histograms.h_trigger.h_Snapshot_hits_per_module)  { histograms.h_trigger.h_Snapshot_hits_per_module -> Write();}
    if (histograms.h_trigger.m_Snapshot_hits_per_pmt) histograms.Write_manager_in_key_dir(*f , histograms.h_trigger.m_Snapshot_hits_per_pmt);
  
  }
};
 
#endif