Files
2016-06-09 16:25:56 +02:00

63 lines
1.8 KiB
C++

//
// Fit of longitudinal profile with Gamma function
// See Review of Particles Physics - Electromagnetic cascades
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Double_t GammaFunction(Double_t* t, Double_t* params ) {
if (ROOT::TMath::Gamma(params[0]) != 0 ) {
return (100 * pow(params[1],params[0]) / ROOT::TMath::Gamma(params[0]))
* pow(t[0],(params[0]-1)) * exp (-params[1]*t[0]);
} else {
return 0;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void GammaFit() {
cout
<< "Try fit [100*pow(b,a)/ROOT::TMath::Gamma(a)] * pow(t,(a-1)) * exp(-b*t)"
<< endl;
TFile* f = new TFile("93ref0.root");
TH1D* h1 = (TH1D*) f->Get("4");
//put error in h1
TH1D* h2 = (TH1D*) f->Get("5");
int nmax = h1->GetNbinsX();
for (int n=0; n<nmax; n++) {
Double_t er = h2->GetBinContent(n);
h1->SetBinError(n,er);
}
TF1* func = new TF1("fit", GammaFunction,-0.5,40.,2);
func->SetParameter(0,1.);
func->SetParameter(1,1.);
func->SetParNames("a", "b");
h1->Fit("fit","r","");
h1->SetMarkerColor(kRed);
h1->SetMarkerStyle(3);
h1->Draw("epl");
Double_t a = func->GetParameter(0);
Double_t b = func->GetParameter(1);
Double_t tmax = (a-1)/b;
// Chi-square distribution
double chisq=func->GetChisquare();
double ndf=func->GetNDF();
double chisqdf=chisq/ndf;
gStyle->SetOptFit(1111);
cout << "----------------------------------------------";
cout << "\n Chisquare: " << chisq << " / " << ndf << " = " << chisqdf;
cout << "\n----------------------------------------------";
cout << "\n a = " << a;
cout << "\n b = " << b;
cout << "\n tmax = " << tmax;
cout << "\n----------------------------------------------" << endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......