Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 langaufun(double *x, double *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 invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
48 double mpshift = -0.22278298; // Landau maximum location
49
50 // Control constants
51 double np = 100.0; // number of convolution steps
52 double sc = 5.0; // convolution extends to +-sc Gaussian sigmas
53
54 // Variables
55 double xx;
56 double mpc;
57 double fland;
58 double sum = 0.0;
59 double xlow,xupp;
60 double step;
61 double 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 *fitrange, double *startvalues, double *parlimitslo, double *parlimitshi, double *fitparams, double *fiterrors, double *ChiSqr, int *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 i;
109 char 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 langaupro(double *params, double &maxx, double &FWHM) {
139
140 // Searches 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 p,x,fy,fxr,fxl;
146 double step;
147 double l,lold;
148 int i = 0;
149 int 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 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 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 fr[2];
256 double 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 chisqr;
265 int ndf;
266 TF1 *fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
267
268 double 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
double fy() const
winID h TVirtualViewer3D TVirtualGLPainter p
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
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
#define gROOT
Definition TROOT.h:406
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
virtual void SetRange(Int_t first=0, Int_t last=0)
Set the viewing range for the axis using bin numbers.
Definition TAxis.cxx:1052
1-Dim function class
Definition TF1.h:233
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:1891
virtual Double_t GetParError(Int_t ipar) const
Return value of parameter number ipar.
Definition TF1.cxx:1932
Double_t GetChisquare() const
Return the Chisquare after fitting. See ROOT::Fit::FitResult::Chi2()
Definition TF1.h:470
virtual void SetParNames(const char *name0="", const char *name1="", const char *name2="", const char *name3="", const char *name4="", const char *name5="", const char *name6="", const char *name7="", const char *name8="", const char *name9="", const char *name10="")
Set up to 10 parameter names.
Definition TF1.cxx:3463
virtual Double_t * GetParameters() const
Definition TF1.h:546
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition TF1.cxx:1335
virtual void SetParLimits(Int_t ipar, Double_t parmin, Double_t parmax)
Set lower and upper limits for parameter ipar.
Definition TF1.cxx:3507
virtual void SetParameters(const Double_t *params)
Definition TF1.h:670
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:621
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:7499
TAxis * GetXaxis()
Definition TH1.h:324
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:3894
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3340
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
virtual TObject * FindObject(const char *name) const
Must be redefined in derived classes.
Definition TObject.cxx:403
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:1636
void SetLabelSize(Float_t size=0.04, Option_t *axis="X")
Set size of axis labels.
Definition TStyle.cxx:1440
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:1589
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