diff --git a/resum/resint.C b/resum/resint.C index 89145ae..be51de3 100644 --- a/resum/resint.C +++ b/resum/resint.C @@ -1,953 +1,1042 @@ #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 "bequad.h" #include #include "intde2_c.h" #include 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 resint::loga; double resint::rloga; complex resint::logmuf2q2; complex resint::logq2muf2; complex 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 jac(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr)); complex 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 jac(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr)); complex 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 jac(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr)); complex 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 jac(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr)); complex 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 jacu(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr)); complex bu = resint::bc + jacu*x; fcomplex argu = fcx(resint::_qt*bu); complex jacd(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr)); complex bd = resint::bc + jacd*x; fcomplex argd = fcx(resint::_qt*bd); complex 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 bint = besselint::bint(bb); - double J0 = fort_besj0_(qtb); - double Y0 = fort_besy0_(qtb); - double res = real(bint)*J0-imag(bint)*Y0; + double JN, YN; + if (resint::_mode == 3) + { + JN = fort_besj1_(qtb); + YN = fort_besy1_(qtb); + } + else + { + JN = fort_besj0_(qtb); + YN = fort_besy0_(qtb); + } + double res = real(bint)*JN-imag(bint)*YN; //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 jacu(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr)); complex bu = resint::bc + jacu*x; fcomplex qtbu = fcx(resint::_qt*bu); complex jacd(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr)); complex bd = resint::bc + jacd*x; fcomplex qtbd = fcx(resint::_qt*bd); complex h1, h2; // //call modified Hankel functions // h1 = cx(h1_(qtbu)); // h2 = cx(h2_(qtbd)); // //cout << bu << " " << h1 << " " << bd << " " << h2 << endl; //call original Hankel functions + // --> Here should use Hankel with J1 and Y1 for qt-integrated mode = 3!!! int n = 1; + int nn; + if (resint::_mode == 3) + nn = 1; + else + nn = 0; int nm; + //fcomplex chf1[n+1], chd1[n+1], chf2[n+1], chd2[n+1]; fcomplex chf1[2], chd1[2], chf2[2], chd2[2]; ch12n_ (n, qtbu, nm, chf1, chd1, chf2, chd2); - h1 = cx(chf1[0]); + h1 = cx(chf1[nn]); ch12n_ (n, qtbd, nm, chf1, chd1, chf2, chd2); - h2 = cx(chf2[0]); + h2 = cx(chf2[nn]); //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 (fort_besj0_(qtbu), fort_besy0_(qtbu)); //h2 = complex (fort_besj0_(qtbd), -fort_besy0_(qtbd)); //h2 = conj(h1); complex bint_up = besselint::bint(bu); complex bint_dn = besselint::bint(bd); //complex 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 jacu(cos(M_PI/opts.phibr),sin(M_PI/opts.phibr)); complex bu = resint::bc + jacu*x; fcomplex argu = fcx(resint::_qt*bu); complex 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 jacd(cos(M_PI/opts.phibr),-sin(M_PI/opts.phibr)); complex bd = resint::bc + jacd*x; fcomplex argd = fcx(resint::_qt*bd); complex 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.g1; //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; alpqren=pdf::alphas(muren)/4./M_PI; alpqfac=pdf::alphas(mufac)/4./M_PI; alpqres=pdf::alphas(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 = 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; +// 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); + + double res1, res2; _qt = qtmx; - intdeo(besselint_dequad, 0.0, qtmx, awinf, &res2, &err2); - - double res1, err1; + res2 = bintegral(qtmx); _qt = qtmn; if (qtmn < opts.qtcutoff) res1 = 0.; else - intdeo(besselint_dequad, 0.0, qtmn, awinf, &res1, &err1); + res1 = bintegral(qtmn); 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)) //if (opts.order == 0 && opts.xspace && (opts.evolmode == 1)) + if (mode < 2) //rapidity differential + res *= mesq::loxs(x1, x2, mufac); + else //rapidity integrated + res *= mesq::loxs(tau, mufac); + + /* 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;; } } } + //cout << "resint " << fun << endl; 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; //Warning!!! 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) //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; + 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); + // cout << endl; + // cout << "dequad result of inverse bessel transform, pt=" << _qt << " m=" << _m << " y=" << _y << " : " << setprecision(16) << res << " +- " << err << endl; + + int bqrule = 20; + complex h1 = 0; + complex h2 = 0; + complex bb; + for (int i = 0; i < bqrule; i++) + { + bb = complex {0,bq::xxx[bqrule-1][i]/qt}; + h1 += 2./M_PI/qt * bq::www[bqrule-1][i] * besselint::bint(bb); + // cout << "bb " << bb << " h1 " << h1 << endl; + bb = complex {0,-bq::xxx[bqrule-1][i]/qt}; + h2 += 2./M_PI/qt * bq::www[bqrule-1][i] * besselint::bint(bb); + // cout << "bb " << bb << " h2 " << h2 << endl; + } + res = real(1./2.*(h1+h2)); + + // cout << "bequad result " << res << endl; + // cout.precision(6); cout.unsetf(ios_base::floatfield); + /* + 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 NU = (_mode != 3)?0:1; //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_; + res = 0; + if (opts.bcf > 0) + { + //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; + //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 = (_mode != 3)?0:1; //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 r = 0.; complex 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 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; + if (errdequad < 0) + cout << "warning: dequad abnormal termination, pt=" << _qt << " m=" << _m << " y=" << _y << " bint: " << resdequad << " err " << errdequad << endl; res += resdequad; //cout << "Total bp = 2 " << res << endl; + + /* + //test alternative integrations + // cout << endl; + // cout << "dequad result of inverse bessel transform, pt=" << _qt << " m=" << _m << " y=" << _y << " : " << setprecision(16) << resdequad << " +- " << errdequad << endl; + + int bqrule = 40; + complex h1 = 0; + complex h2 = 0; + complex bb; + for (int i = 0; i < bqrule; i++) + { + bb = complex {0,bq::xxx[bqrule-1][i]/qt}; + h1 += 2./M_PI/qt * bq::www[bqrule-1][i] * besselint::bint(bb); + // cout << "bb " << bb << " h1 " << h1 << endl; + bb = complex {0,-bq::xxx[bqrule-1][i]/qt}; + h2 += 2./M_PI/qt * bq::www[bqrule-1][i] * besselint::bint(bb); + // cout << "bb " << bb << " h2 " << h2 << endl; + } + res = real(1./2.*(h1+h2)); + + // cout << "bequad result " << res << endl; + // cout.precision(6); cout.unsetf(ios_base::floatfield); + */ + } //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) + res = 0; + if (opts.bcf > 0) + { + bc = opts.bcf*blimit_.rblim_; - 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; + //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 = (_mode != 3)?0:1; //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; + //cout << "Real dequad result " << resdequad << " error " << errdequad << endl; res += resdequad; - cout << "Total bp = 3 " << res << endl; + //cout << "Total bp = 3 " << res << endl; } // //vfn integration --> Not correct!!! // //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; }