00001
00002
00003
00004 #ifndef FROOT_FIN_AN
00005 #define FROOT_FIN_AN
00006
00008
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
00019
00020 namespace FIN_AN {
00021
00022 enum E_FinAn_NMaxPeaks { FinAn_Nmaxpeaks = 500 };
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
00029 static Int_t GetHistogramLimits(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy);
00030
00031
00032
00033
00034 static Double_t FitFCN_GausYield(Double_t *x, Double_t *par);
00035
00036 static void SetFitFCN_ParNames_GausYield(TF1 *fcn);
00037
00038
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
00052 inline Int_t FIN_AN::GetHistogramLimits(TH1F *his,Double_t &minx,Double_t &maxx,Double_t &miny,Double_t &maxy){
00053
00054
00055
00056
00057
00058
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
00070
00071
00072
00073
00074 return 0;
00075 }
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 inline Double_t FIN_AN::FitFCN_GausYield(Double_t *x, Double_t *par){
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
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);
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
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
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183 if(!his){
00184 gROOT->Warning("FIN_AN::FindPeaks","can not call method for NULL histogram.");
00185 return 0;
00186 }
00187
00188
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
00200
00201
00202
00203
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
00228
00229
00230
00231
00232
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
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