Logo ROOT  
Reference Guide
langaus.C File Reference

Detailed Description

View in nbviewer Open in SWAN Convoluted Landau and Gaussian Fitting Function (using ROOT's Landau and Gauss functions)

Based on a Fortran code by R.Fruehwirth (fruhw.nosp@m.irth.nosp@m.@heph.nosp@m.y.oe.nosp@m.aw.ac.nosp@m..at)

to execute this example, do:

root > .x langaus.C

or

root > .x langaus.C++
Fitting...
FCN=5.25252 FROM MIGRAD STATUS=CONVERGED 193 CALLS 194 TOTAL
EDM=3.227e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 Width 1.25740e+00 3.04846e-01 1.37307e-04 1.57021e-04
2 MP 2.08890e+01 1.28244e+00 3.40886e-05 8.22843e-04
3 Area 1.15515e+04 2.42312e+03 1.66644e-05 -1.66635e-03
4 GSigma 4.06350e+00 7.58845e-01 2.30703e-04 -2.48397e-04
Fitting done
Plotting results...
#include "TH1.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TMath.h"
double langaufun(double *x, double *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 invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
double mpshift = -0.22278298; // Landau maximum location
// Control constants
double np = 100.0; // number of convolution steps
double sc = 5.0; // convolution extends to +-sc Gaussian sigmas
// Variables
double xx;
double mpc;
double fland;
double sum = 0.0;
double xlow,xupp;
double step;
double 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(TH1F *his, double *fitrange, double *startvalues, double *parlimitslo, double *parlimitshi, double *fitparams, double *fiterrors, double *ChiSqr, int *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 i;
char 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 langaupro(double *params, double &maxx, double &FWHM) {
// Searches 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 p,x,fy,fxr,fxl;
double step;
double l,lold;
int i = 0;
int 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() {
// Fill Histogram
int data[100] = {0,0,0,0,0,0,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 i=0; i<100; i++) hSNR->Fill(i,data[i]);
// Fitting SNR histo
printf("Fitting...\n");
// Setting fit range and start values
double fr[2];
double sv[4], pllo[4], plhi[4], fp[4], fpe[4];
fr[0]=0.3*hSNR->GetMean();
fr[1]=3.0*hSNR->GetMean();
pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.4;
plhi[0]=5.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=5.0;
sv[0]=1.8; sv[1]=20.0; sv[2]=50000.0; sv[3]=3.0;
double chisqr;
int ndf;
TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
double SNRPeak, SNRFWHM;
langaupro(fp,SNRPeak,SNRFWHM);
printf("Fitting done\nPlotting results...\n");
// Global style settings
gStyle->SetLabelSize(0.03,"x");
gStyle->SetLabelSize(0.03,"y");
hSNR->GetXaxis()->SetRange(0,70);
hSNR->Draw();
fitsnr->Draw("lsame");
}
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
#define gROOT
Definition: TROOT.h:404
R__EXTERN TStyle * gStyle
Definition: TStyle.h:414
virtual void SetRange(Int_t first=0, Int_t last=0)
Set the viewing range for the axis using bin numbers.
Definition: TAxis.cxx:947
1-Dim function class
Definition: TF1.h:213
virtual Int_t GetNDF() const
Return the number of degrees of freedom in the fit the fNDF parameter has been previously computed du...
Definition: TF1.cxx:1909
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1950
Double_t GetChisquare() const
Definition: TF1.h:449
virtual Double_t * GetParameters() const
Definition: TF1.h:525
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition: TF1.cxx:1347
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set lower and upper limits for parameter ipar.
Definition: TF1.cxx:3549
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:649
virtual void SetParNames(const char *name0="p0", const char *name1="p1", const char *name2="p2", const char *name3="p3", const char *name4="p4", const char *name5="p5", const char *name6="p6", const char *name7="p7", const char *name8="p8", const char *name9="p9", const char *name10="p10")
Set up to 10 parameter names.
Definition: TF1.cxx:3510
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:574
virtual Double_t GetMean(Int_t axis=1) const
For axis = 1,2 or 3 returns the mean value of the histogram along X,Y or Z axis.
Definition: TH1.cxx:7456
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:319
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3904
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3348
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3070
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1589
void SetLabelSize(Float_t size=0.04, Option_t *axis="X")
Set size of axis labels.
Definition: TStyle.cxx:1393
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1542
Double_t x[n]
Definition: legend1.C:17
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculates a gaussian function with mean and sigma.
Definition: TMath.cxx:471
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
Definition: TMath.cxx:492
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
TLine l
Definition: textangle.C:4
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345
Authors
H.Pernegger, Markus Friedl

Definition in file langaus.C.