Logo ROOT   6.10/09
Reference Guide
minuit2FitBench.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_fit
3 /// \notebook
4 /// Demonstrate performance and usage of Minuit2 and Fumili2 for monodimensional fits.
5 ///
6 /// \macro_image
7 /// \macro_output
8 /// \macro_code
9 ///
10 /// \author Lorenzo Moneta
11 
12 #include "TH1.h"
13 #include "TF1.h"
14 #include "TCanvas.h"
15 #include "TStopwatch.h"
16 #include "TSystem.h"
17 #include "TRandom3.h"
18 #include "Math/MinimizerOptions.h"
19 #include "TPaveLabel.h"
20 #include "TStyle.h"
21 #include "TMath.h"
22 #include "TROOT.h"
23 #include "TFrame.h"
24 /*#include "Fit/FitConfig.h"*/
25 
26 
27 TF1 *fitFcn;
28 TH1 *histo;
29 
30 // Quadratic background function
32  return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
33 }
34 
35 // Lorenzian Peak function
36 Double_t lorentzianPeak(Double_t *x, Double_t *par) {
37  return (0.5*par[0]*par[1]/TMath::Pi()) /
38  TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]);
39 }
40 
41 // Sum of background and peak function
43  return background(x,par) + lorentzianPeak(x,&par[3]);
44 }
45 
46 bool DoFit(const char* fitter, TVirtualPad *pad, Int_t npass) {
47  printf("\n*********************************************************************************\n");
48  printf("\t %s \n",fitter);
49  printf("*********************************************************************************\n");
50 
51  gRandom = new TRandom3();
53  // timer.Start();
55  pad->SetGrid();
56  pad->SetLogy();
57  fitFcn->SetParameters(1,1,1,6,.03,1);
58  fitFcn->Update();
59  std::string title = std::string(fitter) + " fit bench";
60  histo = new TH1D(fitter,title.c_str(),200,0,3);
61 
62  TString fitterType(fitter);
63 
64  timer.Start();
65  bool ok = true;
66  // fill histogram many times
67  // every time increase its statistics and re-use previous fitted
68  // parameter values as starting point
69  for (Int_t pass=0;pass<npass;pass++) {
70  if (pass%100 == 0) printf("pass : %d\n",pass);
71  else printf(".");
72  if (pass == 0)fitFcn->SetParameters(1,1,1,6,.03,1);
73  for (Int_t i=0;i<5000;i++) {
74  histo->Fill(fitFcn->GetRandom());
75  }
76  int iret = histo->Fit(fitFcn,"Q0");
77  ok &= (iret == 0);
78  if (iret!=0) Error("DoFit","Fit pass %d failed !",pass);
79  }
80  // do last fit computing Minos Errors (except for Fumili)
81  if (!fitterType.Contains("Fumili")) // Fumili does not implement Error options (MINOS)
82  histo->Fit(fitFcn,"E");
83  else
84  histo->Fit(fitFcn,"");
85  timer.Stop();
86 
87  (histo->GetFunction("fitFcn"))->SetLineColor(kRed+3);
88  gPad->SetFillColor(kYellow-10);
89 
90 
91  Double_t cputime = timer.CpuTime();
92  printf("%s, npass=%d : RT=%7.3f s, Cpu=%7.3f s\n",fitter,npass,timer.RealTime(),cputime);
93  TPaveLabel *p = new TPaveLabel(0.45,0.7,0.88,0.8,Form("%s CPU= %g s",fitter,cputime),"brNDC");
94  p->Draw();
95  p->SetTextColor(kRed+3);
96  p->SetFillColor(kYellow-8);
97  pad->Update();
98  return ok;
99 }
100 
101 int minuit2FitBench(Int_t npass=20) {
103  TCanvas *c1 = new TCanvas("FitBench","Fitting Demo",10,10,900,900);
104  c1->Divide(2,2);
105  c1->SetFillColor(kYellow-9);
106  // create a TF1 with the range from 0 to 3 and 6 parameters
107  fitFcn = new TF1("fitFcn",fitFunction,0,3,6);
108  fitFcn->SetNpx(200);
109  gStyle->SetOptFit();
110  gStyle->SetStatY(0.6);
111 
112  bool ok = true;
113  //with Minuit
114  c1->cd(1);
115  ok &= DoFit("Minuit",gPad,npass);
116 
117  //with Fumili
118  c1->cd(2);
119  ok &= DoFit("Fumili",gPad,npass);
120 
121  //with Minuit2
122  c1->cd(3);
123  ok &= DoFit("Minuit2",gPad,npass);
124 
125  //with Fumili2
126  c1->cd(4);
127  ok &= DoFit("Fumili2",gPad,npass);
128 
129  c1->SaveAs("FitBench.root");
130  return (ok) ? 0 : 1;
131 }
double par[1]
Definition: unuranDistr.cxx:38
const int npass
Definition: testPermute.cxx:21
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3126
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:588
virtual Double_t GetRandom()
Return a random number following this function shape.
Definition: TF1.cxx:1815
Random number generator class based on M.
Definition: TRandom3.h:27
virtual void SetLogy(Int_t value=1)=0
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition: TF1.cxx:3202
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition: TH1.cxx:8163
return c1
Definition: legend1.C:41
Definition: Rtypes.h:56
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
virtual void Update()=0
void DoFit(THnSparse *s, TF1 *f, ROOT::Fit::BinData &bd)
Definition: SparseFit4.cxx:162
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:679
Basic string class.
Definition: TString.h:129
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
int Int_t
Definition: RtypesCore.h:41
Definition: Rtypes.h:56
void SetStatY(Float_t y=0)
Definition: TStyle.h:371
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1218
TStopwatch timer
Definition: pirndm.C:37
Double_t fitFunction(Double_t *x, Double_t *par)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
constexpr Double_t Pi()
Definition: TMath.h:40
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition: TVirtualPad.h:49
void Error(const char *location, const char *msgfmt,...)
lv SetLineColor(kBlue)
A Pave (see TPave) with a text centered in the Pave.
Definition: TPaveLabel.h:20
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
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:1219
char * Form(const char *fmt,...)
virtual void Update()
Called by functions such as SetRange, SetNpx, SetParameters to force the deletion of the associated h...
Definition: TF1.cxx:3430
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
const Bool_t kFALSE
Definition: RtypesCore.h:92
The Canvas class.
Definition: TCanvas.h:31
double Double_t
Definition: RtypesCore.h:55
The TH1 histogram class.
Definition: TH1.h:56
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
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:1135
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
1-Dim function class
Definition: TF1.h:150
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:310
#define gPad
Definition: TVirtualPad.h:284
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:5468
static void SetDefaultMinimizer(const char *type, const char *algo=0)
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:3564
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)=0
Stopwatch class.
Definition: TStopwatch.h:28