NBPulse.hh

Go to the documentation of this file.

#ifndef __NBPULSE__
#define __NBPULSE__


#include "RooRealVar.h"
#include "RooWorkspace.h"
#include "RooGaussModel.h"
#include "RooArgList.h"
#include "RooLandau.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooAddPdf.h"
#include "RooGaussian.h"
#include "RooMsgService.h"

#include "TH1D.h"
#include "TH2D.h"
#include "TRandom3.h"



/**
 * \author rgruiz
 */

using namespace RooFit;
using namespace std;


/**
 * Analyzes the signal of a nanobeacon in a PMT
 */
class NBPulse {

  TH2D* time_vs_tot ;

  TH1D *htime_full , *htime_peak , *tot_peak ;

  RooWorkspace w ;

  pair<double , double> baseline ;

  double tlow , thigh , tpeak , nhits , meanToT , sigmaToT , max_ToT ;
  
  int n_bins_right , n_bins_left , n_bins_before_peak , n_bins_afterpulse , max_n_bins_baseline , min_n_bins_baseline , n_sigma_peak , n_sigma_before_peak  ; 

  double fit_region_min , fit_region_max , fit_range_min , fit_range_max , landau_m_min , landau_m_max , landau_s_min , landau_s_max , gaus_m_min , gaus_m_max , gaus_s_min , gaus_s_max , l_g_ratio_min , l_g_ratio_max , KS_threshold;

  bool analyzed  ;

  bool fitted  ;

  bool is_good  ;

  bool saturated_hit  ;

  bool saturated_tot  ; 
  
  bool weak  ;

 public :

  /**
   * Default constructor.
   */
  NBPulse(){}


  /**
   * Constructor.
   * This constructor initializes the nanobeacon pulse from a 2D histogram containing the 
   * time vs ToT hit distribution. 
   *
   * \param  h    time vs ToT histogram
   */
  NBPulse(TH2D* h) :
    max_ToT (26.4) ,
    n_bins_right (5) , 
    n_bins_left  (5) ,
    n_bins_before_peak(2) ,
    n_bins_afterpulse (10),
    max_n_bins_baseline (50),
    min_n_bins_baseline (10),
    n_sigma_peak (5) ,
    n_sigma_before_peak (4) ,
    fit_region_min (20.0), 
    fit_region_max (20.0) , 
    fit_range_min (6.0) , 
    fit_range_max (15.0) , 
    landau_m_min (2.0) , 
    landau_m_max (2.0) , 
    landau_s_min (1.0) , 
    landau_s_max (9.0) , 
    gaus_m_min (2.0) , 
    gaus_m_max (2.0) , 
    gaus_s_min (1.0) , 
    gaus_s_max (9.0) , 
    l_g_ratio_min (0.0) , 
    l_g_ratio_max (1.0) , 
    KS_threshold(1e-3) ,
    analyzed (false) ,
    fitted (false) ,
    is_good (false) ,
    saturated_hit (false) ,
    saturated_tot (false) , 
    weak (false) 

  {

    time_vs_tot = h ;

    initialize() ;

  }


  /**
   * Destructor.
   */
  ~NBPulse(){

    delete(htime_full) ;

    if (analyzed == true){

      delete(tot_peak) ; 

      delete(htime_peak) ;

    }

  } 


  /**
   * Performs a Kolmogorov-Smirnov compatibility test between the hit time distribution and a uniform distribution
   *
   * \param h The hit time distribution
   * \return The probability that the hit time distribution is compatible with a uniform distribution. 
   */
  double Kolmogorov_Smirnov_uniform(TH1D* h){

    TRandom3* r = new TRandom3() ;

    double low_x = 0 ;

    double high_x = h->GetBinCenter(h->GetNbinsX() - 1) ;

    int nentries = h->Integral(1 , h->GetNbinsX() - 1) ;

    if (h->GetBinCenter(1)>=0){
      
      low_x = h->GetBinCenter(1) ;

    }

    TH1D* u = (TH1D*)h->Clone("uniform") ;

    u->Reset() ;

    for (int i = 0 ; i<nentries ; i++){

      u->Fill(r->Uniform(low_x, high_x));

    }
    
    double proba = h->KolmogorovTest(u) ;

    delete(r) ;

    delete(u) ;

    return proba ;

  } ;




  /**
   * Leaves the pulse ready for analysis by projecting the time vs ToT distribution to 
   * the time axis and moving it out of the Root file directory.
   */
  void initialize(){
    
    htime_full = time_vs_tot->ProjectionX(time_vs_tot->GetName() , 1 , time_vs_tot->GetNbinsY());
  
    htime_full->SetDirectory(0) ;

    RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
  
  }


  /**
   * Performs a fast analysis of the hit time distribution.
   */ 
  void analyzeFast(){

    if (analyzed == true) return ;

    computeBaseline(htime_full) ;
    
    int peakbin = htime_full->GetMaximumBin() ;

    double height = htime_full->GetBinContent(peakbin) ;

    int leftbin = peakbin - n_bins_left ;

    int rightbin = peakbin + n_bins_right ;

    tpeak = htime_full->GetBinCenter(peakbin) ;

    tlow = htime_full->GetBinCenter(leftbin) ;

    thigh = htime_full->GetBinCenter(rightbin) ;

    htime_peak = (TH1D*)htime_full->Clone();

    for (int i=0 ; i<htime_peak->GetNbinsX() ; i++){

      if (i<leftbin || i>rightbin){

        htime_peak->SetBinContent(i,0);

      }

    }

    htime_peak->SetName((std::string("htime_")+htime_peak->GetName()).c_str());

    htime_peak->SetDirectory(0) ;
  
    tot_peak = (TH1D*)time_vs_tot->ProjectionY((std::string("ToT_") + time_vs_tot->GetName()).c_str() , leftbin , rightbin) ; 

    tot_peak->SetDirectory(0);

    meanToT = tot_peak->GetMean() ;
  
    if (meanToT > max_ToT){
    
      saturated_tot = true ;
    
    }

    sigmaToT = tot_peak->GetStdDev() ;

    double p_KS = Kolmogorov_Smirnov_uniform(htime_full) ;

    if (baseline.first==0 && baseline.second==0) weak = true ;

    if (p_KS > KS_threshold) weak = true ;

    if (height < (baseline.first + n_sigma_peak * baseline.second)) weak = true ;

    for (int i = 1 ; i<=n_bins_before_peak ; i++){

      if (htime_full->GetBinContent(peakbin - n_bins_before_peak) < (baseline.first + n_sigma_before_peak * baseline.second)) weak = true ;

    }

    if (weak == false && htime_full->GetBinContent(peakbin + n_bins_afterpulse) == 0) saturated_hit = true ;

    if (weak == false && saturated_hit == false) is_good = true ;

    nhits = htime_peak->Integral() ;

    analyzed = true ;

  }



  /**
   * Computes baseline of the hit time distribution measured by the PMT. For this the first bins of the distribution (where the nanobeacon pulse is absent) are used.
   *
   * \param h Histogram with the hit time distribution
   */
  void computeBaseline(TH1* h){

    int N = 0;

    double sum = 0 ;

    vector<double> bins ;

    for (int i=1 ; i<max_n_bins_baseline + 1 ; i++){

      if(h->GetBinContent(i) > 0 ){

        sum += h->GetBinContent(i);

        bins.push_back(h->GetBinContent(i));

        N += 1 ;

      }

    }
 
    if(N < min_n_bins_baseline){

      bins.clear();

      sum = 0 ;

      N = 0 ;

      for (int i=1 ; i<max_n_bins_baseline + 1 ; i++){

        if(h->GetBinContent(h->GetNbinsX() - i)>0){

          sum += h->GetBinContent(h->GetNbinsX() - i);

          bins.push_back(h->GetBinContent(h->GetNbinsX() - i));

          N += 1 ;

        }

      }

    }

    if(N < min_n_bins_baseline){

      baseline.first = 0 ;

      baseline.second = 0 ;

      return ;

    }

    double bg_mean = sum/N ;

    baseline.first = bg_mean ;

    double sdev = 0;

    for (int i=0 ; i<N ; i++){
    
      sdev += pow((bins[i] - bg_mean) , 2) ;
    
    }

    baseline.second = sqrt(sdev/N);

  }


  /**
   * Fits the hit nanobeacon signal to a model composed by a Landau and a Gauss functions.
   */
  void fit(){

    if (fitted==true)   return ;

    if (analyzed==false) analyzeFast() ;

    w.SetName((std::string("W_") + htime_full->GetName()).c_str());

    double tpeak = htime_full->GetBinCenter(htime_full->GetMaximumBin()) ;

    RooRealVar x("time" , "hit time" , tpeak - fit_region_min , tpeak + fit_region_max) ;

    RooRealVar ml("mean landau" , "mean landau" , tpeak - landau_m_min , tpeak + landau_m_max);

    RooRealVar sl("sigma landau" , "sigma landau" , landau_s_min , landau_s_max);
  
    RooLandau L1("L1" , "Main_Landau" , x , ml , sl) ;
    
    RooRealVar mg("mean Gauss" , "mean Gauss" , tpeak - gaus_m_min , tpeak + gaus_m_max);

    RooRealVar sg("sigma Gauss" , "sigma Gauss" , gaus_s_min , gaus_s_max);

    RooGaussian G1("G1" , "Main_Gaussian" , x , mg , sg) ;
  
    RooRealVar c1("c1" , "c1" , l_g_ratio_min , l_g_ratio_max);

    RooAddPdf M("Model" , "Model" , RooArgList(G1,L1) , RooArgList(c1));

    w.import(M) ;

    RooDataHist Data ("data" , "data" , *w.var("time") , Import(*htime_full));
  
    w.pdf("Model")->fitTo(Data , Range(tpeak - fit_range_min  , tpeak + fit_range_max) , RooFit::PrintLevel(-1000) , RooFit::PrintEvalErrors(-1000)) ;

    RooRealVar Tmin ("tmin" , "tmin" , tpeak - fit_range_min) ;
  
    RooRealVar Tmax ("tmax" , "tmax" , tpeak + fit_range_max) ;
  
    RooRealVar Tpeak ("tpeak" , "tpeak" , tpeak) ;

    RooRealVar Mpv("mpv" , "mpv" , Tmin.getVal()) ;

    RooRealVar Diff("peak-mpv" , "peak-mpv" , Tpeak.getVal() - Tmin.getVal()) ;

    double max_pdfval = 0 ;

    RooRealVar *var = (RooRealVar*)w.function("Model")->getObservables(RooArgSet(*w.var("time")))->first();

    for(double x = Tmin.getVal() ; x<Tmax.getVal() ; x+=0.1){

      var->setVal(x) ;

      double pdfval = w.function("Model")->getVal(RooArgSet(*var)) ;

      if (pdfval>max_pdfval){

        max_pdfval = pdfval ;

        Mpv.setVal(x) ;

      }

    }

    w.import(Mpv) ;
  
    w.import(Tmin) ;

    w.import(Tmax) ;

    w.import(Tpeak) ;

    w.import(Diff) ;

    fitted = true ;

  }

 
  /**
   * Check if the peak is good.
   * 
   * \result true if the peak is classified as good.
   */ 
  bool IsGood() {

    return is_good ;

  }


  /**
   * Check if the peak is fitted.
   * 
   * \result true if the peak has been fitted.
   */ 
  bool IsFitted() {

    return fitted ;

  }


  /**
   * Check if the hit time distribution is saturated.
   * 
   * \result true if the hit time distributon is classified as saturated.
   */ 
  bool IsSaturatedHit(){
    
    return saturated_hit ;
  
  }


  /**
   * Check if the hit time distribution is weak.
   * 
   * \result true if the hit time distributon is classified as weak.
   */
  bool IsWeak(){

    return weak ;
    
  }


 
  /**
   * Get the hit time distribution.
   *
   * \return histogram with hit time distribution
   */ 
  TH1D* getHtime_full(){
    
    return htime_full ;
    
  } 



  /**
   * Get the hit time distribution.
   *
   * \return histogram with hit time distribution
   */ 
  TH1D* gettime_peak(){
    
    return htime_peak ;
    
  }
  

  
  /**
   * Get the hit time distribution.
   *
   * \return histogram with hit time distribution
   */ 
  TH1D* gettot_peak(){
    
    return tot_peak ;
    
  }




  /**
   * Get the hit time vs ToT distribution.
   *
   * \return  a 2d histogram with hit time vs ToT distribution
   */ 
  TH2D* gettime_vs_tot(){

    return  time_vs_tot ;
    
  }
  

 
  /**
   * Get all the information from the fit.
   *
   * \return Roofit Workspace with all the information and parameters from the fit 
   */ 
  RooWorkspace getWorkspace(){

    return w ;
    
  }
  

 
  /**
   * Get the arrival time of the nanobeacon pulse.
   *
   * \return pulse arrival time
   */
  double getPeakTime() {

    return tpeak ;
    
  }
  

 
  /**
   * Get the mean of the ToT distribution of the hits produced by the nanobeacon.
   *
   * \return mean ToT value
   */
  double getMeanToT() {
    
    return meanToT ;
    
  }
  

  
  /**
   * Get the standard deviation of the ToT distribution of the hits produced by the nanobeacon.
   *
   * \return standard deviation from ToT distribution
   */
  double getSigmaToT() {

    return sigmaToT ;  
    
  }
  

  
  /**
   * Get name of the hit time vs ToT histogram.
   *
   * \return The name of the histogram.
   */
  string getName() {
    
    return time_vs_tot->GetName();
    
  }
  

};



#endif