Logo ROOT   6.16/01
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
25
27
28namespace 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
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
#define b(i)
Definition: RSha256.hxx:100
int Int_t
Definition: RtypesCore.h:41
void Error(const char *location, const char *msgfmt,...)
void Warning(const char *location, const char *msgfmt,...)
int gDefMaxPts
Definition: TF1Helper.cxx:26
double sqrt(double)
const std::vector< double > & Errors() const
parameter errors (return st::vector)
Definition: FitResult.h:170
void GetCovarianceMatrix(Matrix &mat) const
fill covariance matrix elements using a generic matrix class implementing operator(i,...
Definition: FitResult.h:237
function class representing the derivative with respect a parameter of a given TF1
Definition: TF1Helper.h:27
Backward compatible implementation of TVirtualFitter.
const ROOT::Fit::FitResult & GetFitResult() const
virtual Int_t GetNumberTotalParameters() const
Number of total parameters.
virtual Double_t GetParameter(Int_t ipar) const
Parameter value.
1-Dim function class
Definition: TF1.h:211
virtual Int_t GetNpar() const
Definition: TF1.h:465
virtual Double_t * GetParameters() const
Definition: TF1.h:504
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:2815
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
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:2586
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:496
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
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.
static TVirtualFitter * GetFitter()
static: return the current Fitter
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
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617