ROOT   Reference Guide
Searching...
No Matches
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
31double background(double *x, double *par) {
32 return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
33}
34
35// Lorenzian Peak function
36double lorentzianPeak(double *x, double *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 fitFunction(double *x, double *par) {
43 return background(x,par) + lorentzianPeak(x,&par[3]);
44}
45
47 printf("\n*********************************************************************************\n");
48 printf("\t %s \n",fitter);
49 printf("*********************************************************************************\n");
50
51 gRandom = new TRandom3();
52 TStopwatch timer;
53 // timer.Start();
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 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 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);
89
90
91 double 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);
98 return ok;
99}
100
101int minuit2FitBench(int 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);
116
117 //with Fumili
118 c1->cd(2);
120
121 //with Minuit2
122 c1->cd(3);
124
125 //with Fumili2
126 c1->cd(4);
128
129 c1->SaveAs("FitBench.root");
130 return (ok) ? 0 : 1;
131}
@ kRed
Definition Rtypes.h:66
@ kYellow
Definition Rtypes.h:66
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:185
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t SetLineColor
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2467
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
static void SetDefaultMinimizer(const char *type, const char *algo=nullptr)
Set the default Minimizer type and corresponding algorithms.
The Canvas class.
Definition TCanvas.h:23
1-Dim function class
Definition TF1.h:214
virtual void Update()
Called by functions such as SetRange, SetNpx, SetParameters to force the deletion of the associated h...
Definition TF1.cxx:3613
virtual Double_t GetRandom(TRandom *rng=nullptr, Option_t *opt=nullptr)
Return a random number following this function shape.
Definition TF1.cxx:2192
virtual void SetNpx(Int_t npx=100)
Set the number of points used to draw the function.
Definition TF1.cxx:3433
virtual void SetParameters(const Double_t *params)
Definition TF1.h:650
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:620
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
Sets the flag controlling the automatic add of histograms in memory.
Definition TH1.cxx:1267
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:3901
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3345
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition TH1.cxx:8968
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
A Pave (see TPave) with a text centered in the Pave.
Definition TPaveLabel.h:20
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...
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
void Stop()
Stop the stopwatch.
Basic string class.
Definition TString.h:139
void SetStatY(Float_t y=0)
Definition TStyle.h:395
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:1589
TVirtualPad is an abstract base class for the Pad and Canvas classes.