JPulseFinder_utils.hh

Go to the documentation of this file.

#ifndef __PULSE_UTILS__
#define __PULSE_UTILS__


// c++ standard library
#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
 
// JPP
#include "JDAQ/JDAQ.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JModuleRouter.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitToolkit.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JTrigger/JTimesliceRouter.hh"

// root
#include "TH1D.h"
#include "TH2D.h"

// namespaces
using namespace std;
using namespace KM3NETDAQ ; // for JDAQTimeSlice
using namespace JSUPPORT;   // for JFileScanner
using namespace JLANG;   // for JMultipleFileScanner
using namespace JDETECTOR ;

using namespace JPP ;

/**
 * \author rgruiz
 */



/**
 * Calculates an offset for a combinataion nanobeacon-pmt that allows to define the histogram range in a time window containing the corresponding pulse.
 *
 * \param pmt The target pmt
 * \param NB The source module
 * \param c_w The speed light in water
 * \param options user options containing the stagger and pulse delay information
 * \return The corresponding offset to apply in the definition of the hit time histogram range. 
 */
inline int get_offset(JPMT pmt , JModule NB , double c_w , IO options){

  double distance_pmt_nb = (pmt.getPosition() - NB.getPMT(0).getPosition()).getLength() ;

  double time_of_flight = distance_pmt_nb / c_w ;

  double delta_t0 = pmt.getT0() - NB.getPMT(0).getT0() ;

  double delta_traw = time_of_flight - delta_t0 ;

  double stagger_delay = 16 * options.stagger_16ns * (NB.getFloor() - 1) ;

  double overall_delay = 16 * options.overall_delay_16ns ; 

  int total_offset = (int)round( delta_traw + stagger_delay + overall_delay);
  
  return total_offset ;

}



/**
 * Allocates histograms to store the tot vs time hit distribution for all the combinations nanobeacon - pmt.
 *
 * \param detector The detector containing the DU to be analyzed
 * \param options user options containing the string number and the list of used nanobeacons
 * \param c_w The speed light in water
 * \param top_pmts List of FPGA channels for pmts in the upper hemisphere
 * \param bottom_pmts List of FPGA channels for pmts in the lower hemisphere.
 * \param max_distance Maximum distance between modules
 * \return A table of TH2D histograms defined in the correct hit time ranges. 
 */
inline vector < vector < TH2D* > > allocate_histograms (JDetector detector , IO options , double c_w , vector<int> top_pmts , vector<int> bottom_pmts , int max_distance){
  
  int nFloors = getNumberOfFloors(detector) ;

  int nHistosPerFloor = top_pmts.size() + max_distance * bottom_pmts.size() ;

  vector < vector < TH2D* > > histograms (nFloors , vector<TH2D*> (nHistosPerFloor, NULL)) ;
  
  char name[250] ;
  
  for (int tgt_floor = 1 ; tgt_floor <= nFloors ; ++tgt_floor){

    JModule tgt_module = detector.getModule( JModuleLocation(options.string_number , tgt_floor) ) ;

    int index = 0 ;

    for(auto & top_pmt : top_pmts){

      int total_offset = get_offset(tgt_module.getPMT(top_pmt) , tgt_module , c_w , options) ;

      sprintf(name , "S%iF%iPMT%i_beaconS%iF%i" , tgt_module.getString() , tgt_module.getFloor() , top_pmt , tgt_module.getString() , tgt_module.getFloor());
      
      histograms[tgt_floor - 1][index] = new TH2D( name , NULL , 400 , total_offset - 0.5 , 400 + total_offset - 0.5 , 256 , -0.5 , 255.5 ) ;

      histograms[tgt_floor - 1][index]->SetOption("colz") ;
      
      index += 1 ;
    }

    
    for (int src_floor = tgt_floor - max_distance ; src_floor < tgt_floor ; ++src_floor){

      if(src_floor > 0){
        
        JModule src_module = detector.getModule( JModuleLocation (options.string_number , src_floor) ) ;
        
        for (auto & bottom_pmt : bottom_pmts){
          
          int total_offset = get_offset(tgt_module.getPMT(bottom_pmt) , src_module , c_w , options) ;
          
          sprintf(name , "S%iF%iPMT%i_beaconS%iF%i" , tgt_module.getString() , tgt_module.getFloor() , bottom_pmt , src_module.getString() , src_module.getFloor());
          
          histograms[tgt_floor - 1][index] = new TH2D( name , NULL , 400 , total_offset - 0.5 , 400 + total_offset - 0.5 , 256 , -0.5 , 255.5 ) ;
          
          histograms[tgt_floor - 1][index]->SetOption("colz") ;
          
          index += 1 ;
          
        }
        
      }
      
    }
    
  }

  return histograms ;

}




/**
 * Fills the hit time vs ToT histograms for all the combinations nanobeacon - pmt.
 *
 * \param histograms The vector of histograms to be filled 
 * \param options User options containing the string number and the list files to be read
 * \param detector The detector containing the DU to be analyzed
 * \param top_pmts List of FPGA channels for pmts in the upper hemisphere
 * \param bottom_pmts List of FPGA channels for pmts in the lower hemisphere.
 * \param max_distance Maximum distance between modules
 */
inline void fill_histograms (vector < vector < TH2D* > > &histograms , IO options , JDetector detector , vector<int> top_pmts , vector<int> bottom_pmts , int max_distance){
    
  JTriggerParameters parameters ;

  unsigned int pulse_period = 16 * options.pulse_period_16ns ;

  JMultipleFileScanner<JDAQTimesliceL0> inputFile(options.ifnames) ;

  JModuleRouter moduleRouter(detector) ;

  int slicecounter = 0 ;

  while( inputFile.hasNext() ) {

    slicecounter ++ ;
    
    JDAQTimesliceL0 slice = *(inputFile.next()) ;

    JTimesliceRouter timesliceRouter(parameters.numberOfBins) ;

    timesliceRouter.configure(slice);

    for (int i = 0 ; i < getNumberOfFloors(detector) ; i++){
      
      const JModule& module = detector.getModule(JModuleLocation(options.string_number , i+1)) ;
      
      if(timesliceRouter.hasSuperFrame(module.getID())){

        JDAQSuperFrame SuperFrame = timesliceRouter.getSuperFrame(module.getID()) ;

        vector<JDAQHit> buffer[getNumberOfPMTs(module)];

        for (JDAQSuperFrame::const_iterator hit = SuperFrame.begin(); hit != SuperFrame.end(); ++hit) {
           
          buffer[hit->getPMT()].push_back(*hit) ;
         
        }
         
        int i = 0 ;

        for (auto & top_pmt : top_pmts){
         
          for (int hit = 0 ; hit < (int)buffer[top_pmt].size() - 1 ; hit++){

            double tot = buffer[top_pmt][hit].getToT() ;

            double raw_time = buffer[top_pmt][hit].getT() ;

            //double t0 = module.getPMT(buffer[top_pmt][hit].getPMT()).getT0() ;

            //double calibrated_time = raw_time + t0 ; 
             
            double mod_double = raw_time - pulse_period * int(raw_time/pulse_period) ;

            histograms[module.getFloor()-1][i]->Fill(mod_double , tot) ;

          }

          i++ ;

        }

         
        i = 0 ;

        for (auto & bottom_pmt : bottom_pmts){
         
          for (int hit = 0 ; hit < (int)buffer[bottom_pmt].size() - 1 ; hit++){
             
            double tot = buffer[bottom_pmt][hit].getToT() ;
             
            double raw_time = buffer[bottom_pmt][hit].getT() ;
             
            //double t0 = module.getPMT(buffer[bottom_pmt][hit].getPMT()).getT0() ;
             
            //double calibrated_time = raw_time + t0 ; 
             
            double mod_double = raw_time - pulse_period * int(raw_time/pulse_period) ;
             
            for (int j = 1 ; j <= max_distance ; j++){
               
              if (module.getFloor() - j <=0) continue ;
               
              histograms[module.getFloor() - 1][top_pmts.size() + (j-1)*bottom_pmts.size() + i]->Fill(mod_double , tot) ;
               
            }
             
          }
           
          i++ ;
         
        }
 
      } // end of loop over JDAQTimeSlices
      
    }
    
  }

}





/**
 * Projects the hit time vs ToT histograms for all the combinations nanobeacon - pmt on the Y axis.
 *
 * \param ToT_vs_time The vector of histograms to be projected 
 * \return A vector with ToT distribiutions. 
 */
vector < vector < TH1D* > > project_ToT(vector < vector < TH2D* > > ToT_vs_time){

  vector < vector < TH1D* > > ToT (ToT_vs_time.size() , vector < TH1D* > (ToT_vs_time[0].size() , NULL)) ;

  for (int i=0 ; i<(int)ToT_vs_time.size() ; i++){

    for (int j=0 ; j<(int)ToT_vs_time[0].size() ; j++){

      if (ToT_vs_time[i][j]==NULL) continue ;

      TH1D* p_y = ToT_vs_time[i][j]->ProjectionY((string("ToT_") + ToT_vs_time[i][j]->GetName()).c_str() , 1 , ToT_vs_time[i][j]->GetNbinsX()-1) ;
      
      ToT[i][j] = p_y ;
      
    }
    
  }

  return ToT ;

}



/**
 * Projects the hit time vs ToT histograms for all the combinations nanobeacon - pmt on the X axis.
 *
 * \param ToT_vs_time The vector of histograms to be projected 
 * \return A vector with hit time distribiutions. 
 */
inline vector < vector < TH1D* > > project_time(vector < vector < TH2D* > > ToT_vs_time){

  vector < vector < TH1D* > > time (ToT_vs_time.size() , vector < TH1D* > (ToT_vs_time[0].size() , NULL)) ;

  for (int i=0 ; i<(int)ToT_vs_time.size() ; i++){

    for (int j=0 ; j<(int)ToT_vs_time[0].size() ; j++){

      if (ToT_vs_time[i][j]==NULL) continue ;
      
      TH1D* p_x = ToT_vs_time[i][j]->ProjectionX((string("htime_") + ToT_vs_time[i][j]->GetName()).c_str() , 1 , ToT_vs_time[i][j]->GetNbinsY()-1) ;
      
      time[i][j] = p_x ;
      
    }
    
  }

  return time ;

}


#endif