Logo ROOT   6.08/07
Reference Guide
testFit.cxx
Go to the documentation of this file.
1 #include "TH1.h"
2 #include "TH2.h"
3 #include "TF1.h"
4 #include "TF2.h"
5 #include "TGraphErrors.h"
6 #include "TGraphAsymmErrors.h"
7 #include "TGraph2DErrors.h"
8 #include "TSystem.h"
9 #include "TRandom3.h"
10 #include "TROOT.h"
11 #include "TVirtualFitter.h"
12 
13 #include "Fit/BinData.h"
14 #include "Fit/UnBinData.h"
15 #include "HFitInterface.h"
16 #include "Fit/Fitter.h"
17 
18 #include "Math/WrappedMultiTF1.h"
20 #include "Math/WrappedTF1.h"
21 //#include "Math/Polynomial.h"
22 #include "RConfigure.h"
23 
24 #include <string>
25 #include <iostream>
26 #include <cmath>
27 
28 // print the data
30  std::cout << "Bin data, point size is " << data.PointSize() << " data dimension is " << data.NDim() << " type is " << data.GetErrorType() << std::endl;
31  for (unsigned int i = 0; i < data.Size(); ++i) {
33  std::cout << data.Coords(i)[0] << " " << data.Value(i) << std::endl;
35  std::cout << data.Coords(i)[0] << " " << data.Value(i) << " +/- " << data.Error(i) << std::endl;
36  else if (data.GetErrorType() == ROOT::Fit::BinData::kCoordError) {
37  double ey = 0;
38  const double * ex = data.GetPointError(i,ey);
39  std::cout << data.Coords(i)[0] << " +/- " << ex[0] << " " << data.Value(i) << " +/- " << ey << std::endl;
40  }
41  else if (data.GetErrorType() == ROOT::Fit::BinData::kAsymError) {
42  double eyl = 0; double eyh = 0;
43  const double * ex = data.GetPointError(i,eyl,eyh);
44  std::cout << data.Coords(i)[0] << " +/- " << ex[0] << " " << data.Value(i) << " - " << eyl << " + " << eyh << std::endl;
45  }
46 
47  }
48  std::cout << "\ndata size is " << data.Size() << std::endl;
49 }
51  for (unsigned int i = 0; i < data.Size(); ++i) {
52  std::cout << data.Coords(i)[0] << "\t";
53  }
54  std::cout << "\ndata size is " << data.Size() << std::endl;
55 }
56 
57 int compareResult(double v1, double v2, std::string s = "", double tol = 0.01) {
58  // compare v1 with reference v2
59  // give 1% tolerance
60  if (std::abs(v1-v2) < tol * std::abs(v2) ) return 0;
61  std::cerr << s << " Failed comparison of fit results \t chi2 = " << v1 << " it should be = " << v2 << std::endl;
62  return -1;
63 }
64 
65 double chi2FromFit(const TF1 * func ) {
66  // return last chi2 obtained from Fit method function
67  assert(TVirtualFitter::GetFitter() != 0 );
68  return (TVirtualFitter::GetFitter()->Chisquare(func->GetNpar(), func->GetParameters() ) );
69 }
70 
72 
73 
74  std::string fname("gaus");
75  TF1 * func = (TF1*)gROOT->GetFunction(fname.c_str());
76  func->SetParameter(0,10);
77  func->SetParameter(1,0);
78  func->SetParameter(2,3.0);
79 
80  TRandom3 rndm;
81  int iret = 0;
82  double chi2ref = 0;
83 
84  // fill an histogram
85  TH1D * h1 = new TH1D("h1","h1",30,-5.,5.);
86 // h1->FillRandom(fname.c_str(),100);
87  for (int i = 0; i <1000; ++i)
88  h1->Fill( rndm.Gaus(0,1) );
89 
90  h1->Print();
91  //h1->Draw();
92 
93 // gSystem->Load("libMinuit2");
94 // gSystem->Load("libFit");
95 
96 
97  // ROOT::Fit::DataVector<ROOT::Fit::BinPoint> dv;
98 
100  ROOT::Fit::FillData(d,h1,func);
101 
102 
103  printData(d);
104 
105  // create the function
106  ROOT::Math::WrappedMultiTF1 wf(*func);
107  // need to do that to avoid gradient calculation
109 
110  double p[3] = {100,0,3.};
111  f.SetParameters(p);
112 
113  // create the fitter
114 
115  ROOT::Fit::Fitter fitter;
116 
117  bool ret = fitter.Fit(d, f);
118  if (ret)
119  fitter.Result().Print(std::cout);
120  else {
121  std::cout << "Chi2 Fit Failed " << std::endl;
122  return -1;
123  }
124  chi2ref = fitter.Result().Chi2();
125 
126  // compare with TH1::Fit
128  std::cout << "\n******************************\n\t TH1::Fit Result \n" << std::endl;
129  func->SetParameters(p);
130  h1->Fit(func);
131 
132  iret |= compareResult( chi2FromFit(func), chi2ref,"1D histogram chi2 fit");
133 
134 
135  // test using binned likelihood
136  std::cout << "\n\nTest Binned Likelihood Fit" << std::endl;
137 
139  opt.fUseEmpty = true;
140  ROOT::Fit::BinData dl(opt);
141  ROOT::Fit::FillData(dl,h1,func);
142 
143  ret = fitter.LikelihoodFit(dl, f, true);
144  f.SetParameters(p);
145  if (ret)
146  fitter.Result().Print(std::cout);
147  else {
148  std::cout << "Binned Likelihood Fit Failed " << std::endl;
149  return -1;
150  }
151  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D histogram likelihood fit",0.3);
152 
153  // compare with TH1::Fit
154  std::cout << "\n******************************\n\t TH1::Fit Result \n" << std::endl;
155  func->SetParameters(p);
156  h1->Fit(func,"L");
157 
158  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TH1::Fit likelihood ",0.001);
159  //std::cout << "Equivalent Chi2 from TF1::Fit " << func->GetChisquare() << std::endl;
160 
161  std::cout << "\n\nTest Chi2 Fit using integral option" << std::endl;
162 
163  // need to re-create data
164  opt = ROOT::Fit::DataOptions();
165  opt.fIntegral = true;
166  ROOT::Fit::BinData d2(opt);
167  ROOT::Fit::FillData(d2,h1,func);
168 
169  f.SetParameters(p);
170  ret = fitter.Fit(d2, f);
171  if (ret)
172  fitter.Result().Print(std::cout);
173  else {
174  std::cout << "Integral Chi2 Fit Failed " << std::endl;
175  return -1;
176  }
177  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D histogram integral chi2 fit",0.2);
178 
179  // compare with TH1::Fit
180  std::cout << "\n******************************\n\t TH1::Fit Result \n" << std::endl;
181  func->SetParameters(p);
182  h1->Fit(func,"I");
183  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TH1::Fit integral ",0.001);
184 
185  // test integral likelihood
186  opt = ROOT::Fit::DataOptions();
187  opt.fIntegral = true;
188  opt.fUseEmpty = true;
189  ROOT::Fit::BinData dl2(opt);
190  ROOT::Fit::FillData(dl2,h1,func);
191  f.SetParameters(p);
192  ret = fitter.LikelihoodFit(dl2, f, true);
193  if (ret)
194  fitter.Result().Print(std::cout);
195  else {
196  std::cout << "Integral Likelihood Fit Failed " << std::endl;
197  return -1;
198  }
199  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D histogram integral likelihood fit",0.3);
200 
201  // compare with TH1::Fit
202  std::cout << "\n******************************\n\t TH1::Fit Result \n" << std::endl;
203  func->SetParameters(p);
204  h1->Fit(func,"IL");
205  //std::cout << "Equivalent Chi2 from TF1::Fit " << func->GetChisquare() << std::endl;
206  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TH1::Fit likelihood integral ",0.001);
207 
208 
209 
210  // redo chi2fit
211  std::cout << "\n\nRedo Chi2 Hist Fit" << std::endl;
212  f.SetParameters(p);
213  ret = fitter.Fit(d, f);
214  if (ret)
215  fitter.Result().Print(std::cout);
216  else {
217  std::cout << "Chi2 Fit Failed " << std::endl;
218  return -1;
219  }
220  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D histogram chi2 fit (2)",0.001);
221 
222 
223 
224  // test grapherrors fit
225  std::cout << "\n\nTest Same Fit from a TGraphErrors - no coord errors" << std::endl;
226  TGraphErrors gr(h1);
228  dg.Opt().fCoordErrors = false; // do not use coordinate errors (default is using )
229  ROOT::Fit::FillData(dg,&gr);
230 
231  f.SetParameters(p);
232  ret = fitter.Fit(dg, f);
233  if (ret)
234  fitter.Result().Print(std::cout);
235  else {
236  std::cout << "Chi2 Graph Errors Fit Failed " << std::endl;
237  return -1;
238  }
239  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"TGraphErrors chi2 fit",0.001);
240 
241 
242  // fit using error on X
243  std::cout << "\n\nTest Same Fit from a TGraphErrors - use coord errors" << std::endl;
244  ROOT::Fit::BinData dger;
245  // not needed since they are used by default
246  //dger.Opt().fCoordErrors = true; // use coordinate errors
247  dger.Opt().fUseEmpty = true; // this will set error 1 for the empty bins
248  ROOT::Fit::FillData(dger,&gr);
249 
250  f.SetParameters(p);
251  ret = fitter.Fit(dger, f);
252  if (ret)
253  fitter.Result().Print(std::cout);
254  else {
255  std::cout << "Chi2 Graph Errors Fit Failed " << std::endl;
256  return -1;
257  }
258  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"TGraphErrors effective chi2 fit ",0.7);
259 
260  // compare with TGraphErrors::Fit
261  std::cout << "\n******************************\n\t TGraphErrors::Fit Result \n" << std::endl;
262  func->SetParameters(p);
263  // set error = 1 for empty bins
264  for (int ip = 0; ip < gr.GetN(); ++ip)
265  if (gr.GetErrorY(ip) <= 0) gr.SetPointError(ip, gr.GetErrorX(ip), 1.);
266 
267  gr.Fit(func);
268  std::cout << "Ndf of TGraphErrors::Fit = " << func->GetNDF() << std::endl;
269  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TGraphErrors::Fit ",0.001);
270 
271 
272  // test graph fit (errors are 1) do a re-normalization
273  std::cout << "\n\nTest Same Fit from a TGraph" << std::endl;
274  fitter.Config().SetNormErrors(true);
275  TGraph gr2(h1);
276  ROOT::Fit::BinData dg2;
277  ROOT::Fit::FillData(dg2,&gr2);
278 
279  f.SetParameters(p);
280  ret = fitter.Fit(dg2, f);
281  if (ret)
282  fitter.Result().Print(std::cout);
283  else {
284  std::cout << "Chi2 Graph Fit Failed " << std::endl;
285  return -1;
286  }
287  //iret |= compareResult(fitter.Result().Chi2(), chi2ref,"TGraph fit (no errors) ",0.3);
288 
289 
290  // compare with TGraph::Fit (no errors)
291  std::cout << "\n******************************\n\t TGraph::Fit Result \n" << std::endl;
292  func->SetParameters(p);
293  gr2.Fit(func);
294  std::cout << "Ndf of TGraph::Fit = " << func->GetNDF() << std::endl;
295 
296  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TGraph::Fit ",0.001);
297 
298 
299  // reddo chi2fit using Fumili2
300  std::cout << "\n\nRedo Chi2 Hist Fit using FUMILI2" << std::endl;
301  f.SetParameters(p);
302  fitter.Config().SetMinimizer("Minuit2","Fumili");
303  ret = fitter.Fit(d, f);
304  if (ret)
305  fitter.Result().Print(std::cout);
306  else {
307  std::cout << "Chi2 Fit Failed " << std::endl;
308  return -1;
309  }
310  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D Histo Fumili2 fit ");
311 
312  // reddo chi2fit using old Fumili
313  std::cout << "\n\nRedo Chi2 Hist Fit using FUMILI" << std::endl;
314  f.SetParameters(p);
315  fitter.Config().SetMinimizer("Fumili");
316  ret = fitter.Fit(d, f);
317  if (ret)
318  fitter.Result().Print(std::cout);
319  else {
320  std::cout << "Chi2 Fit Failed " << std::endl;
321  return -1;
322  }
323  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D Histo Fumili fit ");
324 
325  // test using GSL multi fit (L.M. method)
326  std::cout << "\n\nRedo Chi2 Hist Fit using GSLMultiFit" << std::endl;
327  f.SetParameters(p);
328  fitter.Config().SetMinimizer("GSLMultiFit");
329  ret = fitter.Fit(d, f);
330  if (ret)
331  fitter.Result().Print(std::cout);
332  else {
333  std::cout << "Chi2 Fit Failed " << std::endl;
334  return -1;
335  }
336  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D Histo GSL NLS fit ");
337 
338  // test using GSL multi min method
339  std::cout << "\n\nRedo Chi2 Hist Fit using GSLMultiMin" << std::endl;
340  f.SetParameters(p);
341  fitter.Config().SetMinimizer("GSLMultiMin","BFGS2");
342  ret = fitter.Fit(d, f);
343  if (ret)
344  fitter.Result().Print(std::cout);
345  else {
346  std::cout << "Chi2 Fit Failed " << std::endl;
347  return -1;
348  }
349  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"1D Histo GSL Minimizer fit ");
350 
351  return iret;
352 }
353 
354 
355 class Func1D : public ROOT::Math::IParamFunction {
356 public:
357  void SetParameters(const double *p) { std::copy(p,p+NPar(),fp);}
358  const double * Parameters() const { return fp; }
359  ROOT::Math::IGenFunction * Clone() const {
360  Func1D * f = new Func1D();
361  f->SetParameters(fp);
362  return f;
363  };
364  unsigned int NPar() const { return 3; }
365 private:
366  double DoEvalPar( double x, const double *p) const {
367  return p[0]*x*x + p[1]*x + p[2];
368  }
369  double fp[3];
370 
371 };
372 
373 // gradient 2D function
374 class GradFunc2D : public ROOT::Math::IParamMultiGradFunction {
375 public:
376  void SetParameters(const double *p) { std::copy(p,p+NPar(),fp);}
377  const double * Parameters() const { return fp; }
378  ROOT::Math::IMultiGenFunction * Clone() const {
379  GradFunc2D * f = new GradFunc2D();
380  f->SetParameters(fp);
381  return f;
382  };
383  unsigned int NDim() const { return 2; }
384  unsigned int NPar() const { return 5; }
385 
386  void ParameterGradient( const double * x, const double * , double * grad) const {
387  grad[0] = x[0]*x[0];
388  grad[1] = x[0];
389  grad[2] = x[1]*x[1];
390  grad[3] = x[1];
391  grad[4] = 1;
392  }
393 
394 private:
395 
396  double DoEvalPar( const double *x, const double * p) const {
397  return p[0]*x[0]*x[0] + p[1]*x[0] + p[2]*x[1]*x[1] + p[3]*x[1] + p[4];
398  }
399 // double DoDerivative(const double *x, unsigned int icoord = 0) const {
400 // assert(icoord <= 1);
401 // if (icoord == 0)
402 // return 2. * fp[0] * x[0] + fp[1];
403 // else
404 // return 2. * fp[2] * x[1] + fp[3];
405 // }
406 
407  double DoParameterDerivative(const double * x, const double * p, unsigned int ipar) const {
408  std::vector<double> grad(NPar());
409  ParameterGradient(x, p, &grad[0] );
410  return grad[ipar];
411  }
412 
413  double fp[5];
414 
415 };
416 
418 
419 
420 
421  std::string fname("pol2");
422  TF1 * func = (TF1*)gROOT->GetFunction(fname.c_str());
423  func->SetParameter(0,1.);
424  func->SetParameter(1,2.);
425  func->SetParameter(2,3.0);
426 
427  TRandom3 rndm;
428  int iret = 0;
429  double chi2ref = 0;
430 
431  // fill an histogram
432  TH1D * h2 = new TH1D("h2","h2",30,-5.,5.);
433 // h1->FillRandom(fname.c_str(),100);
434  for (int i = 0; i <1000; ++i)
435  h2->Fill( func->GetRandom() );
436 
437  // fill fit data
439  ROOT::Fit::FillData(d,h2,func);
440 
441 
442  printData(d);
443 
444  // create the function
445  Func1D f;
446 
447  double p[3] = {100,0,3.};
448  f.SetParameters(p);
449 
450 
451  // create the fitter
452  //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl;
453  std::cout << "\n\nTest histo polynomial fit using Fitter" << std::endl;
454 
455  ROOT::Fit::Fitter fitter;
456  bool ret = fitter.Fit(d, f);
457  if (ret)
458  fitter.Result().Print(std::cout,true);
459  else {
460  std::cout << " Fit Failed " << std::endl;
461  return -1;
462  }
463  chi2ref = fitter.Result().Chi2();
464 
465  // compare with TH1::Fit
466  std::cout << "\n******************************\n\t TH1::Fit(pol2) Result \n" << std::endl;
467  func->SetParameters(p);
468  h2->Fit(func,"F");
469  iret |= compareResult(func->GetChisquare(),chi2ref,"TH1::Fit ",0.001);
470 
471 
472  std::cout << "\n\nTest histo polynomial linear fit " << std::endl;
473 
474  ROOT::Math::WrappedTF1 pf(*func);
475  //ROOT::Math::Polynomial pf(2);
476  pf.SetParameters(p);
477 
478  fitter.Config().SetMinimizer("Linear");
479  ret = fitter.Fit(d, pf);
480  if (ret) {
481  fitter.Result().Print(std::cout);
482  fitter.Result().PrintCovMatrix(std::cout);
483  }
484  else {
485  std::cout << " Fit Failed " << std::endl;
486  return -1;
487  }
488  iret |= compareResult(fitter.Result().Chi2(),chi2ref,"1D histo linear Fit ");
489 
490  // compare with TH1::Fit
491  std::cout << "\n******************************\n\t TH1::Fit(pol2) Result with TLinearFitter \n" << std::endl;
492  func->SetParameters(p);
493  h2->Fit(func);
494  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TH1::Fit linear",0.001);
495 
496 
497 
498  return iret;
499 }
500 
502 
503  // fit using a 2d parabola (test also gradient)
504 
505 
506  std::string fname("pol2");
507  TF2 * func = new TF2("f2d",ROOT::Math::ParamFunctor(GradFunc2D() ), 0.,10.,0,10,5);
508  double p0[5] = { 1.,2.,0.5,1.,3. };
509  func->SetParameters(p0);
510  assert(func->GetNpar() == 5);
511 
512  TRandom3 rndm;
513  double chi2ref = 0;
514  int iret = 0;
515 
516  // fill an histogram
517  TH2D * h2 = new TH2D("h2d","h2d",30,0,10.,30,0.,10.);
518 // h1->FillRandom(fname.c_str(),100);
519  for (int i = 0; i <10000; ++i) {
520  double x,y = 0;
521  func->GetRandom2(x,y);
522  h2->Fill(x,y);
523  }
524  // fill fit data
526  ROOT::Fit::FillData(d,h2,func);
527 
528 
529  //printData(d);
530 
531  // create the function
532  GradFunc2D f;
533 
534  double p[5] = { 2.,1.,1,2.,100. };
535  f.SetParameters(p);
536 
537 
538  // create the fitter
539  ROOT::Fit::Fitter fitter;
540  //fitter.Config().MinimizerOptions().SetPrintLevel(3);
541  fitter.Config().SetMinimizer("Minuit2");
542 
543  std::cout <<"\ntest 2D histo fit using gradient" << std::endl;
544  bool ret = fitter.Fit(d, f);
545  if (ret)
546  fitter.Result().Print(std::cout);
547  else {
548  std::cout << "Gradient Fit Failed " << std::endl;
549  return -1;
550  }
551  chi2ref = fitter.Result().Chi2();
552 
553  // test without gradient
554  std::cout <<"\ntest result without using gradient" << std::endl;
556  ret = fitter.Fit(d, f2);
557  if (ret)
558  fitter.Result().Print(std::cout);
559  else {
560  std::cout << " Chi2 Fit Failed " << std::endl;
561  return -1;
562  }
563  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"2D histogram chi2 fit");
564 
565  // test Poisson bin likelihood fit (no gradient)
566  std::cout <<"\ntest result without gradient and binned likelihood" << std::endl;
567  f.SetParameters(p);
568  fitter.SetFunction(static_cast<const ROOT::Math::IParamMultiFunction &>(f) );
569  fitter.Config().ParSettings(0).SetLimits(0,100);
570  fitter.Config().ParSettings(1).SetLimits(0,100);
571  fitter.Config().ParSettings(2).SetLimits(0,100);
572  fitter.Config().ParSettings(3).SetLimits(0,100);
573  fitter.Config().ParSettings(4).SetLowerLimit(0);
574  //fitter.Config().MinimizerOptions().SetPrintLevel(3);
575  ret = fitter.LikelihoodFit(d);
576  if (ret)
577  fitter.Result().Print(std::cout);
578  else {
579  std::cout << "Poisson 2D Bin Likelihood Fit Failed " << std::endl;
580  return -1;
581  }
582 
583  // test binned likelihood gradient
584  std::cout <<"\ntest result using gradient and binned likelihood" << std::endl;
585  f.SetParameters(p);
586  fitter.SetFunction(f);
587  //fitter.Config().MinimizerOptions().SetPrintLevel(3);
588  fitter.Config().ParSettings(0).SetLimits(0,100);
589  fitter.Config().ParSettings(1).SetLimits(0,100);
590  fitter.Config().ParSettings(2).SetLimits(0,100);
591  fitter.Config().ParSettings(3).SetLimits(0,100);
592  fitter.Config().ParSettings(4).SetLowerLimit(0);
593  ret = fitter.LikelihoodFit(d);
594  if (ret) {
595  // redo fit releasing the parameters
596  f.SetParameters(&(fitter.Result().Parameters().front()) );
597  ret = fitter.LikelihoodFit(d,f, true);
598  }
599  if (ret)
600  fitter.Result().Print(std::cout);
601  else {
602  std::cout << "Gradient Bin Likelihood Fit Failed " << std::endl;
603  return -1;
604  }
605 
606  // test with linear fitter
607  std::cout <<"\ntest result using linear fitter" << std::endl;
608  fitter.Config().SetMinimizer("Linear");
609  f.SetParameters(p);
610  ret = fitter.Fit(d, f);
611  if (ret)
612  fitter.Result().Print(std::cout);
613  else {
614  std::cout << "Linear 2D Fit Failed " << std::endl;
615  return -1;
616  }
617  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"2D histogram linear fit");
618 
619  // test fitting using TGraph2D ( chi2 will be larger since errors are 1)
620  // should test with a TGraph2DErrors
621  TGraph2D g2(h2);
622 
623  std::cout <<"\ntest using TGraph2D" << std::endl;
625  ROOT::Fit::FillData(d2,&g2,func);
626  //g2.Dump();
627  std::cout << "data size from graph " << d2.Size() << std::endl;
628 
629  f2.SetParameters(p);
630  fitter.Config().SetMinimizer("Minuit2");
631  ret = fitter.Fit(d2, f2);
632  if (ret)
633  fitter.Result().Print(std::cout);
634  else {
635  std::cout << " TGraph2D Fit Failed " << std::endl;
636  return -1;
637  }
638  double chi2ref2 = fitter.Result().Chi2();
639 
640  // compare with TGraph2D::Fit
641  std::cout << "\n******************************\n\t TGraph::Fit Result \n" << std::endl;
642  func->SetParameters(p);
643  g2.Fit(func);
644 
645  std::cout <<"\ntest using TGraph2D and gradient function" << std::endl;
646  f.SetParameters(p);
647  ret = fitter.Fit(d2, f);
648  if (ret)
649  fitter.Result().Print(std::cout);
650  else {
651  std::cout << " TGraph2D Grad Fit Failed " << std::endl;
652  return -1;
653  }
654  iret |= compareResult(fitter.Result().Chi2(), chi2ref2,"TGraph2D chi2 fit");
655 
656  std::cout <<"\ntest using TGraph2DErrors - error only in Z" << std::endl;
657  TGraph2DErrors g2err(g2.GetN() );
658  // need to set error by hand since constructor from TH2 does not exist
659  for (int i = 0; i < g2.GetN(); ++i) {
660  double x = g2.GetX()[i];
661  double y = g2.GetY()[i];
662  g2err.SetPoint(i,x,y,g2.GetZ()[i]);
663  g2err.SetPointError(i,0,0,h2->GetBinError(h2->FindBin(x,y) ) );
664  }
665  func->SetParameters(p);
666  // g2err.Fit(func);
667  f.SetParameters(p);
669  ROOT::Fit::FillData(d3,&g2err,func);
670  ret = fitter.Fit(d3, f);
671  if (ret)
672  fitter.Result().Print(std::cout);
673  else {
674  std::cout << " TGraph2DErrors Fit Failed " << std::endl;
675  return -1;
676  }
677 
678  iret |= compareResult(fitter.Result().Chi2(), chi2ref,"TGraph2DErrors chi2 fit");
679 
680 
681 
682  std::cout <<"\ntest using TGraph2DErrors - with error in X,Y,Z" << std::endl;
683  for (int i = 0; i < g2err.GetN(); ++i) {
684  double x = g2.GetX()[i];
685  double y = g2.GetY()[i];
686  g2err.SetPointError(i,0.5* h2->GetXaxis()->GetBinWidth(1),0.5*h2->GetXaxis()->GetBinWidth(1),h2->GetBinError(h2->FindBin(x,y) ) );
687  }
688  std::cout << "\n******************************\n\t TGraph2DErrors::Fit Result \n" << std::endl;
689  func->SetParameters(p);
690  g2err.Fit(func);
691 
692 
693  return iret;
694 }
695 
696 
698 
699  int iret = 0;
700 
701  std::string fname("gausn");
702  TF1 * func = (TF1*)gROOT->GetFunction(fname.c_str());
703 
704  TRandom3 rndm;
705 
706  int n = 100;
708 
709  for (int i = 0; i <n; ++i)
710  d.Add( rndm.Gaus(0,1) );
711 
712 
713  // printData(d);
714 
715  // create the function
716  ROOT::Math::WrappedMultiTF1 wf(*func);
717  // need to do that to avoid gradient calculation
719  double p[3] = {1,2,10.};
720  f.SetParameters(p);
721 
722  // create the fitter
723  //std::cout << "Fit parameters " << f.Parameters()[2] << std::endl;
724 
725  ROOT::Fit::Fitter fitter;
726  fitter.SetFunction(f);
727  std::cout << "fix parameter 0 " << " to value " << f.Parameters()[0] << std::endl;
728  fitter.Config().ParSettings(0).Fix();
729  // set range in sigma sigma > 0
730  std::cout << "set lower range to 0 for sigma " << std::endl;
731  fitter.Config().ParSettings(2).SetLowerLimit(0);
732 
733 #ifdef DEBUG
734  fitter.Config().MinimizerOptions().SetPrintLevel(3);
735 #endif
736 
737 // double x[1]; x[0] = 0.;
738 // std::cout << "fval " << f(x) << std::endl;
739 // x[0] = 1.; std::cout << "fval " << f(x) << std::endl;
740 
741  //fitter.Config().SetMinimizer("Minuit2");
742  // fails with minuit (t.b. investigate)
743  fitter.Config().SetMinimizer("Minuit2");
744 
745 
746  bool ret = fitter.Fit(d);
747  if (ret)
748  fitter.Result().Print(std::cout);
749  else {
750  std::cout << "Unbinned Likelihood Fit Failed " << std::endl;
751  iret |= 1;
752  }
753  double lref = fitter.Result().MinFcnValue();
754 
755  std::cout << "\n\nRedo Fit using FUMILI2" << std::endl;
756  f.SetParameters(p);
757  fitter.Config().SetMinimizer("Fumili2");
758  // need to set function first (need to change this)
759  fitter.SetFunction(f);
760  fitter.Config().ParSettings(0).Fix(); //need to re-do it
761  // set range in sigma sigma > 0
762  fitter.Config().ParSettings(2).SetLowerLimit(0);
763 
764  ret = fitter.Fit(d);
765  if (ret)
766  fitter.Result().Print(std::cout);
767  else {
768  std::cout << "Unbinned Likelihood Fit using FUMILI2 Failed " << std::endl;
769  iret |= 1;
770  }
771 
772  iret |= compareResult(fitter.Result().MinFcnValue(), lref,"1D unbin FUMILI2 fit");
773 
774  std::cout << "\n\nRedo Fit using FUMILI" << std::endl;
775  f.SetParameters(p);
776  fitter.Config().SetMinimizer("Fumili");
777  // fitter.Config().MinimizerOptions().SetPrintLevel(3);
778  // need to set function first (need to change this)
779  fitter.SetFunction(f);
780  fitter.Config().ParSettings(0).Fix(); //need to re-do it
781  // set range in sigma sigma > 0
782  fitter.Config().ParSettings(2).SetLowerLimit(0);
783 
784  ret = fitter.Fit(d);
785  if (ret)
786  fitter.Result().Print(std::cout);
787  else {
788  std::cout << "Unbinned Likelihood Fit using FUMILI Failed " << std::endl;
789  iret |= 1;
790  }
791 
792  iret |= compareResult(fitter.Result().MinFcnValue(), lref,"1D unbin FUMILI fit");
793 
794 
795  return iret;
796 }
797 
799 
800  int iret = 0;
801 
802  // simple test of fitting a Tgraph
803 
804  double x[5] = {1,2,3,4,5};
805  double y[5] = {2.1, 3.5, 6.5, 8.8, 9.5};
806  double ex[5] = {.3,.3,.3,.3,.3};
807  double ey[5] = {.5,.5,.5,.5,.5};
808  double eyl[5] = {.2,.2,.2,.2,.2};
809  double eyh[5] = {.8,.8,.8,.8,.8};
810 
811  std::cout << "\n********************************************************\n";
812  std::cout << "Test simple fit of Tgraph of 5 points" << std::endl;
813  std::cout << "\n********************************************************\n";
814 
815 
816  double p[2] = {1,1};
817  TF1 * func = new TF1("f","pol1",0,10);
818  func->SetParameters(p);
819 
820  ROOT::Math::WrappedMultiTF1 wf(*func);
821  // need to do that to avoid gradient calculation
823  f.SetParameters(p);
824 
825  ROOT::Fit::Fitter fitter;
826  fitter.SetFunction(f);
827 
828 
829  std::cout <<"\ntest TGraph (no errors) " << std::endl;
830  TGraph gr(5, x,y);
831 
832  ROOT::Fit::BinData dgr;
833  ROOT::Fit::FillData(dgr,&gr);
834 
835  //printData(dgr);
836 
837  f.SetParameters(p);
838  bool ret = fitter.Fit(dgr, f);
839  if (ret)
840  fitter.Result().Print(std::cout);
841  else {
842  std::cout << "Chi2 Graph Fit Failed " << std::endl;
843  return -1;
844  }
845  double chi2ref = fitter.Result().Chi2();
846 
847 
848  // compare with TGraph::Fit
849  std::cout << "\n******************************\n\t TGraph::Fit Result \n" << std::endl;
850  func->SetParameters(p);
851  gr.Fit(func,"F"); // use Minuit
852 
853  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TGraph::Fit ",0.001);
854 
855 
856  std::cout <<"\ntest TGraphErrors " << std::endl;
857  TGraphErrors grer(5, x,y,ex,ey);
858 
859  ROOT::Fit::BinData dgrer;
860  dgrer.Opt().fCoordErrors = true;
861  ROOT::Fit::FillData(dgrer,&grer);
862 
863  //printData(dgrer);
864 
865  f.SetParameters(p);
866  ret = fitter.Fit(dgrer, f);
867  if (ret)
868  fitter.Result().Print(std::cout);
869  else {
870  std::cout << "Chi2 Graph Fit Failed " << std::endl;
871  return -1;
872  }
873 
874  iret |= compareResult(fitter.Result().Chi2(),chi2ref,"TGraphErrors fit with coord errors",0.8);
875 
876 
877  // compare with TGraph::Fit
878  std::cout << "\n******************************\n\t TGraphErrors::Fit Result \n" << std::endl;
879  func->SetParameters(p);
880  grer.Fit(func,"F"); // use Minuit
881  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TGraphErrors::Fit ",0.001);
882 
883  chi2ref = fitter.Result().Chi2();
884 
885  std::cout <<"\ntest TGraphAsymmErrors " << std::endl;
886  TGraphAsymmErrors graer(5, x,y,ex,ex,eyl, eyh);
887 
888  ROOT::Fit::BinData dgraer;
889  // option error on coordinate and asymmetric on values
890  dgraer.Opt().fCoordErrors = true;
891  dgraer.Opt().fAsymErrors = true;
892  ROOT::Fit::FillData(dgraer,&graer);
893  //printData(dgraer);
894 
895  f.SetParameters(p);
896  ret = fitter.Fit(dgraer, f);
897  if (ret)
898  fitter.Result().Print(std::cout);
899  else {
900  std::cout << "Chi2 Graph Fit Failed " << std::endl;
901  return -1;
902  }
903  //iret |= compareResult(fitter.Result().Chi2(),chi2ref,"TGraphAsymmErrors fit ",0.5);
904 
905  // compare with TGraph::Fit
906  std::cout << "\n******************************\n\t TGraphAsymmErrors::Fit Result \n" << std::endl;
907  func->SetParameters(p);
908  graer.Fit(func,"F"); // use Minuit
909  iret |= compareResult(func->GetChisquare(),fitter.Result().Chi2(),"TGraphAsymmErrors::Fit ",0.001);
910 
911 
912 
913  return iret;
914 }
915 
916 
917 template<typename Test>
918 int testFit(Test t, std::string name) {
919  std::cout << name << "\n\t\t";
920  int iret = t();
921  std::cout << "\n" << name << ":\t\t";
922  if (iret == 0)
923  std::cout << "OK" << std::endl;
924  else
925  std::cout << "Failed" << std::endl;
926  return iret;
927 }
928 
929 int main() {
930 
931  int iret = 0;
932  iret |= testFit( testHisto1DFit, "Histogram1D Fit");
933  iret |= testFit( testHisto1DPolFit, "Histogram1D Polynomial Fit");
934  iret |= testFit( testHisto2DFit, "Histogram2D Gradient Fit");
935  iret |= testFit( testUnBin1DFit, "Unbin 1D Fit");
936  iret |= testFit( testGraphFit, "Graph 1D Fit");
937 
938  std::cout << "\n******************************\n";
939  if (iret) std::cerr << "\n\t testFit FAILED !!!!!!!!!!!!!!!! \n";
940  else std::cout << "\n\t testFit all OK \n";
941  return iret;
942 }
943 
unsigned int Size() const
return number of fit points
Definition: BinData.h:447
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Axis_t xmin=0, Axis_t xmax=0)
Fit this graph with function with name fname.
Definition: TGraph.cxx:1047
int testHisto2DFit()
Definition: testFit.cxx:501
virtual void Print(Option_t *option="") const
Print some global quantities for this histogram.
Definition: TH1.cxx:6324
virtual const double * Parameters() const =0
Access the parameter values.
virtual Int_t FindBin(Double_t x, Double_t y=0, Double_t z=0)
Return Global bin number corresponding to x,y,z.
Definition: TH1.cxx:3440
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3125
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
virtual Double_t GetRandom()
Return a random number following this function shape.
Definition: TF1.cxx:1798
void SetPrintLevel(int level)
set print level
Random number generator class based on M.
Definition: TRandom3.h:29
int testHisto1DPolFit()
Definition: testFit.cxx:417
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:133
double Chisquare(const TH1 &h1, TF1 &f1, bool useRange, bool usePL=false)
compute the chi2 value for an histogram given a function (see TH1::Chisquare for the documentation) ...
Definition: HFitImpl.cxx:997
ErrorType GetErrorType() const
Definition: BinData.h:75
int compareResult(double v1, double v2, std::string s="", double tol=0.01)
Definition: testFit.cxx:57
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
Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface of multi-dimensions...
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
#define gROOT
Definition: TROOT.h:364
Class to Wrap a ROOT Function class (like TF1) in a IParamFunction interface of one dimensions to be ...
Definition: WrappedTF1.h:41
int testUnBin1DFit()
Definition: testFit.cxx:697
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:47
double MinFcnValue() const
Return value of the objective function (chi2 or likelihood) used in the fit.
Definition: FitResult.h:125
ROOT::Math::MinimizerOptions & MinimizerOptions()
access to the minimizer control parameter (non const method)
Definition: FitConfig.h:138
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition: FitConfig.h:80
const FitResult & Result() const
get fit result
Definition: Fitter.h:354
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:1602
TGraph with asymmetric error bars.
void SetParameters(const double *p)
set parameter values
Definition: WrappedTF1.h:90
void SetParameters(TFitEditor::FuncParams_t &pars, TF1 *func)
Restore the parameters from pars into the function.
Definition: TFitEditor.cxx:287
Double_t x[n]
Definition: legend1.C:17
void Fix()
fix the parameter
int main()
Definition: testFit.cxx:929
int testGraphFit()
Definition: testFit.cxx:798
TH1F * h1
Definition: legend1.C:5
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:382
double Error(unsigned int ipoint) const
return error on the value for the given fit point Safe (but slower) method returning correctly the er...
Definition: BinData.h:249
int testHisto1DFit()
Definition: testFit.cxx:71
Double_t GetErrorX(Int_t bin) const
This function is called by GraphFitChisquare.
double chi2FromFit(const TF1 *func)
Definition: testFit.cxx:65
void SetLowerLimit(double low)
set a single lower limit
virtual unsigned int NPar() const =0
Return the number of Parameters.
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
const double tol
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
IParamFunction interface (abstract class) describing multi-dimensional parameteric functions It is a ...
RooCmdArg Parameters(const RooArgSet &params)
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:94
const double * Coords(unsigned int ipoint) const
return a pointer to the coordinates data for the given fit point
Definition: BinData.h:226
static TVirtualFitter * GetFitter()
static: return the current Fitter
const DataOptions & Opt() const
access to options
Definition: DataVector.h:97
Int_t GetN() const
Definition: TGraph.h:133
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:61
unsigned int PointSize() const
return the size of a fit point (is the coordinate dimension + 1 for the value and eventually the numb...
Definition: BinData.h:150
A 2-Dim function with parameters.
Definition: TF2.h:33
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
void SetFunction(const IModelFunction &func, bool useGradient=false)
Set the fitted function (model function) from a parametric function interface.
Definition: Fitter.cxx:104
TGraphErrors * gr
Definition: legend1.C:25
bool LikelihoodFit(const BinData &data, bool extended=true)
Binned Likelihood fit.
Definition: Fitter.h:181
Double_t GetChisquare() const
Definition: TF1.h:336
virtual void SetParameters(const double *p)=0
Set the parameter values.
const double * Coords(unsigned int ipoint) const
return pointer to coordinate data
Definition: UnBinData.h:282
void PrintCovMatrix(std::ostream &os) const
print error matrix and correlations
Definition: FitResult.cxx:501
Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
double f(double x)
Specialized IParamFunction interface (abstract class) for one-dimensional parametric functions It is ...
virtual void GetRandom2(Double_t &xrandom, Double_t &yrandom)
Return 2 random numbers following this function shape.
Definition: TF2.cxx:504
void Add(double x)
add one dim coordinate data (unweighted)
Definition: UnBinData.h:199
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
const std::vector< double > & Parameters() const
parameter values (return std::vector)
Definition: FitResult.h:179
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Double_t ey[n]
Definition: legend1.C:17
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:307
int testFit(Test t, std::string name)
Definition: testFit.cxx:918
double Value(unsigned int ipoint) const
return the value for the given fit point
Definition: BinData.h:236
unsigned int NDim() const
return coordinate data dimension
Definition: BinData.h:452
virtual Int_t GetNpar() const
Definition: TF1.h:349
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width.
Definition: TAxis.cxx:526
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
double f2(const double *x)
1-Dim function class
Definition: TF1.h:149
const double * GetPointError(unsigned int ipoint, double &errvalue) const
Retrieve the errors on the point (coordinate and value) for the given fit point It must be called onl...
Definition: BinData.h:349
void SetNormErrors(bool on=true)
set the option to normalize the error on the result according to chi2/ndf
Definition: FitConfig.h:192
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
Param Functor class for Multidimensional functions.
Definition: ParamFunctor.h:209
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:28
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:301
virtual IBaseFunctionOneDim * Clone() const =0
Clone a function.
virtual Double_t * GetParameters() const
Definition: TF1.h:365
void SetLimits(double low, double up)
set a double side limit, if low == up the parameter is fixed if low > up the limits are removed ...
void SetParameters(const double *p)
Set the parameter values.
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:431
const int NPar
virtual void SetPointError(Double_t ex, Double_t ey)
Set ex and ey values for point pointed by the mouse.
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:50
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:63
virtual double DoEvalPar(double x, const double *p) const =0
Implementation of the evaluation function using the x value and the parameters.
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
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:3563
Double_t ex[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
double Chi2() const
Chi2 fit value in case of likelihood must be computed ?
Definition: FitResult.h:165
unsigned int Size() const
return number of contained points
Definition: UnBinData.h:316
void printData(const ROOT::Fit::BinData &data)
Definition: testFit.cxx:29
char name[80]
Definition: TGX11.cxx:109
TAxis * GetXaxis()
Definition: TH1.h:324
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8173
Graph 2D class with errors.
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:296