LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GammaFit.C
Go to the documentation of this file.
1 //
2 // Fit of longitudinal profile with Gamma function
3 // See Review of Particles Physics - Electromagnetic cascades
4 //
5 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
6 
7 Double_t GammaFunction(Double_t* t, Double_t* params ) {
8  if (ROOT::TMath::Gamma(params[0]) != 0 ) {
9  return (100 * pow(params[1],params[0]) / ROOT::TMath::Gamma(params[0]))
10  * pow(t[0],(params[0]-1)) * exp (-params[1]*t[0]);
11  } else {
12  return 0;
13  }
14 }
15 
16 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
17 
18 void GammaFit() {
19  cout
20  << "Try fit [100*pow(b,a)/ROOT::TMath::Gamma(a)] * pow(t,(a-1)) * exp(-b*t)"
21  << endl;
22 TFile* f = new TFile("93ref0.root");
23 TH1D* h1 = (TH1D*) f->Get("4");
24 //put error in h1
25 TH1D* h2 = (TH1D*) f->Get("5");
26 int nmax = h1->GetNbinsX();
27 for (int n=0; n<nmax; n++) {
28  Double_t er = h2->GetBinContent(n);
29  h1->SetBinError(n,er);
30 }
31 
32 TF1* func = new TF1("fit", GammaFunction,-0.5,40.,2);
33 func->SetParameter(0,1.);
34 func->SetParameter(1,1.);
35 func->SetParNames("a", "b");
36 
37 h1->Fit("fit","r","");
38 h1->SetMarkerColor(kRed);
39 h1->SetMarkerStyle(3);
40 h1->Draw("epl");
41 
42 Double_t a = func->GetParameter(0);
43 Double_t b = func->GetParameter(1);
44 Double_t tmax = (a-1)/b;
45 
46 // Chi-square distribution
47 double chisq=func->GetChisquare();
48 double ndf=func->GetNDF();
49 double chisqdf=chisq/ndf;
50 
51 gStyle->SetOptFit(1111);
52 
53 cout << "----------------------------------------------";
54 cout << "\n Chisquare: " << chisq << " / " << ndf << " = " << chisqdf;
55 cout << "\n----------------------------------------------";
56 cout << "\n a = " << a;
57 cout << "\n b = " << b;
58 cout << "\n tmax = " << tmax;
59 cout << "\n----------------------------------------------" << endl;
60 }
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void GammaFit()
Definition: GammaFit.C:18
TFile f
Definition: plotHisto.C:6
TH1F * h2
Definition: plot.C:44
Double_t GammaFunction(Double_t *t, Double_t *params)
Definition: GammaFit.C:7
TH1F * h1
Definition: plot.C:41
Char_t n[5]