ANALYSIS/FIN_AN.h

00001 // @(#)fROOT/ANALYSIS:$Name:  $:$Id: FIN_AN.h,v 1.10 2007/09/24 07:32:41 Diego_Faso Exp $
00002 // Author: Diego Faso <mailto:faso@to.infn.it>, 2007/05/10
00003 
00004 #ifndef FROOT_FIN_AN
00005 #define FROOT_FIN_AN
00006 
00008 // The FINUDA analysys nemespace //
00010 
00011 #include <TString.h>
00012 #include <TThread.h>
00013 #include <TMath.h>
00014 #include "TSpectrum.h"
00015 #include "TH1F.h"
00016 #include "TF1.h"
00017 #include "TROOT.h"
00018 //#include <Riostream.h>
00019 
00020 namespace FIN_AN {
00021   // --- enum and constants
00022   enum E_FinAn_NMaxPeaks { FinAn_Nmaxpeaks = 500 }; // software limit (could be useful)
00023 
00024   static const Int_t    K_FinAn_int    = 12345;
00025   static const Float_t  K_FinAn_float  = 12345.;
00026   static const Double_t K_FinAn_double = 12345.;
00027 
00028   // --- histograms information/manipulation
00029   static Int_t GetHistogramLimits(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy);
00030   //  static Int_t GetHistogramRanges(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy);
00031 
00032   // --- static methods used as fit functions
00033   //      (and functions attributes setting)
00034   static Double_t FitFCN_GausYield(Double_t *x, Double_t *par); // par[0,1,2]: yield,mean,sigma
00035   //
00036   static void SetFitFCN_ParNames_GausYield(TF1 *fcn);
00037 
00038   // --- histograms evaluation (background and peaks estimation)
00039   static TH1F *EvalBackground(TH1F *his,Int_t n_max_peaks=10,Int_t niter=40,Option_t* option = "same");
00040 
00041   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");
00042 
00043   static Int_t PrintFoundPeaks(TSpectrum *spec);
00044   static Int_t LocatePeak(TSpectrum *spec,Double_t low_lim,Double_t high_lim);
00045 
00046 };
00047 
00049 
00050 
00051 // --- histograms informationmanipulation
00052 inline Int_t FIN_AN::GetHistogramLimits(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy){
00053   // required value are passed to arguments 
00054   // for reference
00055   //
00056   // return value:
00057   //              0: ok
00058   //             -1: error
00059   
00060   if(!his){
00061     gROOT->Warning("TFndGraph::GetHistogramLimits","Pointer to the given histogram is NULL!");
00062     return -1;
00063   }
00064   minx = static_cast<Double_t>( his->GetBinCenter(1) );
00065   maxx = static_cast<Double_t>( his->GetBinCenter(his->GetNbinsX()) + his->GetBinWidth(1) );
00066   miny = static_cast<Double_t>( his->GetMinimum() );
00067   maxy = static_cast<Double_t>( his->GetMaximum() );
00068   //
00069   /*   cout << " -> minx: " << minx << endl; */
00070   /*   cout << " -> maxx: " << maxx << endl; */
00071   /*   cout << " -> miny: " << miny << endl; */
00072   /*   cout << " -> maxy: " << maxy << endl; */
00073   //
00074   return 0;  
00075 }
00076 
00077 //inline Int_t FIN_AN::GetHistogramRanges(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy){
00078 //  // required value are passed to arguments for reference
00079 //  //
00080 //  // return value:
00081 //  //              0: ok
00082 //  //             -1: error
00083 //  
00084 //  if(!his){
00085 //    gROOT->Warning("TFndGraph::GetHistogramRanges","Pointer to the given histogram is NULL!");
00086 //    return -1;
00087 //  }
00088 //  minx = static_cast<Double_t>( his->GetXaxis()->GetXmin() );
00089 //  maxx = static_cast<Double_t>( his->GetXaxis()->GetXmax() );
00090 //  miny = static_cast<Double_t>( his->GetYaxis()->GetXmin() );
00091 //  maxy = static_cast<Double_t>( his->GetYaxis()->GetXmax() );
00092 //  //
00093 //  /*   cout << " -> minx: " << minx << endl; */
00094 //  /*   cout << " -> maxx: " << maxx << endl; */
00095 //  /*   cout << " -> miny: " << miny << endl; */
00096 //  /*   cout << " -> maxy: " << maxy << endl; */
00097 //  //
00098 //  return 0;  
00099 //}
00100 
00101 // --- static methods used as fit functions
00102 //      (and functions attributes setting)
00103 
00104 //_______________________
00105 inline Double_t FIN_AN::FitFCN_GausYield(Double_t *x, Double_t *par){
00106   // simple gaussian peak with the yield as parameter
00107   //                       (instead of normalization)
00108   // par[0,1,2]: yield,mean,sigma
00109   //
00110   // -------- IMPORTANT NOTE --------
00111   //   Remember that the number of events under the
00112   //   fitted gaussian peak (area) depends on the bin-width:
00113   //              Area = yield / bin-width
00114   //
00115   //     remember (confidence levels):
00116   //               +- 1 * sigma ===> 68.30 % of the total area
00117   //               +- 2 * sigma ===> 95.45 % of the total area
00118   //               +- 3 * sigma ===> 99.70 % of the total area
00119   // --------------------------------
00120 
00121   Double_t result = 0;
00122   
00123   Double_t DenNorm = (2*TMath::Pi())*TMath::Power(par[2],2);
00124   Double_t Norm =  par[0]/TMath::Sqrt(DenNorm); // par[0]: yield
00125   Double_t N1 = TMath::Power((x[0]-par[1]),2);
00126   Double_t N2 = 2 * TMath::Power(par[2],2);
00127   
00128   result = Norm * TMath::Exp(-N1/N2);
00129   
00130   return result;
00131 }
00132 
00133 //_______________________
00134 inline void FIN_AN::SetFitFCN_ParNames_GausYield(TF1 *fcn){
00135 
00136   fcn->SetParName(0,"yield");
00137   fcn->SetParName(1,"mean");
00138   fcn->SetParName(2,"sigma");
00139 
00140   fcn->SetLineWidth(2);
00141   fcn->SetLineColor(4);
00142 }
00143 
00144 // --- histograms evaluation (background and peaks estimation)
00145 
00146 //_______________________
00147 inline TH1F *FIN_AN::EvalBackground(TH1F *his,Int_t n_max_peaks,Int_t niter,Option_t* option){
00148   
00149   if(!his){
00150     gROOT->Warning("FIN_AN::EvalBackground","can not call method for NULL histogram.");
00151     return 0;
00152   }
00153 
00154   TSpectrum TmpSpec(n_max_peaks);
00155   return (TH1F *) ( TmpSpec.Background(his,niter,option) );
00156   
00157 }
00158 
00159 //_______________________
00160 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){
00161   // peaks and markers are added into the his->ListOfFinctions
00162   // return value: TSpectrum used
00163   //
00164   // threshold:  [ 0<threshold<1 ]  (default: 0.25):
00165   // peaks with amplitude less than threshold*highest_peak are discarded.              
00166   //
00167   // resol_fact [resolution factor] (default: 1):
00168   //              determines resolution of the neighboring peaks
00169   //              default value is 1 correspond to 3 sigma distance
00170   //              between peaks. Higher values allow higher resolution
00171   //              (smaller distance between peaks.
00172   //
00173   // option: (default:"")
00174   //  if option is not equal to "goff" , then          
00175   //   a polymarker object is created and added to the list of functions of  
00176   //   the histogram. The histogram is drawn with the specified option and   
00177   //   the polymarker object drawn on top of the histogram.                  
00178   //   The polymarker coordinates correspond to the npeaks peaks found in    
00179   //   the histogram.                                                        
00180   //   A pointer to the polymarker object can be retrieved later via:        
00181   //    TList *functions = hin->GetListOfFunctions();                        
00182   //    TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker")                                                        
00183   if(!his){
00184     gROOT->Warning("FIN_AN::FindPeaks","can not call method for NULL histogram.");
00185     return 0;
00186   }
00187   
00188   // Printf("Finding peaks for \"%s\"",his->GetName());
00189                
00190   TSpectrum *TmpSpectrum = new TSpectrum(n_max_peaks);
00191   TmpSpectrum->SetResolution(resol_fact);
00192   Int_t nfound = TmpSpectrum->Search(his,sigma,option,threshold);
00193   if(nfound > FinAn_Nmaxpeaks) Printf("too many peaks found: setting to %d",nfound);
00194   return TmpSpectrum;
00195 }
00196 
00197 //_____________________
00198 inline Int_t FIN_AN::PrintFoundPeaks(TSpectrum *spec){
00199   // Print all peaks found by the given "TSpectrum"
00200   //
00201   // return value:
00202   //              >=0: number of peaks found
00203   //               -1: errors
00204 
00205   if(!spec){
00206     gROOT->Warning("FIN_AN::PrintFoundPeaks","can not call method for NULL spectrum.");
00207     return -1;
00208   }
00209 
00210   Int_t n_found_peaks = spec->GetNPeaks();
00211   TString msg = "\nPeaks evaluated position: \n";
00212   TString linestr ="";
00213   for(Int_t pi=0;pi<n_found_peaks;pi++){
00214     linestr.Form("Peak number %d: position: %.4f ; height: %.4f \n",
00215                  pi+1,
00216                  (*(spec->GetPositionX() + pi)),
00217                  (*(spec->GetPositionY() + pi))
00218                  );
00219     msg+=linestr;
00220   }
00221   Printf(msg.Data());
00222   return n_found_peaks;
00223 }
00224 
00225 //_____________________
00226 inline Int_t FIN_AN::LocatePeak(TSpectrum *spec,Double_t low_lim,Double_t high_lim){
00227   // Print all peaks found by the given "TSpectrum"
00228   //
00229   // return value:
00230   //               > 0: number of the searched peak (starting from 1)
00231   //                -1: peak not found
00232   //                -2: error
00233 
00234   if(!spec){
00235     gROOT->Warning("FIN_AN::LocatePeak","can not call method for NULL spectrum.");
00236     return -2;
00237   }
00238 
00239   Int_t n_found_peaks = spec->GetNPeaks();
00240   for(Int_t pi=0;pi<n_found_peaks;pi++){
00241     Double_t cval = (*(spec->GetPositionX() + pi));
00242     if( cval >= low_lim && cval <= high_lim){
00243       //      Printf("Peak found: %.3f",cval);
00244       return pi+1;
00245     }
00246   }
00247   Printf("FIN_AN::LocatePeak---> Peak not found in range [%.4f ; %.4f]",low_lim,high_lim);
00248   return -1;
00249 }
00250 
00251 #endif // FROOT_FIN_AN
00253 

Generated on Tue Oct 16 15:40:46 2007 by  doxygen 1.5.2