JRunAnalyzer Class Reference

Class dedicated to the analysis of KM3NeT runs. More...

#include <JRunAnalyzer.hh>

Classes
struct	dom_type

Public Member Functions
	JRunAnalyzer (const JSingleFileScanner<> &file, const JDetector &detector, TFile *out, JLimit_t nTimeslices, JLimit_t nSummaryslices, JLimit_t nEvents, int analysislevel)
	Constructor. More...
	~JRunAnalyzer ()
	Destructor. More...
void	iterateEventTree (JTreeScanner< JDAQEvent > &scanner)
void	iterateRecoEventTree (JParallelFileScanner< JTypeList< JDAQEvent, JFIT::JEvt > > &scanner, TFile &f)
void	iterateSummarysliceTree (JTreeScanner< JDAQSummaryslice > &scanner)
template<class T >
void	iterateTimesliceTree (JTreeScanner< T > &scanner)
void	readSummaryData ()
template<class T >
void	readTimesliceData ()
void	readEvents ()
void	readRecoEvents ()
void	writeToFile (TFile *f, int analysis_level)

Private Attributes
JSingleFileScanner	inputFile
JModuleRouter	router
TFile *	outputFile
JRA_Histograms	histograms
int	highest_floor
JLimit_t	numberOfTimeslices
JLimit_t	numberOfSummaryslices
JLimit_t	numberOfEvents
int	analysis_level

Detailed Description

Class dedicated to the analysis of KM3NeT runs.

Definition at line 32 of file JRunAnalyzer.hh.

Constructor & Destructor Documentation

JRunAnalyzer::JRunAnalyzer ( const JSingleFileScanner<> &  file, const JDetector &  detector, TFile *  out, JLimit_t  nTimeslices, JLimit_t  nSummaryslices, JLimit_t  nEvents, int  analysislevel  )

inline

Constructor.

Parameters

file	file name
detector	detector
out	pointer to output file
nTimeslices	number of timeslices to be read
nSummaryslices	number of summary slices to be read
nEvents	number of events to be read
analysislevel	option for analysis of PMT or reco data

Definition at line 57 of file JRunAnalyzer.hh.

                                                                                                   :
    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);
    highest_floor  = getRangeOfFloors(detector).getUpperLimit();
  } 

JRunAnalyzer::~JRunAnalyzer ( )

inline

Destructor.

Definition at line 77 of file JRunAnalyzer.hh.

{}

Member Function Documentation

void JRunAnalyzer::iterateEventTree ( JTreeScanner< JDAQEvent > &  scanner )

inline

Definition at line 85 of file JRunAnalyzer.hh.

                                                               {

    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());
        }
      }
    }
  }

void JRunAnalyzer::iterateRecoEventTree ( JParallelFileScanner< JTypeList< JDAQEvent, JFIT::JEvt > > &  scanner, TFile &  f  )

inline

Definition at line 186 of file JRunAnalyzer.hh.

                                                                                                          {

    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();
  }

void JRunAnalyzer::iterateSummarysliceTree ( JTreeScanner< JDAQSummaryslice > &  scanner )

inline

Definition at line 257 of file JRunAnalyzer.hh.

                                                                             {

    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(auto& i : DOM_rate_quantiles){
      
      for (int j=1; j <= highest_floor; j++) {
      
        if (i.second[j].getCount() > 0) histograms.h_summary.h_rate_summary->Fill(i.first, j, i.second[j].getMean());
      }
    }

    if(analysis_level == 1){

      for(auto& i : PMT_rate_quantiles){

        for (int j = 1; j <= highest_floor; 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, 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());
            }
          }
        }
      }
    }
  }

template<class T >

void JRunAnalyzer::iterateTimesliceTree ( JTreeScanner< T > &  scanner )

inline

Definition at line 372 of file JRunAnalyzer.hh.

                                                           {
    
    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() ) ;
      }
    }
  }

void JRunAnalyzer::readSummaryData ( )

inline

Definition at line 436 of file JRunAnalyzer.hh.

                         {

    JTreeScanner <JDAQSummaryslice> scanner(inputFile, numberOfSummaryslices);
    if (scanner.hasNext()) {
      histograms.initialize_summary_histograms();
      scanner.rewind();
      iterateSummarysliceTree(scanner);
    }
  }

template<class T >

void JRunAnalyzer::readTimesliceData ( )

inline

Definition at line 452 of file JRunAnalyzer.hh.

                           {

    JTreeScanner <T> scanner(inputFile, numberOfTimeslices);
    if(scanner.hasNext()) {
      
      histograms.initialize_timeslice_histograms <T>();
      scanner.rewind();
      iterateTimesliceTree(scanner);
    }
  }

void JRunAnalyzer::readEvents ( )

inline

Definition at line 468 of file JRunAnalyzer.hh.

                    {

    JTreeScanner <JDAQEvent> scanner(inputFile, numberOfEvents);
        
    if(scanner.hasNext()) {

      histograms.initialize_trigger_histograms();
      scanner.rewind();
      iterateEventTree(scanner);
    }
  }

void JRunAnalyzer::readRecoEvents ( )

inline

Definition at line 484 of file JRunAnalyzer.hh.

                        {

    JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> >  scanner(inputFile, numberOfEvents);
        
    if(scanner.hasNext()) {

      scanner.rewind();
      iterateRecoEventTree(scanner, *outputFile);
    }
  }

void JRunAnalyzer::writeToFile ( TFile *  f, int  analysis_level  )

inline

Definition at line 499 of file JRunAnalyzer.hh.

                                               {

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

Member Data Documentation

JSingleFileScanner JRunAnalyzer::inputFile

private

Definition at line 34 of file JRunAnalyzer.hh.

JModuleRouter JRunAnalyzer::router

private

Definition at line 35 of file JRunAnalyzer.hh.

TFile* JRunAnalyzer::outputFile

private

Definition at line 36 of file JRunAnalyzer.hh.

JRA_Histograms JRunAnalyzer::histograms

private

Definition at line 37 of file JRunAnalyzer.hh.

int JRunAnalyzer::highest_floor

private

Definition at line 38 of file JRunAnalyzer.hh.

JLimit_t JRunAnalyzer::numberOfTimeslices

private

Definition at line 39 of file JRunAnalyzer.hh.

JLimit_t JRunAnalyzer::numberOfSummaryslices

private

Definition at line 40 of file JRunAnalyzer.hh.

JLimit_t JRunAnalyzer::numberOfEvents

private

Definition at line 41 of file JRunAnalyzer.hh.

int JRunAnalyzer::analysis_level

private

Definition at line 42 of file JRunAnalyzer.hh.

---

The documentation for this class was generated from the following file:

• JRunAnalyzer.hh