diff --git a/resum/resint.C b/resum/resint.C
index bde9fe3..7905c54 100644
--- a/resum/resint.C
+++ b/resum/resint.C
@@ -1,955 +1,955 @@
#include "resint.h"
#include "phasespace.h"
#include "settings.h"
#include "dyres_interface.h"
#include "resconst.h"
#include "mesq.h"
#include "omegaintegr.h"
#include "gaussrules.h"
#include "hcoefficients.h"
#include "hcoeff.h"
#include "pdfevol.h"
#include "besselint.h"
#include "cubature.h"
#include "pdf.h"
#include "scales.h"
#include "sudakovff.h"
#include "besche_interface.h"
#include "hankel.h"
#include "minprescription_interface.h"
#include "specialfunctions_interface.h"
#include "588_interface.h"
#include <iostream>
#include "intde2_c.h"
#include <LHAPDF/LHAPDF.h>
double tiny = 1.0e-307;
double awinf[lenaw];
double awdei[lenaw];
//double awfin[lenaw];
double resint::_qt;
double resint::_m;
double resint::_y;
double resint::_costh;
double resint::tau;
double resint::x1;
double resint::x2;
int resint::_mode;
double resint::muren;
double resint::mufac;
double resint::mures;
double resint::muren2;
double resint::mufac2;
double resint::mures2;
double resint::a;
complex <double> resint::loga;
double resint::rloga;
complex <double> resint::logmuf2q2;
complex <double> resint::logq2muf2;
complex <double> resint::logq2mur2;
double resint::rlogq2mur2;
double resint::alpqr;
double resint::alpqf;
double resint::aass;
double resint::alpqfac;
double resint::alpqren;
double resint::alpqres;
double resint::bc;
//using namespace resint;
//#pragma omp threadprivate(_qt,_m,_y,_costh,_mode,muren,mufac,mures,muren2,mufac2,mures2,a,loga,rloga,logmuf2q2,logq2muf2,logq2mur2,rlogq2mur2,alpqr,alpqf,aass,alpqfac,alpqren,alpqres)
extern "C"
{
double besselint_besche_(double &x)
{
return real(besselint::bint(x));
}
double bint_mp_up_real_(double &x)
{
complex <double> jac(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr));
complex <double> b = resint::bc + jac*x;
fcomplex arg = fcx(resint::_qt*b);
return real(besselint::bint(b)*jac*cx(h1_(arg)));
}
double bint_mp_up_imag_(double &x)
{
complex <double> jac(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr));
complex <double> b = resint::bc + jac*x;
fcomplex arg = fcx(resint::_qt*b);
return imag(besselint::bint(b)*jac*cx(h1_(arg)));
}
double bint_mp_dn_real_(double &x)
{
complex <double> jac(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr));
complex <double> b = resint::bc + jac*x;
fcomplex arg = fcx(resint::_qt*b);
return real(besselint::bint(b)*jac*cx(h2_(arg)));
}
double bint_mp_dn_imag_(double &x)
{
complex <double> jac(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr));
complex <double> b = resint::bc + jac*x;
fcomplex arg = fcx(resint::_qt*b);
return imag(besselint::bint(b)*jac*cx(h2_(arg)));
}
// double besselint_mp_real_(double &x)
// {
// return besselint_mp_real_dequad(x);
// }
}
int besselint_cubature(unsigned ndim, const double x[], void *data, unsigned ncomp, double f[])
{
f[0] = real(besselint::bint(x[0]));
return 0;
}
double besselint_dequad(double x)
{
return real(besselint::bint(x));
}
double besselint_hankel(double x)
{
double qtb = x*resint::_qt;
return real(besselint::bint(x))/x/fort_besj0_(qtb);
}
double besselint_588(double &x)
{
double qtb = x*resint::_qt;
return real(besselint::bint(x))/fort_besj0_(qtb);
}
double besselint_mp_real_dequad(double x)
{
/*
complex <double> jacu(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr));
complex <double> bu = resint::bc + jacu*x;
fcomplex argu = fcx(resint::_qt*bu);
complex <double> jacd(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr));
complex <double> bd = resint::bc + jacd*x;
fcomplex argd = fcx(resint::_qt*bd);
complex <double> h1, h2;
// //call modified Hankel functions
// h1 = cx(h1_(argu));
// h2 = cx(h2_(argd));
// //cout << bu << " " << h1 << " " << bd << " " << h2 << endl;
//call original Hankel functions
int n = 1;
int nm;
fcomplex chf1[2], chd1[2], chf2[2], chd2[2];
ch12n_ (n, argu, nm, chf1, chd1, chf2, chd2);
h1 = cx(chf1[0]);
ch12n_ (n, argd, nm, chf1, chd1, chf2, chd2);
h2 = cx(chf2[0]);
//cout << endl;
//cout << " b " << bu << " " << bd << " bint " << besselint::bint(bu) << " " << besselint::bint(bd) << " h " << h1 << " " << h2 << " jac " << jacu << " " << jacd << endl;
//return (real(besselint::bint(bu)*jacu*h1) + real(besselint::bint(bd)*jacd*h2))/2.;
*/
//Real axis integration (phib = 0)
double bb = resint::bc + x;
double qtb = bb*resint::_qt;
complex <double> bint = besselint::bint(bb);
double J0 = fort_besj0_(qtb);
double Y0 = fort_besy0_(qtb);
double res = real(bint)*J0-imag(bint)*Y0;
//cout << " b " << bb << " " << " bint " << bint << " J0 " << J0 << " " << Y0 << endl;
//cout << " b " << bb << " h " << (real(besselint::bint(bu)*jacu*h1) + real(besselint::bint(bd)*jacd*h2))/2. << " j " << res << endl;
return res;
}
double besselint_mp_complex_dequad(double x)
{
complex <double> jacu(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr));
complex <double> bu = resint::bc + jacu*x;
fcomplex qtbu = fcx(resint::_qt*bu);
complex <double> jacd(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr));
complex <double> bd = resint::bc + jacd*x;
fcomplex qtbd = fcx(resint::_qt*bd);
complex <double> h1, h2;
// //call modified Hankel functions
// h1 = cx(h1_(qtbu));
// h2 = cx(h2_(qtbd));
// //cout << bu << " " << h1 << " " << bd << " " << h2 << endl;
//call original Hankel functions
int n = 1;
int nm;
fcomplex chf1[2], chd1[2], chf2[2], chd2[2];
ch12n_ (n, qtbu, nm, chf1, chd1, chf2, chd2);
h1 = cx(chf1[0]);
ch12n_ (n, qtbd, nm, chf1, chd1, chf2, chd2);
h2 = cx(chf2[0]);
//cout << endl;
//cout << " b " << bu << " " << bd << " bint " << besselint::bint(bu) << " " << besselint::bint(bd) << " h " << h1 << " " << h2 << " jac " << jacu << " " << jacd << endl;
//Compute Hankel functions from fortran bessel functions (does not work, because they accept only real values...)
//h1 = complex <double> (fort_besj0_(qtbu), fort_besy0_(qtbu));
//h2 = complex <double> (fort_besj0_(qtbd), -fort_besy0_(qtbd));
//h2 = conj(h1);
complex <double> bint_up = besselint::bint(bu);
complex <double> bint_dn = besselint::bint(bd);
//complex <double> bint_dn = conj(besselint::bint(bu));
return (real(bint_up*jacu*h1) + real(bint_dn*jacd*h2))/2.;
//return real(bint_up*jacu*h1);
}
double frealu_dequad(double x)
{
complex <double> jacu(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr));
complex <double> bu = resint::bc + jacu*x;
fcomplex argu = fcx(resint::_qt*bu);
complex <double> h1;
//call modified Hankel functions
/*
h1 = cx(h1_(argu));
//cout << bu << " " << h1 << endl;
*/
//call original Hankel functions
int n = 1;
int nm;
fcomplex chf1[2], chd1[2], chf2[2], chd2[2];
ch12n_ (n, argu, nm, chf1, chd1, chf2, chd2);
h1 = cx(chf1[0]);
return real(besselint::bint(bu)*jacu*h1);
}
double freald_dequad(double x)
{
complex <double> jacd(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr));
complex <double> bd = resint::bc + jacd*x;
fcomplex argd = fcx(resint::_qt*bd);
complex <double> h2;
//call modified Hankel functions
/*
h2 = cx(h2_(argd));
//cout << bd << " " << h2 << endl;
*/
//call original Hankel functions
int n = 1;
int nm;
fcomplex chf1[2], chd1[2], chf2[2], chd2[2];
ch12n_ (n, argd, nm, chf1, chd1, chf2, chd2);
h2 = cx(chf2[0]);
return real(besselint::bint(bd)*jacd*h2);
}
void resint::init()
{
//Set values in the common block for the sudakov
flag1_.flag1_ = opts.order; //order for the sudakov
iorder_.iord_ = opts.order - 1; //order for LL/NLL running of alphas
//modified_.imod_ = 1; // normal (imod=0) or modified (imod=1) sudakov
modified_.imod_ = opts.modlog; // canonical (imod=0) or modified (imod=1) logarithms in the Sudakov
// choose real axis (complex plane) integration of bstar (b) (always 0)
if (opts.bprescription == 2 || opts.bprescription == 3)
flagrealcomplex_.flagrealcomplex_ = 1;
else
flagrealcomplex_.flagrealcomplex_ = 0;
//v parameter of the modified Hankel functions
v_.v_=1.;
//g-parameter of the non perturbative form factor
np_.g_ = opts.g_param;
//a-parameter of the resummation scale, set it for the dynamic case
if (opts.fmures > 0)
{
a = 1./opts.kmures;
a_param_.a_param_ = a;
a_param_.b0p_ = resconst::b0*a;
// precompute log of scales
loga = log(a);
rloga = log(a);
rlogs_.rloga_ = rloga;
// clogs_.loga_ = loga; //complex loga is not used
}
// define points for quadratures integration
intdeoini(lenaw, tiny, opts.bintaccuracy, awinf);
//intdeoini(lenaw, tiny, 1.0e-6, awinf);
//intdeiini(lenaw, tiny, 1.0e-6, awdei);
//intdeini(lenaw, tiny, 1.0e-2, awfin);
}
double resint::rint(double costh, double m, double qt, double y, int mode)
{
/*
//limit the integration to qt = m
double jac = 1;
double qtp = qt;
if (qtp >= m*0.999)
return 0.;
else
qt = qtp/sqrt(1-pow(qtp/m,2));
jac = pow(m/sqrt(m*m-qtp*qtp),3);
*/
//point in phase space (should use the phasespace namespace instead of storing them in resint)
_qt = qt;
_y = y;
_m = m;
_costh = costh;
//compute x1 and x2 for the NP ff
double exppy = exp(y);
double expmy = 1./exppy;
tau = m/opts.sroot;
x1 = tau*exppy;
x2 = tau*expmy;
//mode 0: differential
//mode 1: integrated in costh
//mode 2: integrated in costh and y
//mode 3: integrated in costh, y, and pt
_mode = mode;
// Kinematical limit in rapidity
double q2 = m*m;
double x = q2/pow(opts.sroot,2);
double ax = log(x);
double ylim=-0.5*ax;
if (abs(y) > ylim)
return 0;
//move the scale setting to an independent namelist, so it can be called by all the calculations
scales::set(m);
muren = scales::ren;
mufac = scales::fac;
mures = scales::res;
//for fixed resummation scale need to recompute a_param
if (opts.fmures == 0)
{
a = m/mures;
loga = log(a);
rloga = log(a);
a_param_.a_param_ = a;
a_param_.b0p_ = resconst::b0*a;
rlogs_.rloga_ = rloga;
}
/*
//Set factorization and renormalization scales (set dynamic scale here)
if (opts.dynamicscale)
{
muren = m*opts.kmuren;
mufac = m*opts.kmufac;
}
else
{
muren = opts.rmass*opts.kmuren;
mufac = opts.rmass*opts.kmufac;
}
if (opts.dynamicresscale)
mures = m*opts.kmures;
else
//for fixed resummation scale need to recompute a_param
{
mures = opts.rmass*opts.kmures;
a = m/mures;
loga = log(a);
rloga = log(a);
a_param_.a_param_ = a;
a_param_.b0p_ = resconst::b0*a;
rlogs_.rloga_ = rloga;
}
*/
muren2=pow(muren,2);
mufac2=pow(mufac,2);
mures2=pow(mures,2);
//update PDFs in Mellin space at the starting scale, if the factorisation scale is proportional to mll
if (opts.fmufac > 0)
{
pdfevol::allocate();
pdfevol::update();
}
//alphas at various scales (alphas convention is alphas/4/pi)
alpqren=LHAPDF::alphasPDF(muren)/4./M_PI;
alpqfac=LHAPDF::alphasPDF(mufac)/4./M_PI;
alpqres=LHAPDF::alphasPDF(mures)/4./M_PI;
//alpqr is alphas at the renormalization scale
if (opts_.approxpdf_ == 1)
{
int nloop = 3;
double scale = muren;
alpqr=dyalphas_mcfm_(scale,couple_.amz_,nloop)/4./M_PI;
}
else
alpqr=alpqren;
//alpqf is alphas at the resummation scale. It is used as starting scale for the PDF evolution
if (opts_.approxpdf_ == 1)
{
int nloop = 3;
double scale = mures;
alpqf=dyalphas_mcfm_(scale,couple_.amz_,nloop)/4./M_PI;
}
else
alpqf=alpqres;
/***************************************************/
//set values in the common blocks for the sudakov
//aaas is alphas/pi at renormalization scale
aass = alpqr*4.;
aass_.aass_ = aass;
double q = m;
double blim;
//dynamic blim
if (opts.blim < 0)
{
//if (opts.order == 2 && q2 > 2.2*muren2) //avoid large values of b in f2(y) when q2>mur2
if (q2 > 2.2*muren2) //avoid large values of b in f2(y) when q2>mur2
blim = a_param_.b0p_*(1./q)*exp(1./(2.*aass*resconst::beta0))*sqrt(muren2/q2);
else
blim = a_param_.b0p_*(1./q)*exp(1./(2.*aass*resconst::beta0)); // avoid Landau pole
/*
double lambdaqcd = muren/(exp(1./(2.*aass*resconst::beta0))); //--> corresponds to a divergence in alphas
double lambdasudakov = mures/(exp(1./(2.*aass*resconst::beta0))); //--> correspond to a divergence in the Sudakov
//should blim depend or not on a_param? In principle yes because PDFs are evolved to b0/bstar*a
blim = min(resconst::b0/lambdaqcd,resconst::b0/lambdasudakov);
//cout << _m << " " << resconst::b0/lambdaqcd << " " << resconst::b0/lambdasudakov << " " << a_param_.b0p_*(1./q)*exp(1./(2.*aass*resconst::beta0)) << endl;
*/
blim = blim/(opts.blim*(-1));
}
//fixed blim
if (opts.blim > 0)
blim = opts.blim;
// blim = 100;
// blim = 4.0;
// blim = 1.5;
blimit_.rblim_ = blim;
blimit_.cblim_ = fcx(blim);
// cout << q << " " << mur << " " << blimit_.rblim_ << " " << a_param_.b0p_/lambdaqcd << " " << lambdaqcd <<" " << a_param_.b0p_/blim << endl;
//fill in common block for scales (used in the sudakov)
scaleh_.qt_ = _qt;
scaleh_.qt2_ = _qt*_qt;
scaleh_.xtau_ = 0.;
scaleh_.q_ = _m;
scaleh_.q2_ = _m*_m;
scaleh_.shad_ = pow(opts.sroot,2);
scaleh_.sroot_ = opts.sroot;
scaleh_.mur_ = muren;
scaleh_.mur2_ = muren2;
scaleh_.muf_ = mufac;
scaleh_.muf2_ = mufac2;
/***************************************************/
//******************************************
// Initialise the born-level matrix elements mesqij
// ---> mesqij depends on mass and costh
// when costh is integrated and costh moments are provided, mesqij_expy
// are convoluted with the rapidity depent expression
// and they depend also on rapidity
if (!(opts.order == 0 && opts.xspace))
{
mesq::allocate();
mesq::setmesq_expy(mode, m, costh, y);
}
else //No Mellin transform case
{
mesq::setpropagators(m);
if (mode == 0) //costh differential
mesq::setmesq(1., costh, pow(costh,2));
else //costh integrated
{
double cthmom0, cthmom1, cthmom2;
omegaintegr::cthmoments(cthmom0, cthmom1, cthmom2);
mesq::setmesq(cthmom0, cthmom1, cthmom2);
}
}
//******************************************
//*****************************************
//precompute scales: there is mass dependence in logq2muf2 logq2mur2, only with fixed factorization and renormalization scales
logmuf2q2=log(mufac2/q2);
logq2muf2=log(q2/mufac2);
logq2mur2=log(q2/muren2);
rlogq2mur2=log(q2/muren2);
rlogs_.rlogq2mur2_ = real(logq2mur2);
//no need to recompute coefficients with variable scales
//if (!dynamicscale || fixedresummationscale)
if (opts.mellin1d)
{
hcoeff::allocate();
hcoeff::calc(aass,logmuf2q2,logq2muf2,logq2mur2,loga);
}
else
{
hcoefficients::allocate();
hcoefficients::calc(aass,logmuf2q2,logq2muf2,logq2mur2,loga);
}
pdfevol::allocate_fx();
//*****************************************
double res;
if (mode == 3)
{
/*
Since qt appears only in the Bessel function J_0(bqt),
the integration in qt could be factorised by using
the following result: int_a^b J0(x) dx = Lambda0(b) - Lambda0(a)
where Lambda0(x) = xJ0(x) + pi x/2 [J1(x)H0(x) - J0(x)H1(x)], with Hn Struve functions.
!!! Things are actually much simpler because there is a qt factor (res *= qt/2./pow(opts.sroot,2);) !!!
int_a^b x J0(x) dx = b J1(b) - a J1(a)
See http://fisica.ciens.ucv.ve/~svincenz/TISPISGIMR.pdf and http://web.eah-jena.de/~rsh/Forschung/Stoer/besint.pdf (1.1.1)
!!! The problem with this is the switching function, which depends on qt !!!
--> The approach would be valid only for qt < m*k
Now the integration call will automatically switch to mode = 2 for bins where qtmax > mmin*k
*/
//double qtmn = max(opts.qtcutoff,phasespace::qtmin);
double qtmn = phasespace::qtmin;
//double kinqtlim = sqrt(max(0.,pow(pow(opts.sroot,2)+phasespace::m2,2)/(4*pow(opts.sroot,2))-phasespace::m2)); //introduced max to avoid neqative argument of sqrt
double kinqtlim = 1e10;
- double switchqtlim = phasespace::m*opts.dampk; //here use the value of qt where the switching start
+ double switchqtlim = opts.damp?phasespace::m*opts.dampk:1e10; //here use the value of qt where the switching start
double qtlim = min(kinqtlim, switchqtlim);
double qtmx = min(qtlim, phasespace::qtmax);
if (qtmn >= qtmx || qtmx < opts.qtcutoff)
res = 0;
else
{
double res2, err2;
_qt = qtmx;
intdeo(besselint_dequad, 0.0, qtmx, awinf, &res2, &err2);
double res1, err1;
_qt = qtmn;
if (qtmn < opts.qtcutoff)
res1 = 0.;
else
intdeo(besselint_dequad, 0.0, qtmn, awinf, &res1, &err1);
res = qtmx*res2 - qtmn*res1;
// Normalization
res *= 1./2./pow(opts.sroot,2);
//cout << " res " << res << endl;
}
}
else
{
// //if (opts.fast && opts.xspace) {
// //Fast function
// res = exp(-qt);
//
// // Normalization
// double shad = pow(opts.sroot,2);
// res *= qt/shad;
// //}
// //else
// //{
//*****************************************
//dependence on qt, m, also y, and costh unless integrated
//perform b integration (int_0^inf db)
res = bintegral(qt);
// cout << "phase space point in resumm" << " " << _qt << " " << _y << " " << _m << " " << _costh << " " << res << endl;
//*****************************************
// Normalization
res *= qt/2./pow(opts.sroot,2);
}
//Do not need the Mellin transform for the LL case, because the HN coefficient is 1 --> Use PDFs in x space
if (opts.order == 0 && opts.xspace && (opts.evolmode == 0 || opts.evolmode == 1))
{
double fun = 0.;
if (mode < 2) //rapidity differential
{
//PDFs
double fx1[2*MAXNF+1],fx2[2*MAXNF+1];
fdist_(opts.ih1,x1,mufac,fx1);
fdist_(opts.ih2,x2,mufac,fx2);
for (int sp = 0; sp < mesq::totpch; sp++)
fun += fx1[mesq::pid1[sp]]*fx2[mesq::pid2[sp]]*real(mesq::mesqij[sp]);
}
else //rapidity integration
{
for (int i=0; i < opts.yintervals; i++)
{
double ya = phasespace::ymin+(phasespace::ymax-phasespace::ymin)*i/opts.yintervals;
double yb = phasespace::ymin+(phasespace::ymax-phasespace::ymin)*(i+1)/opts.yintervals;
double xc = 0.5*(ya+yb);
double xm = 0.5*(yb-ya);
for (int j=0; j < opts.yrule; j++)
{
double y = xc+xm*gr::xxx[opts.yrule-1][j];
double exppy = exp(y);
double expmy = 1./exppy;
x1 = tau*exppy;
x2 = tau*expmy;
//PDFs
double fx1[2*MAXNF+1],fx2[2*MAXNF+1];
fdist_(opts.ih1,x1,mufac,fx1);
fdist_(opts.ih2,x2,mufac,fx2);
for (int sp = 0; sp < mesq::totpch; sp++)
fun += fx1[mesq::pid1[sp]]*fx2[mesq::pid2[sp]]*real(mesq::mesqij[sp])*gr::www[opts.yrule-1][j]*xm;;
}
}
}
res *= fun;
}
//free allocated Local Thread Storage (LTS) memory
if (opts.mellin1d)
hcoeff::free();
else
hcoefficients::free();
if (!(opts.order == 0 && opts.xspace))
mesq::free();
if (opts.fmufac > 0)
pdfevol::free();
pdfevol::free_fx();
//res *= jac;//jacobian for the change of variable qt=qtp/sqrt(1-qtp^2/m^2)
return res;
}
double resint::bintegral(double qt)
{
double res, err;
if (!opts.vfnsudakov)
{
//Warning!!! even if the sudakov is not VFN, when evolmode = 4, the integrand of bessel transform is discontinous at mb and mc
//--> Need to split the integral in this case
// Compute integral using double exponential quadratures for oscillatory fuctions (intde2.c)
sudakov::setnf(5);
//bstar prescription
if (opts.bprescription == 0)
{
//intdeo(besselint::bint, 0.0, qt, awinf, &res, &err);
intdeo(besselint_dequad, 0.0, qt, awinf, &res, &err);
if (err < 0)
cout << "warning: dequad abnormal termination, pt=" << _qt << " m=" << _m << " y=" << _y << " bint: " << res << " err " << err << endl;
/*
//test alternative integrations
cout << endl;
cout << "dequad result of inverse bessel transform, pt=" << _qt << " m=" << _m << " y=" << _y << " : " << setprecision(16) << res << " +- " << err << endl;
hankel::init(0,0,0.1);
hankel::transform(besselint_hankel, qt, res, err);
hankel::free();
cout << "hankel result " << res << " +- " << err << endl;
cout.precision(6); cout.unsetf(ios_base::floatfield);
//input
int nb = 1; //lagged convolutions
int nrel = 1; //related convolutions
int ntol = 1;
int nord[1] = {0}; //order of the transform
int ijrel[2*1]; //exponents of the related convolutions (not used)
double dwork [1*801]; //work array
//output
double dans[1*1]; //result
double arg[1]; //arguments of the lagged convolutions
int nofun1; //number of function evalutions
int ierr; //error code
dhankl_(qt, nb, nrel, opts.bintaccuracy, ntol, nord, besselint_588, ijrel, dwork, dans, arg, nofun1, ierr);
res = dans[0];
err = 0.;
cout << "588 result " << res << " +- " << err << " " << nofun1 << endl;
cout.precision(6); cout.unsetf(ios_base::floatfield);
*/
}
//Integrate up to blim with besche
else if (opts.bprescription == 1)
{
//cout << endl;
//cout << "pt " << _qt << " m " << _m << endl;
double C = blimit_.rblim_; //upper limit of integration
double ALFA[1] = {_qt}; //oscillation frequency
int NUM = 1; //number of integrals to be computed
int NU = 0; //order of the bessel function
int N = 30; //degree of the chebyshev approximation
double RESULT[1]; //output result
int INFO[1]; //output status code
besche_(besselint_besche_,C,ALFA,NUM,NU,N,RESULT,INFO);
res = RESULT[0];
err = 0;
//cout << "besche result " << res << " status " << INFO[0] << endl;
}
//Minimal prescription to avoid Landau singularity (see hep-ph/9604351 and hep-ph/0002078)
else if (opts.bprescription == 2)
{
//cout << endl;
//cout << "pt " << _qt << " m " << _m << endl;
bc = opts.bcf*blimit_.rblim_;
//besche from 0 to bc
double C = bc; //upper limit of integration
double ALFA[1] = {_qt}; //oscillation frequency
int NUM = 1; //number of integrals to be computed
int NU = 0; //order of the bessel function
int N = 30; //degree of the chebyshev approximation
double RESULT[1]; //output result
int INFO[1]; //output status code
besche_(besselint_besche_,C,ALFA,NUM,NU,N,RESULT,INFO);
res = RESULT[0];
err = 0;
//cout << "besche result up to bc " << res << " status " << INFO[0] << endl;
/*
double min = 0.;
double max = bc+blimit_.rblim_/_qt;
//adzint settings
double abserr = 1e-10;
double relerr = 1e-6;
double errest;
double ier;
double iacta = 0;
double iactb = 0;
//res = 0.;
complex <double> r = 0.;
complex <double> i(0.,1.);
//upper real integral
//r += adzint_(fureal_,min,max,abserr,relerr,errest,ier,iacta,iactb);
//lower real integral
//r += adzint_(fdreal_,min,max,abserr,relerr,errest,ier,iacta,iactb);
r += adzint_(freal_,min,max,abserr,relerr,errest,ier,iacta,iactb);
//res += real(r)/2.;
cout << "adzint result " << real(r)/2. << " error " << errest << endl;
*/
double resdequad;
double errdequad;
/*
complex <double> r = 0.;
intdeo(frealu_dequad, 0, qt*cos(M_PI/opts.phibr), awinf, &resdequad, &errdequad);
r += resdequad;
cout << "dequad up result " << resdequad << " error " << errdequad << endl;
intdeo(freald_dequad, 0, qt*cos(M_PI/opts.phibr), awinf, &resdequad, &errdequad);
r += resdequad;
cout << "dequad dn result " << resdequad << " error " << errdequad << endl;
res += real(r)/2.;
cout << "Split dequad result " << real(r)/2. << " error " << errdequad << endl;
*/
intdeo(besselint_mp_complex_dequad, 0, qt*cos(M_PI/opts.phibr), awinf, &resdequad, &errdequad);
//cout << "Complex dequad result " << resdequad << " error " << errdequad << endl;
res += resdequad;
//cout << "Total bp = 2 " << res << endl;
}
//Minimal prescription along the real axis, crossing the Landau singularity
else if (opts.bprescription == 3)
{
//minimal prescription with real b integration (not always more efficient, but save the calculation of PDFs at complex scales)
bc = opts.bcf*blimit_.rblim_;
//besche from 0 to bc
double C = bc; //upper limit of integration
double ALFA[1] = {_qt}; //oscillation frequency
int NUM = 1; //number of integrals to be computed
int NU = 0; //order of the bessel function
int N = 30; //degree of the chebyshev approximation
double RESULT[1]; //output result
int INFO[1]; //output status code
besche_(besselint_besche_,C,ALFA,NUM,NU,N,RESULT,INFO);
res = RESULT[0];
err = 0;
double resdequad;
double errdequad;
intdeo(besselint_mp_real_dequad, 0, qt, awinf, &resdequad, &errdequad);
cout << "Real dequad result " << resdequad << " error " << errdequad << endl;
res += resdequad;
cout << "Total bp = 3 " << res << endl;
}
}
else
{
//split the integral above and below the b and c masses
//scale b0/b without a_param
double bstar_mb = resconst::b0/(LHAPDF::getThreshold(5)*opts.kmub);
//convert bstar to b
double b_mb = bstar_mb / sqrt(1-(bstar_mb*bstar_mb)/(blimit_.rblim_*blimit_.rblim_));
//same for charm mass
double bstar_mc = resconst::b0/(LHAPDF::getThreshold(4)*opts.kmuc);
double b_mc = bstar_mc / sqrt(1-(bstar_mc*bstar_mc)/(blimit_.rblim_*blimit_.rblim_));
//Integrate from pt=0 to pt=mc (from b_mc to infinity)
//Set NF = 3
sudakov::setnf(3);
double res1, err1;
//intdeo(besselint::bint, b_mc, qt, awinf, &res1, &err1);
intdeo(besselint_dequad, b_mc, qt, awinf, &res1, &err1);
//Integrate from pt=mb to pt=infinity (from 0 to b_mb)
//Set NF = 5
sudakov::setnf(5);
double res2, err2;
//dequad seems not to be appropriate for this integral
// intde(besselint::bint, 0.0, b_mb, awfin, &res2, &err2);
//try simple pcubature
const int ndim = 1; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata = NULL;
double integral[1];
double error[1];
const int eval = 0;
const double epsrel = 1e-2;//opts.bintaccuracy;
const double epsabs = 0.;
double xmin[1] = {0.000001};
double xmax[1] = {b_mb};
/*
pcubature(ncomp, besselint_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
res2 = integral[0];
err2 = error[0];
cout << endl;
cout << "pcubature result " << res2 << " error " << err2 << endl;
*/
double C = b_mb; //upper limit of integration
double ALFA[1] = {_qt}; //oscillation frequency
int NUM = 1; //number of integrals to be computed
int NU = 0; //order of the bessel function
int N = 30; //degree of the chebyshev approximation
double RESULT[1]; //output result
int INFO[1]; //output status code
besche_(besselint_besche_,C,ALFA,NUM,NU,N,RESULT,INFO);
res2 = RESULT[0];
err2 = 0;
//cout << "besche result " << res2 << " status " << INFO[0] << endl;
//Integrate from pt=mc to pt=mb (from b_mb to b_mc)
//Set NF = 4
sudakov::setnf(4);
double res3, err3;
xmin[0] = {b_mb};
xmax[0] = {b_mc};
pcubature(ncomp, besselint_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
res3 = integral[0];
err3 = error[0];
res = res1+res2+res3;
/*
cout << "split: dequad result of inverse bessel transform, " << endl;
cout << "total: "
//<< setprecision(16)
<< res << " " << (err1+err2+err3)/res*100 << endl;
cout << " b(mc)-inf: " << res1 << " " << err1/res1*100 << "%"
<< " 0-b(mb): " << res2 << " " << err2/res2*100 << "%"
<< " b(mb)-b(mc): " << res3 << " " << err3/res3*100 << "%"
<< endl;
//cout.precision(6); cout.unsetf(ios_base::floatfield);
//cout << "b_mb " << b_mb << " b_mc " << b_mc << " blim " << blimit_.rblim_ << endl;
cout << endl;
*/
}
/*
//plot the b-integrand
cout << "{" << endl;
cout << "TGraph *g = new TGraph();" << endl;
for (int i = 0; i < 1000; i++)
{
// double b = b_mb*0.95+i*(b_mb*1.05 -b_mb*0.95)/100.;
double b = 0.+i*(b_mb*5 -0)/1000.;
// cout << b << " " << b_mb << " " << besselint::bint(b) << endl;;
cout << "g->SetPoint(g->GetN(), " << b << ", " << besselint::bint(b) << ");" << endl;
}
cout << "g->Draw();" << endl;
cout << "}" << endl;
*/
return res;
}
diff --git a/src/cubacall.C b/src/cubacall.C
index 7550a59..4f97cc1 100644
--- a/src/cubacall.C
+++ b/src/cubacall.C
@@ -1,994 +1,994 @@
#include "cubacall.h"
#include "settings.h"
#include "resintegr.h"
#include "ctintegr.h"
#include "finintegr.h"
#include "bornintegr.h"
#include "cubature.h"
#include "smolpack.h"
#include "phasespace.h"
#include "HistoHandler.h"
#include <cuba.h>
#include <iostream>
//flags for cuba Vegas integration:
//flags += 0 or 4; //collect only weights from final iteration (4) or from all iterations (0)
//flags = 8; //smoothing of importance sampling (0) or not (8)
/***************************************************************/
//resummation
void resintegr1d(vector <double> &res, double &err)
{
//Force qt differential mode when crossing the value of qt where the switching start
- if (phasespace::qtmax > phasespace::mmin*opts.dampk)
+ if (opts.damp && (phasespace::qtmax > phasespace::mmin*opts.dampk))
{
resintegr2d(res, err);
return;
}
const int ndim = 1; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata = NULL;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 65+2*65*opts.niterBORN;
const int maxeval = 65+2*65*opts.niterBORN;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) resintegrand2d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[1] = {0};
double xmax[1] = {1};
if (opts.cubacores == 0)
pcubature(ncomp, resintegrand1d_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, resintegrand1d_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
}
res.clear();
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
return;
}
void resintegr2d(vector <double> &res, double &err)
{
const int ndim = 2; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata = NULL;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 65+2*65*opts.niterBORN;
const int maxeval = 65+2*65*opts.niterBORN;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) resintegrand2d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[2] = {0, 0};
double xmax[2] = {1, 1};
if (opts.cubacores == 0)
pcubature(ncomp, resintegrand2d_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, resintegrand2d_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
}
res.clear();
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
return;
}
void resintegr3d(vector <double> &res, double &err)
{
const int ndim = 3; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 127+2*127*opts.niterBORN;
const int maxeval = 127+2*127*opts.niterBORN;
const int key = 13;
int nregions;
Cuhre(ndim, ncomp,
(integrand_t) resintegrand3d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
res.clear();
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
return;
}
//original integration, can use this to sample phase space and fill histograms
void resintegrMC(vector <double> &res, double &err)
{
const int ndim = 6; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = (4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsBORN;
const int maxeval = opts.vegasncallsBORN;
const int nstart = max(10, int(opts.vegasncallsBORN/10));
const int nincrease = max(10, int(opts.vegasncallsBORN/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)resintegrandMC, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
return;
}
/***************************************************************/
/***************************************************************/
// fixed order
//Cuba integration of analytical V+j
void vjintegr3d(vector <double> &res, double &err)
{
const int ndim = 3; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 127+2*127*opts.niterVJ;
const int maxeval = 127+2*127*opts.niterVJ;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) vjintegrand, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[3] = {0, 0., 0.};
double xmax[3] = {1, 1., 1.};
if (opts.cubacores == 0)
pcubature(ncomp, vjintegrand_cubature, userdata,
//hcubature(ncomp, vjintegrand_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, vjintegrand_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
/*
//smolyak
int print_stats = 0;
int dim = ndim;
int l = 10;
integral[0] = int_smolyak (ndim, ndim+l, vjintegrand_smolyak, print_stats );
error[0] = 0.000000001;
*/
}
res.clear();
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
return;
}
//Cuba integration of V+j LO
void vjlointegr(vector <double> &res, double &err)
{
const int ndim = 7; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsVJLO;
const int maxeval = opts.vegasncallsVJLO;
const int nstart = max(10, int(opts.vegasncallsVJLO/10));
const int nincrease = max(10, int(opts.vegasncallsVJLO/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)lowintegrand, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
void vjlointegr7d(vector <double> &res, double &err)
{
const int ndim = 7; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsVJLO;
const int maxeval = opts.vegasncallsVJLO;
const int nstart = max(10, int(opts.vegasncallsVJLO/10));
const int nincrease = max(10, int(opts.vegasncallsVJLO/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)vjlointegrandMC, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
void vjlointegr5d(vector <double> &res, double &err)
{
const int ndim = 4;//5; //dimensions of the integral
const int ncomp = (opts.helicity >= 0 ? 2 : 1); //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 153+2*153*opts.niterVJ;
const int maxeval = 153+2*153*opts.niterVJ;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) vjlointegrand, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double tiny = 0.;//1e-6;
/*
double xmin[5] = {0., 0., tiny, 0., 0.};
double xmax[5] = {1., 1., 1.-tiny, 1., 1.};
//if (opts.cubacores == 0)
//pcubature(ncomp, vjlointegrand_cubature, userdata, //--> pcubature has an issue when phi is symmetric, it is resonant for nested quadrature rules, and pcubature misses the substructure of the phi distribution
hcubature(ncomp, vjlointegrand_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
*/
// /*
//4d integration (phi_lep integrated inside) works better in full phase space, and to calculate moments
double xmin[4] = {0., 0., tiny, 0.};
double xmax[4] = {1., 1., 1.-tiny, 1.};
if (opts.cubacores == 0)
//pcubature(ncomp, vjlointegrand_cubature, userdata, //--> pcubature has an issue when phi is symmetric, it is resonant for nested quadrature rules, and pcubature misses the substructure of the phi distribution
hcubature(ncomp, vjlointegrand_cubature, userdata,
4, xmin, xmax,
eval, epsabs, epsrel, ERROR_LINF, integral, error);
else
//pcubature_v(ncomp, vjlointegrand_cubature_v, userdata,
hcubature_v(ncomp, vjlointegrand_cubature_v, userdata,
4, xmin, xmax,
eval, epsabs, epsrel, ERROR_LINF, integral, error);
// */
/*
//smolyak
int print_stats = 0;
int dim = 4;
int l = 20;
integral[0] = int_smolyak (ndim, ndim+l, vjlointegrand_smolyak, print_stats );
error[0] = 0.000001;
*/
}
res.clear();
if (opts.helicity >= 0)
res.push_back(integral[0] != 0? integral[1]/integral[0] : 0);
else
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
if (opts.helicity >= 0)
err = integral[0] != 0? error[0]/integral[0]*integral[1]/integral[0] : 0; // error[1]/integral[1]*integral[1]/integral[0]
else
err = error[0];
return;
}
//Cuba integration of the V+J NLO real part
void vjrealintegr(vector <double> &res, double &err)
{
const int ndim = 10; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const long long int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
long long int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const long long int mineval = opts.vegasncallsVJREAL;
const long long int maxeval = opts.vegasncallsVJREAL;
const long long int nstart = max(10, int(opts.vegasncallsVJREAL/10));
const long long int nincrease = max(10, int(opts.vegasncallsVJREAL/10));
const long long int nbatch = opts.cubanbatch;
const int gridno = 0;
llVegas(ndim, ncomp, (integrand_t)realintegrand, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
/*
//smolyak
int print_stats = 0;
int dim = 10;
int l = 12;
integral[0] = int_smolyak (ndim, ndim+l, realintegrand_smolyak, print_stats );
error[0] = 0.000001;
*/
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
//Cuba integration of the V+J NLO virtual part
void vjvirtintegr(vector <double> &res, double &err)
{
const int ndim = 8; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsVJVIRT;
const int maxeval = opts.vegasncallsVJVIRT;
const int nstart = max(10, int(opts.vegasncallsVJVIRT/10));
const int nincrease = max(10, int(opts.vegasncallsVJVIRT/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)virtintegrand, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
//Cuba integration of double virtual born configuration
void bornintegrMC6d(vector <double> &res, double &err)
{
const int ndim = 6; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsBORN;
const int maxeval = opts.vegasncallsBORN;
const int nstart = max(10, int(opts.vegasncallsBORN/10));
const int nincrease = max(10, int(opts.vegasncallsBORN/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)doublevirtintegrand, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
//Cuba integration of LO
void bornintegrMC4d(vector <double> &res, double &err)
{
const int ndim = 4; //dimensions of the integral
const int ncomp = 1; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsBORN;
const int maxeval = opts.vegasncallsBORN;
const int nstart = max(10, int(opts.vegasncallsBORN/10));
const int nincrease = max(10, int(opts.vegasncallsBORN/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)lointegrandMC, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
return;
}
//Cuba integration of LO
void bornintegr2d(vector <double> &res, double &err)
{
const int ndim = 2; //dimensions of the integral
const int ncomp = (opts.helicity >= 0 ? 2 : 1); //components of the integrand
void *userdata = NULL;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 65+2*65*opts.niterBORN;
const int maxeval = 65+2*65*opts.niterBORN;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) lointegrand2d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[2] = {0, 0};
double xmax[2] = {1, 1};
if (opts.cubacores == 0)
pcubature(ncomp, lointegrand2d_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, lointegrand2d_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
/*
//smolyak
int print_stats = 0;
int dim = 2;
int l = 10;
integral[0] = int_smolyak (ndim, ndim+l, lointegrand2d_smolyak, print_stats );
error[0] = 0.000001;
*/
}
res.clear();
/*
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
err = error[0];
*/
if (opts.helicity >= 0)
res.push_back(integral[0] != 0? integral[1]/integral[0] : 0);
else
res.push_back(integral[0]);
for (int i = 1; i < opts.totpdf; i++)
res.push_back(0);
if (opts.helicity >= 0)
err = integral[0] != 0? error[0]/integral[0]*integral[1]/integral[0] : 0; // error[1]/integral[1]*integral[1]/integral[0]
else
err = error[0];
return;
return;
}
//Cuba integration of V + 2 jets
void v2jintegr(vector <double> &res, double &err)
{
const int ndim = 10; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsVJREAL;
const int maxeval = opts.vegasncallsVJREAL;
const int nstart = max(10, int(opts.vegasncallsVJREAL/10));
const int nincrease = max(10, int(opts.vegasncallsVJREAL/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)v2jintegrand, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
/***************************************************************/
/***************************************************************/
// Counter term
//Cuba integration of the counterterm
void ctintegr(vector <double> &res, double &err)
{
const int ndim = 8; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsCT;
const int maxeval = opts.vegasncallsCT;
const int nstart = max(10, int(opts.vegasncallsCT/10));
const int nincrease = max(10, int(opts.vegasncallsCT/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)ctintegrand, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
//Cuba integration of the counterterm
void ctintegrMC(vector <double> &res, double &err)
{
const int ndim = 6; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 8+(4*!opts.vegascollect)+opts.cubaverbosity;
const int seed = opts.rseed;
const int mineval = opts.vegasncallsCT;
const int maxeval = opts.vegasncallsCT;
const int nstart = max(10, int(opts.vegasncallsCT/10));
const int nincrease = max(10, int(opts.vegasncallsCT/10));
const int nbatch = opts.cubanbatch;
const int gridno = 0;
Vegas(ndim, ncomp, (integrand_t)ctintegrandMC, userdata, nvec,
epsrel, epsabs,
flags, seed,
mineval, maxeval,
nstart, nincrease, nbatch,
gridno, statefile, spin,
&neval, &fail,
integral, error, prob);
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
return;
}
void ctintegr3d(vector <double> &res, double &err)
{
const int ndim = 3; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 127+2*127*opts.niterCT;
const int maxeval = 127+2*127*opts.niterCT;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) ctintegrand3d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[3] = {0, 0., 0.};
double xmax[3] = {1, 1., 1.};
if (opts.cubacores == 0)
pcubature(ncomp, ctintegrand3d_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, ctintegrand3d_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
}
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
//hists.FillQuadrature(res[0],err);
return;
}
void ctintegr2d(vector <double> &res, double &err)
{
const int ndim = 2; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata = NULL;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 65+2*65*opts.niterCT;
const int maxeval = 65+2*65*opts.niterCT;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) ctintegrand2d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[2] = {0, 0.};
double xmax[2] = {1, 1.};
if (opts.cubacores == 0)
pcubature(ncomp, ctintegrand2d_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, ctintegrand2d_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
}
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
//hists.FillQuadrature(res[0],err);
return;
}
void ctintegr1d(vector <double> &res, double &err)
{
const int ndim = 1; //dimensions of the integral
const int ncomp = opts.totpdf; //components of the integrand
void *userdata = NULL;
const int nvec = 1;
const double epsrel = 0.;
const double epsabs = 0.;
const char *statefile = "";
void *spin=NULL;
int neval;
int fail;
double integral[ncomp];
double error[ncomp];
double prob[ncomp];
const int flags = 0+opts.cubaverbosity;
const int mineval = 65+2*65*opts.niterCT;
const int maxeval = 65+2*65*opts.niterCT;
const int key = 13;
int nregions;
if (!opts.pcubature)
Cuhre(ndim, ncomp,
(integrand_t) ctintegrand1d, userdata, nvec,
epsrel, epsabs,
flags,
mineval, maxeval,
key, statefile, NULL,
&nregions, &neval, &fail,
integral, error, prob);
else
{
const int eval = 0;
const double epsrel = opts.relaccuracy;
const double epsabs = opts.absaccuracy;
double xmin[1] = {0.};
double xmax[1] = {1.};
if (opts.cubacores == 0)
pcubature(ncomp, ctintegrand1d_cubature, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
else
pcubature_v(ncomp, ctintegrand1d_cubature_v, userdata,
ndim, xmin, xmax,
eval, epsabs, epsrel, ERROR_INDIVIDUAL, integral, error);
}
res.clear();
for (int i = 0; i < opts.totpdf; i++)
res.push_back(integral[i]);
err = error[0];
//hists.FillQuadrature(res[0],err);
return;
}
/***************************************************************/
void initfun(void * input, const int &core){
if(core != PARENT_PROC) HistoHandler::Reset();
}
void exitfun(void * input, const int &core){
HistoHandler::Save(core);
}
void tell_to_grid_we_are_alive(){
if(opts.gridverbose && ICALL % 100000==0)
printf (" Hi Grid, we are sitll alive! Look, our event is %d\n",ICALL);
ICALL++;
}