RunAnalyzer.hh

Go to the documentation of this file.

#ifndef __RUNANALYZER__
#define __RUNANALYZER__
 
#include "JSupport/JSupport.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JGizmo/JManager.hh"
#include "JTools/JQuantile.hh"
#include "TFile.h"
#include "TH1D.h"
#include "JRunHistograms.hh"
 
using namespace std ;
using namespace KM3NETDAQ ;
using namespace JDETECTOR ;
using namespace JLANG ;   
using namespace JPP ;
 
/**
 * \author rgruiz
 */
 
/**
 * Class dedicated to the analysis of KM3NeT runs. 
 *
 */
class RunAnalyzer {
 
  string          inputFile;
  JDetector       detector;
  JModuleRouter*  router;
  JRA_Histograms  histograms;
  int             modulesPerString;
  JLimit_t        numberOfTimeslices;
  JLimit_t        numberOfSummaryslices;
  JLimit_t        numberOfEvents;
 
public :
 
  /**
   * Default constructor.
   */
  RunAnalyzer(){}
 
  /**
   * Constructor.
   *
   * \param file Path to the file to be analyzed
   * \param detectorFile path to a detector file
   * \param nTimeslices Number of timeslices to be read
   * \param nSummaryslices Number of summary slices to be read
   * \param nEvents Number of frames to be read
   */
  RunAnalyzer (string file , string detectorFile , JLimit_t nTimeslices , JLimit_t nSummaryslices , JLimit_t nEvents):
    inputFile                       (file),
    router                          (NULL),
    numberOfTimeslices       (nTimeslices),
    numberOfSummaryslices (nSummaryslices),
    numberOfEvents               (nEvents)
  {
    try {
      load (detectorFile, detector);
    }
    catch(const JException & error) {
      cerr << "FATAL ERROR. Could not open detector file '" << detectorFile << "'." << endl;
      exit(1);
    }
    
    histograms       = JRA_Histograms(detector);
    router           = new JModuleRouter(detector);
    modulesPerString = JDETECTOR::getNumberOfModules (detector) / getNumberOfStrings (detector);
  } 
 
  /**
   * Destructor.
   */
  ~RunAnalyzer(){}
 
  /*
   * Function template used to read time ordered events from a tree.
   *
   * \param scanner A JTreeScanner
   * \param 
   */
  void iterateEventTree (JTreeScanner < JDAQEvent , JDAQEvaluator > & 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> ());
           
      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);
        }
      }
      
      for (JDAQEvent::const_iterator<JDAQTriggeredHit> hit = event.begin<JDAQTriggeredHit>(); hit != event.end<JDAQTriggeredHit>(); ++hit) {
    
        histograms.h_trigger.h_pmt_distribution_triggered_hits->Fill(hit->getPMT());
        histograms.h_trigger.h_tot_distribution_triggered_hits->Fill(hit->getToT());
 
        int String = router->getModule (hit->getModuleID()).getString();
        int Floor  = router->getModule (hit->getModuleID()).getFloor();
 
        histograms.h_trigger.h_Triggered_hits_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);
          }
        }
      }
 
      for (JDAQEvent::const_iterator<JDAQSnapshotHit> hit = event.begin<JDAQSnapshotHit>() ; hit != event.end<JDAQSnapshotHit>() ; ++hit){
 
        int String = router -> getModule(hit -> getModuleID()).getString();
        int Floor  = router -> getModule(hit -> getModuleID()).getFloor();
        int pmt    = hit    -> getPMT() ;
 
        (*histograms.h_trigger.m_Snapshot_hits_per_pmt)[MAKE_STRING("Detector/DU" + to_string(String))] -> Fill(pmt, Floor);
        histograms.h_trigger.h_Snapshot_hits_per_module                                                 -> Fill(String, Floor);
        histograms.h_trigger.h_pmt_distribution_snapshot_hits                                           -> Fill(hit->getPMT()); 
        histograms.h_trigger.h_tot_distribution_snapshot_hits                                           -> Fill(hit->getToT());
      }
    }
  }
 
 
   /*
   * 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 , JDAQEvaluator > & scanner) {
    
    std::map <int,vector<vector <JQuantile>>> PMT_rate_quantiles;
    std::map <int,vector<JQuantile>>          DOM_rate_quantiles;
    
    while (scanner.hasNext()){
 
      JDAQSummaryslice slice = *(scanner.next());
 
      for (JDAQSummaryslice::const_iterator frame = slice.begin() ; frame != slice.end() ; ++ frame) {
 
        PMT_rate_quantiles[router->getModule(frame->getModuleID()).getString()].resize(modulesPerString , vector <JQuantile>(NUMBER_OF_PMTS)); 
 
        JDAQFrameStatus status = frame -> getDAQFrameStatus();
        int nFIFOcount         = status.countFIFOStatus();
        int nHRVcount          = status.countHighRateVeto();
 
        int string = router->getModule(frame->getModuleID()).getString();
        int floor  = router->getModule(frame->getModuleID()).getFloor ();
 
        if (nHRVcount > 0) {
          histograms.h_summary.h_hrv_per_dom->Fill(router -> getModule (frame->getModuleID()).getString() , router->getModule (frame->getModuleID()).getFloor());
        }
        
        if (status.testDAQStatus() == false) {
          histograms.h_summary.h_daq_status_per_dom->Fill(router -> getModule (frame->getModuleID()).getString() , router->getModule (frame->getModuleID()).getFloor());          
        }
        
        histograms.h_summary.h_fifo_per_dom->Fill(router -> getModule (frame->getModuleID()).getString() , router->getModule (frame->getModuleID()).getFloor() , nFIFOcount); 
 
        TH2D*  h2   = (*histograms.h_summary.m_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))];
        double rate = 0;        
 
        for (int i = 0 ; i < NUMBER_OF_PMTS ; i++){
          
          h2->Fill(i , frame->getRate(i) * 1e-3);
          
          rate += frame->getRate(i);
          histograms.h_summary.h_pmt_rate_distribution->Fill(frame->getRate(i)*1e-3);
          PMT_rate_quantiles[router->getModule(frame->getModuleID()).getString()][router->getModule(frame->getModuleID()).getFloor()-1][i].put(frame->getRate(i)*1e-3); 
        }
 
        histograms.h_summary.h_dom_rate_distribution->Fill(rate*1e-3);
 
        DOM_rate_quantiles[router->getModule (frame->getModuleID()).getString()].resize(modulesPerString);
        DOM_rate_quantiles[router->getModule (frame->getModuleID()).getString()][router->getModule(frame->getModuleID()).getFloor()-1].put(rate * 1e-3);
      }      
    }
        
    for (std::map<int , vector<JQuantile> >::const_iterator i = DOM_rate_quantiles.begin() ; i!= DOM_rate_quantiles.end() ; ++i) {
      
      for (int j=1 ; j < modulesPerString + 1 ; j++) {
      
        if (i->second[j-1].getCount() > 0) histograms.h_summary.h_rate_summary->Fill(i->first, j, i->second[j-1].getMean());
      }
    }
    
    for (std::map<int , vector < vector<JQuantile> > >::const_iterator i = PMT_rate_quantiles.begin() ; i!= PMT_rate_quantiles.end() ; ++i) {
    
      for (int j=0 ; j < modulesPerString ; j++){
      
        for (int k=0 ; k < NUMBER_OF_PMTS ; k++){
 
          if (i -> second[j][k].getCount() > 0){
 
            (*histograms.h_summary.m_mean_summary_rate)             [MAKE_STRING("Detector/DU" + to_string(i->first))]->Fill(k, j+1, i->second[j][k].getMean());
            (*histograms.h_summary.m_mean_summary_rate_distribution)[MAKE_STRING("Detector/DU" + to_string(i->first))]->Fill(        i->second[j][k].getMean());
          }
        }
      }
    }
  }
 
  
  /*
   * Function template to read and process time ordered timeslices.
   * 
   * \param scanner A JTreeScanner from which time ordered timeslices are accessed. 
   * \param frame_index_range the range of the indices of the frames to be accessed.
   */
  template <class T>
  void iterateTimesliceTree (JTreeScanner < T , JDAQEvaluator > & 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) {
          
        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 < vector < int > >  hit_time_buffers ; 
 
        hit_time_buffers.resize(NUMBER_OF_PMTS);
 
        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()) ;
 
          hit_time_buffers [hit->getPMT()].push_back (hit->getT()) ;
          
          if (hit->getToT() == 255) {
 
            (*histograms.h_timeslice.m_ToT_255[ts_type])[MAKE_STRING("Detector/DU" + to_string(string))] -> Fill(floor , hit->getPMT()) ;
 
            histograms.h_timeslice.h_ToT_255_vs_time        [ts_type] -> Fill(hit -> getT());
            histograms.h_timeslice.h_ToT_255_Floor_vs_time  [ts_type] -> Fill(hit -> getT() , floor);
            histograms.h_timeslice.h_ToT_255_Floor_vs_time_2[ts_type] -> Fill(hit -> getT() , floor);       
          }
        }
        
        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 * hit_time_buffers[pmt].size() * inverseFrameTimeSec);
                  
          TH2D* h2 = (*histograms.h_timeslice.m_pmt_dt_consecutive_hits[ts_type])[MAKE_STRING("Detector/DU" + to_string(string) + "/F" + to_string(floor))];
 
          if (hit_time_buffers[pmt].size() > 1) {
            for (vector<int>::const_iterator q = hit_time_buffers[pmt].begin(), p = q++; q != hit_time_buffers[pmt].end(); ++p, ++q) {
              h2 -> Fill ( log10(*q - *p) , pmt);
            }
          }
        }
      }      
    }
 
    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, JDAQEvaluator> 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 , JDAQEvaluator> 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 , JDAQEvaluator> scanner(inputFile, numberOfEvents);
        
    if(scanner.hasNext()) {
 
      histograms.initialize_trigger_histograms();
      scanner.rewind();
      iterateEventTree(scanner);
    }
  }
 
  /*
   * Returns the histograms produced from the data in the file.
   */
  JRA_Histograms getHistograms(){
 
    return histograms ;
  }
};
 
#endif