ROOT logo
// 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
// Author: Lorenzo Moneta

#include "TF1.h"
#include "TH1D.h"
#include "TVirtualFitter.h"
#include "TMath.h"
#include <assert.h>
#include <iostream>
#include <cmath>

//#define HAVE_OLD_ROOT_VERSION

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] ); 
}

// when TF1::IntegralError was not available 

#ifdef HAVE_OLD_ROOT_VERSION
//____________________________________________________________________
double df_dPar(double * x, double * p) { 
   // derivative of the function w.r..t parameters
   // use calculated derivatives from TF1::GradientPar

   double grad[NPAR]; 
   // p is used to specify for which parameter the derivative is computed 
   int ipar = int(p[0] ); 
   assert (ipar >=0 && ipar < NPAR );

   assert(fitFunc);
   fitFunc->GradientPar(x, grad);

   return grad[ipar]; 
}

//____________________________________________________________________
double IntegralError(int npar, double * c, double * errPar, 
   double * covMatrix = 0) {   
// calculate the error on the integral given the parameter error and 
// the integrals of the gradient functions c[] 

   double err2 = 0; 
   for (int i = 0; i < npar; ++i) { 
      if (covMatrix == 0) { // assume error are uncorrelated
         err2 += c[i] * c[i] * errPar[i] * errPar[i]; 
      } else {
         double s = 0; 
         for (int j = 0; j < npar; ++j) {
            s += covMatrix[i*npar + j] * c[j]; 
         }
         err2 += c[i] * s; 
      }
   }

   return TMath::Sqrt(err2);
}
#endif

//____________________________________________________________________
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
   h1->Fit(fitFunc);             // fit the histogram 

   h1->Draw();

   // calculate the integral 
   double integral = fitFunc->Integral(0,1);

   TVirtualFitter * fitter = TVirtualFitter::GetFitter();
   assert(fitter != 0);
   double * covMatrix = fitter->GetCovarianceMatrix(); 

#ifdef HAVE_OLD_ROOT_VERSION

   // calculate now the error (needs the derivatives of the function 
   // w..r.t the parameters)
   TF1 * deriv_par0 = new TF1("dfdp0",df_dPar,0,1,1);
   deriv_par0->SetParameter(0,0);

   TF1 * deriv_par1 = new TF1("dfdp1",df_dPar,0,1,1);
   deriv_par1->SetParameter(0,1.);

   double c[2]; 

   c[0] = deriv_par0->Integral(0,1); 
   c[1] = deriv_par1->Integral(0,1); 

   double * epar = fitFunc->GetParErrors();

   // without correlations
   double sigma_integral_0 = IntegralError(2,c,epar);



   // with correlations
   double sigma_integral = IntegralError(2,c,epar,covMatrix);

#else 
   
   // using new function in TF1 (from 12/6/2007) 
   double sigma_integral = fitFunc->IntegralError(0,1);

#endif

   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] + c1c*c1c * covMatrix[3] 
      + 2.* c0c*c1c * covMatrix[1]); 

   if ( std::fabs(sigma_integral-sic) > 1.E-6*sic ) 
      std::cout << " ERROR: test failed : different analytical  integral : " 
                << ic << " +/- " << sic << std::endl;
}
 ErrorIntegral.C:1
 ErrorIntegral.C:2
 ErrorIntegral.C:3
 ErrorIntegral.C:4
 ErrorIntegral.C:5
 ErrorIntegral.C:6
 ErrorIntegral.C:7
 ErrorIntegral.C:8
 ErrorIntegral.C:9
 ErrorIntegral.C:10
 ErrorIntegral.C:11
 ErrorIntegral.C:12
 ErrorIntegral.C:13
 ErrorIntegral.C:14
 ErrorIntegral.C:15
 ErrorIntegral.C:16
 ErrorIntegral.C:17
 ErrorIntegral.C:18
 ErrorIntegral.C:19
 ErrorIntegral.C:20
 ErrorIntegral.C:21
 ErrorIntegral.C:22
 ErrorIntegral.C:23
 ErrorIntegral.C:24
 ErrorIntegral.C:25
 ErrorIntegral.C:26
 ErrorIntegral.C:27
 ErrorIntegral.C:28
 ErrorIntegral.C:29
 ErrorIntegral.C:30
 ErrorIntegral.C:31
 ErrorIntegral.C:32
 ErrorIntegral.C:33
 ErrorIntegral.C:34
 ErrorIntegral.C:35
 ErrorIntegral.C:36
 ErrorIntegral.C:37
 ErrorIntegral.C:38
 ErrorIntegral.C:39
 ErrorIntegral.C:40
 ErrorIntegral.C:41
 ErrorIntegral.C:42
 ErrorIntegral.C:43
 ErrorIntegral.C:44
 ErrorIntegral.C:45
 ErrorIntegral.C:46
 ErrorIntegral.C:47
 ErrorIntegral.C:48
 ErrorIntegral.C:49
 ErrorIntegral.C:50
 ErrorIntegral.C:51
 ErrorIntegral.C:52
 ErrorIntegral.C:53
 ErrorIntegral.C:54
 ErrorIntegral.C:55
 ErrorIntegral.C:56
 ErrorIntegral.C:57
 ErrorIntegral.C:58
 ErrorIntegral.C:59
 ErrorIntegral.C:60
 ErrorIntegral.C:61
 ErrorIntegral.C:62
 ErrorIntegral.C:63
 ErrorIntegral.C:64
 ErrorIntegral.C:65
 ErrorIntegral.C:66
 ErrorIntegral.C:67
 ErrorIntegral.C:68
 ErrorIntegral.C:69
 ErrorIntegral.C:70
 ErrorIntegral.C:71
 ErrorIntegral.C:72
 ErrorIntegral.C:73
 ErrorIntegral.C:74
 ErrorIntegral.C:75
 ErrorIntegral.C:76
 ErrorIntegral.C:77
 ErrorIntegral.C:78
 ErrorIntegral.C:79
 ErrorIntegral.C:80
 ErrorIntegral.C:81
 ErrorIntegral.C:82
 ErrorIntegral.C:83
 ErrorIntegral.C:84
 ErrorIntegral.C:85
 ErrorIntegral.C:86
 ErrorIntegral.C:87
 ErrorIntegral.C:88
 ErrorIntegral.C:89
 ErrorIntegral.C:90
 ErrorIntegral.C:91
 ErrorIntegral.C:92
 ErrorIntegral.C:93
 ErrorIntegral.C:94
 ErrorIntegral.C:95
 ErrorIntegral.C:96
 ErrorIntegral.C:97
 ErrorIntegral.C:98
 ErrorIntegral.C:99
 ErrorIntegral.C:100
 ErrorIntegral.C:101
 ErrorIntegral.C:102
 ErrorIntegral.C:103
 ErrorIntegral.C:104
 ErrorIntegral.C:105
 ErrorIntegral.C:106
 ErrorIntegral.C:107
 ErrorIntegral.C:108
 ErrorIntegral.C:109
 ErrorIntegral.C:110
 ErrorIntegral.C:111
 ErrorIntegral.C:112
 ErrorIntegral.C:113
 ErrorIntegral.C:114
 ErrorIntegral.C:115
 ErrorIntegral.C:116
 ErrorIntegral.C:117
 ErrorIntegral.C:118
 ErrorIntegral.C:119
 ErrorIntegral.C:120
 ErrorIntegral.C:121
 ErrorIntegral.C:122
 ErrorIntegral.C:123
 ErrorIntegral.C:124
 ErrorIntegral.C:125
 ErrorIntegral.C:126
 ErrorIntegral.C:127
 ErrorIntegral.C:128
 ErrorIntegral.C:129
 ErrorIntegral.C:130
 ErrorIntegral.C:131
 ErrorIntegral.C:132
 ErrorIntegral.C:133
 ErrorIntegral.C:134
 ErrorIntegral.C:135
 ErrorIntegral.C:136
 ErrorIntegral.C:137
 ErrorIntegral.C:138
 ErrorIntegral.C:139
 ErrorIntegral.C:140
 ErrorIntegral.C:141
 ErrorIntegral.C:142
 ErrorIntegral.C:143