Logo ROOT   6.12/07
Reference Guide
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 
52 double 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 }
57 double 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
63 TRandom3 rndm;
64 TH2D *h1, *h2;
65 Int_t npfits;
66 
67 void 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 
108 void 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 
144 int 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");
193  else TVirtualFitter::SetDefaultFitter("Minuit");
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 }
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6063
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
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:481
Random number generator class based on M.
Definition: TRandom3.h:27
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF2.cxx:217
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:100
return c1
Definition: legend1.C:41
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
TH1 * h
Definition: legend2.C:5
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
virtual Double_t GetParameter(Int_t ipar) const =0
int Int_t
Definition: RtypesCore.h:41
void SetStatY(Float_t y=0)
Definition: TStyle.h:371
Double_t x[n]
Definition: legend1.C:17
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
TH1F * h1
Definition: legend1.C:5
virtual Double_t GetParError(Int_t ipar) const =0
virtual const char * GetParName(Int_t ipar) const
Definition: TF1.h:513
virtual void SetChisquare(Double_t chi2)
Definition: TF1.h:596
ROOT::R::TRInterface & r
Definition: Object.C:4
Class to manage histogram axis.
Definition: TAxis.h:30
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
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:1218
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
TAxis * GetYaxis()
Definition: TH1.h:316
A 2-Dim function with parameters.
Definition: TF2.h:29
virtual Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const =0
The Canvas class.
Definition: TCanvas.h:31
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF2.h:150
double Double_t
Definition: RtypesCore.h:55
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes-*.
Definition: TF2.cxx:244
Double_t y[n]
Definition: legend1.C:17
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...
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:290
Abstract Base Class for Fitting.
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1153
virtual void Add(TObject *obj)
Definition: TList.h:87
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:3371
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:496
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8276
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:84
#define gPad
Definition: TVirtualPad.h:285
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:3306
virtual Int_t GetNbinsX() const
Definition: TH1.h:291
double exp(double)
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
TList * GetListOfFunctions() const
Definition: TH1.h:238
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:3688
const Int_t n
Definition: legend1.C:16
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:315
virtual Int_t GetNbinsY() const
Definition: TH1.h:292
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8319
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290