Example to fit two histograms at the same time.
Processing /mnt/build/workspace/root-makedoc-v610/rootspi/rdoc/src/v6-10-00-patches/tutorials/fit/TwoHistoFit2D.C...
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.6486e-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.29664e-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.33229e-02
4 p3 7.09821e+00 3.32869e-02 2.39025e-04 2.42668e-02
5 p4 2.94271e+00 2.42010e-02 -1.88552e-04 2.78492e-03
6 p5 2.63145e+01 2.69272e-01 -2.31447e-03 -2.60065e-03
7 p6 1.19850e+01 3.51596e-02 4.24095e-04 -3.93611e-02
8 p7 2.90086e+00 2.64547e-02 8.06256e-05 -5.19640e-03
9 p8 1.09762e+01 1.47334e-02 -6.74373e-05 -1.09628e-02
10 p9 1.95760e+00 1.14466e-02 2.85422e-05 -1.15592e-01
Chi2 Fit = 4015.63 ndf = 3921 3921
(int) 0
#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);
}
std::vector<std::pair<double, double> > coords;
std::vector<double > values;
std::vector<double > errors;
{
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;
}
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 = 0.5;
for (
int i = 0; i <
n; ++i) {
if (r < w1) {
x = x*sx1 + mx1;
y = y*sy1 + my1;
}
else {
x = x*sx2 + mx2;
y = y*sy2 + my2;
}
}
}
int TwoHistoFit2D(bool global = true) {
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. };
TF2 * func =
new TF2(
"func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10);
int n1 = 50000;
int n2 = 50000;
FillHisto(h1,n1,iniParams);
FillHisto(h2,n2,iniParams);
double dx1 = (xup1-xlow1)/double(nbx1);
double dy1 = (yup1-ylow1)/double(nby1);
double dx2 = (xup2-xlow2)/double(nbx2);
double dy2 = (yup2-ylow2)/double(nby2);
h2->
Scale( (
double(n1) * dx1 * dy1 ) / (
double(n2) * dx2 * dy2 ) );
if (global) {
std::cout << "Do global fit" << std::endl;
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) {
}
}
}
for (int ix = 1; ix <= nbinX2; ++ix) {
for (int iy = 1; iy <= nbinY2; ++iy) {
}
}
}
for (int i = 0; i < 10; ++i) {
}
double arglist[100];
arglist[0] = 0;
arglist[0] = 5000;
arglist[1] = 0.01;
double minParams[10];
double parErrors[10];
for (int i = 0; i < 10; ++i) {
}
double chi2, edm, errdef;
int nvpar, nparx;
minuit->
GetStats(chi2,edm,errdef,nvpar,nparx);
int ndf = coords.size()-nvpar;
std::cout <<
"Chi2 Fit = " << chi2 <<
" ndf = " << ndf <<
" " << func->
GetNDF() << std::endl;
}
else {
}
TCanvas * c1 =
new TCanvas(
"c1",
"Two HIstogram Fit example",100,10,900,800);
func->
Draw(
"surf1 same");
return 0;
}
- Author
- Rene Brun
Definition in file TwoHistoFit2D.C.