Logo ROOT   6.12/07
Reference Guide
TF1Helper.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Lorenzo Moneta 12/06/07
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7  * *
8  * *
9  **********************************************************************/
10 
11 // helper functions and classes used internally by TF1
12 
13 #include "TF1Helper.h"
14 #include "TError.h"
15 #include "TMath.h"
16 #include <vector>
17 #include <cmath>
18 #include <cassert>
19 
20 #include "TBackCompFitter.h"
21 #include "TVectorD.h"
22 #include "TMatrixD.h"
23 
24 #include "Math/IntegratorOptions.h"
25 
27 
28 namespace ROOT {
29 
30 
31 
32  namespace TF1Helper{
33 
34 
35 
36 
37 
38  double IntegralError(TF1 * func, Int_t ndim, const double * a, const double * b, const double * params, const double * covmat, double epsilon) {
39 
40  // calculate the eror on an integral from a to b of a parametetric function f when the parameters
41  // are estimated from a fit and have an error represented by the covariance matrix of the fit.
42  // The latest fit result is used
43 
44  // need to create the gradient functions w.r.t to the parameters
45 
46 
47  // loop on all parameters
48  bool onedim = ndim == 1;
49  int npar = func->GetNpar();
50  if (npar == 0) {
51  Error("TF1Helper","Function has no parameters");
52  return 0;
53  }
54 
55  std::vector<double> oldParams;
56  if (params) {
57  // when using an external set of parameters
58  oldParams.resize(npar);
59  std::copy(func->GetParameters(), func->GetParameters()+npar, oldParams.begin());
60  func->SetParameters(params);
61  }
62 
63 
64  TMatrixDSym covMatrix(npar);
65  if (covmat == 0) {
66  // use matrix from last fit (needs to be a TBackCompFitter)
68  TBackCompFitter * fitter = dynamic_cast<TBackCompFitter*> (vfitter);
69  if (fitter == 0) {
70  Error("TF1Helper::IntegralError","No existing fitter can be used for computing the integral error");
71  return 0;
72  }
73  // check that fitter and function are in sync
74  if (fitter->GetNumberTotalParameters() != npar) {
75  Error("TF1Helper::IntegralError","Last used fitter is not compatible with the current TF1");
76  return 0;
77  }
78  // check that errors are provided
79  if (int(fitter->GetFitResult().Errors().size()) != npar) {
80  Warning("TF1Helper::INtegralError","Last used fitter does no provide parameter errors and a covariance matrix");
81  return 0;
82  }
83 
84  // check also the parameter values
85  for (int i = 0; i < npar; ++i) {
86  if (fitter->GetParameter(i) != func->GetParameter(i) ) {
87  Error("TF1Helper::IntegralError","Last used Fitter has different parameter values");
88  return 0;
89  }
90  }
91 
92  // fill the covariance matrix
93  fitter->GetFitResult().GetCovarianceMatrix(covMatrix);
94  }
95  else {
96  covMatrix.Use(npar,covmat);
97  }
98 
99 
100 
101  // loop on the parameter and calculate the errors
102  TVectorD ig(npar);
103 
104  if (epsilon <= 0) epsilon = 1.e7*ROOT::Math::IntegratorMultiDimOptions::DefaultRelTolerance();
105 
106  double numError2 = 0;
107  for (int i=0; i < npar; ++i) {
108 
109  // use tolerance factor of 10 smaller than parameter errors
110  double epsrel = TMath::Min(0.1, epsilon / std::abs(func->GetParameter(i)) );
111  double epsabs = epsrel;
112 
113  // check that parameter error is not zero - otherwise skip it
114  // should check the limits
115  double integral = 0;
116  double error = 0;
117  if (covMatrix(i,i) > 0 ) {
118  TF1 gradFunc("gradFunc",TGradientParFunction(i,func),0,0,0);
119  if (onedim) {
120  integral = gradFunc.IntegralOneDim(*a,*b,epsrel,epsabs, error);
121  }
122  else {
123  double relerr;
124  int ifail = 0;
125  int nfnevl = 0;
127  if (maxpts == gDefMaxPts) maxpts = 10*gDefMaxPts; // increate points
128  integral = gradFunc.IntegralMultiple(ndim,a,b,maxpts,epsrel,epsabs,relerr,nfnevl,ifail);
129  //if (ifail) Warning("TF1Helper::IntegralError","n-dim integration failed code=%d I = %g, relerr =%g, ncall = %d, maxpts = %d, epsrel = %g, epsabs = %g, ",ifail,integral,relerr,nfnevl,maxpts,epsrel,epsabs);
130  error = relerr*std::abs(integral);
131  }
132  }
133  ig[i] = integral;
134  // std::cout << " ipar " << i << " sigma " << sqrt(covMatrix(i,i)) << " rel " << sqrt(covMatrix(i,i))/std::abs(func->GetParameter(i)) << " integral " << integral << " +/- " << error << " " <<
135  // error/std::abs(integral) << std::endl;
136 
137  // estimate numerical error (neglect correlations)
138  numError2 += covMatrix(i,i)*covMatrix(i,i) * integral * integral * error * error;
139  }
140  double err2 = covMatrix.Similarity(ig);
141  double result = sqrt(err2);
142  double numError = sqrt( numError2/sqrt(err2) );
143 
144  //std::cout << "integral error is " << result << " num error is " << numError << std::endl;
145  if ( result > epsilon && numError > 2*epsilon*result )
146  Warning("TF1Helper::IntegralError","numerical error from integration is too large. Integral error = %g +/- %g - eps = %g",result,numError,epsilon);
147 
148  // restore old parameters in TF1
149  if (!oldParams.empty()) {
150  func->SetParameters(&oldParams.front());
151  }
152 
153 
154 
155  return std::sqrt(err2);
156 
157 }
158 
159 
160 } // end namespace TF1Helper
161 
162 
163 } // end namespace ROOT
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const std::vector< double > & Errors() const
parameter errors (return st::vector)
Definition: FitResult.h:170
virtual Double_t IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Int_t maxpts, Double_t epsrel, Double_t epsabs, Double_t &relerr, Int_t &nfnevl, Int_t &ifail)
This function computes, to an attempted specified accuracy, the value of the integral.
Definition: TF1.cxx:2719
int gDefMaxPts
Definition: TF1Helper.cxx:26
Backward compatible implementation of TVirtualFitter.
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
double sqrt(double)
virtual Int_t GetNumberTotalParameters() const
Number of total parameters.
const ROOT::Fit::FitResult & GetFitResult() const
function class representing the derivative with respect a parameter of a given TF1 ...
Definition: TF1Helper.h:27
auto * a
Definition: textangle.C:12
double IntegralError(TF1 *func, Int_t ndim, const double *a, const double *b, const double *params, const double *covmat, double epsilon)
Definition: TF1Helper.cxx:38
static TVirtualFitter * GetFitter()
static: return the current Fitter
void Warning(const char *location, const char *msgfmt,...)
REAL epsilon
Definition: triangle.c:617
void GetCovarianceMatrix(Matrix &mat) const
fill covariance matrix elements using a generic matrix class implementing operator(i,j) the matrix must be previously allocates with right size (npar * npar)
Definition: FitResult.h:237
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
Abstract Base Class for Fitting.
virtual Int_t GetNpar() const
Definition: TF1.h:465
virtual Double_t IntegralOneDim(Double_t a, Double_t b, Double_t epsrel, Double_t epsabs, Double_t &err)
Return Integral of function between a and b using the given parameter values and relative and absolut...
Definition: TF1.cxx:2492
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:496
1-Dim function class
Definition: TF1.h:211
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual Double_t * GetParameters() const
Definition: TF1.h:504
void Error(ErrorHandler_t func, int code, const char *va_(fmt),...)
Write error message and call a handler, if required.
virtual Double_t GetParameter(Int_t ipar) const
Parameter value.