Logo ROOT  

Reference Guide
Loading...
Searching...
No Matches
minuit2FitBench.C File Reference
Tutorials » Math tutorials » Fit Tutorials

Detailed Description

View in nbviewer Open in SWAN
Demonstrate performance and usage of Minuit2 and Fumili2 for monodimensional fits.

*********************************************************************************
Minuit
*********************************************************************************
pass : 0
................... FCN=205.276 FROM MINOS STATUS=SUCCESSFUL 44 CALLS 429 TOTAL
EDM=3.83274e-10 STRATEGY= 1 ERROR MATRIX ACCURATE
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 5.13639e+01 2.01330e+00 -2.79355e-04 -2.05444e-06
2 p1 5.57813e+01 4.80586e+00 3.09152e-03 -9.98323e-07
3 p2 7.42112e+01 1.87043e+00 -1.20323e-03 -1.92352e-07
4 p3 4.27344e+02 2.93232e+00 -1.66250e-02 -7.79026e-07
5 p4 3.58604e-02 3.47006e-04 1.73862e-07 9.80669e-02
6 p5 1.00001e+00 1.64203e-04 1.64203e-04 3.19367e-02
Minuit, npass=20 : RT= 0.238 s, Cpu= 0.230 s
*********************************************************************************
Fumili
*********************************************************************************
pass : 0
...................****************************************
Minimizer is Fumili
Chi2 = 206.284
NDf = 194
NCalls = 4
p0 = 51.4325 +/- 2.01397
p1 = 55.5412 +/- 4.81253
p2 = 74.2976 +/- 1.87298
p3 = 427.425 +/- 2.93868
p4 = 0.0358559 +/- 0.000357243
p5 = 1.00001 +/- 0.00016009
Fumili, npass=20 : RT= 0.104 s, Cpu= 0.060 s
*********************************************************************************
Minuit2
*********************************************************************************
pass : 0
...................****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 205.34
NDf = 194
Edm = 1.91396e-10
NCalls = 85
p0 = 51.3576 +/- 2.0133 -2.01329 +2.0133 (Minos)
p1 = 55.8172 +/- 4.80582 -4.80605 +4.80561 (Minos)
p2 = 74.1981 +/- 1.8704 -1.87032 +1.8705 (Minos)
p3 = 427.317 +/- 2.93218 -2.93197 +2.93246 (Minos)
p4 = 0.0358574 +/- 0.000346961 -0.000345228 +0.000348733 (Minos)
p5 = 1.00001 +/- 0.000164204 -0.000164437 +0.000163973 (Minos)
Minuit2, npass=20 : RT= 0.080 s, Cpu= 0.090 s
*********************************************************************************
Fumili2
*********************************************************************************
pass : 0
...................****************************************
Minimizer is Minuit2 / Fumili
Chi2 = 207.495
NDf = 194
Edm = 1.09727e-10
NCalls = 2
p0 = 51.4278 +/- 2.01392
p1 = 55.5389 +/- 4.81246
p2 = 74.3004 +/- 1.87297
p3 = 427.406 +/- 2.93858
p4 = 0.0358584 +/- 0.000357283
p5 = 1.00001 +/- 0.000160099
Fumili2, npass=20 : RT= 0.036 s, Cpu= 0.030 s
(int) 0
#include "TH1.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "TRandom3.h"
#include "Math/MinimizerOptions.h"
#include "TPaveLabel.h"
#include "TStyle.h"
#include "TMath.h"
#include "TROOT.h"
#include "TFrame.h"
/*#include "Fit/FitConfig.h"*/
TF1 *fitFcn;
TH1 *histo;
// Quadratic background function
double background(double *x, double *par) {
return par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
}
// Lorenzian Peak function
double lorentzianPeak(double *x, double *par) {
return (0.5*par[0]*par[1]/TMath::Pi()) /
TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]);
}
// Sum of background and peak function
double fitFunction(double *x, double *par) {
return background(x,par) + lorentzianPeak(x,&par[3]);
}
bool DoFit(const char* fitter, TVirtualPad *pad, int npass) {
printf("\n*********************************************************************************\n");
printf("\t %s \n",fitter);
printf("*********************************************************************************\n");
gRandom = new TRandom3();
TStopwatch timer;
// timer.Start();
ROOT::Math::MinimizerOptions::SetDefaultMinimizer(fitter);
pad->SetGrid();
pad->SetLogy();
fitFcn->SetParameters(1,1,1,6,.03,1);
fitFcn->Update();
std::string title = std::string(fitter) + " fit bench";
histo = new TH1D(fitter,title.c_str(),200,0,3);
TString fitterType(fitter);
timer.Start();
bool ok = true;
// fill histogram many times
// every time increase its statistics and re-use previous fitted
// parameter values as starting point
for (int pass=0;pass<npass;pass++) {
if (pass%100 == 0) printf("pass : %d\n",pass);
else printf(".");
if (pass == 0)fitFcn->SetParameters(1,1,1,6,.03,1);
for (int i=0;i<5000;i++) {
histo->Fill(fitFcn->GetRandom());
}
int iret = histo->Fit(fitFcn,"Q0");
ok &= (iret == 0);
if (iret!=0) Error("DoFit","Fit pass %d failed !",pass);
}
// do last fit computing Minos Errors (except for Fumili)
if (!fitterType.Contains("Fumili")) // Fumili does not implement Error options (MINOS)
histo->Fit(fitFcn,"E");
else
histo->Fit(fitFcn,"");
timer.Stop();
(histo->GetFunction("fitFcn"))->SetLineColor(kRed+3);
gPad->SetFillColor(kYellow-10);
double cputime = timer.CpuTime();
printf("%s, npass=%d : RT=%7.3f s, Cpu=%7.3f s\n",fitter,npass,timer.RealTime(),cputime);
TPaveLabel *p = new TPaveLabel(0.45,0.7,0.88,0.8,Form("%s CPU= %g s",fitter,cputime),"brNDC");
p->Draw();
p->SetTextColor(kRed+3);
p->SetFillColor(kYellow-8);
pad->Update();
return ok;
}
int minuit2FitBench(int npass=20) {
TDirectory::TContext ctx{nullptr}; // Don't register histograms to the current directory
TCanvas *c1 = new TCanvas("FitBench","Fitting Demo",10,10,900,900);
c1->Divide(2,2);
c1->SetFillColor(kYellow-9);
// create a TF1 with the range from 0 to 3 and 6 parameters
fitFcn = new TF1("fitFcn",fitFunction,0,3,6);
fitFcn->SetNpx(200);
gStyle->SetOptFit();
gStyle->SetStatY(0.6);
bool ok = true;
//with Minuit
c1->cd(1);
ok &= DoFit("Minuit",gPad,npass);
//with Fumili
c1->cd(2);
ok &= DoFit("Fumili",gPad,npass);
//with Minuit2
c1->cd(3);
ok &= DoFit("Minuit2",gPad,npass);
//with Fumili2
c1->cd(4);
ok &= DoFit("Fumili2",gPad,npass);
c1->SaveAs("FitBench.root");
return (ok) ? 0 : 1;
}
kRed
@ kRed
Definition Rtypes.h:67
kYellow
@ kYellow
Definition Rtypes.h:67
Error
Error("WriteTObject","The current directory (%s) is not associated with a file. The object (%s) has not been written.", GetName(), objname)
gRandom
externTRandom * gRandom
Definition TRandom.h:62
Form
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2496
gStyle
externTStyle * gStyle
Definition TStyle.h:442
gPad
#define gPad
Definition TVirtualPad.h:322
ROOT::Math::MinimizerOptions::SetDefaultMinimizer
static void SetDefaultMinimizer(const char *type, const char *algo=nullptr)
Set the default Minimizer type and corresponding algorithms.
Definition MinimizerOptions.cxx:43
TAttFill::SetFillColor
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:40
TAttText::SetTextColor
virtual void SetTextColor(Color_t tcolor=1)
Set the text color.
Definition TAttText.h:50
TCanvas
The Canvas class.
Definition TCanvas.h:23
TDirectory::TContext
TDirectory::TContext keeps track and restore the current directory.
Definition TDirectory.h:89
TF1
Definition TF1.h:182
TF1::SetNpx
virtual void SetNpx(Int_t npx=100)
TF1::GetRandom
virtual Double_t GetRandom(TRandom *rng=nullptr, Option_t *opt=nullptr)
TF1::Update
virtual void Update()
TF1::SetParameters
virtual void SetParameters(const Double_t *params)
Definition TF1.h:618
TH1D
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:926
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
TH1::Fit
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:3954
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3393
TH1::GetFunction
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition TH1.cxx:9185
TPaveLabel
A Pave (see TPave) with a text centered in the Pave.
Definition TPaveLabel.h:20
TPaveLabel::Draw
void Draw(Option_t *option="") override
Draw this pavelabel with its current attributes.
Definition TPaveLabel.cxx:87
TRandom3
Random number generator class based on M.
Definition TRandom3.h:27
TStopwatch
Stopwatch class.
Definition TStopwatch.h:28
TStopwatch::RealTime
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition TStopwatch.cxx:109
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition TStopwatch.cxx:57
TStopwatch::CpuTime
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition TStopwatch.cxx:124
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition TStopwatch.cxx:76
TString
Basic string class.
Definition TString.h:138
TVirtualPad
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition TVirtualPad.h:51
TVirtualPad::SetGrid
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)=0
TVirtualPad::Update
virtual void Update()=0
TVirtualPad::SetLogy
virtual void SetLogy(Int_t value=1)=0
c1
return c1
Definition legend1.C:41
x
Double_t x[n]
Definition legend1.C:17
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:249
TMath::Pi
constexpr Double_t Pi()
Definition TMath.h:40
Author
Lorenzo Moneta

Definition in file minuit2FitBench.C.