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