Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 cling, 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)"
#define b(i)
Definition RSha256.hxx:100
float * q

or execute interactively from the shell

root fit2dHist.C+
root "fit2dHist.C+(12)"
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 2613.61
NDf = 2478
Edm = 2.04928e-05
NCalls = 1133
p0 = 534.106 +/- 2.22027
p1 = 6.00013 +/- 0.00562856
p2 = 1.98723 +/- 0.00361514
p3 = 7.02975 +/- 0.0260609
p4 = 2.9968 +/- 0.0137019
p5 = 519.282 +/- 50.0512
p6 = 11.5486 +/- 0.456627
p7 = 2.72866 +/- 0.245596
p8 = 11.1983 +/- 0.233757
p9 = 2.08449 +/- 0.0982633
****************************************
Minimizer is Minuit2 / Migrad
Chi2 = 2220.46
NDf = 2231
Edm = 2.08188e-06
NCalls = 311
p0 = 530.876 +/- 1.55657
p1 = 6.0121 +/- 0.00596885
p2 = 1.99427 +/- 0.0055483
p3 = 6.98634 +/- 0.0177191
p4 = 2.98763 +/- 0.0113697
p5 = 532.75 +/- 1.15393
p6 = 11.9894 +/- 0.00876067
p7 = 2.99537 +/- 0.00613849
p8 = 10.9975 +/- 0.00338646
p9 = 1.98879 +/- 0.00240993
(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, "0");
h2->Fit(func, "0");
}
// 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("cont3 same");
c1->cd(2);
h1->Draw("lego");
func->DrawCopy("surf1 same");
c1->cd(3);
func->SetRange(xlow2,ylow2,xup2,yup2);
h2->Draw();
func->DrawCopy("cont3 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:185
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:433
#define gPad
void SetMinimizer(const char *type, const char *algo=nullptr)
set minimizer type
Definition FitConfig.h:179
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:179
double MinFcnValue() const
Return value of the objective function (chi2 or likelihood) used in the fit.
Definition FitResult.h:111
double Edm() const
Expected distance from minimum.
Definition FitResult.h:117
unsigned int NFreeParameters() const
get total number of free parameters
Definition FitResult.h:125
double Parameter(unsigned int i) const
parameter value by index
Definition FitResult.h:174
Fitter class, entry point for performing all type of fits.
Definition Fitter.h:77
const FitResult & Result() const
get fit result
Definition Fitter.h:394
bool FitFCN(unsigned int npar, Function &fcn, const double *params=nullptr, unsigned int dataSize=0, int fitType=0)
Fit using the a generic FCN function as a C++ callable object implementing double () (const double *)...
Definition Fitter.h:649
const FitConfig & Config() const
access to the fit configuration (const method)
Definition Fitter.h:422
bool SetFCN(unsigned int npar, Function &fcn, const double *params=nullptr, unsigned int dataSize=0, int fitType=0)
Set a generic FCN function as a C++ callable object implementing double () (const double *) Note that...
Definition Fitter.h:656
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
Class to manage histogram axis.
Definition TAxis.h:31
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:3419
virtual void SetChisquare(Double_t chi2)
Definition TF1.h:638
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:3490
virtual const char * GetParName(Int_t ipar) const
Definition TF1.h:555
virtual void SetParameters(const Double_t *params)
Definition TF1.h:670
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:538
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:146
virtual Int_t GetNbinsY() const
Definition TH1.h:298
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9027
TAxis * GetXaxis()
Definition TH1.h:324
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:297
TAxis * GetYaxis()
Definition TH1.h:325
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
TList * GetListOfFunctions() const
Definition TH1.h:244
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5025
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:338
void Add(TObject *obj) override
Definition TList.h:81
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
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:507
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
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.