ROOT   Reference Guide
Searching...
No Matches
ErrorIntegral.C File Reference

## Detailed Description

Estimate the error in the integral of a fitted function taking into account the errors in the parameters resulting from the fit.

The error is estimated also using the correlations values obtained from the fit

run the macro doing:

.x ErrorIntegral.C

After having computed the integral and its error using the integral and the integral error using the generic functions TF1::Integral and TF1::IntegralError, we compute the integrals and its error analytically using the fact that the fitting function is $$f(x) = p[1]* sin(p[0]*x)$$.

Therefore we have:

• integral in [0,1] : ic = p[1]* (1-std::cos(p[0]) )/p[0]
• derivative of integral with respect to p0: c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0]
• derivative of integral with respect to p1: c1c = (1-std::cos(p[0]) )/p[0]

and then we can compute the integral error using error propagation and the covariance matrix for the parameters p obtained from the fit.

integral error : sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1) + 2.* c0c*c1c * covMatrix(0,1))

Note that, if possible, one should fit directly the function integral, which are the number of events of the different components (e.g. signal and background). In this way one obtains a better and more correct estimate of the integrals uncertainties, since they are obtained directly from the fit without using the approximation of error propagation. This is possible in ROOT. when using the TF1NormSum class, see the tutorial fitNormSum.C

****************************************
Chi2 = 49.5952
NDf = 48
Edm = 1.61787e-06
NCalls = 58
p0 = 3.13205 +/- 0.0312726
p1 = 29.7625 +/- 1.00876
Covariance matrix from the fit
2x2 matrix is as follows
| 0 | 1 |
-------------------------------
0 | 0.000978 0.009147
1 | 0.009147 1.018
Integral = 19.0047 +/- 0.616472
#include "TF1.h"
#include "TH1D.h"
#include "TFitResult.h"
#include "TMath.h"
#include <cassert>
#include <iostream>
#include <cmath>
TF1 * fitFunc; // fit function pointer
const int NPAR = 2; // number of function parameters;
//____________________________________________________________________
double f(double * x, double * p) {
// function used to fit the data
return p[1]*TMath::Sin( p[0] * x[0] );
}
//____________________________________________________________________
void ErrorIntegral() {
fitFunc = new TF1("f",f,0,1,NPAR);
TH1D * h1 = new TH1D("h1","h1",50,0,1);
double par[NPAR] = { 3.14, 1.};
fitFunc->SetParameters(par);
h1->FillRandom("f",1000); // fill histogram sampling fitFunc
fitFunc->SetParameter(0,3.); // vary a little the parameters
auto fitResult = h1->Fit(fitFunc,"S"); // fit the histogram and get fit result pointer
h1->Draw();
/* calculate the integral*/
double integral = fitFunc->Integral(0,1);
auto covMatrix = fitResult->GetCovarianceMatrix();
std::cout << "Covariance matrix from the fit ";
covMatrix.Print();
// need to pass covariance matrix to fit result.
// Parameters values are are stored inside the function but we can also retrieve from TFitResult
double sigma_integral = fitFunc->IntegralError(0,1, fitResult->GetParams() , covMatrix.GetMatrixArray());
std::cout << "Integral = " << integral << " +/- " << sigma_integral
<< std::endl;
// estimated integral and error analytically
double * p = fitFunc->GetParameters();
double ic = p[1]* (1-std::cos(p[0]) )/p[0];
double c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0];
double c1c = (1-std::cos(p[0]) )/p[0];
// estimated error with correlations
double sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1)
+ 2.* c0c*c1c * covMatrix(0,1));
if ( std::fabs(sigma_integral-sic) > 1.E-6*sic )
std::cout << " ERROR: test failed : different analytical integral : "
<< ic << " +/- " << sic << std::endl;
}
#define f(i)
Definition RSha256.hxx:104
winID h TVirtualViewer3D TVirtualGLPainter p
1-Dim function class
Definition TF1.h:233
virtual Double_t IntegralError(Double_t a, Double_t b, const Double_t *params=nullptr, const Double_t *covmat=nullptr, Double_t epsilon=1.E-2)
Return Error on Integral of a parametric function between a and b due to the parameter uncertainties ...
Definition TF1.cxx:2708
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition TF1.cxx:2531
virtual Double_t * GetParameters() const
Definition TF1.h:546
virtual void SetParameters(const Double_t *params)
Definition TF1.h:670
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:660
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:669
virtual void FillRandom(const char *fname, Int_t ntimes=5000, TRandom *rng=nullptr)
Fill histogram following distribution in function fname.
Definition TH1.cxx:3519
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
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3066
Double_t x[n]
Definition legend1.C:17
TH1F * h1
Definition legend1.C:5
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:588

Definition in file ErrorIntegral.C.