Logo ROOT  
Reference Guide
TwoHistoFit2D.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.
5///
6/// \macro_image
7/// \macro_output
8/// \macro_code
9///
10/// \author Rene Brun
11
12#include "TH2D.h"
13#include "TF2.h"
14#include "TCanvas.h"
15#include "TStyle.h"
16#include "TRandom3.h"
17#include "TVirtualFitter.h"
18#include "TList.h"
19
20#include <vector>
21#include <map>
22#include <iostream>
23
24double gauss2D(double *x, double *par) {
25 double z1 = double((x[0]-par[1])/par[2]);
26 double z2 = double((x[1]-par[3])/par[4]);
27 return par[0]*exp(-0.5*(z1*z1+z2*z2));
28}
29double my2Dfunc(double *x, double *par) {
30 double *p1 = &par[0];
31 double *p2 = &par[5];
32 return gauss2D(x,p1) + gauss2D(x,p2);
33}
34
35
36
37// data need to be globals to be visible by fcn
38
39std::vector<std::pair<double, double> > coords;
40std::vector<double > values;
41std::vector<double > errors;
42
43void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */ )
44{
45 int n = coords.size();
46 double chi2 = 0;
47 double tmp,x[2];
48 for (int i = 0; i <n; ++i ) {
49 x[0] = coords[i].first;
50 x[1] = coords[i].second;
51 tmp = ( values[i] - my2Dfunc(x,p))/errors[i];
52 chi2 += tmp*tmp;
53 }
54 fval = chi2;
55}
56TRandom3 rndm;
57void FillHisto(TH2D * h, int n, double * p) {
58
59
60 const double mx1 = p[1];
61 const double my1 = p[3];
62 const double sx1 = p[2];
63 const double sy1 = p[4];
64 const double mx2 = p[6];
65 const double my2 = p[8];
66 const double sx2 = p[7];
67 const double sy2 = p[9];
68 //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2);
69 const double w1 = 0.5;
70
71 double x, y;
72 for (int i = 0; i < n; ++i) {
73 // generate randoms with larger gaussians
74 rndm.Rannor(x,y);
75
76 double r = rndm.Rndm(1);
77 if (r < w1) {
78 x = x*sx1 + mx1;
79 y = y*sy1 + my1;
80 }
81 else {
82 x = x*sx2 + mx2;
83 y = y*sy2 + my2;
84 }
85 h->Fill(x,y);
86
87 }
88}
89
90
91
92
93int TwoHistoFit2D(bool global = true) {
94
95 // create two histograms
96
97 int nbx1 = 50;
98 int nby1 = 50;
99 int nbx2 = 50;
100 int nby2 = 50;
101 double xlow1 = 0.;
102 double ylow1 = 0.;
103 double xup1 = 10.;
104 double yup1 = 10.;
105 double xlow2 = 5.;
106 double ylow2 = 5.;
107 double xup2 = 20.;
108 double yup2 = 20.;
109
110 TH2D * h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
111 TH2D * h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2);
112
113 double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. };
114 // create fit function
115 TF2 * func = new TF2("func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);
116 func->SetParameters(iniParams);
117
118 // fill Histos
119 int n1 = 50000;
120 int n2 = 50000;
121 // h1->FillRandom("func", n1);
122 //h2->FillRandom("func",n2);
123 FillHisto(h1,n1,iniParams);
124 FillHisto(h2,n2,iniParams);
125
126 // scale histograms to same heights (for fitting)
127 double dx1 = (xup1-xlow1)/double(nbx1);
128 double dy1 = (yup1-ylow1)/double(nby1);
129 double dx2 = (xup2-xlow2)/double(nbx2);
130 double dy2 = (yup2-ylow2)/double(nby2);
131// h1->Sumw2();
132// h1->Scale( 1.0 / ( n1 * dx1 * dy1 ) );
133 // scale histo 2 to scale of 1
134 h2->Sumw2();
135 h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) );
136
137
138 if (global) {
139 // fill data structure for fit (coordinates + values + errors)
140 std::cout << "Do global fit" << std::endl;
141 // fit now all the function together
142
143 // fill data structure for fit (coordinates + values + errors)
144 TAxis *xaxis1 = h1->GetXaxis();
145 TAxis *yaxis1 = h1->GetYaxis();
146 TAxis *xaxis2 = h2->GetXaxis();
147 TAxis *yaxis2 = h2->GetYaxis();
148
149 int nbinX1 = h1->GetNbinsX();
150 int nbinY1 = h1->GetNbinsY();
151 int nbinX2 = h2->GetNbinsX();
152 int nbinY2 = h2->GetNbinsY();
153
154 /// reset data structure
155 coords = std::vector<std::pair<double,double> >();
156 values = std::vector<double>();
157 errors = std::vector<double>();
158
159
160 for (int ix = 1; ix <= nbinX1; ++ix) {
161 for (int iy = 1; iy <= nbinY1; ++iy) {
162 if ( h1->GetBinContent(ix,iy) > 0 ) {
163 coords.push_back( std::make_pair(xaxis1->GetBinCenter(ix), yaxis1->GetBinCenter(iy) ) );
164 values.push_back( h1->GetBinContent(ix,iy) );
165 errors.push_back( h1->GetBinError(ix,iy) );
166 }
167 }
168 }
169 for (int ix = 1; ix <= nbinX2; ++ix) {
170 for (int iy = 1; iy <= nbinY2; ++iy) {
171 if ( h2->GetBinContent(ix,iy) > 0 ) {
172 coords.push_back( std::make_pair(xaxis2->GetBinCenter(ix), yaxis2->GetBinCenter(iy) ) );
173 values.push_back( h2->GetBinContent(ix,iy) );
174 errors.push_back( h2->GetBinError(ix,iy) );
175 }
176 }
177 }
178
180 TVirtualFitter * minuit = TVirtualFitter::Fitter(0,10);
181 for (int i = 0; i < 10; ++i) {
182 minuit->SetParameter(i, func->GetParName(i), func->GetParameter(i), 0.01, 0,0);
183 }
184 minuit->SetFCN(myFcn);
185
186 double arglist[100];
187 arglist[0] = 0;
188 // set print level
189 minuit->ExecuteCommand("SET PRINT",arglist,2);
190
191 // minimize
192 arglist[0] = 5000; // number of function calls
193 arglist[1] = 0.01; // tolerance
194 minuit->ExecuteCommand("MIGRAD",arglist,2);
195
196 //get result
197 double minParams[10];
198 double parErrors[10];
199 for (int i = 0; i < 10; ++i) {
200 minParams[i] = minuit->GetParameter(i);
201 parErrors[i] = minuit->GetParError(i);
202 }
203 double chi2, edm, errdef;
204 int nvpar, nparx;
205 minuit->GetStats(chi2,edm,errdef,nvpar,nparx);
206
207 func->SetParameters(minParams);
208 func->SetParErrors(parErrors);
209 func->SetChisquare(chi2);
210 int ndf = coords.size()-nvpar;
211 func->SetNDF(ndf);
212
213 std::cout << "Chi2 Fit = " << chi2 << " ndf = " << ndf << " " << func->GetNDF() << std::endl;
214
215 // add to list of functions
216 h1->GetListOfFunctions()->Add(func);
217 h2->GetListOfFunctions()->Add(func);
218 }
219 else {
220 // fit independently
221 h1->Fit(func);
222 h2->Fit(func);
223 }
224
225
226
227 // Create a new canvas.
228 TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800);
229 c1->Divide(2,2);
230 gStyle->SetOptFit();
231 gStyle->SetStatY(0.6);
232
233 c1->cd(1);
234 h1->Draw();
235 func->SetRange(xlow1,ylow1,xup1,yup1);
236 func->DrawCopy("cont1 same");
237 c1->cd(2);
238 h1->Draw("lego");
239 func->DrawCopy("surf1 same");
240 c1->cd(3);
241 func->SetRange(xlow2,ylow2,xup2,yup2);
242 h2->Draw();
243 func->DrawCopy("cont1 same");
244 c1->cd(4);
245 h2->Draw("lego");
246 gPad->SetLogz();
247 func->Draw("surf1 same");
248
249 return 0;
250}
double
Definition: Converters.cxx:921
ROOT::R::TRInterface & r
Definition: Object.C:4
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:43
double Double_t
Definition: RtypesCore.h:57
double exp(double)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
#define gPad
Definition: TVirtualPad.h:287
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:475
The Canvas class.
Definition: TCanvas.h:27
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 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:3421
virtual void SetChisquare(Double_t chi2)
Definition: TF1.h:606
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:3486
virtual const char * GetParName(Int_t ipar) const
Definition: TF1.h:523
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:506
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:268
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF2.h:150
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF2.cxx:241
virtual Int_t GetNbinsY() const
Definition: TH1.h:293
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
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 GetNbinsX() const
Definition: TH1.h:292
TAxis * GetYaxis()
Definition: TH1.h:317
TList * GetListOfFunctions() const
Definition: TH1.h:239
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6246
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8476
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:87
Random number generator class based on M.
Definition: TRandom3.h:27
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:100
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:489
void SetStatY(Float_t y=0)
Definition: TStyle.h:379
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
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