Hello,
I am using the ROOT::Math::GSLMCIntegrator with kVEGAS to integrate a multi-dimensional function. The function is a polynomial with negative exponentials, defined as ff6 in the code. The integration of ff6 gives the matrix elements of row number 6. This is a 10x10 correction matrix. I would like to calculate the integrals with an absolute tolerance of at least 1e-4 for all elements of the matrix. However, not even for an absolute tolerance of 1e-1 the integration converges. ROOT keeps running it forever.
The same integration if performed using ROOT::Math::AdaptiveIntegratorMultiDim with absolute tolerance of 1e-3 takes about 4 hours to complete. When I calculate the elements in rows 8, 9 and 10, each element takes around 10 hours with an absolute tolerance of 1e-3. With an absolute tolerance of 1e-4 it would take days to run a few elements. I tried substituting the exponentials by a series expansion, but the series only converge after 20 terms which would complicate even more the function being integrated. I did this because I have the impression that the exponential terms are the ones responsible for slowing down the integration.
Do you have an idea of why the MC integration is not converging? Is there another method I could try?
I appreciate your help,
Kind regards,
Jeremias
double f(double x) {
TF1 *tf1 = new TF1("tf1","([0]**2*[1])/([1]-[0])**2*(exp(-[0]*x)*(([1]-[0])*x-1)+exp(-[1]*x))*[2]+[3]",0.0,100.0); //[0] is lambda and [1] stands for beta.
tf1->SetParameters(0.154,0.0722,100.0,0.0001,0.0);
Double_t top=tf1->Eval(tf1->GetMaximumX())+100.0;
TF1 *tf2 = new TF1("tf2","tf1-[3]",0.0,100.0);
TF1 *tf4 = new TF1("tf4","tf2/([0])",0.0,60.0);
tf4->SetParameter(0,tf2->Integral(0.0,60.0));
Double_t p0=tf1->GetParameter(0), p1=tf1->GetParameter(1), p2=tf1->GetParameter(2), p3=tf2->Integral(0.0,60.0);
delete tf1; delete tf2; delete tf4;
return ((pow(p0,2)*p1/pow(p1-p0,2))*(exp(-p0*x)*((p1-p0)*x-1.0)+exp(-p1*x))*p2)/p3;
}
//____________________________
double fint1(double t)
{
ROOT::Math::Functor1D wf(&f);
ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR);
ig.SetFunction(wf);
Double_t tau=0.9;
Double_t fintval = ig.Integral(t,t+tau);
return fintval;
}
//_______________________________i=6-------------------------------------------------------
double ff6(const double *x, const double *par)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t t2min=x[0]+tau;
Double_t t2max=T-(i-2.0)*tau;
Double_t t2=(t2max-t2min)*x[1]+t2min;
Double_t dt2=t2max-t2min;
Double_t t3min=t2+tau;
Double_t t3max=T-(i-3.0)*tau;
Double_t t3=(t3max-t3min)*x[2]+t3min;
Double_t dt3=t3max-t3min;
Double_t t4min=t3+tau;
Double_t t4max=T-(i-4.0)*tau;
Double_t t4=(t4max-t4min)*x[3]+t4min;
Double_t dt4=t4max-t4min;
Double_t t5min=t4+tau;
Double_t t5max=T-(i-5.0)*tau;
Double_t t5=(t5max-t5min)*x[4]+t5min;
Double_t dt5=t5max-t5min;
Double_t t6min=t5+tau;
Double_t t6max=T-(i-6.0)*tau;
Double_t t6=(t6max-t6min)*x[5]+t6min;
Double_t dt6=t6max-t6min;
Double_t ffval = (TMath::Factorial(par[0])/TMath::Factorial(par[0]-i)) *f(x[0]) *f(t2) *f(t3) *f(t4) *f(t5) *f(t6) *dt2 *dt3 *dt4 *dt5 *dt6 *pow(fint1(x[0]) +fint1(t2) +fint1(t3) +fint1(t4) +fint1(t5) +fint1(t6),par[0]-i);
return ffval;
}
//____________________________
double S6j(double j)
{
const int i=6;
Double_t tau=0.9;
Double_t T=60.0;
Double_t epsabs=1e-1;
ROOT::Math::WrappedParamFunction<> wf( &ff6, i, 1);
Double_t p[1] = {j};
wf.SetParameters(p);
// Create the Integrator
ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kVEGAS);
// Set parameters of the integration
ig.SetFunction(wf);
ig.SetAbsTolerance(epsabs);
Double_t xmin[i] = {0.0,0.0,0.0,0.0,0.0,0.0};
Double_t xmax[i] = {T-(i-1.0)*tau,1.0,1.0,1.0,1.0,1.0};
Double_t Sjvalue =ig.Integral(xmin, xmax);
return Sjvalue;
}
void spencerIntegrals() {
ROOT::EnableImplicitMT(0);
cout << "Matrix element S66: " << S6j(6) << endl;
cout << "Matrix element S67: " << S6j(7) << endl;
}
ROOT Version: v6-28-06
Platform: Red Hat Enterprise Linux 8.8; Memory 125.1 GiB; Intel Core i7-9800x CPU 3.8 GHz x16.
Compiler: C++ (GCC) 8.5.0 20210514 (Red Hat 8.5.0-18)