Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 16.85
**********
**********
** 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 & /*nPar*/, double * /*grad*/ , double &fval, double *p, int /*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) );
}
}
}
TVirtualFitter * minuit = TVirtualFitter::Fitter(nullptr,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 = 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;
}
#define h(i)
Definition RSha256.hxx:106
winID h TVirtualViewer3D TVirtualGLPainter p
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
R__EXTERN TStyle * gStyle
Definition TStyle.h:436
#define gPad
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 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:1889
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:640
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:557
virtual void SetParameters(const Double_t *params)
Definition TF1.h:677
virtual Double_t GetParameter(Int_t ipar) const
Definition TF1.h:540
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:9063
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:3898
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:3066
TList * GetListOfFunctions() const
Definition TH1.h:244
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5061
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:357
void Add(TObject *obj) override
Definition TList.h:83
Random number generator class based on M.
Definition TRandom3.h:27
Double_t Rndm() override
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:507
void SetStatY(Float_t y=0)
Definition TStyle.h:398
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:1593
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.
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
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.