Logo ROOT  
Reference Guide
fit2dHist.C File Reference

Detailed Description

View in nbviewer Open in SWAN

Example to fit two histograms at the same time via the Fitter class

To execute this tutorial, you can do:

root > .x fit2dHist.C (executing via CINT, slow)

or

root > .x fit2dHist.C+ (executing via ACLIC , fast, with Minuit)
root > .x fit2dHist.C+(2) (executing via ACLIC , fast, with Minuit2)
Double_t x[n]
Definition: legend1.C:17

or using the option to fit independently the 2 histos

root > .x fit2dHist.C+(10) (via ACLIC, fast, independent fits with Minuit)
root > .x fit2dHist.C+(12) (via ACLIC, fast, independent fits with Minuit2)

Note that you can also execute this script in batch with eg,

root -b -q "fit2dHist.C+(12)"
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
float * q
Definition: THbookFile.cxx:89

or execute interactively from the shell

root fit2dHist.C+
root "fit2dHist.C+(12)"
FCN=2613.61 FROM MIGRAD STATUS=CONVERGED 1090 CALLS 1091 TOTAL
EDM=1.5599e-08 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 3.7 per cent
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 5.34104e+02 2.25626e+00 -1.45532e-03 2.18435e-05
2 p1 6.00014e+00 5.67524e-03 -2.38292e-06 -3.27949e-02
3 p2 1.98724e+00 3.63694e-03 -2.58780e-06 -5.50559e-03
4 p3 7.02973e+00 2.65118e-02 -2.44324e-05 -1.18794e-02
5 p4 2.99679e+00 1.39392e-02 -1.13807e-05 1.50299e-02
6 p5 5.19346e+02 5.08272e+01 3.87334e-02 -2.41684e-05
7 p6 1.15499e+01 4.81865e-01 5.78545e-04 4.63910e-03
8 p7 2.72921e+00 2.57821e-01 2.95923e-04 -5.50344e-03
9 p8 1.11977e+01 2.40323e-01 -9.65097e-05 7.16022e-03
10 p9 2.08422e+00 1.01013e-01 -2.43739e-05 -1.10303e-02
FCN=2220.46 FROM MIGRAD STATUS=CONVERGED 333 CALLS 334 TOTAL
EDM=6.12528e-07 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 1.1 per cent
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 5.30875e+02 1.56318e+00 -8.45958e-04 2.13120e-04
2 p1 6.01215e+00 1.39029e-02 5.56975e-05 1.20143e-01
3 p2 1.99424e+00 1.02676e-02 -3.65137e-05 1.99143e-01
4 p3 6.98634e+00 1.77537e-02 4.59097e-06 1.88058e-02
5 p4 2.98764e+00 1.14564e-02 5.16037e-06 2.41585e-02
6 p5 5.32751e+02 1.16044e+00 9.85153e-04 1.59279e-03
7 p6 1.19894e+01 8.92126e-03 7.18458e-06 8.54804e-02
8 p7 2.99536e+00 6.32688e-03 -5.28473e-06 3.40897e-01
9 p8 1.09975e+01 3.41959e-03 -1.79221e-06 2.14024e-02
10 p9 1.98880e+00 2.41489e-03 5.35898e-07 3.99846e-01
(int) 0
#include "TH2D.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "Fit/Fitter.h"
#include "TList.h"
#include <iostream>
double gauss2D(double *x, double *par) {
double z1 = double((x[0]-par[1])/par[2]);
double z2 = double((x[1]-par[3])/par[4]);
return par[0]*exp(-0.5*(z1*z1+z2*z2));
}
double my2Dfunc(double *x, double *par) {
return gauss2D(x,&par[0]) + gauss2D(x,&par[5]);
}
class MyFcn {
public:
TH2D *h1 = nullptr;
TH2D *h2 = nullptr;
int npfits = 0;
MyFcn(TH2D * _h1, TH2D * _h2) :
h1(_h1), h2(_h2) {}
double operator()(const double *p) {
TAxis *xaxis1 = h1->GetXaxis();
TAxis *yaxis1 = h1->GetYaxis();
TAxis *xaxis2 = h2->GetXaxis();
TAxis *yaxis2 = h2->GetYaxis();
int nbinX1 = h1->GetNbinsX();
int nbinY1 = h1->GetNbinsY();
int nbinX2 = h2->GetNbinsX();
int nbinY2 = h2->GetNbinsY();
double chi2 = 0;
double x[2];
double tmp;
npfits = 0;
for (int ix = 1; ix <= nbinX1; ++ix) {
x[0] = xaxis1->GetBinCenter(ix);
for (int iy = 1; iy <= nbinY1; ++iy) {
if ( h1->GetBinError(ix,iy) > 0 ) {
x[1] = yaxis1->GetBinCenter(iy);
tmp = (h1->GetBinContent(ix,iy) - my2Dfunc(x,(double*) p))/h1->GetBinError(ix,iy);
chi2 += tmp*tmp;
npfits++;
}
}
}
for (int ix = 1; ix <= nbinX2; ++ix) {
x[0] = xaxis2->GetBinCenter(ix);
for (int iy = 1; iy <= nbinY2; ++iy) {
if ( h2->GetBinError(ix,iy) > 0 ) {
x[1] = yaxis2->GetBinCenter(iy);
tmp = (h2->GetBinContent(ix,iy) - my2Dfunc(x,(double *) p))/h2->GetBinError(ix,iy);
chi2 += tmp*tmp;
npfits++;
}
}
}
return chi2;
}
};
void FillHisto(TH2D * h, int n, double * p) {
const double mx1 = p[1];
const double my1 = p[3];
const double sx1 = p[2];
const double sy1 = p[4];
const double mx2 = p[6];
const double my2 = p[8];
const double sx2 = p[7];
const double sy2 = p[9];
//const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2);
const double w1 = 0.5;
double x, y;
for (int i = 0; i < n; ++i) {
// generate randoms with larger Gaussians
double r = gRandom->Rndm(1);
if (r < w1) {
x = x*sx1 + mx1;
y = y*sy1 + my1;
}
else {
x = x*sx2 + mx2;
y = y*sy2 + my2;
}
h->Fill(x,y);
}
}
int fit2dHist(int option=1) {
// create two histograms
int nbx1 = 50;
int nby1 = 50;
int nbx2 = 50;
int nby2 = 50;
double xlow1 = 0.;
double ylow1 = 0.;
double xup1 = 10.;
double yup1 = 10.;
double xlow2 = 5.;
double ylow2 = 5.;
double xup2 = 20.;
double yup2 = 20.;
auto h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
auto h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2);
double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. };
// create fit function
TF2 * func = new TF2("func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);
func->SetParameters(iniParams);
// fill Histos
int n1 = 1000000;
int n2 = 1000000;
FillHisto(h1,n1,iniParams);
FillHisto(h2,n2,iniParams);
// scale histograms to same heights (for fitting)
double dx1 = (xup1-xlow1)/double(nbx1);
double dy1 = (yup1-ylow1)/double(nby1);
double dx2 = (xup2-xlow2)/double(nbx2);
double dy2 = (yup2-ylow2)/double(nby2);
// scale histo 2 to scale of 1
h2->Sumw2();
h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) );
bool global = false;
if (option > 10) global = true;
// do global combined fit
if (global) {
// fill data structure for fit (coordinates + values + errors)
std::cout << "Do global fit" << std::endl;
// fit now all the function together
//The default minimizer is Minuit, you can also try Minuit2
MyFcn myFcn(h1,h2);
fitter.SetFCN(10, myFcn);
if (option%10 == 2) fitter.Config().SetMinimizer("Minuit2");
// set parameter initial value, name and step size
for (int i = 0; i < 10; ++i) {
fitter.Config().ParSettings(i) = ROOT::Fit::ParameterSettings(func->GetParName(i), func->GetParameter(i), 0.01);
}
bool ret = fitter.FitFCN();
if (!ret) {
Error("fit2DHist","Fit Failed to converge");
return -1;
}
//get result
double minParams[10];
double parErrors[10];
for (int i = 0; i < 10; ++i) {
minParams[i] = fitter.Result().Parameter(i);
parErrors[i] = fitter.Result().Error(i);
}
double chi2 = fitter.Result().MinFcnValue();
double edm = fitter.Result().Edm();
int npfits = myFcn.npfits;
func->SetParameters(minParams);
func->SetParErrors(parErrors);
func->SetChisquare(chi2);
int ndf = npfits - fitter.Result().NFreeParameters();
func->SetNDF(ndf);
// add to list of functions
h2->GetListOfFunctions()->Add(func);
}
else {
// fit independently
h1->Fit(func);
h2->Fit(func);
}
// Create a new canvas.
TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800);
c1->Divide(2,2);
c1->cd(1);
h1->Draw();
func->SetRange(xlow1,ylow1,xup1,yup1);
func->DrawCopy("cont1 same");
c1->cd(2);
h1->Draw("lego");
func->DrawCopy("surf1 same");
c1->cd(3);
func->SetRange(xlow2,ylow2,xup2,yup2);
h2->Draw();
func->DrawCopy("cont1 same");
c1->cd(4);
h2->Draw("lego");
gPad->SetLogz();
func->Draw("surf1 same");
return 0;
}
#define h(i)
Definition: RSha256.hxx:106
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
TRObject operator()(const T1 &t1) const
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
R__EXTERN TStyle * gStyle
Definition: TStyle.h:414
#define gPad
Definition: TVirtualPad.h:288
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:181
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition: FitConfig.h:76
double Error(unsigned int i) const
parameter error by index
Definition: FitResult.h:186
double MinFcnValue() const
Return value of the objective function (chi2 or likelihood) used in the fit.
Definition: FitResult.h:120
double Edm() const
Expected distance from minimum.
Definition: FitResult.h:126
unsigned int NFreeParameters() const
get total number of free parameters
Definition: FitResult.h:134
double Parameter(unsigned int i) const
parameter value by index
Definition: FitResult.h:181
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:77
bool FitFCN(unsigned int npar, Function &fcn, const double *params=0, unsigned int dataSize=0, bool chi2fit=false)
Fit using the a generic FCN function as a C++ callable object implementing double () (const double *)...
Definition: Fitter.h:625
bool SetFCN(unsigned int npar, Function &fcn, const double *params=0, unsigned int dataSize=0, bool chi2fit=false)
Set a generic FCN function as a C++ callable object implementing double () (const double *) Note that...
Definition: Fitter.h:632
const FitResult & Result() const
get fit result
Definition: Fitter.h:382
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:410
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:478
The Canvas class.
Definition: TCanvas.h:23
virtual void SetNDF(Int_t ndf)
Set the number of degrees of freedom ndf should be the number of points used in a fit - the number of...
Definition: TF1.cxx:3463
virtual void SetChisquare(Double_t chi2)
Definition: TF1.h:617
virtual void SetParErrors(const Double_t *errors)
Set errors for all active parameters when calling this function, the array errors must have at least ...
Definition: TF1.cxx:3528
virtual const char * GetParName(Int_t ipar) const
Definition: TF1.h:534
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:649
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:517
A 2-Dim function with parameters.
Definition: TF2.h:29
TF1 * DrawCopy(Option_t *option="") const override
Draw a copy of this function with its current attributes-*.
Definition: TF2.cxx:286
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition: TF2.cxx:259
void SetRange(Double_t xmin, Double_t xmax) override
Initialize the upper and lower bounds to draw the function.
Definition: TF2.h:148
virtual Int_t GetNbinsY() const
Definition: TH1.h:296
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8930
TAxis * GetXaxis()
Definition: TH1.h:319
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:3894
virtual Int_t GetNbinsX() const
Definition: TH1.h:295
TAxis * GetYaxis()
Definition: TH1.h:320
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3060
TList * GetListOfFunctions() const
Definition: TH1.h:242
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5025
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6586
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8887
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:300
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition: TH2.h:89
void Add(TObject *obj) override
Definition: TList.h:81
Double_t Rndm() override
Machine independent random number generator.
Definition: TRandom.cxx:552
virtual void Rannor(Float_t &a, Float_t &b)
Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
Definition: TRandom.cxx:500
void SetStatY(Float_t y=0)
Definition: TStyle.h:382
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:1542
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1761
Double_t y[n]
Definition: legend1.C:17
return c1
Definition: legend1.C:41
const Int_t n
Definition: legend1.C:16
TH1F * h1
Definition: legend1.C:5
Authors
: Lorenzo Moneta, Rene Brun 18/01/2006

Definition in file fit2dHist.C.