ErrorIntegral.C: Estimate the error in the integral of a fitted function

// 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;
}