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 TVirtualFitter

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)"
#define b(i)
Definition: RSha256.hxx:100
float * q
Definition: THbookFile.cxx:87

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 "TVirtualFitter.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]);
}
// data need to be globals to be visible by fcn
TRandom3 rndm;
TH2D *h1, *h2;
Int_t npfits;
void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */ )
{
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,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,p))/h2->GetBinError(ix,iy);
chi2 += tmp*tmp;
npfits++;
}
}
}
fval = 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
rndm.Rannor(x,y);
double r = rndm.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.;
h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
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;
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
if (option%10 == 2) TVirtualFitter::SetDefaultFitter("Minuit2");
for (int i = 0; i < 10; ++i) {
minuit->SetParameter(i, func->GetParName(i), func->GetParameter(i), 0.01, 0,0);
}
minuit->SetFCN(myFcn);
double arglist[100];
arglist[0] = 0;
// set print level
minuit->ExecuteCommand("SET PRINT",arglist,2);
// minimize
arglist[0] = 5000; // number of function calls
arglist[1] = 0.01; // tolerance
minuit->ExecuteCommand("MIGRAD",arglist,2);
//get result
double minParams[10];
double parErrors[10];
for (int i = 0; i < 10; ++i) {
minParams[i] = minuit->GetParameter(i);
parErrors[i] = minuit->GetParError(i);
}
double chi2, edm, errdef;
int nvpar, nparx;
minuit->GetStats(chi2,edm,errdef,nvpar,nparx);
func->SetParameters(minParams);
func->SetParErrors(parErrors);
func->SetChisquare(chi2);
int ndf = npfits-nvpar;
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;
}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
double exp(double)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:407
#define gPad
Definition: TVirtualPad.h:286
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
The Canvas class.
Definition: TCanvas.h:31
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:3418
virtual void SetChisquare(Double_t chi2)
Definition: TF1.h:606
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:3483
virtual const char * GetParName(Int_t ipar) const
Definition: TF1.h:523
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:638
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:506
A 2-Dim function with parameters.
Definition: TF2.h:29
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes-*.
Definition: TF2.cxx:268
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF2.h:150
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF2.cxx:241
virtual Int_t GetNbinsY() const
Definition: TH1.h:293
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8507
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
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 GetNbinsX() const
Definition: TH1.h:292
TAxis * GetYaxis()
Definition: TH1.h:317
TList * GetListOfFunctions() const
Definition: TH1.h:239
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4899
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6234
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8464
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:88
virtual void Add(TObject *obj)
Definition: TList.h:87
Random number generator class based on M.
Definition: TRandom3.h:27
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:100
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:489
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
Abstract Base Class for Fitting.
static void SetDefaultFitter(const char *name="")
static: set name of default fitter
virtual Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const =0
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
virtual Double_t GetParError(Int_t ipar) const =0
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
virtual Double_t GetParameter(Int_t ipar) const =0
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
return c1
Definition: legend1.C:41
Double_t y[n]
Definition: legend1.C:17
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.