Logo ROOT  
Reference Guide
TwoHistoFit2D.C File Reference

Detailed Description

View in nbviewer Open in SWAN Example to fit two histograms at the same time.

Do global fit
PARAMETER DEFINITIONS:
NO. NAME VALUE STEP SIZE LIMITS
1 p0 1.00000e+02 1.00000e-02 no limits
2 p1 6.00000e+00 1.00000e-02 no limits
3 p2 2.00000e+00 1.00000e-02 no limits
4 p3 7.00000e+00 1.00000e-02 no limits
5 p4 3.00000e+00 1.00000e-02 no limits
6 p5 1.00000e+02 1.00000e-02 no limits
7 p6 1.20000e+01 1.00000e-02 no limits
8 p7 3.00000e+00 1.00000e-02 no limits
9 p8 1.10000e+01 1.00000e-02 no limits
10 p9 2.00000e+00 1.00000e-02 no limits
**********
** 1 **SET PRINT 0 3.096e-314
**********
**********
** 2 **MIGRAD 5000 0.01
**********
MIGRAD MINIMIZATION HAS CONVERGED.
FCN=4015.63 FROM MIGRAD STATUS=CONVERGED 525 CALLS 526 TOTAL
EDM=7.64858e-07 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 4.8 per cent
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 2.55114e+01 2.22488e-01 1.18177e-03 1.29669e-03
2 p1 6.03551e+00 1.56999e-02 1.78147e-04 5.19787e-02
3 p2 1.95953e+00 1.34972e-02 1.02338e-04 -2.33217e-02
4 p3 7.09821e+00 3.32869e-02 2.39024e-04 2.42669e-02
5 p4 2.94271e+00 2.42010e-02 -1.88552e-04 2.78529e-03
6 p5 2.63145e+01 2.69272e-01 -2.31447e-03 -2.60063e-03
7 p6 1.19850e+01 3.51596e-02 4.24094e-04 -3.93616e-02
8 p7 2.90086e+00 2.64547e-02 8.06260e-05 -5.19624e-03
9 p8 1.09762e+01 1.47334e-02 -6.74372e-05 -1.09627e-02
10 p9 1.95760e+00 1.14466e-02 2.85422e-05 -1.15591e-01
Chi2 Fit = 4015.63 ndf = 3921 3921
(int) 0
#include "TH2D.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "TVirtualFitter.h"
#include "TList.h"
#include <vector>
#include <map>
#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) {
double *p1 = &par[0];
double *p2 = &par[5];
return gauss2D(x,p1) + gauss2D(x,p2);
}
// data need to be globals to be visible by fcn
std::vector<std::pair<double, double> > coords;
std::vector<double > values;
std::vector<double > errors;
void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */ )
{
int n = coords.size();
double chi2 = 0;
double tmp,x[2];
for (int i = 0; i <n; ++i ) {
x[0] = coords[i].first;
x[1] = coords[i].second;
tmp = ( values[i] - my2Dfunc(x,p))/errors[i];
chi2 += tmp*tmp;
}
fval = chi2;
}
TRandom3 rndm;
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 TwoHistoFit2D(bool global = true) {
// 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.;
TH2D * h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1);
TH2D * 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 = 50000;
int n2 = 50000;
// h1->FillRandom("func", n1);
//h2->FillRandom("func",n2);
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);
// h1->Sumw2();
// h1->Scale( 1.0 / ( n1 * dx1 * dy1 ) );
// scale histo 2 to scale of 1
h2->Sumw2();
h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) );
if (global) {
// fill data structure for fit (coordinates + values + errors)
std::cout << "Do global fit" << std::endl;
// fit now all the function together
// fill data structure for fit (coordinates + values + errors)
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();
/// reset data structure
coords = std::vector<std::pair<double,double> >();
values = std::vector<double>();
errors = std::vector<double>();
for (int ix = 1; ix <= nbinX1; ++ix) {
for (int iy = 1; iy <= nbinY1; ++iy) {
if ( h1->GetBinContent(ix,iy) > 0 ) {
coords.push_back( std::make_pair(xaxis1->GetBinCenter(ix), yaxis1->GetBinCenter(iy) ) );
values.push_back( h1->GetBinContent(ix,iy) );
errors.push_back( h1->GetBinError(ix,iy) );
}
}
}
for (int ix = 1; ix <= nbinX2; ++ix) {
for (int iy = 1; iy <= nbinY2; ++iy) {
if ( h2->GetBinContent(ix,iy) > 0 ) {
coords.push_back( std::make_pair(xaxis2->GetBinCenter(ix), yaxis2->GetBinCenter(iy) ) );
values.push_back( h2->GetBinContent(ix,iy) );
errors.push_back( h2->GetBinError(ix,iy) );
}
}
}
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 = coords.size()-nvpar;
func->SetNDF(ndf);
std::cout << "Chi2 Fit = " << chi2 << " ndf = " << ndf << " " << func->GetNDF() << std::endl;
// 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
Definition: Converters.cxx:921
ROOT::R::TRInterface & r
Definition: Object.C:4
#define h(i)
Definition: RSha256.hxx:106
int Int_t
Definition: RtypesCore.h:43
double Double_t
Definition: RtypesCore.h:57
double exp(double)
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
#define gPad
Definition: TVirtualPad.h:287
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:475
The Canvas class.
Definition: TCanvas.h:27
virtual Int_t GetNDF() const
Return the number of degrees of freedom in the fit the fNDF parameter has been previously computed du...
Definition: TF1.cxx:1872
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:3421
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:3486
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:8519
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:4907
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6246
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8476
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:379
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
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
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TH1F * h1
Definition: legend1.C:5
Author
Rene Brun

Definition in file TwoHistoFit2D.C.