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///
5/// Example to fit two histograms at the same time via TVirtualFitter
6///
7/// To execute this tutorial, you can do:
8///
9/// ~~~{.cpp}
10/// root > .x fit2dHist.C (executing via CINT, slow)
11/// ~~~
12///
13/// or
14/// ~~~{.cpp}
15/// root > .x fit2dHist.C+ (executing via ACLIC , fast, with Minuit)
16/// root > .x fit2dHist.C+(2) (executing via ACLIC , fast, with Minuit2)
17/// ~~~
18///
19/// or using the option to fit independently the 2 histos
20/// ~~~{.cpp}
21/// root > .x fit2dHist.C+(10) (via ACLIC, fast, independent fits with Minuit)
22/// root > .x fit2dHist.C+(12) (via ACLIC, fast, independent fits with Minuit2)
23/// ~~~
24///
25/// Note that you can also execute this script in batch with eg,
26/// ~~~{.cpp}
27/// root -b -q "fit2dHist.C+(12)"
28/// ~~~
29///
30/// or execute interactively from the shell
31/// ~~~{.cpp}
32/// root fit2dHist.C+
33/// root "fit2dHist.C+(12)"
34/// ~~~
35///
36/// \macro_image
37/// \macro_output
38/// \macro_code
39///
40/// \authors: Lorenzo Moneta, Rene Brun 18/01/2006
41
42#include "TH2D.h"
43#include "TF2.h"
44#include "TCanvas.h"
45#include "TStyle.h"
46#include "TRandom3.h"
47#include "TVirtualFitter.h"
48#include "TList.h"
49
50#include <iostream>
51
52double gauss2D(double *x, double *par) {
53 double z1 = double((x[0]-par[1])/par[2]);
54 double z2 = double((x[1]-par[3])/par[4]);
55 return par[0]*exp(-0.5*(z1*z1+z2*z2));
56}
57double my2Dfunc(double *x, double *par) {
58 return gauss2D(x,&par[0]) + gauss2D(x,&par[5]);
59}
60
61
62// data need to be globals to be visible by fcn
63TRandom3 rndm;
64TH2D *h1, *h2;
65Int_t npfits;
66
67void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */ )
68{
69 TAxis *xaxis1 = h1->GetXaxis();
70 TAxis *yaxis1 = h1->GetYaxis();
71 TAxis *xaxis2 = h2->GetXaxis();
72 TAxis *yaxis2 = h2->GetYaxis();
73
74 int nbinX1 = h1->GetNbinsX();
75 int nbinY1 = h1->GetNbinsY();
76 int nbinX2 = h2->GetNbinsX();
77 int nbinY2 = h2->GetNbinsY();
78
79 double chi2 = 0;
80 double x[2];
81 double tmp;
82 npfits = 0;
83 for (int ix = 1; ix <= nbinX1; ++ix) {
84 x[0] = xaxis1->GetBinCenter(ix);
85 for (int iy = 1; iy <= nbinY1; ++iy) {
86 if ( h1->GetBinError(ix,iy) > 0 ) {
87 x[1] = yaxis1->GetBinCenter(iy);
88 tmp = (h1->GetBinContent(ix,iy) - my2Dfunc(x,p))/h1->GetBinError(ix,iy);
89 chi2 += tmp*tmp;
90 npfits++;
91 }
92 }
93 }
94 for (int ix = 1; ix <= nbinX2; ++ix) {
95 x[0] = xaxis2->GetBinCenter(ix);
96 for (int iy = 1; iy <= nbinY2; ++iy) {
97 if ( h2->GetBinError(ix,iy) > 0 ) {
98 x[1] = yaxis2->GetBinCenter(iy);
99 tmp = (h2->GetBinContent(ix,iy) - my2Dfunc(x,p))/h2->GetBinError(ix,iy);
100 chi2 += tmp*tmp;
101 npfits++;
102 }
103 }
104 }
105 fval = chi2;
106}
107
108void FillHisto(TH2D * h, int n, double * p) {
109
110
111 const double mx1 = p[1];
112 const double my1 = p[3];
113 const double sx1 = p[2];
114 const double sy1 = p[4];
115 const double mx2 = p[6];
116 const double my2 = p[8];
117 const double sx2 = p[7];
118 const double sy2 = p[9];
119 //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2);
120 const double w1 = 0.5;
121
122 double x, y;
123 for (int i = 0; i < n; ++i) {
124 // generate randoms with larger gaussians
125 rndm.Rannor(x,y);
126
127 double r = rndm.Rndm(1);
128 if (r < w1) {
129 x = x*sx1 + mx1;
130 y = y*sy1 + my1;
131 }
132 else {
133 x = x*sx2 + mx2;
134 y = y*sy2 + my2;
135 }
136 h->Fill(x,y);
137
138 }
139}
140
141
142
143
144int fit2dHist(int option=1) {
145
146 // create two histograms
147
148 int nbx1 = 50;
149 int nby1 = 50;
150 int nbx2 = 50;
151 int nby2 = 50;
152 double xlow1 = 0.;
153 double ylow1 = 0.;
154 double xup1 = 10.;
155 double yup1 = 10.;
156 double xlow2 = 5.;
157 double ylow2 = 5.;
158 double xup2 = 20.;
159 double yup2 = 20.;
160
161 h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
162 h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2);
163
164 double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. };
165 // create fit function
166 TF2 * func = new TF2("func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);
167 func->SetParameters(iniParams);
168
169 // fill Histos
170 int n1 = 1000000;
171 int n2 = 1000000;
172 FillHisto(h1,n1,iniParams);
173 FillHisto(h2,n2,iniParams);
174
175 // scale histograms to same heights (for fitting)
176 double dx1 = (xup1-xlow1)/double(nbx1);
177 double dy1 = (yup1-ylow1)/double(nby1);
178 double dx2 = (xup2-xlow2)/double(nbx2);
179 double dy2 = (yup2-ylow2)/double(nby2);
180 // scale histo 2 to scale of 1
181 h2->Sumw2();
182 h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) );
183
184 bool global = false;
185 if (option > 10) global = true;
186 if (global) {
187 // fill data structure for fit (coordinates + values + errors)
188 std::cout << "Do global fit" << std::endl;
189 // fit now all the function together
190
191 //The default minimizer is Minuit, you can also try Minuit2
192 if (option%10 == 2) TVirtualFitter::SetDefaultFitter("Minuit2");
194 TVirtualFitter * minuit = TVirtualFitter::Fitter(0,10);
195 for (int i = 0; i < 10; ++i) {
196 minuit->SetParameter(i, func->GetParName(i), func->GetParameter(i), 0.01, 0,0);
197 }
198 minuit->SetFCN(myFcn);
199
200 double arglist[100];
201 arglist[0] = 0;
202 // set print level
203 minuit->ExecuteCommand("SET PRINT",arglist,2);
204
205 // minimize
206 arglist[0] = 5000; // number of function calls
207 arglist[1] = 0.01; // tolerance
208 minuit->ExecuteCommand("MIGRAD",arglist,2);
209
210 //get result
211 double minParams[10];
212 double parErrors[10];
213 for (int i = 0; i < 10; ++i) {
214 minParams[i] = minuit->GetParameter(i);
215 parErrors[i] = minuit->GetParError(i);
216 }
217 double chi2, edm, errdef;
218 int nvpar, nparx;
219 minuit->GetStats(chi2,edm,errdef,nvpar,nparx);
220
221 func->SetParameters(minParams);
222 func->SetParErrors(parErrors);
223 func->SetChisquare(chi2);
224 int ndf = npfits-nvpar;
225 func->SetNDF(ndf);
226
227 // add to list of functions
228 h1->GetListOfFunctions()->Add(func);
229 h2->GetListOfFunctions()->Add(func);
230 }
231 else {
232 // fit independently
233 h1->Fit(func);
234 h2->Fit(func);
235 }
236
237 // Create a new canvas.
238 TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800);
239 c1->Divide(2,2);
240 gStyle->SetOptFit();
241 gStyle->SetStatY(0.6);
242
243 c1->cd(1);
244 h1->Draw();
245 func->SetRange(xlow1,ylow1,xup1,yup1);
246 func->DrawCopy("cont1 same");
247 c1->cd(2);
248 h1->Draw("lego");
249 func->DrawCopy("surf1 same");
250 c1->cd(3);
251 func->SetRange(xlow2,ylow2,xup2,yup2);
252 h2->Draw();
253 func->DrawCopy("cont1 same");
254 c1->cd(4);
255 h2->Draw("lego");
256 gPad->SetLogz();
257 func->Draw("surf1 same");
258
259 return 0;
260}
double
ROOT::R::TRInterface & r
Definition Object.C:4
#define h(i)
Definition RSha256.hxx:106
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
R__EXTERN TStyle * gStyle
Definition TStyle.h:413
#define gPad
Class to manage histogram axis.
Definition TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition TAxis.cxx:478
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:3439
virtual void SetChisquare(Double_t chi2)
Definition TF1.h:612
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:3504
virtual const char * GetParName(Int_t ipar) const
Definition TF1.h:529
virtual void SetParameters(const Double_t *params)
Definition TF1.h:644
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:512
A 2-Dim function with parameters.
Definition TF2.h:29
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes-*.
Definition TF2.cxx:287
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition TF2.h:148
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition TF2.cxx:260
virtual Int_t GetNbinsY() const
Definition TH1.h:297
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:8893
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320
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:3893
virtual Int_t GetNbinsX() const
Definition TH1.h:296
TAxis * GetYaxis()
Definition TH1.h:321
TList * GetListOfFunctions() const
Definition TH1.h:243
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3074
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4994
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6553
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition TH1.cxx:8850
2-D histogram with a double per channel (see TH1 documentation)}
Definition TH2.h:292
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH2.h:88
virtual void Add(TObject *obj)
Definition TList.h:81
Random number generator class based on M.
Definition TRandom3.h:27
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom3.cxx:99
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:500
void SetStatY(Float_t y=0)
Definition TStyle.h:381
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:1541
Abstract Base Class for Fitting.
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
virtual Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const =0
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
virtual Double_t GetParError(Int_t ipar) const =0
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
virtual Double_t GetParameter(Int_t ipar) const =0
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
return c1
Definition legend1.C:41
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5