Logo ROOT  
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"
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
27TF1 *fitFcn;
28TH1 *histo;
29
30// Quadratic background function
31Double_t background(Double_t *x, Double_t *par) {
32 return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
33}
34
35// Lorenzian Peak function
36Double_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
42Double_t fitFunction(Double_t *x, Double_t *par) {
43 return background(x,par) + lorentzianPeak(x,&par[3]);
44}
45
46bool 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();
52 TStopwatch timer;
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);
97 pad->Update();
98 return ok;
99}
100
101int 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}
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
double Double_t
Definition: RtypesCore.h:55
@ kRed
Definition: Rtypes.h:64
@ kYellow
Definition: Rtypes.h:64
void Error(const char *location, const char *msgfmt,...)
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
char * Form(const char *fmt,...)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:407
#define gPad
Definition: TVirtualPad.h:286
static void SetDefaultMinimizer(const char *type, const char *algo=0)
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:37
virtual void SetTextColor(Color_t tcolor=1)
Set the text color.
Definition: TAttText.h:43
The Canvas class.
Definition: TCanvas.h:31
1-Dim function class
Definition: TF1.h:211
virtual void Update()
Called by functions such as SetRange, SetNpx, SetParameters to force the deletion of the associated h...
Definition: TF1.cxx:3615
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition: TF1.cxx:3432
virtual Double_t GetRandom()
Return a random number following this function shape.
Definition: TF1.cxx:2074
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
The TH1 histogram class.
Definition: TH1.h:56
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1226
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 Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition: TH1.cxx:8495
A Pave (see TPave) with a text centered in the Pave.
Definition: TPaveLabel.h:20
virtual void Draw(Option_t *option="")
Draw this pavelabel with its current attributes.
Definition: TPaveLabel.cxx:77
Random number generator class based on M.
Definition: TRandom3.h:27
Stopwatch class.
Definition: TStopwatch.h:28
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
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
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
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Basic string class.
Definition: TString.h:131
void SetStatY(Float_t y=0)
Definition: TStyle.h:376
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:1402
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition: TVirtualPad.h:50
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)=0
virtual void Update()=0
virtual void SetLogy(Int_t value=1)=0
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
constexpr Double_t Pi()
Definition: TMath.h:38
lv SetLineColor(kBlue)