//pdfevol::bscale = C3*b0/b is the "soft" factorisation scale at which the Wilson coefficient functions are evaluated (see arxiv:1309.1393 for details)
//double bstar = real(b / sqrt(1.+(b*b)/(1.1229190*1.1229190))); //// -->hard-coded!! --> should allow a different blim in the PDF evolution as a setting
*****************************************************/ // Moved to pdfevol::scales (end)
//The integration from b to qt space is done with the bstar prescription (real axis in the b space), and use the bessel function
//The (complex) integration in the minimal prescription would require hankel functions
//********************
//qt and b dependence (bessel function) (is xj0 a jacobian? Probably yes, the Jacobian for the change of variable from cartesian to radial coordinate, which translates from Fourier to Hankel transform)
complex <double> xj0;
double qtb = qt*real(b);
if (resint::_mode == 3 || resint::_mode == 4) //qt-integrated mode --> Needs to be implemented also for the other prescriptions
//alphasl gives the LL/NLL/NNLL evolution of alphas from Qres=Q/a_param to q2=b0^2/b^2
pdfevol::alphasl(b);
pdfevol::scales(b);
/********************************************
//pdfevol::alpqf is used as starting scale of the PDF evolution in pdfevol::evolution
//according to Eq. 42 of arXiv:hep-ph/0508068 it should be the factorisation scale,
//but in Eq. 98 in the resummation scale is used
double alpqf = resint::alpqres; //alphas at resummation scale (as in dyres)
//double alpqf = resint::alpqfac; //alphas at factorisation scale (as in Eq. 42 of arXiv:hep-ph/0508068, but see also Eq. 98)
********************************************/
/********************************************
//alpq is used in hcoefficients::calcb, it is alpq = alphas(b0^2/b^2)
//it is used only at NLL, at NNLL instead aexp and aexpb are used (aexp is the same as alphasl, but with a different blim)
complex <double> alpq;
alpq = resint::alpqres * cx(alphasl_(fscale2_mub)); //--> Attention! in DYRES it is alphas(qres), in DyQT it is alphas(mur)
if (opts.evolmode == 2 || opts.evolmode == 3)
//In order to have variable flavour evolution, use here a VFN definition of alpq
//There is possibly an issue here when the renormalisation scale is (very) different from mll, since aass=alphas(muren) is used in xlambda = beta0*aass*blog
// //In order to have variable flavour evolution, use here a VFN definition of aexp (aexpB instead, should be independent from NF)
// //It seems that the pt-integrated cross section is invariant for variations of a_param if aexp is the ratio of alpq/as(muren) (and not alpq/as(Qres))
// //Indeed, aexp is always multiplied by aass, which is alphas at the renormalisation scale
// if (opts.evolmode == 2 || opts.evolmode == 3) //account for VFN evolution
// Cache the positive and negative branch of coefficients which depend only on one I index
expc::reset();
expc::calc(b);
if (opts.mufvar)
{
muf::reset();
muf::calc(b);
}
/*
if (opts.mellin1d)
{
//if (opts.mufevol)
//expc::reset();
//hcoeff::calcb(resint::aass,resint::logmuf2q2,resint::loga,pdfevol::alpq,expc::aexp,expc::aexpB); // --> Need to access aass,logmuf2q2,loga,alpq,aexp,aexpb
}
else
//--> Implement 2D expc here
{
//complex <double> aexpb = cx(exponent_.aexpb_);
//complex <double> aexp = cx(exponent_.aexp_); //aexp is actually the ratio alphas(a*b0/b)/alphas(muren)
//hcoefficients::calcb(resint::aass,resint::logmuf2q2,resint::loga,pdfevol::alpq,expc::aexp,expc::aexpB); // --> Need to access aass,logmuf2q2,loga,alpq,aexp,aexpb
//expc::reset();
//hcoefficients::calcb(resint::aass,resint::logmuf2q2,resint::loga,pdfevol::alpq,expc::aexp,expc::aexpB); // --> Need to access aass,logmuf2q2,loga,alpq,aexp,aexpb
}
*/
//Inverting the HN coefficients from N to z space does not work, because it would miss the convolution with PDFs
//double q2 = pow(phasespace::m,2);
//if (opts.mellin1d)
//hcoeff::invert(q2);
//Truncate HN coefficients, to implement mellin1d also for rapdity intervals (currently not working, not sure if it is eventually possible. Need to figure out how to truncate the moments)
//if (opts.mellin1d)
//hcoeff::truncate();
//Apply flavour dependent non-perturbative form factors
if (opts.flavour_kt)
pdfevol::flavour_kt();
complex <double> fun = 0.;
mellinint::reset();
//Do not need the Mellin transform for the LL case, because the HN coefficient is 1 --> Use PDFs in x space