ROOT   Reference Guide
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/// \macro_image
16/// \macro_output
17/// \macro_code
18///
19/// \author Lorenzo Moneta
20
21#include "TF1.h"
22#include "TH1D.h"
23#include "TFitResult.h"
24#include "TMath.h"
25#include <assert.h>
26#include <iostream>
27#include <cmath>
28
29TF1 * fitFunc; // fit function pointer
30
31const int NPAR = 2; // number of function parameters;
32
33//____________________________________________________________________
34double f(double * x, double * p) {
35 // function used to fit the data
36 return p[1]*TMath::Sin( p[0] * x[0] );
37}
38
39//____________________________________________________________________
40void ErrorIntegral() {
41 fitFunc = new TF1("f",f,0,1,NPAR);
42 TH1D * h1 = new TH1D("h1","h1",50,0,1);
43
44 double par[NPAR] = { 3.14, 1.};
45 fitFunc->SetParameters(par);
46
47 h1->FillRandom("f",1000); // fill histogram sampling fitFunc
48 fitFunc->SetParameter(0,3.); // vary a little the parameters
49 auto fitResult = h1->Fit(fitFunc,"S"); // fit the histogram and get fit result pointer
50
51 h1->Draw();
52
53 /* calculate the integral*/
54 double integral = fitFunc->Integral(0,1);
55
56 auto covMatrix = fitResult->GetCovarianceMatrix();
57 std::cout << "Covariance matrix from the fit ";
58 covMatrix.Print();
59
60 // need to pass covariance matrix to fit result.
61 // Parameters values are are stored inside the function but we can also retrieve from TFitResult
62 double sigma_integral = fitFunc->IntegralError(0,1, fitResult->GetParams() , covMatrix.GetMatrixArray());
63
64 std::cout << "Integral = " << integral << " +/- " << sigma_integral
65 << std::endl;
66
67 // estimated integral and error analytically
68
69 double * p = fitFunc->GetParameters();
70 double ic = p[1]* (1-std::cos(p[0]) )/p[0];
71 double c0c = p[1] * (std::cos(p[0]) + p[0]*std::sin(p[0]) -1.)/p[0]/p[0];
72 double c1c = (1-std::cos(p[0]) )/p[0];
73
74 // estimated error with correlations
75 double sic = std::sqrt( c0c*c0c * covMatrix(0,0) + c1c*c1c * covMatrix(1,1)
76 + 2.* c0c*c1c * covMatrix(0,1));
77
78 if ( std::fabs(sigma_integral-sic) > 1.E-6*sic )
79 std::cout << " ERROR: test failed : different analytical integral : "
80 << ic << " +/- " << sic << std::endl;
81}
#define f(i)
Definition RSha256.hxx:104
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
Double_t Sin(Double_t)
Definition TMath.h:639