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

Processing /mnt/build/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/fit/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
Integral = 19.005 +/- 0.6159

#include "TF1.h"
#include "TH1D.h"
#include "TVirtualFitter.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
   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();

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

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

Author

Lorenzo Moneta

Definition in file ErrorIntegral.C.