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 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

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;
}
double
ROOT::R::TRInterface & r
Definition Object.C:4
#define h(i)
Definition RSha256.hxx:106
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
R__EXTERN TStyle * gStyle
Definition TStyle.h:413
#define gPad
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:3439
virtual void SetChisquare(Double_t chi2)
Definition TF1.h:612
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:3504
virtual const char * GetParName(Int_t ipar) const
Definition TF1.h:529
virtual void SetParameters(const Double_t *params)
Definition TF1.h:644
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:512
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:287
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition TF2.h:148
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition TF2.cxx:260
virtual Int_t GetNbinsY() const
Definition TH1.h:297
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:8893
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320
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:3893
virtual Int_t GetNbinsX() const
Definition TH1.h:296
TAxis * GetYaxis()
Definition TH1.h:321
TList * GetListOfFunctions() const
Definition TH1.h:243
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3074
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4994
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6553
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition TH1.cxx:8850
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:81
Random number generator class based on M.
Definition TRandom3.h:27
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom3.cxx:99
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:381
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:1541
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.