TestEm2/egs4/GammaFit.C

0001 //
0002 // Fit of longitudinal profile with Gamma function
0003 // See Review of Particles Physics - Electromagnetic cascades
0004 //
0005 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0006
0007 Double_t GammaFunction(Double_t* t, Double_t* params ) {
0008   if (ROOT::TMath::Gamma(params[0]) != 0 ) {
0009    return  (100 * pow(params[1],params[0]) / ROOT::TMath::Gamma(params[0]))
0010            * pow(t[0],(params[0]-1)) * exp (-params[1]*t[0]);
0011   } else {
0012    return 0;
0013   }
0014 }
0015
0016 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0017
0018 void GammaFit() {
0019  cout
0020   << "Try fit [100*pow(b,a)/ROOT::TMath::Gamma(a)] * pow(t,(a-1)) * exp(-b*t)"
0021   << endl;
0022 TFile* f = new  TFile("93ref0.root");
0023 TH1D* h1 = (TH1D*) f->Get("4");
0024 //put error in h1
0025 TH1D* h2 = (TH1D*) f->Get("5");
0026 int nmax = h1->GetNbinsX();
0027 for (int n=0; n<nmax; n++) {
0028   Double_t er = h2->GetBinContent(n);
0029   h1->SetBinError(n,er);
0030 }
0031
0032 TF1* func = new TF1("fit", GammaFunction,-0.5,40.,2);
0033 func->SetParameter(0,1.);
0034 func->SetParameter(1,1.);
0035 func->SetParNames("a", "b");
0036
0037 h1->Fit("fit","r","");
0038 h1->SetMarkerColor(kRed);
0039 h1->SetMarkerStyle(3);
0040 h1->Draw("epl");
0041
0042 Double_t a = func->GetParameter(0);
0043 Double_t b = func->GetParameter(1);
0044 Double_t tmax = (a-1)/b;
0045
0046 // Chi-square distribution
0047 double chisq=func->GetChisquare();
0048 double ndf=func->GetNDF();
0049 double chisqdf=chisq/ndf;
0050
0051 gStyle->SetOptFit(1111);
0052
0053 cout << "----------------------------------------------";
0054 cout << "\n Chisquare: " << chisq << " / " << ndf << " = " << chisqdf;
0055 cout << "\n----------------------------------------------";
0056 cout << "\n a = " << a;
0057 cout << "\n b = " << b;
0058 cout << "\n tmax = " << tmax;
0059 cout << "\n----------------------------------------------" << endl;
0060 }
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......