Logo ROOT  
Reference Guide
langaus.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_fit
3/// \notebook
4/// Convoluted Landau and Gaussian Fitting Function
5/// (using ROOT's Landau and Gauss functions)
6///
7/// Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
8///
9/// to execute this example, do:
10///
11/// ~~~{.cpp}
12/// root > .x langaus.C
13/// ~~~
14///
15/// or
16///
17/// ~~~{.cpp}
18/// root > .x langaus.C++
19/// ~~~
20///
21/// \macro_image
22/// \macro_output
23/// \macro_code
24///
25/// \authors H.Pernegger, Markus Friedl
26
27#include "TH1.h"
28#include "TF1.h"
29#include "TROOT.h"
30#include "TStyle.h"
31#include "TMath.h"
32
33Double_t langaufun(Double_t *x, Double_t *par) {
34
35 //Fit parameters:
36 //par[0]=Width (scale) parameter of Landau density
37 //par[1]=Most Probable (MP, location) parameter of Landau density
38 //par[2]=Total area (integral -inf to inf, normalization constant)
39 //par[3]=Width (sigma) of convoluted Gaussian function
40 //
41 //In the Landau distribution (represented by the CERNLIB approximation),
42 //the maximum is located at x=-0.22278298 with the location parameter=0.
43 //This shift is corrected within this function, so that the actual
44 //maximum is identical to the MP parameter.
45
46 // Numeric constants
47 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
48 Double_t mpshift = -0.22278298; // Landau maximum location
49
50 // Control constants
51 Double_t np = 100.0; // number of convolution steps
52 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
53
54 // Variables
55 Double_t xx;
56 Double_t mpc;
57 Double_t fland;
58 Double_t sum = 0.0;
59 Double_t xlow,xupp;
60 Double_t step;
61 Double_t i;
62
63
64 // MP shift correction
65 mpc = par[1] - mpshift * par[0];
66
67 // Range of convolution integral
68 xlow = x[0] - sc * par[3];
69 xupp = x[0] + sc * par[3];
70
71 step = (xupp-xlow) / np;
72
73 // Convolution integral of Landau and Gaussian by sum
74 for(i=1.0; i<=np/2; i++) {
75 xx = xlow + (i-.5) * step;
76 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
77 sum += fland * TMath::Gaus(x[0],xx,par[3]);
78
79 xx = xupp - (i-.5) * step;
80 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
81 sum += fland * TMath::Gaus(x[0],xx,par[3]);
82 }
83
84 return (par[2] * step * sum * invsq2pi / par[3]);
85}
86
87
88
89TF1 *langaufit(TH1F *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)
90{
91 // Once again, here are the Landau * Gaussian parameters:
92 // par[0]=Width (scale) parameter of Landau density
93 // par[1]=Most Probable (MP, location) parameter of Landau density
94 // par[2]=Total area (integral -inf to inf, normalization constant)
95 // par[3]=Width (sigma) of convoluted Gaussian function
96 //
97 // Variables for langaufit call:
98 // his histogram to fit
99 // fitrange[2] lo and hi boundaries of fit range
100 // startvalues[4] reasonable start values for the fit
101 // parlimitslo[4] lower parameter limits
102 // parlimitshi[4] upper parameter limits
103 // fitparams[4] returns the final fit parameters
104 // fiterrors[4] returns the final fit errors
105 // ChiSqr returns the chi square
106 // NDF returns ndf
107
108 Int_t i;
109 Char_t FunName[100];
110
111 sprintf(FunName,"Fitfcn_%s",his->GetName());
112
113 TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
114 if (ffitold) delete ffitold;
115
116 TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
117 ffit->SetParameters(startvalues);
118 ffit->SetParNames("Width","MP","Area","GSigma");
119
120 for (i=0; i<4; i++) {
121 ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
122 }
123
124 his->Fit(FunName,"RB0"); // fit within specified range, use ParLimits, do not plot
125
126 ffit->GetParameters(fitparams); // obtain fit parameters
127 for (i=0; i<4; i++) {
128 fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors
129 }
130 ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2
131 NDF[0] = ffit->GetNDF(); // obtain ndf
132
133 return (ffit); // return fit function
134
135}
136
137
138Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
139
140 // Seaches for the location (x value) at the maximum of the
141 // Landau-Gaussian convolute and its full width at half-maximum.
142 //
143 // The search is probably not very efficient, but it's a first try.
144
145 Double_t p,x,fy,fxr,fxl;
146 Double_t step;
147 Double_t l,lold;
148 Int_t i = 0;
149 Int_t MAXCALLS = 10000;
150
151
152 // Search for maximum
153
154 p = params[1] - 0.1 * params[0];
155 step = 0.05 * params[0];
156 lold = -2.0;
157 l = -1.0;
158
159
160 while ( (l != lold) && (i < MAXCALLS) ) {
161 i++;
162
163 lold = l;
164 x = p + step;
165 l = langaufun(&x,params);
166
167 if (l < lold)
168 step = -step/10;
169
170 p += step;
171 }
172
173 if (i == MAXCALLS)
174 return (-1);
175
176 maxx = x;
177
178 fy = l/2;
179
180
181 // Search for right x location of fy
182
183 p = maxx + params[0];
184 step = params[0];
185 lold = -2.0;
186 l = -1e300;
187 i = 0;
188
189
190 while ( (l != lold) && (i < MAXCALLS) ) {
191 i++;
192
193 lold = l;
194 x = p + step;
195 l = TMath::Abs(langaufun(&x,params) - fy);
196
197 if (l > lold)
198 step = -step/10;
199
200 p += step;
201 }
202
203 if (i == MAXCALLS)
204 return (-2);
205
206 fxr = x;
207
208
209 // Search for left x location of fy
210
211 p = maxx - 0.5 * params[0];
212 step = -params[0];
213 lold = -2.0;
214 l = -1e300;
215 i = 0;
216
217 while ( (l != lold) && (i < MAXCALLS) ) {
218 i++;
219
220 lold = l;
221 x = p + step;
222 l = TMath::Abs(langaufun(&x,params) - fy);
223
224 if (l > lold)
225 step = -step/10;
226
227 p += step;
228 }
229
230 if (i == MAXCALLS)
231 return (-3);
232
233
234 fxl = x;
235
236 FWHM = fxr - fxl;
237 return (0);
238}
239
240void langaus() {
241 // Fill Histogram
242 Int_t data[100] = {0,0,0,0,0,0,2,6,11,18,18,55,90,141,255,323,454,563,681,
243 737,821,796,832,720,637,558,519,460,357,291,279,241,212,
244 153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
245 22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
246 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};
247 TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
248
249 for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
250
251 // Fitting SNR histo
252 printf("Fitting...\n");
253
254 // Setting fit range and start values
255 Double_t fr[2];
256 Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
257 fr[0]=0.3*hSNR->GetMean();
258 fr[1]=3.0*hSNR->GetMean();
259
260 pllo[0]=0.5; pllo[1]=5.0; pllo[2]=1.0; pllo[3]=0.4;
261 plhi[0]=5.0; plhi[1]=50.0; plhi[2]=1000000.0; plhi[3]=5.0;
262 sv[0]=1.8; sv[1]=20.0; sv[2]=50000.0; sv[3]=3.0;
263
264 Double_t chisqr;
265 Int_t ndf;
266 TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
267
268 Double_t SNRPeak, SNRFWHM;
269 langaupro(fp,SNRPeak,SNRFWHM);
270
271 printf("Fitting done\nPlotting results...\n");
272
273 // Global style settings
274 gStyle->SetOptStat(1111);
275 gStyle->SetOptFit(111);
276 gStyle->SetLabelSize(0.03,"x");
277 gStyle->SetLabelSize(0.03,"y");
278
279 hSNR->GetXaxis()->SetRange(0,70);
280 hSNR->Draw();
281 fitsnr->Draw("lsame");
282}
283
int Int_t
Definition: RtypesCore.h:43
char Char_t
Definition: RtypesCore.h:31
double Double_t
Definition: RtypesCore.h:57
#define gROOT
Definition: TROOT.h:406
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
virtual void SetRange(Int_t first=0, Int_t last=0)
Set the viewing range for the axis from bin first to last.
Definition: TAxis.cxx:914
1-Dim function class
Definition: TF1.h:210
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:1872
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition: TF1.cxx:1913
Double_t GetChisquare() const
Definition: TF1.h:438
virtual Double_t * GetParameters() const
Definition: TF1.h:514
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set limits for parameter ipar.
Definition: TF1.cxx:3500
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1320
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
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:3464
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
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:7086
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
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:3808
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual const char * GetName() const
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:1590
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)
Calculate a gaussian function with mean and sigma.
Definition: TMath.cxx:448
Double_t Landau(Double_t x, Double_t mpv=0, Double_t sigma=1, Bool_t norm=kFALSE)
The LANDAU function.
Definition: TMath.cxx:469
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * l
Definition: textangle.C:4
static long int sum(long int i)
Definition: Factory.cxx:2275