ROOT  6.06/09
Reference Guide
testRooFit.cxx
Go to the documentation of this file.
1 // test fitting using also RooFit and new Fitter
2 #include <RooAbsPdf.h>
3 #include <RooRealVar.h>
4 #include <RooArgSet.h>
5 #include <RooGaussian.h>
6 
7 #include "RooDataSet.h"
8 #include "RooRealVar.h"
9 #include "RooGaussian.h"
10 #include "RooGlobalFunc.h"
11 #include "RooFitResult.h"
12 #include "RooProdPdf.h"
13 
14 #include <TF1.h>
15 #include <TTree.h>
16 #include <TRandom3.h>
17 
18 #include "TStopwatch.h"
19 
20 #include "Math/DistFunc.h"
21 
22 #include "Fit/UnBinData.h"
23 #include "Fit/BinData.h"
24 #include "Fit/Fitter.h"
25 
26 
27 
28 #include "WrapperRooPdf.h"
29 
30 #include <string>
31 #include <iostream>
32 #include <vector>
33 #include <memory>
34 
35 #include "MinimizerTypes.h"
36 
38 #include <cmath>
39 
40 const int N = 3; // n must be greater than 1
41 const int nfit = 1;
42 const int nEvents = 10000;
43 double iniPar[2*N];
44 
45 //#define DEBUG
46 
48 
49 void fillTree(TTree & t2) {
50 
51 
52  double x[N];
53  Int_t ev;
54  for (int j = 0; j < N; ++j) {
55  std::string xname = "x_" + ROOT::Math::Util::ToString(j);
56  std::string xname2 = "x_" + ROOT::Math::Util::ToString(j) + "/D";
57  t2.Branch(xname.c_str(),&x[j],xname2.c_str());
58  }
59  t2.Branch("ev",&ev,"ev/I");
60  //fill the tree
61  TRandom3 r;
62  for (Int_t i=0;i<nEvents;i++) {
63  for (int j = 0; j < N; ++j) {
64  double mu = double(j)/10.;
65  double s = 1.0 + double(j)/10.;
66  x[j] = r.Gaus(mu,s);
67  }
68 
69  ev = i;
70  t2.Fill();
71  }
72  // t2.Print();
73  // std::cout << "number of branches " << t2.GetNbranches() << std::endl;
74 
76 }
77 
79 
80  // fill the unbin data set from a TTree
81 
82  // large tree
83  unsigned int n = tree->GetEntries();
84  std::cout << "number of unbin data is " << n << " of dim " << N << std::endl;
85  d.Initialize(n,N);
86 
87  double vx[N];
88  for (int j = 0; j <N; ++j) {
89  std::string bname = "x_" + ROOT::Math::Util::ToString(j);
90  TBranch * bx = tree->GetBranch(bname.c_str());
91  bx->SetAddress(&vx[j]);
92  }
93 
94  std::vector<double> m(N);
95  for (int unsigned i = 0; i < n; ++i) {
96  tree->GetEntry(i);
97  d.Add(vx);
98  for (int j = 0; j < N; ++j)
99  m[j] += vx[j];
100  }
101 
102 #ifdef DEBUG
103  std::cout << "average values of means :\n";
104  for (int j = 0; j < N; ++j)
105  std::cout << m[j]/n << " ";
106  std::cout << "\n";
107 #endif
108 
109  tree->ResetBranchAddresses() ;
110 
111  return;
112 
113 }
114 
115 
116 
117 // class describing product of gaussian pdf
118 class MultiGaussRooPdf {
119 
120  public:
121 
122  // create all pdf
123  MultiGaussRooPdf(unsigned int n) :
124  x(n), m(n), s(n), g(n), pdf(n)
125  {
126  //assert(n >= 2);
127 
128  // create the gaussians
129  for (unsigned int j = 0; j < n; ++j) {
130 
131  std::string xname = "x_" + ROOT::Math::Util::ToString(j);
132  x[j] = new RooRealVar(xname.c_str(),xname.c_str(),-10000,10000) ;
133 
134  std::string mname = "m_" + ROOT::Math::Util::ToString(j);
135  std::string sname = "s_" + ROOT::Math::Util::ToString(j);
136 
137 
138 // m[j] = new RooRealVar(mname.c_str(),mname.c_str(),iniPar[2*j],-10000,10000) ;
139 // s[j] = new RooRealVar(sname.c_str(),sname.c_str(),iniPar[2*j+1],-10000,10000) ;
140  m[j] = new RooRealVar(mname.c_str(),mname.c_str(),iniPar[2*j],-10,10) ;
141  s[j] = new RooRealVar(sname.c_str(),sname.c_str(),iniPar[2*j+1],-10,10) ;
142 
143  std::string gname = "g_" + ROOT::Math::Util::ToString(j);
144  g[j] = new RooGaussian(gname.c_str(),"gauss(x,mean,sigma)",*x[j],*m[j],*s[j]);
145 
146  std::string pname = "prod_" + ROOT::Math::Util::ToString(j);
147  if (j == 0)
148  pdf[0] = g[0];
149  else if (j == 1) {
150  pdf[1] = new RooProdPdf(pname.c_str(),pname.c_str(),RooArgSet(*g[1],*g[0]) );
151  }
152  else if (j > 1) {
153  pdf[j] = new RooProdPdf(pname.c_str(),pname.c_str(),RooArgSet(*g[j],*pdf[j-1]) );
154  }
155  }
156 
157 
158  }
159 
160  RooAbsPdf & getPdf() { return *pdf.back(); }
161 
162  std::auto_ptr<RooArgSet> getVars() {
163  std::auto_ptr<RooArgSet> vars(new RooArgSet() );
164  for (unsigned int i = 0; i < x.size(); ++i)
165  vars->add(*x[i]);
166  return vars;
167  }
168 
169  ~MultiGaussRooPdf() {
170  // free
171  int n = x.size();
172  for (int j = 0; j < n; ++j) {
173  delete x[j];
174  delete m[j];
175  delete s[j];
176  delete g[j];
177  if (j> 0) delete pdf[j]; // no pdf allocated for j = 0
178  }
179  }
180 
181  private:
182 
183  std::vector<RooRealVar *> x;
184  std::vector<RooRealVar *> m;
185  std::vector<RooRealVar *> s;
186 
187  std::vector<RooAbsPdf *> g;
188  std::vector<RooAbsPdf *> pdf;
189 
190 };
191 
192 
193 //unbinned roo fit (large tree)
194 int FitUsingRooFit(TTree & tree, RooAbsPdf & pdf, RooArgSet & xvars) {
195 
196  int iret = 0;
197  std::cout << "\n************************************************************\n";
198  std::cout << "\tFit using RooFit (Likelihood Fit) on " << pdf.GetName() << std::endl;
199 
200 
201 
202  TStopwatch w;
203  w.Start();
204 
205  for (int i = 0; i < nfit; ++i) {
206 
207  RooDataSet data("unbindata","unbin dataset with x",&tree,xvars) ;
208 
209 
210 
211 #ifdef DEBUG
212  int level = 2;
213  std::cout << "num entries = " << data.numEntries() << std::endl;
214  bool save = true;
215  pdf.getVariables()->Print("v"); // print the parameters
216  std::cout << "\n\nDo the fit now \n\n";
217 #else
218  int level = -1;
219  bool save = false;
220 #endif
221 
222 #ifndef _WIN32 // until a bug 30762 is fixed
223  RooFitResult * result = pdf.fitTo(data, RooFit::Minos(0), RooFit::Hesse(0) , RooFit::PrintLevel(level), RooFit::Save(save) );
224 #else
225  RooFitResult * result = pdf.fitTo(data );
226 #endif
227 
228 #ifdef DEBUG
229  assert(result != 0);
230  std::cout << " Roofit status " << result->status() << std::endl;
231  result->Print();
232 #endif
233  iret |= (result == 0);
234 
235  }
236 
237  w.Stop();
238 
239  std::cout << "RooFit result " << std::endl;
240  RooArgSet * params = pdf.getParameters(xvars);
241  params->Print("v");
242  delete params;
243 
244 
245  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
246  std::cout << "\n************************************************************\n";
247  return iret;
248 }
249 
250 
251 // unbin fit
252 template <class MinType>
253 int DoFit(TTree * tree, Func & func, bool debug = false, bool = false ) {
254 
256  // need to have done Tree->Draw() before fit
257  FillUnBinData(d,tree);
258 
259  //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl;
260 
261  ROOT::Fit::Fitter fitter;
262  fitter.Config().MinimizerOptions().SetPrintLevel(0);
263  fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str());
264  fitter.Config().MinimizerOptions().SetTolerance(1.); // to be consistent with RooFit
265 
266  if (debug)
267  fitter.Config().MinimizerOptions().SetPrintLevel(3);
268 
269 
270  // create the function
271 
272  fitter.SetFunction(func);
273  // need to fix param 0 , normalization in the unbinned fits
274  //fitter.Config().ParSettings(0).Fix();
275 
276  bool ret = fitter.Fit(d);
277  if (!ret) {
278  std::cout << " Fit Failed " << std::endl;
279  return -1;
280  }
281  fitter.Result().Print(std::cout);
282  return 0;
283 
284 }
285 
286 template <class MinType, class FitObj>
287 int FitUsingNewFitter(FitObj * fitobj, Func & func, bool useGrad=false) {
288 
289  std::cout << "\n************************************************************\n";
290  std::cout << "\tFit using new Fit::Fitter\n";
291  std::cout << "\tMinimizer is " << MinType::name() << " " << MinType::name2() << std::endl;
292 
293  int iret = 0;
294  TStopwatch w; w.Start();
295 
296 #ifdef DEBUG
297  func.SetParameters(iniPar);
298  iret |= DoFit<MinType>(fitobj,func,true, useGrad);
299 
300 #else
301  for (int i = 0; i < nfit; ++i) {
302  func.SetParameters(iniPar);
303  iret = DoFit<MinType>(fitobj,func, false, useGrad);
304  if (iret != 0) break;
305  }
306 #endif
307  w.Stop();
308  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
309  std::cout << "\n************************************************************\n";
310 
311  return iret;
312 }
313 
314 template <int N>
315 struct GausNorm {
316  static double F(const double *x, const double *p) {
317  return ROOT::Math::normal_pdf(x[N-1], p[2*N-1], p[2*N-2] ) * GausNorm<N-1>::F(x,p);
318  }
319 };
320 template <>
321 struct GausNorm<1> {
322 
323  static double F(const double *x, const double *p) {
324  return ROOT::Math::normal_pdf(x[0], p[1], p[0] );
325  }
326 };
327 
328 
329 // double gausnorm(
330 // return ROOT::Math::normal_pdf(x[0], p[1], p[0] );
331 // // //return p[0]*TMath::Gaus(x[0],p[1],p[2]);
332 // // double invsig = 1./p[1];
333 // // double tmp = (x[0]-p[0]) * invsig;
334 // // const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 );
335 // // return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig;
336 // }
337 
338 
339 
340 int main() {
341 
342  TTree tree("t","a large Tree with many gaussian variables");
343  fillTree(tree);
344 
345  // set initial parameters
346  for (int i = 0; i < N; ++i) {
347  iniPar[2*i] = 1;
348  iniPar[2*i+1] = 1;
349  }
350 
351  // start counting the time
352  MultiGaussRooPdf multipdf(N);
353  RooAbsPdf & pdf = multipdf.getPdf();
354 
355  std::auto_ptr<RooArgSet> xvars = multipdf.getVars();
356 
357  WrapperRooPdf wpdf( &pdf, *xvars );
358 
359 
360  std::cout << "ndim " << wpdf.NDim() << std::endl;
361  std::cout << "npar " << wpdf.NPar() << std::endl;
362  for (unsigned int i = 0; i < wpdf.NPar(); ++i)
363  std::cout << " par " << i << " is " << wpdf.ParameterName(i) << " value " << wpdf.Parameters()[i] << std::endl;
364 
365 
366  FitUsingNewFitter<TMINUIT>(&tree,wpdf);
367 
368  // reset pdf original values
369  wpdf.SetParameters(iniPar);
370 
371  FitUsingRooFit(tree,pdf, *xvars);
372 
373  // in case of N = 1 do also a simple gauss fit
374  // using TF1 gausN
375 // if (N == 1) {
376  ROOT::Math::WrappedParamFunction<> gausn(&GausNorm<N>::F,2*N,iniPar,iniPar+2*N);
377  FitUsingNewFitter<TMINUIT>(&tree,gausn);
378  //}
379 
380 
381 
382 }
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:382
virtual void Print(Option_t *options=0) const
Print TNamed name and title.
Definition: RooFitResult.h:65
void SetPrintLevel(int level)
set print level
Random number generator class based on M.
Definition: TRandom3.h:29
const double * Parameters() const
Access the parameter values.
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
virtual void SetAddress(void *add)
Set address of this branch.
Definition: TBranch.cxx:2048
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
void SetTolerance(double tol)
set the tolerance
unsigned int NPar() const
Return the number of Parameters.
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
RooCmdArg PrintLevel(Int_t code)
RooCmdArg Minos(Bool_t flag=kTRUE)
#define assert(cond)
Definition: unittest.h:542
int FitUsingNewFitter(FitObj *fitobj, Func &func, bool useGrad=false)
Definition: testRooFit.cxx:287
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4328
int main()
Definition: testRooFit.cxx:340
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition: TTree.cxx:5163
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:47
int Int_t
Definition: RtypesCore.h:41
ROOT::Math::MinimizerOptions & MinimizerOptions()
access to the minimizer control parameter (non const method)
Definition: FitConfig.h:138
std::string ParameterName(unsigned int i) const
Return the name of the i-th parameter (starting from zero) Overwrite if want to avoid the default nam...
double normal_pdf(double x, double sigma=1, double x0=0)
Probability density function of the normal (Gaussian) distribution.
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Double_t x[n]
Definition: legend1.C:17
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2081
void SetParameters(const double *p)
set parameters.
const FitResult & Result() const
get fit result
Definition: Fitter.h:354
const int nfit
Definition: testRooFit.cxx:41
ROOT::Math::IParamMultiFunction Func
Definition: testRooFit.cxx:47
const int N
Definition: testRooFit.cxx:40
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:559
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition: TTree.cxx:4822
const int nEvents
Definition: testRooFit.cxx:42
#define F(x, y, z)
unsigned int NDim() const
Retrieve the dimension of the function.
bool Fit(const Data &data, const Function &func)
fit a data set using any generic model function If data set is binned a least square fit is performed...
Definition: Fitter.h:146
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:152
ROOT::R::TRInterface & r
Definition: Object.C:4
IParamFunction interface (abstract class) describing multi-dimensional parameteric functions It is a ...
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:94
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:291
void fillTree(TTree &t2)
Definition: testRooFit.cxx:49
TMarker * m
Definition: textangle.C:8
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
void SetFunction(const IModelFunction &func, bool useGradient=false)
Set the fitted function (model function) from a parametric function interface.
Definition: Fitter.cxx:99
virtual void SetParameters(const double *p)=0
Set the parameter values.
Int_t status() const
Definition: RooFitResult.h:76
RooCmdArg Hesse(Bool_t flag=kTRUE)
int FitUsingRooFit(TTree &tree, RooAbsPdf &pdf, RooArgSet &xvars)
Definition: testRooFit.cxx:194
void Add(double x)
add one dim coordinate data (unweighted)
Definition: UnBinData.h:199
virtual void ResetBranchAddresses()
Tell all of our branches to drop their current objects and allocate new ones.
Definition: TTree.cxx:7285
double func(double *x, double *p)
Definition: stressTF1.cxx:213
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
#define name(a, b)
Definition: linkTestLib0.cpp:5
virtual Int_t Branch(TCollection *list, Int_t bufsize=32000, Int_t splitlevel=99, const char *name="")
Create one branch for each element in the collection.
Definition: TTree.cxx:1624
RooCmdArg Save(Bool_t flag=kTRUE)
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:42
bool debug
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1056
int DoFit(TTree *tree, Func &func, bool debug=false, bool=false)
Definition: testRooFit.cxx:253
virtual Long64_t GetEntries() const
Definition: TTree.h:382
A TTree object has a header with a name and a title.
Definition: TTree.h:94
double result[121]
A TTree is a list of TBranches.
Definition: TBranch.h:58
const Int_t n
Definition: legend1.C:16
void FillUnBinData(ROOT::Fit::UnBinData &d, TTree *tree)
Definition: testRooFit.cxx:78
double iniPar[2 *N]
Definition: testRooFit.cxx:43
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
void Initialize(unsigned int maxpoints, unsigned int dim=1, bool isWeighted=false)
preallocate a data set given size and dimension of the coordinates if a vector already exists with co...
Definition: UnBinData.cxx:177
Stopwatch class.
Definition: TStopwatch.h:30