ROOT   Reference Guide
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
FCN=49.5952 FROM MIGRAD STATUS=CONVERGED 52 CALLS 53 TOTAL
EDM=1.22682e-09 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.5 per cent
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 3.13201e+00 3.12699e-02 -3.64656e-05 2.15221e-03
2 p1 2.97626e+01 1.00773e+00 6.67621e-05 -4.02033e-06
Covariance matrix from the fit
2x2 matrix is as follows
| 0 | 1 |
-------------------------------
0 | 0.0009778 0.009128
1 | 0.009128 1.016
Integral = 19.005 +/- 0.6159
#include "TF1.h"
#include "TH1D.h"
#include "TFitResult.h"
#include "TMath.h"
#include <assert.h>
#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
double cos(double)
double sqrt(double)
double sin(double)
1-Dim function class
Definition: TF1.h:213
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2515
virtual Double_t * GetParameters() const
Definition: TF1.h:520
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:644
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:634
virtual Double_t IntegralError(Double_t a, Double_t b, const Double_t *params=0, const Double_t *covmat=0, 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:2692
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
virtual void FillRandom(const char *fname, Int_t ntimes=5000, TRandom *rng=nullptr)
Fill histogram following distribution in function fname.
Definition: TH1.cxx:3525
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:3892
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3073
Double_t x[n]
Definition: legend1.C:17
TH1F * h1
Definition: legend1.C:5
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Double_t Sin(Double_t)
Definition: TMath.h:639

Definition in file ErrorIntegral.C.