//Example for fitting signal/background. // This example can be executed with: // root > .x FittingDemo.C (using the CINT interpreter) // root > .x FittingDemo.C+ (using the native complier via ACLIC) //Author: Rene Brun #include "TH1.h" #include "TMath.h" #include "TF1.h" #include "TLegend.h" #include "TCanvas.h" // Quadratic background function Double_t background(Double_t *x, Double_t *par) { return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; } // Lorenzian Peak function Double_t lorentzianPeak(Double_t *x, Double_t *par) { return (0.5*par[0]*par[1]/TMath::Pi()) / TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]); } // Sum of background and peak function Double_t fitFunction(Double_t *x, Double_t *par) { return background(x,par) + lorentzianPeak(x,&par[3]); } void NewFittingDemo() { //Bevington Exercise by Peter Malzacher, modified by Rene Brun ifstream in; in.open("./output.txt"); Int_t x,y, data[1000000][2], Hash[1000000]; Double_t xAry[1000000], yAry[1000000]; Long64_t index[1000000], n=0; Double_t yMean=0,xMean=0, XShift; //Choose the mean 0:the first mean, 1: the second and so on int Mean=1; for(Int_t i=0; i<1000000; i++) Hash[i]=0; while (1) { in >> x; if (!in.good()) break; Hash[x]++; } for(i=0; i<1000000; i++) { if (Hash[i]!=0) { data[n][0]=i-150; data[n][1]=Hash[i]; n++; } } in.close(); Int_t p1=0; //the peak on intervals- set the size of interval on x axis XShift=1000; for(Int_t i=0; (i+XShift)=i)&&(data[r][0]<=i+XShift)) { if(data[r][1]>yMean)) { yMean=data[r][1]; xMean=data[r][0]; } } } index[p1]=p1; yAry[p1]=yMean; xAry[p1]=xMean; p1++; } TMath::Sort(p1, yAry, index, 1); const int nBins = n; TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,700,500); //c1->SetFillColor(30); //c1->SetFrameFillColor(41); c1->SetGrid(); TH1F *histo = new TH1F("histo"," Peak Fitting",8192,0,8192); histo->SetMarkerStyle(20); histo->SetMarkerSize(0.8); histo->SetStats(0); for(int i=0; i < nBins; i++) histo->SetBinContent(data[i][0],data[i][1]); // create a TF1 with the range from 0 to 3 and 6 parameters TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,8192,6); fitFcn->SetNpx(500); fitFcn->SetLineWidth(4); fitFcn->SetLineColor(kMagenta); // first try without starting values for the parameters // This defaults to 1 for each param. // this results in an ok fit for the polynomial function // however the non-linear part (lorenzian) does not // respond well. //fitFcn->SetParameters(1,1,1,1,1,1); histo->Fit("fitFcn","0"); // second try: set start values for some parameters fitFcn->SetParameter(4,80); // width fitFcn->SetParameter(5,xAry[index[Mean]]); // peak histo->Fit("fitFcn","V+","ep"); // improve the picture: TF1 *backFcn = new TF1("backFcn",background,0,8192,3); backFcn->SetLineColor(kRed); TF1 *signalFcn = new TF1("signalFcn",lorentzianPeak,0,8192,3); signalFcn->SetLineColor(kBlue); signalFcn->SetNpx(500); Double_t par[6]; // writes the fit results into the par array fitFcn->GetParameters(par); backFcn->SetParameters(par); backFcn->Draw("same"); signalFcn->SetParameters(&par[3]); signalFcn->Draw("same"); // draw the legend TLegend *legend=new TLegend(0.6,0.65,0.88,0.85); legend->SetTextFont(72); legend->SetTextSize(0.04); legend->AddEntry(histo,"Data","lpe"); legend->AddEntry(backFcn,"Background fit","l"); legend->AddEntry(signalFcn,"Signal fit","l"); legend->AddEntry(fitFcn,"Global Fit","l"); legend->Draw(); printf("=============display peaks==============\n"); for(i=0; iIntegral(0,n)); cout<<"\Integration="<