Logo ROOT   6.10/09
Reference Guide
fit2dHist.C File Reference
Tutorials » Fit Tutorials

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)

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

or execute interactively from the shell

root fit2dHist.C+
root "fit2dHist.C+(12)"
pict1_fit2dHist.C.png
Processing /mnt/build/workspace/root-makedoc-v610/rootspi/rdoc/src/v6-10-00-patches/tutorials/fit/fit2dHist.C...
FCN=2613.61 FROM MIGRAD STATUS=CONVERGED 1090 CALLS 1091 TOTAL
EDM=1.55989e-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.45533e-03 2.18412e-05
2 p1 6.00014e+00 5.67524e-03 -2.38293e-06 -3.27949e-02
3 p2 1.98724e+00 3.63694e-03 -2.58781e-06 -5.50660e-03
4 p3 7.02973e+00 2.65118e-02 -2.44325e-05 -1.18793e-02
5 p4 2.99679e+00 1.39392e-02 -1.13807e-05 1.50297e-02
6 p5 5.19346e+02 5.08272e+01 3.87335e-02 -2.41686e-05
7 p6 1.15499e+01 4.81865e-01 5.78546e-04 4.63913e-03
8 p7 2.72921e+00 2.57821e-01 2.95924e-04 -5.50348e-03
9 p8 1.11977e+01 2.40323e-01 -9.65099e-05 7.16026e-03
10 p9 2.08422e+00 1.01013e-01 -2.43739e-05 -1.10304e-02
FCN=2220.46 FROM MIGRAD STATUS=CONVERGED 333 CALLS 334 TOTAL
EDM=6.12518e-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.45948e-04 2.13124e-04
2 p1 6.01215e+00 1.39029e-02 5.56970e-05 1.20140e-01
3 p2 1.99424e+00 1.02676e-02 -3.65134e-05 1.99145e-01
4 p3 6.98634e+00 1.77537e-02 4.59085e-06 1.88054e-02
5 p4 2.98764e+00 1.14564e-02 5.16037e-06 2.41589e-02
6 p5 5.32751e+02 1.16044e+00 9.85146e-04 1.59278e-03
7 p6 1.19894e+01 8.92126e-03 7.18450e-06 8.54783e-02
8 p7 2.99536e+00 6.32688e-03 -5.28469e-06 3.40895e-01
9 p8 1.09975e+01 3.41959e-03 -1.79220e-06 2.14026e-02
10 p9 1.98880e+00 2.41489e-03 5.35908e-07 3.99849e-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");
else TVirtualFitter::SetDefaultFitter("Minuit");
TVirtualFitter * minuit = TVirtualFitter::Fitter(0,10);
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
h1->GetListOfFunctions()->Add(func);
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);
gStyle->SetOptFit();
gStyle->SetStatY(0.6);
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;
}
Authors
: Lorenzo Moneta, Rene Brun 18/01/2006

Definition in file fit2dHist.C.