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