ANALYSIS/FIN_AN.h

// @(#)fROOT/ANALYSIS:$Name:  $:$Id: FIN_AN.h,v 1.10 2007/09/24 07:32:41 Diego_Faso Exp $
// Author: Diego Faso <mailto:faso@to.infn.it>, 2007/05/10

#ifndef FROOT_FIN_AN
#define FROOT_FIN_AN

// The FINUDA analysys nemespace //

#include <TString.h>
#include <TThread.h>
#include <TMath.h>
#include "TSpectrum.h"
#include "TH1F.h"
#include "TF1.h"
#include "TROOT.h"
//#include <Riostream.h>

namespace FIN_AN {
  // --- enum and constants
  enum E_FinAn_NMaxPeaks { FinAn_Nmaxpeaks = 500 }; // software limit (could be useful)

  static const Int_t    K_FinAn_int    = 12345;
  static const Float_t  K_FinAn_float  = 12345.;
  static const Double_t K_FinAn_double = 12345.;

  // --- histograms information/manipulation
  static Int_t GetHistogramLimits(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy);
  //  static Int_t GetHistogramRanges(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy);

  // --- static methods used as fit functions
  //      (and functions attributes setting)
  static Double_t FitFCN_GausYield(Double_t *x, Double_t *par); // par[0,1,2]: yield,mean,sigma
  //
  static void SetFitFCN_ParNames_GausYield(TF1 *fcn);

  // --- histograms evaluation (background and peaks estimation)
  static TH1F *EvalBackground(TH1F *his,Int_t n_max_peaks=10,Int_t niter=40,Option_t* option = "same");

  static TSpectrum *FindPeaks(TH1F *his,Int_t n_max_peaks=10,Double_t sigma=0.02,Double_t threshold=0.25,Double_t resol_fact=1,Option_t* option="goff");

  static Int_t PrintFoundPeaks(TSpectrum *spec);
  static Int_t LocatePeak(TSpectrum *spec,Double_t low_lim,Double_t high_lim);

};



// --- histograms informationmanipulation
inline Int_t FIN_AN::GetHistogramLimits(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy){
  // required value are passed to arguments 
  // for reference
  //
  // return value:
  //              0: ok
  //             -1: error
  
  if(!his){
    gROOT->Warning("TFndGraph::GetHistogramLimits","Pointer to the given histogram is NULL!");
    return -1;
  }
  minx = static_cast<Double_t>( his->GetBinCenter(1) );
  maxx = static_cast<Double_t>( his->GetBinCenter(his->GetNbinsX()) + his->GetBinWidth(1) );
  miny = static_cast<Double_t>( his->GetMinimum() );
  maxy = static_cast<Double_t>( his->GetMaximum() );
  //
  /*   cout << " -> minx: " << minx << endl; */
  /*   cout << " -> maxx: " << maxx << endl; */
  /*   cout << " -> miny: " << miny << endl; */
  /*   cout << " -> maxy: " << maxy << endl; */
  //
  return 0;  
}

//inline Int_t FIN_AN::GetHistogramRanges(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy){
//  // required value are passed to arguments for reference
//  //
//  // return value:
//  //              0: ok
//  //             -1: error
//  
//  if(!his){
//    gROOT->Warning("TFndGraph::GetHistogramRanges","Pointer to the given histogram is NULL!");
//    return -1;
//  }
//  minx = static_cast<Double_t>( his->GetXaxis()->GetXmin() );
//  maxx = static_cast<Double_t>( his->GetXaxis()->GetXmax() );
//  miny = static_cast<Double_t>( his->GetYaxis()->GetXmin() );
//  maxy = static_cast<Double_t>( his->GetYaxis()->GetXmax() );
//  //
//  /*   cout << " -> minx: " << minx << endl; */
//  /*   cout << " -> maxx: " << maxx << endl; */
//  /*   cout << " -> miny: " << miny << endl; */
//  /*   cout << " -> maxy: " << maxy << endl; */
//  //
//  return 0;  
//}

// --- static methods used as fit functions
//      (and functions attributes setting)

//_______________________
inline Double_t FIN_AN::FitFCN_GausYield(Double_t *x, Double_t *par){
  // simple gaussian peak with the yield as parameter
  //                       (instead of normalization)
  // par[0,1,2]: yield,mean,sigma
  //
  // -------- IMPORTANT NOTE --------
  //   Remember that the number of events under the
  //   fitted gaussian peak (area) depends on the bin-width:
  //              Area = yield / bin-width
  //
  //     remember (confidence levels):
  //               +- 1 * sigma ===> 68.30 % of the total area
  //               +- 2 * sigma ===> 95.45 % of the total area
  //               +- 3 * sigma ===> 99.70 % of the total area
  // --------------------------------

  Double_t result = 0;
  
  Double_t DenNorm = (2*TMath::Pi())*TMath::Power(par[2],2);
  Double_t Norm =  par[0]/TMath::Sqrt(DenNorm); // par[0]: yield
  Double_t N1 = TMath::Power((x[0]-par[1]),2);
  Double_t N2 = 2 * TMath::Power(par[2],2);
  
  result = Norm * TMath::Exp(-N1/N2);
  
  return result;
}

//_______________________
inline void FIN_AN::SetFitFCN_ParNames_GausYield(TF1 *fcn){

  fcn->SetParName(0,"yield");
  fcn->SetParName(1,"mean");
  fcn->SetParName(2,"sigma");

  fcn->SetLineWidth(2);
  fcn->SetLineColor(4);
}

// --- histograms evaluation (background and peaks estimation)

//_______________________
inline TH1F *FIN_AN::EvalBackground(TH1F *his,Int_t n_max_peaks,Int_t niter,Option_t* option){
  
  if(!his){
    gROOT->Warning("FIN_AN::EvalBackground","can not call method for NULL histogram.");
    return 0;
  }

  TSpectrum TmpSpec(n_max_peaks);
  return (TH1F *) ( TmpSpec.Background(his,niter,option) );
  
}

//_______________________
inline TSpectrum *FIN_AN::FindPeaks(TH1F *his,Int_t n_max_peaks,Double_t sigma,Double_t threshold,Double_t resol_fact,Option_t* option){
  // peaks and markers are added into the his->ListOfFinctions
  // return value: TSpectrum used
  //
  // threshold:  [ 0<threshold<1 ]  (default: 0.25):
  // peaks with amplitude less than threshold*highest_peak are discarded.              
  //
  // resol_fact [resolution factor] (default: 1):
  //              determines resolution of the neighboring peaks
  //              default value is 1 correspond to 3 sigma distance
  //              between peaks. Higher values allow higher resolution
  //              (smaller distance between peaks.
  //
  // option: (default:"")
  //  if option is not equal to "goff" , then          
  //   a polymarker object is created and added to the list of functions of  
  //   the histogram. The histogram is drawn with the specified option and   
  //   the polymarker object drawn on top of the histogram.                  
  //   The polymarker coordinates correspond to the npeaks peaks found in    
  //   the histogram.                                                        
  //   A pointer to the polymarker object can be retrieved later via:        
  //    TList *functions = hin->GetListOfFunctions();                        
  //    TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker")                                                        
  if(!his){
    gROOT->Warning("FIN_AN::FindPeaks","can not call method for NULL histogram.");
    return 0;
  }
  
  // Printf("Finding peaks for \"%s\"",his->GetName());
               
  TSpectrum *TmpSpectrum = new TSpectrum(n_max_peaks);
  TmpSpectrum->SetResolution(resol_fact);
  Int_t nfound = TmpSpectrum->Search(his,sigma,option,threshold);
  if(nfound > FinAn_Nmaxpeaks) Printf("too many peaks found: setting to %d",nfound);
  return TmpSpectrum;
}

//_____________________
inline Int_t FIN_AN::PrintFoundPeaks(TSpectrum *spec){
  // Print all peaks found by the given "TSpectrum"
  //
  // return value:
  //              >=0: number of peaks found
  //               -1: errors

  if(!spec){
    gROOT->Warning("FIN_AN::PrintFoundPeaks","can not call method for NULL spectrum.");
    return -1;
  }

  Int_t n_found_peaks = spec->GetNPeaks();
  TString msg = "\nPeaks evaluated position: \n";
  TString linestr ="";
  for(Int_t pi=0;pi<n_found_peaks;pi++){
    linestr.Form("Peak number %d: position: %.4f ; height: %.4f \n",
                 pi+1,
                 (*(spec->GetPositionX() + pi)),
                 (*(spec->GetPositionY() + pi))
                 );
    msg+=linestr;
  }
  Printf(msg.Data());
  return n_found_peaks;
}

//_____________________
inline Int_t FIN_AN::LocatePeak(TSpectrum *spec,Double_t low_lim,Double_t high_lim){
  // Print all peaks found by the given "TSpectrum"
  //
  // return value:
  //               > 0: number of the searched peak (starting from 1)
  //                -1: peak not found
  //                -2: error

  if(!spec){
    gROOT->Warning("FIN_AN::LocatePeak","can not call method for NULL spectrum.");
    return -2;
  }

  Int_t n_found_peaks = spec->GetNPeaks();
  for(Int_t pi=0;pi<n_found_peaks;pi++){
    Double_t cval = (*(spec->GetPositionX() + pi));
    if( cval >= low_lim && cval <= high_lim){
      //      Printf("Peak found: %.3f",cval);
      return pi+1;
    }
  }
  Printf("FIN_AN::LocatePeak---> Peak not found in range [%.4f ; %.4f]",low_lim,high_lim);
  return -1;
}

#endif // FROOT_FIN_AN

---