Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
fit2dHist.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_fit
3/// \notebook
4/// Example to fit two histograms at the same time via the Fitter class.
5///
6/// To execute this tutorial, you can do:
7///
8/// ~~~{.cpp}
9/// root > .x fit2dHist.C (executing via cling, slow)
10/// ~~~
11///
12/// or
13/// ~~~{.cpp}
14/// root > .x fit2dHist.C+ (executing via ACLIC , fast, with Minuit)
15/// root > .x fit2dHist.C+(2) (executing via ACLIC , fast, with Minuit2)
16/// ~~~
17///
18/// or using the option to fit independently the 2 histos
19/// ~~~{.cpp}
20/// root > .x fit2dHist.C+(10) (via ACLIC, fast, independent fits with Minuit)
21/// root > .x fit2dHist.C+(12) (via ACLIC, fast, independent fits with Minuit2)
22/// ~~~
23///
24/// Note that you can also execute this script in batch with eg,
25/// ~~~{.cpp}
26/// root -b -q "fit2dHist.C+(12)"
27/// ~~~
28///
29/// or execute interactively from the shell
30/// ~~~{.cpp}
31/// root fit2dHist.C+
32/// root "fit2dHist.C+(12)"
33/// ~~~
34///
35/// \macro_image
36/// \macro_output
37/// \macro_code
38///
39/// \authors: Lorenzo Moneta, Rene Brun 18/01/2006
40
41#include "TH2D.h"
42#include "TF2.h"
43#include "TCanvas.h"
44#include "TStyle.h"
45#include "TRandom3.h"
46#include "Fit/Fitter.h"
47#include "TList.h"
48
49#include <iostream>
50
51double gauss2D(double *x, double *par) {
52 double z1 = double((x[0]-par[1])/par[2]);
53 double z2 = double((x[1]-par[3])/par[4]);
54 return par[0]*exp(-0.5*(z1*z1+z2*z2));
55}
56double my2Dfunc(double *x, double *par) {
57 return gauss2D(x,&par[0]) + gauss2D(x,&par[5]);
58}
59
60
61
62class MyFcn {
63 public:
64 TH2D *h1 = nullptr;
65 TH2D *h2 = nullptr;
66 int npfits = 0;
67
68 MyFcn(TH2D * _h1, TH2D * _h2) :
69 h1(_h1), h2(_h2) {}
70
71 double operator()(const double *p) {
72
73 TAxis *xaxis1 = h1->GetXaxis();
74 TAxis *yaxis1 = h1->GetYaxis();
75 TAxis *xaxis2 = h2->GetXaxis();
76 TAxis *yaxis2 = h2->GetYaxis();
77
78 int nbinX1 = h1->GetNbinsX();
79 int nbinY1 = h1->GetNbinsY();
80 int nbinX2 = h2->GetNbinsX();
81 int nbinY2 = h2->GetNbinsY();
82
83 double chi2 = 0;
84 double x[2];
85 double tmp;
86 npfits = 0;
87 for (int ix = 1; ix <= nbinX1; ++ix) {
88 x[0] = xaxis1->GetBinCenter(ix);
89 for (int iy = 1; iy <= nbinY1; ++iy) {
90 if ( h1->GetBinError(ix,iy) > 0 ) {
91 x[1] = yaxis1->GetBinCenter(iy);
92 tmp = (h1->GetBinContent(ix,iy) - my2Dfunc(x,(double*) p))/h1->GetBinError(ix,iy);
93 chi2 += tmp*tmp;
94 npfits++;
95 }
96 }
97 }
98 for (int ix = 1; ix <= nbinX2; ++ix) {
99 x[0] = xaxis2->GetBinCenter(ix);
100 for (int iy = 1; iy <= nbinY2; ++iy) {
101 if ( h2->GetBinError(ix,iy) > 0 ) {
102 x[1] = yaxis2->GetBinCenter(iy);
103 tmp = (h2->GetBinContent(ix,iy) - my2Dfunc(x,(double *) p))/h2->GetBinError(ix,iy);
104 chi2 += tmp*tmp;
105 npfits++;
106 }
107 }
108 }
109 return chi2;
110}
111};
112
113void FillHisto(TH2D * h, int n, double * p) {
114
115
116 const double mx1 = p[1];
117 const double my1 = p[3];
118 const double sx1 = p[2];
119 const double sy1 = p[4];
120 const double mx2 = p[6];
121 const double my2 = p[8];
122 const double sx2 = p[7];
123 const double sy2 = p[9];
124 //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2);
125 const double w1 = 0.5;
126
127 double x, y;
128 for (int i = 0; i < n; ++i) {
129 // generate randoms with larger Gaussians
130 gRandom->Rannor(x,y);
131
132 double r = gRandom->Rndm(1);
133 if (r < w1) {
134 x = x*sx1 + mx1;
135 y = y*sy1 + my1;
136 }
137 else {
138 x = x*sx2 + mx2;
139 y = y*sy2 + my2;
140 }
141 h->Fill(x,y);
142
143 }
144}
145
146
147
148
149int fit2dHist(int option=1) {
150
151 // create two histograms
152
153 int nbx1 = 50;
154 int nby1 = 50;
155 int nbx2 = 50;
156 int nby2 = 50;
157 double xlow1 = 0.;
158 double ylow1 = 0.;
159 double xup1 = 10.;
160 double yup1 = 10.;
161 double xlow2 = 5.;
162 double ylow2 = 5.;
163 double xup2 = 20.;
164 double yup2 = 20.;
165
166 auto h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
167 auto h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2);
168
169 double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. };
170 // create fit function
171 TF2 * func = new TF2("func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);
172 func->SetParameters(iniParams);
173
174 // fill Histos
175 int n1 = 1000000;
176 int n2 = 1000000;
177 FillHisto(h1,n1,iniParams);
178 FillHisto(h2,n2,iniParams);
179
180 // scale histograms to same heights (for fitting)
181 double dx1 = (xup1-xlow1)/double(nbx1);
182 double dy1 = (yup1-ylow1)/double(nby1);
183 double dx2 = (xup2-xlow2)/double(nbx2);
184 double dy2 = (yup2-ylow2)/double(nby2);
185 // scale histo 2 to scale of 1
186 h2->Sumw2();
187 h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) );
188
189 bool global = false;
190 if (option > 10) global = true;
191 // do global combined fit
192 if (global) {
193 // fill data structure for fit (coordinates + values + errors)
194 std::cout << "Do global fit" << std::endl;
195 // fit now all the function together
196
197 ROOT::Fit::Fitter fitter;
198 //The default minimizer is Minuit, you can also try Minuit2
199
200 MyFcn myFcn(h1,h2);
201 fitter.SetFCN(10, myFcn);
202 if (option%10 == 2) fitter.Config().SetMinimizer("Minuit2");
203
204 // set parameter initial value, name and step size
205 for (int i = 0; i < 10; ++i) {
206 fitter.Config().ParSettings(i) = ROOT::Fit::ParameterSettings(func->GetParName(i), func->GetParameter(i), 0.01);
207 }
208
209 bool ret = fitter.FitFCN();
210 if (!ret) {
211 Error("fit2DHist","Fit Failed to converge");
212 return -1;
213 }
214
215 //get result
216 double minParams[10];
217 double parErrors[10];
218 for (int i = 0; i < 10; ++i) {
219 minParams[i] = fitter.Result().Parameter(i);
220 parErrors[i] = fitter.Result().Error(i);
221 }
222 double chi2 = fitter.Result().MinFcnValue();
223 double edm = fitter.Result().Edm();
224 int npfits = myFcn.npfits;
225
226 func->SetParameters(minParams);
227 func->SetParErrors(parErrors);
228 func->SetChisquare(chi2);
229 int ndf = npfits - fitter.Result().NFreeParameters();
230 func->SetNDF(ndf);
231
232 // add to list of functions
233 h1->GetListOfFunctions()->Add(func);
234 h2->GetListOfFunctions()->Add(func);
235 }
236 else {
237 // fit independently
238 h1->Fit(func, "0");
239 h2->Fit(func, "0");
240 }
241
242 // Create a new canvas.
243 TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800);
244 c1->Divide(2,2);
245 gStyle->SetOptFit();
246 gStyle->SetStatY(0.6);
247
248 c1->cd(1);
249 h1->Draw();
250 func->SetRange(xlow1,ylow1,xup1,yup1);
251 func->DrawCopy("cont3 same");
252 c1->cd(2);
253 h1->Draw("lego");
254 func->DrawCopy("surf1 same");
255 c1->cd(3);
256 func->SetRange(xlow2,ylow2,xup2,yup2);
257 h2->Draw();
258 func->DrawCopy("cont3 same");
259 c1->cd(4);
260 h2->Draw("lego");
261 gPad->SetLogz();
262 func->Draw("surf1 same");
263
264 return 0;
265}
#define h(i)
Definition RSha256.hxx:106
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
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 r
TRObject operator()(const T1 &t1) const
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
#define gPad
void SetMinimizer(const char *type, const char *algo=nullptr)
set minimizer type
Definition FitConfig.h:179
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition FitConfig.h:76
double Error(unsigned int i) const
parameter error by index
Definition FitResult.h:179
double MinFcnValue() const
Return value of the objective function (chi2 or likelihood) used in the fit.
Definition FitResult.h:111
double Edm() const
Expected distance from minimum.
Definition FitResult.h:117
unsigned int NFreeParameters() const
get total number of free parameters
Definition FitResult.h:125
double Parameter(unsigned int i) const
parameter value by index
Definition FitResult.h:174
Fitter class, entry point for performing all type of fits.
Definition Fitter.h:77
const FitResult & Result() const
get fit result
Definition Fitter.h:394
bool FitFCN(unsigned int npar, Function &fcn, const double *params=nullptr, unsigned int dataSize=0, int fitType=0)
Fit using the a generic FCN function as a C++ callable object implementing double () (const double *)...
Definition Fitter.h:649
const FitConfig & Config() const
access to the fit configuration (const method)
Definition Fitter.h:422
bool SetFCN(unsigned int npar, Function &fcn, const double *params=nullptr, unsigned int dataSize=0, int fitType=0)
Set a generic FCN function as a C++ callable object implementing double () (const double *) Note that...
Definition Fitter.h:656
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
Class to manage histogram axis.
Definition TAxis.h:32
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition TAxis.cxx:473
The Canvas class.
Definition TCanvas.h:23
virtual void SetNDF(Int_t ndf)
Set the number of degrees of freedom ndf should be the number of points used in a fit - the number of...
Definition TF1.cxx:3419
virtual void SetChisquare(Double_t chi2)
Definition TF1.h:640
virtual void SetParErrors(const Double_t *errors)
Set errors for all active parameters when calling this function, the array errors must have at least ...
Definition TF1.cxx:3490
virtual const char * GetParName(Int_t ipar) const
Definition TF1.h:557
virtual void SetParameters(const Double_t *params)
Definition TF1.h:677
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:540
A 2-Dim function with parameters.
Definition TF2.h:29
TF1 * DrawCopy(Option_t *option="") const override
Draw a copy of this function with its current attributes-*.
Definition TF2.cxx:286
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition TF2.cxx:259
void SetRange(Double_t xmin, Double_t xmax) override
Initialize the upper and lower bounds to draw the function.
Definition TF2.h:146
virtual Int_t GetNbinsY() const
Definition TH1.h:299
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9096
TAxis * GetXaxis()
Definition TH1.h:325
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:3906
virtual Int_t GetNbinsX() const
Definition TH1.h:298
TAxis * GetYaxis()
Definition TH1.h:326
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3068
TList * GetListOfFunctions() const
Definition TH1.h:245
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5090
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:358
void Add(TObject *obj) override
Definition TList.h:81
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
virtual void Rannor(Float_t &a, Float_t &b)
Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
Definition TRandom.cxx:507
void SetStatY(Float_t y=0)
Definition TStyle.h:398
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:1593
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition RVec.hxx:1837
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5