|
|
/// \ingroup tutorial_fit
/// \notebook
/// Convoluted Landau and Gaussian Fitting Function
/// (using ROOT's Landau and Gauss functions)
///
/// Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
///
/// to execute this example, do:
///
/// ~~~{.cpp}
/// root > .x langaus.C
/// ~~~
///
/// or
///
/// ~~~{.cpp}
/// root > .x langaus.C++
/// ~~~
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \authors H.Pernegger, Markus Friedl
#include "TH1.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"
Double_t langaufun(Double_t *x, Double_t *par) { //Fit parameters:
//par[0]=Width (scale) parameter of Landau density
//par[1]=Most Probable (MP, location) parameter of Landau density
//par[2]=Total area (integral -inf to inf, normalization constant)
//par[3]=Width (sigma) of convoluted Gaussian function
//
//In the Landau distribution (represented by the CERNLIB approximation),
//the maximum is located at x=-0.22278298 with the location parameter=0.
//This shift is corrected within this function, so that the actual
//maximum is identical to the MP parameter.
// Numeric constants
Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
Double_t mpshift = -0.22278298; // Landau maximum location
// Control constants
Double_t np = 100.0; // number of convolution steps
Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
// Variables
Double_t xx; Double_t mpc; Double_t fland; Double_t sum = 0.0; Double_t xlow,xupp; Double_t step; Double_t i; // MP shift correction
mpc = par[1] - mpshift * par[0]; // Range of convolution integral
xlow = x[0] - sc * par[3]; xupp = x[0] + sc * par[3]; step = (xupp-xlow) / np; // Convolution integral of Landau and Gaussian by sum
for(i=1.0; i<=np/2; i++) { xx = xlow + (i-.5) * step; fland = TMath::Landau(xx,mpc,par[0]) / par[0]; sum += fland * TMath::Gaus(x[0],xx,par[3]); xx = xupp - (i-.5) * step; fland = TMath::Landau(xx,mpc,par[0]) / par[0]; sum += fland * TMath::Gaus(x[0],xx,par[3]); } return (par[2] * step * sum * invsq2pi / par[3]); } TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF) { // Once again, here are the Landau * Gaussian parameters:
// par[0]=Width (scale) parameter of Landau density
// par[1]=Most Probable (MP, location) parameter of Landau density
// par[2]=Total area (integral -inf to inf, normalization constant)
// par[3]=Width (sigma) of convoluted Gaussian function
//
// Variables for langaufit call:
// his histogram to fit
// fitrange[2] lo and hi boundaries of fit range
// startvalues[4] reasonable start values for the fit
// parlimitslo[4] lower parameter limits
// parlimitshi[4] upper parameter limits
// fitparams[4] returns the final fit parameters
// fiterrors[4] returns the final fit errors
// ChiSqr returns the chi square
// NDF returns ndf
Int_t i; Char_t FunName[100]; sprintf(FunName,"Fitfcn_%s",his->GetName()); TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName); if (ffitold) delete ffitold; TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4); ffit->SetParameters(startvalues); ffit->SetParNames("Width","MP","Area","GSigma"); for (i=0; i<4; i++) { ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]); } his->Fit(FunName,"RB0"); // fit within specified range, use ParLimits, do not plot
ffit->GetParameters(fitparams); // obtain fit parameters
for (i=0; i<4; i++) { fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors
} ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2
NDF[0] = ffit->GetNDF(); // obtain ndf
return (ffit); // return fit function
} Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) { // Seaches for the location (x value) at the maximum of the
// Landau-Gaussian convolute and its full width at half-maximum.
//
// The search is probably not very efficient, but it's a first try.
Double_t p,x,fy,fxr,fxl; Double_t step; Double_t l,lold; Int_t i = 0; Int_t MAXCALLS = 10000; // Search for maximum
p = params[1] - 0.1 * params[0]; step = 0.05 * params[0]; lold = -2.0; l = -1.0; while ( (l != lold) && (i < MAXCALLS) ) { i++; lold = l; x = p + step; l = langaufun(&x,params); if (l < lold) step = -step/10; p += step; } if (i == MAXCALLS) return (-1); maxx = x; fy = l/2; // Search for right x location of fy
p = maxx + params[0]; step = params[0]; lold = -2.0; l = -1e300; i = 0; while ( (l != lold) && (i < MAXCALLS) ) { i++; lold = l; x = p + step; l = TMath::Abs(langaufun(&x,params) - fy); if (l > lold) step = -step/10; p += step; } if (i == MAXCALLS) return (-2); fxr = x; // Search for left x location of fy
p = maxx - 0.5 * params[0]; step = -params[0]; lold = -2.0; l = -1e300; i = 0; while ( (l != lold) && (i < MAXCALLS) ) { i++; lold = l; x = p + step; l = TMath::Abs(langaufun(&x,params) - fy); if (l > lold) step = -step/10; p += step; } if (i == MAXCALLS) return (-3); fxl = x; FWHM = fxr - fxl; return (0); } void langaus(string rootfilename, int ionsort) { // Fill Histogram
/* Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
737,821,796,832,720,637,558,519,460,357,291,279,241,212, 153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23, 22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4, 9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4}; TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400); */ // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
TFile *rootFile = new TFile(rootfilename.c_str(),"OPEN"); TFile *outFile = new TFile("langausout.root","RECREATE"); TH1D * hSNR; TH1D* h_beamSignal_b0[30]; TH1D* h_beamSignal_b1[30]; TH1D* h_beamSignal_b2[30]; TH1D* h_beamSignal_b3[30]; int energy_indexbin = 0; char histname[30] = ""; TH1D* h_beamSignal_ic1[30];
Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30]; TGraphErrors * graph_1; Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169}; Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313}; Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602}; Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424}; for (int i = 0;i<23;i++){ if (ionsort ==3 ) graph_x[i] = beta_proton[i]; if (ionsort ==1 ) graph_x[i] = beta_helium[i]; if (ionsort ==4 ) graph_x[i] = beta_carbon[i]; if (ionsort ==2 ) graph_x[i] = beta_oxygen[i]; sprintf(histname,"h_beamSignal_b0_%d",i); h_beamSignal_b0[i] = (TH1D*)rootFile->Get(histname); sprintf(histname,"h_beamSignal_b1_%d",i); h_beamSignal_b1[i] = (TH1D*)rootFile->Get(histname); sprintf(histname,"h_beamSignal_b2_%d",i); h_beamSignal_b2[i] = (TH1D*)rootFile->Get(histname); sprintf(histname,"h_beamSignal_b3_%d",i); h_beamSignal_b3[i] = (TH1D*)rootFile->Get(histname); sprintf(histname,"h_beamSignal_ic1_%d",i); h_beamSignal_ic1[i] = (TH1D*)rootFile->Get(histname); } // hSNR = h_beamSignal_b0[13];
// Fitting SNR histo
printf("Fitting...\n"); // Setting fit range and start values
Double_t fr[2]; Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4]; Double_t chisqr; Int_t ndf; TF1 *fitsnr[24]; Double_t SNRPeak, SNRFWHM; TCanvas * c1 = new TCanvas("c1","c1", 800, 600); Double_t norm = 1; for (int i = 0; i<23;i++){ //hSNR = h_beamSignal_b0[i];
h_beamSignal_ic1[i]->Scale(norm/h_beamSignal_ic1[i]->Integral("width")); hSNR = h_beamSignal_ic1[i]; fr[0]=hSNR->GetMean() - 3*hSNR->GetRMS(); fr[1]=hSNR->GetMean() + 3*hSNR->GetRMS(); pllo[0]=0.0075; pllo[1]=1.0; pllo[2]=0.01; pllo[3]=0.07; plhi[0]=1.0E3; plhi[1]=300.0; plhi[2]=2E2; plhi[3]=0.65; sv[0]=20.2; //sv[1]=2.1;
sv[2]=1.0; sv[3]=hSNR->GetRMS(); sv[1] = hSNR->GetMean(); fitsnr[i] = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
langaupro(fp,SNRPeak,SNRFWHM); graph_y[i] = fp[1]; graph_yerr[i] = sqrt(fitsnr[i]->GetParError(1)*fitsnr[i]->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]); printf("Fitting done\nPlotting results...\n"); // Global style settings
gStyle->SetOptStat(1111); gStyle->SetOptFit(111); //gStyle->SetLabelSize(0.03,"x");
//gStyle->SetLabelSize(0.03,"y");
// hSNR->GetXaxis()->SetRange(0,70);
// hSNR->Draw();
// fitsnr[i]->Draw("lsame");
// c1->Update();
// c1->WaitPrimitive();
} graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr); TCanvas * c2 = new TCanvas("c2","c2", 800, 600); c2->cd(); graph_1->Draw("A*"); outFile->Write(); outFile->Close(); }
|
