Logo ROOT  
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#include "TROOT.h"
24
26
28
29namespace ROOT {
30
31
32
33 namespace TF1Helper{
34
35
36
37
38
39 double IntegralError(TF1 * func, Int_t ndim, const double * a, const double * b, const double * params, const double * covmat, double epsilon) {
40
41 // calculate the eror on an integral from a to b of a parametetric function f when the parameters
42 // are estimated from a fit and have an error represented by the covariance matrix of the fit.
43 // The latest fit result is used
44
45 // need to create the gradient functions w.r.t to the parameters
46
47
48 // loop on all parameters
49 bool onedim = ndim == 1;
50 int npar = func->GetNpar();
51 if (npar == 0) {
52 Error("TF1Helper","Function has no parameters");
53 return 0;
54 }
55
56 std::vector<double> oldParams;
57 if (params) {
58 // when using an external set of parameters
59 oldParams.resize(npar);
60 std::copy(func->GetParameters(), func->GetParameters()+npar, oldParams.begin());
61 func->SetParameters(params);
62 }
63
64
65 TMatrixDSym covMatrix(npar);
66 if (covmat == 0) {
67 // with ROOT implicit MT there is no global TVirtualFitter
69 Error("TF1Helper::IntegralError", "ROOT has enabled implicit MT. There is no existing lobal fitter, as shown in the documentation a pointer to the covariance matrix"
70 "from the TFitResult must be passed to TF1::IntegralError");
71 return 0;
72 }
73 // use matrix from last fit (needs to be a TBackCompFitter)
75 TBackCompFitter * fitter = dynamic_cast<TBackCompFitter*> (vfitter);
76 if (fitter == 0) {
77 Error("TF1Helper::IntegralError","No existing fitter can be used for computing the integral error");
78 return 0;
79 }
80 // check that fitter and function are in sync
81 if (fitter->GetNumberTotalParameters() != npar) {
82 Error("TF1Helper::IntegralError","Last used fitter is not compatible with the current TF1");
83 return 0;
84 }
85 // check that errors are provided
86 if (int(fitter->GetFitResult().Errors().size()) != npar) {
87 Warning("TF1Helper::INtegralError","Last used fitter does no provide parameter errors and a covariance matrix");
88 return 0;
89 }
90
91 // check also the parameter values
92 for (int i = 0; i < npar; ++i) {
93 if (fitter->GetParameter(i) != func->GetParameter(i) ) {
94 Error("TF1Helper::IntegralError","Last used Fitter has different parameter values");
95 return 0;
96 }
97 }
98
99 // fill the covariance matrix
100 fitter->GetFitResult().GetCovarianceMatrix(covMatrix);
101 }
102 else {
103 covMatrix.Use(npar,covmat);
104 }
105
106
107
108 // loop on the parameter and calculate the errors
109 TVectorD ig(npar);
110
112
113 double numError2 = 0;
114 for (int i=0; i < npar; ++i) {
115
116 // use tolerance factor of 10 smaller than parameter errors
117 double epsrel = TMath::Min(0.1, epsilon / std::abs(func->GetParameter(i)) );
118 double epsabs = epsrel;
119
120 // check that parameter error is not zero - otherwise skip it
121 // should check the limits
122 double integral = 0;
123 double error = 0;
124 if (covMatrix(i,i) > 0 ) {
125 TF1 gradFunc("gradFunc",TGradientParFunction(i,func),0,0,0);
126 if (onedim) {
127 integral = gradFunc.IntegralOneDim(*a,*b,epsrel,epsabs, error);
128 }
129 else {
130 double relerr;
131 int ifail = 0;
132 int nfnevl = 0;
134 if (maxpts == gDefMaxPts) maxpts = 10*gDefMaxPts; // increate points
135 integral = gradFunc.IntegralMultiple(ndim,a,b,maxpts,epsrel,epsabs,relerr,nfnevl,ifail);
136 //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);
137 error = relerr*std::abs(integral);
138 }
139 }
140 ig[i] = integral;
141 // std::cout << " ipar " << i << " sigma " << sqrt(covMatrix(i,i)) << " rel " << sqrt(covMatrix(i,i))/std::abs(func->GetParameter(i)) << " integral " << integral << " +/- " << error << " " <<
142 // error/std::abs(integral) << std::endl;
143
144 // estimate numerical error (neglect correlations)
145 numError2 += covMatrix(i,i)*covMatrix(i,i) * integral * integral * error * error;
146 }
147 double err2 = covMatrix.Similarity(ig);
148 double result = sqrt(err2);
149 double numError = sqrt( numError2/sqrt(err2) );
150
151 //std::cout << "integral error is " << result << " num error is " << numError << std::endl;
152 if ( result > epsilon && numError > 2*epsilon*result )
153 Warning("TF1Helper::IntegralError","numerical error from integration is too large. Integral error = %g +/- %g - eps = %g",result,numError,epsilon);
154
155 // restore old parameters in TF1
156 if (!oldParams.empty()) {
157 func->SetParameters(&oldParams.front());
158 }
159
160
161
162 return std::sqrt(err2);
163
164}
165
166
167} // end namespace TF1Helper
168
169
170} // end namespace ROOT
#define b(i)
Definition: RSha256.hxx:100
int Int_t
Definition: RtypesCore.h:45
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition: TError.cxx:231
int gDefMaxPts
Definition: TF1Helper.cxx:27
const std::vector< double > & Errors() const
parameter errors (return st::vector)
Definition: FitResult.h:169
void GetCovarianceMatrix(Matrix &mat) const
fill covariance matrix elements using a generic matrix class implementing operator(i,...
Definition: FitResult.h:236
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
Get reference to Fit Result object (NOTE: it will be invalid when class is deleted)
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:213
virtual Int_t GetNpar() const
Definition: TF1.h:481
virtual Double_t * GetParameters() const
Definition: TF1.h:520
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:2841
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:644
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:2611
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:512
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
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
double IntegralError(TF1 *func, Int_t ndim, const double *a, const double *b, const double *params, const double *covmat, double epsilon)
Definition: TF1Helper.cxx:39
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Bool_t IsImplicitMTEnabled()
Returns true if the implicit multi-threading in ROOT is enabled.
Definition: TROOT.cxx:558
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:176
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:618