Hello,
These days, I have been working with multiple integrals using different routines, let’s say: ROOT::Math::GSLMCIntegrator, ROOT::Math::AdaptiveIntegratorMultiDim and the method TF1::IntegralMultiple. I noticed that the last one was not always returning the same error code than ROOT::Math::AdaptiveIntegratorMultiDim although this is actually the integrator that internally is used in that member function.
I looked at the source code at http://root.cern.ch/root/html/src/TF1.cxx and I found that the documentation for Double_t TF1::IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Int_t minpts, Int_t maxpts, Double_t eps, Double_t &relerr,Int_t &nfnevl, Int_t &ifail) tells that:
[quote]
Output parameters:
__ relerr : Contains, on exit, an estimation of the relative accuracy of the result.
__ nfnevl : number of function evaluations performed.
__ ifail :
_____ 0 Normal exit. . At least minpts and at most maxpts calls to the function were performed.
_____ 1 maxpts is too small for the specified accuracy eps.
_______ The result and relerr contain the values obtainable for the
_______ specified value of maxpts.
_____ 3 n<2 or n>15[/quote]
But the code is doing:
[code] ROOT::Math::WrappedMultiFunction<TF1&> wf1(*this, n);
ROOT::Math::AdaptiveIntegratorMultiDim aimd(wf1, eps, eps, maxpts);
aimd.SetMinPts(minpts);
double result = aimd.Integral(a,b);
relerr = aimd.RelError();
nfnevl = aimd.NEval();
ifail = 0;
return result;
[/code]
So ifail is always returning 0, what I think is not the expected behavior. Even more, the member function with the same name but different prototype, which was kept for backwards compatibility, has the code:
[code]Double_t TF1::IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Double_t eps, Double_t &relerr)
{
// See more general prototype below.
// This interface kept for back compatibility
Int_t nfnevl,ifail;
Int_t minpts = 2+2n(n+1)+1; //ie 7 for n=1
Int_t maxpts = 1000;
Double_t result = IntegralMultiple(n,a,b,minpts, maxpts,eps,relerr,nfnevl,ifail);
if (ifail > 0) {
Warning(“IntegralMultiple”,"failed code=%d, ",ifail);
}
return result;[/code]
So with the mentioned coding about ifail in the general prototype, the ‘if’ branch condition is never met here and the warning never issued.
Maybe I have missed something but if not, please let me suggest substituting ifail = 0; with ifail = aimd.Status(); in Double_t TF1::IntegralMultiple(Int_t n, const Double_t *a, const Double_t *b, Int_t minpts, Int_t maxpts, Double_t eps, Double_t &relerr,Int_t &nfnevl, Int_t &ifail) to solve both issues.
Many thanks,
Jose