Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
ErrorIntegral.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_fit
3/// \notebook -js
4/// Estimate the error in the integral of a fitted function
5/// taking into account the errors in the parameters resulting from the fit.
6/// The error is estimated also using the correlations values obtained from
7/// the fit
8///
9/// run the macro doing:
10///
11/// ~~~{.cpp}
12/// .x ErrorIntegral.C
13/// ~~~
14///
15/// After having computed the integral and its error using the integral and the integral
16/// error using the generic functions TF1::Integral and TF1::IntegralError, we compute
17/// the integrals and its error analytically using the fact that the fitting function is
18/// \f$ f(x) = p[1]* sin(p[0]*x) \f$.
19///
20/// Therefore we have:
21/// - integral in [0,1] : `ic = p[1]* (1-std::cos(p[0]) )/p[0]`
22/// - derivative of integral with respect to p0:
23/// `c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0]`
24/// - derivative of integral with respect to p1:
25/// `c1c = (1-std::cos(p[0]) )/p[0]`
26///
27/// and then we can compute the integral error using error propagation and the covariance
28/// matrix for the parameters p obtained from the fit.
29///
30/// integral error : `sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1) + 2.* c0c*c1c * covMatrix(0,1))`
31///
32/// Note that, if possible, one should fit directly the function integral, which are the
33/// number of events of the different components (e.g. signal and background).
34/// In this way one obtains a better and more correct estimate of the integrals uncertainties,
35/// since they are obtained directly from the fit without using the approximation of error propagation.
36/// This is possible in ROOT. when using the TF1NormSum class, see the tutorial fitNormSum.C
37///
38/// \macro_image
39/// \macro_output
40/// \macro_code
41///
42/// \author Lorenzo Moneta
43
44#include "TF1.h"
45#include "TH1D.h"
46#include "TFitResult.h"
47#include "TMath.h"
48#include <cassert>
49#include <iostream>
50#include <cmath>
51
52TF1 * fitFunc; // fit function pointer
53
54const int NPAR = 2; // number of function parameters;
55
56//____________________________________________________________________
57double f(double * x, double * p) {
58 // function used to fit the data
59 return p[1]*TMath::Sin( p[0] * x[0] );
60}
61
62//____________________________________________________________________
63void ErrorIntegral() {
64 fitFunc = new TF1("f",f,0,1,NPAR);
65 TH1D * h1 = new TH1D("h1","h1",50,0,1);
66
67 double par[NPAR] = { 3.14, 1.};
68 fitFunc->SetParameters(par);
69
70 h1->FillRandom("f",1000); // fill histogram sampling fitFunc
71 fitFunc->SetParameter(0,3.); // vary a little the parameters
72 auto fitResult = h1->Fit(fitFunc,"S"); // fit the histogram and get fit result pointer
73
74 h1->Draw();
75
76 /* calculate the integral*/
77 double integral = fitFunc->Integral(0,1);
78
79 auto covMatrix = fitResult->GetCovarianceMatrix();
80 std::cout << "Covariance matrix from the fit ";
81 covMatrix.Print();
82
83 // need to pass covariance matrix to fit result.
84 // Parameters values are are stored inside the function but we can also retrieve from TFitResult
85 double sigma_integral = fitFunc->IntegralError(0,1, fitResult->GetParams() , covMatrix.GetMatrixArray());
86
87 std::cout << "Integral = " << integral << " +/- " << sigma_integral
88 << std::endl;
89
90 // estimated integral and error analytically
91
92 double * p = fitFunc->GetParameters();
93 double ic = p[1]* (1-std::cos(p[0]) )/p[0];
94 double c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0];
95 double c1c = (1-std::cos(p[0]) )/p[0];
96
97 // estimated error with correlations
98 double sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1)
99 + 2.* c0c*c1c * covMatrix(0,1));
100
101 if ( std::fabs(sigma_integral-sic) > 1.E-6*sic )
102 std::cout << " ERROR: test failed : different analytical integral : "
103 << ic << " +/- " << sic << std::endl;
104}
#define f(i)
Definition RSha256.hxx:104
winID h TVirtualViewer3D TVirtualGLPainter p
1-Dim function class
Definition TF1.h:214
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:526
virtual void SetParameters(const Double_t *params)
Definition TF1.h:650
virtual void SetParameter(Int_t param, Double_t value)
Definition TF1.h:640
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:620
virtual void FillRandom(const char *fname, Int_t ntimes=5000, TRandom *rng=nullptr)
Fill histogram following distribution in function fname.
Definition TH1.cxx:3520
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:3901
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3067
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