diff --git a/resum/pdfevol.C b/resum/pdfevol.C index f3bf49a..3da4a92 100644 --- a/resum/pdfevol.C +++ b/resum/pdfevol.C @@ -1,1388 +1,1513 @@ #include "pdfevol.h" #include "interface.h" #include "settings.h" #include "mesq.h" #include "anomalous.h" #include "resconst.h" #include "chebyshev.h" #include "phasespace.h" #include "mellinpdf.h" #include "scales.h" #include "parton.h" #include "npff.h" #include "string.h" #include "clock_real.h" #include #include #include //PDFs mellin moments at the factorisation scale complex *pdfevol::UVP; complex *pdfevol::DVP; complex *pdfevol::USP; complex *pdfevol::DSP; complex *pdfevol::SSP; complex *pdfevol::GLP; complex *pdfevol::CHP; complex *pdfevol::BOP; +complex *pdfevol::SIP; +complex *pdfevol::NS3P; +complex *pdfevol::NS8P; +complex *pdfevol::NS15P; +complex *pdfevol::NS24P; +complex *pdfevol::NS35P; + complex *pdfevol::fx1; complex *pdfevol::fx2; complex pdfevol::fn1[2*MAXNF+1]; complex pdfevol::fn2[2*MAXNF+1]; //scales complex pdfevol::bscale; complex pdfevol::bstarscale; complex pdfevol::bstartilde; complex pdfevol::qbstar; complex pdfevol::bcomplex; complex pdfevol::XL; complex pdfevol::XL1; complex pdfevol::SALP; complex pdfevol::alpr; using namespace parton; //fortran interface void pdfevol_(int& i1, int& i2, int& sign) { pdfevol::retrieve(i1-1, i2-1, sign-1); }; //moments at the starting scale void pdfevol::allocate() { UVP = new complex [mellinint::mdim]; DVP = new complex [mellinint::mdim]; USP = new complex [mellinint::mdim]; DSP = new complex [mellinint::mdim]; SSP = new complex [mellinint::mdim]; GLP = new complex [mellinint::mdim]; CHP = new complex [mellinint::mdim]; BOP = new complex [mellinint::mdim]; + + SIP = new complex [mellinint::mdim]; + NS3P = new complex [mellinint::mdim]; + NS8P = new complex [mellinint::mdim]; + NS15P = new complex [mellinint::mdim]; + NS24P = new complex [mellinint::mdim]; + NS35P = new complex [mellinint::mdim]; } void pdfevol::free() { delete[] UVP; delete[] DVP; delete[] USP; delete[] DSP; delete[] SSP; delete[] GLP; delete[] CHP; delete[] BOP; + + delete[] SIP; + delete[] NS3P; + delete[] NS8P; + delete[] NS15P; + delete[] NS24P; + delete[] NS35P; } //evolved moments void pdfevol::allocate_fx() { fx1 = new complex [mellinint::mdim*(2*MAXNF+1)*2]; fx2 = new complex [mellinint::mdim*(2*MAXNF+1)*2]; } void pdfevol::free_fx() { delete[] fx1; delete[] fx2; } void pdfevol::init() { //calculate Mellin moments of PDFs if (opts.fmufac == 0) if (!opts.silent) cout << "Initialise PDF moments with numerical integration... " << flush; //double xmin = 1e-8; double xmin = pow(bins.mbins.front()/opts.sroot,2); //Restrict the integration of moments to xmin = m/sqrt(s)*exp(-ymax) = (m/sqrt(s))^2 mellinpdf::init(xmin); clock_t begin_time, end_time; //Old Fortran code for the moments /* fcomplex uval,dval,usea,dsea,splus,ssea,glu,charm,bot; begin_time = clock(); for (int k = 0; k < mellinint::mdim; k++) { int hadron = 1; //opts.ih1; fcomplex XN = fcx(mellinint::Np[k]); //compute positive branch only, the negative branch is obtained by complex conjugation double facscale = opts.kmufac*opts.rmass; pdfmoments_(hadron,facscale,XN,uval,dval,usea,dsea,splus,ssea,glu,charm,bot,xmin); // cout << "moment " << k << " " << cx(XN) << " "; // cout << "uval " << cx(uval) << endl; // cout << "dval " << dval << endl; // cout << "usea " << usea << endl; // cout << "dsea " << dsea << endl; // cout << "gluon " << cx(glu) << endl; // cout << "charm " << charm << endl; // cout << "bottom" << bot << endl; UVP[k] = cx(uval); DVP[k] = cx(dval); USP[k] = cx(usea); DSP[k] = cx(dsea); SSP[k] = cx(ssea);//SSP[k] = cx(splus); !! issue for PDFs with s-sbar asymmetry !! GLP[k] = cx(glu); CHP[k] = cx(charm); BOP[k] = cx(bot); } end_time = clock(); cout << "Done " << float(end_time - begin_time)/CLOCKS_PER_SEC << endl; */ /* //polynomial interpolation for analytical continuation int order = 500; complex a = opts.cpoint+1-2*opts.zmax*1i; complex b = opts.cpoint+1+2*opts.zmax*1i; complex c = 0.5*(a+b); complex m = 0.5*(b-a); complex f[order]; for (int i = 1; i <= order; i++) { complex x = c+m*cheb::xxx[order-1][i-1]; int hadron = 1; //opts.ih1; fcomplex XN = fcx(x); double facscale = 1.0;//opts.kmufac*opts.rmass; pdfmoments_(hadron,facscale,XN,uval,dval,usea,dsea,splus,ssea,glu,charm,bot); f[i-1] = cx(uval); } for (int i = 0; i < 100; i++) { complex x = a+(b-a)*double(i)/100.; int hadron = 1; //opts.ih1; fcomplex XN = fcx(x); double facscale = 1.0;//opts.kmufac*opts.rmass; pdfmoments_(hadron,facscale,XN,uval,dval,usea,dsea,splus,ssea,glu,charm,bot); cout << setw(30) << x << setw(30) << cx(uval) << setw(30) << cheb::ipol(a,b,order, f, x) << setw(30) << cx(uval)-cheb::ipol(a,b,order, f, x) << endl; } for (int i = 1; i <= order; i++) { complex x = c+m*cheb::xxx[order-1][i-1]; cout << x << " " << f[i-1] << " " << cheb::ipol(a,b,order, f, x) << endl; } //verify for (int k = 0; k < mellinint::mdim; k++) { complex z = mellinint::Np[k]; cout << setw(30) << z << setw(30) << UVP[k] << setw(30) << cheb::ipol(a,b,order, f, z) << endl; } */ if (opts.fmufac == 0) { allocate(); scales::set(opts.rmass); begin_time = clock(); update(); end_time = clock(); if (!opts.silent) cout << "Done in " << float(end_time - begin_time)/CLOCKS_PER_SEC*1000. << "ms" << endl; } } void pdfevol::release() { if (opts.fmufac == 0) free(); } void pdfevol::update() { clock_t begin_time, end_time; //Restrict the integration of moments to xmin = m/sqrt(s)*exp(-ymax) = (m/sqrt(s))^2 //double xmin = pow(phasepace::m/opts.sroot,2); //mellinpdf::init(xmin); //scales::set(phasespace::m); //assume scales were already set mellinpdf::allocate(); begin_time = clock(); mellinpdf::evalpdfs(scales::fac); if (opts.mellininv == 1 || opts.phi > 0.5) mellinpdf::laguerre_ipol(); else mellinpdf::gauss_quad(); end_time = clock(); //cout << "mll is " << scales::fac << " x to N done in " << float(end_time - begin_time)/CLOCKS_PER_SEC*1000. << "ms" << endl; for (int k = 0; k < mellinint::mdim; k++) { UVP[k] = mellinpdf::UV[k]; DVP[k] = mellinpdf::DV[k]; USP[k] = mellinpdf::US[k]; DSP[k] = mellinpdf::DS[k]; SSP[k] = mellinpdf::SM[k]; GLP[k] = mellinpdf::GL[k]; CHP[k] = mellinpdf::CP[k]; BOP[k] = mellinpdf::BP[k]; //cout << "moment " << k << " " << mellinint::Np[k] << " "; //cout << "uval " << mellinpdf::UV[k] << endl; //cout << "dval " << mellinpdf::DV[k] << endl; //cout << "usea " << mellinpdf::US[k] << endl; //cout << "dsea " << mellinpdf::DS[k] << endl; //cout << "ssea " << mellinpdf::SM[k] << endl; //cout << "gluon " << mellinpdf::GL[k] << endl; //cout << "charm " << mellinpdf::CP[k] << endl; //cout << "bottom" << mellinpdf::BP[k] << endl; } mellinpdf::free(); + + //Precompute singlet/non-singlet decomposition + if (opts.evolmode == 0) + { + //N flavour dependence + int nf = resconst::NF; + for (int k = 0; k < mellinint::mdim; k++) + { + NS3P[k] = UVP[k] + 2.*USP[k] - DVP[k] - 2.*DSP[k]; //(u+ub-d-db) u-d + NS8P[k] = UVP[k] + 2.*USP[k] + DVP[k] + 2.*DSP[k] - 4.*SSP[k]; //(u+ub+d+db-s-sb) u+d-s + + if (nf == 5) + { + SIP[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k] + 2.*CHP[k] + 2.*BOP[k]; //(u+ub+d+db+s+sb+c+cb+b+bb) -> all quarks u+d+s+c+b + NS15P[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k] - 6.*CHP[k]; //(u+ub+d+db+s+sb-3(c+cb)) u+d+s-3c + NS24P[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k] + 2.*CHP[k] - 8.*BOP[k]; //(u+ub+d+db+s+sb+c+cb-4(b+bb)) u+d+s+c-4b + NS35P[k] = SIP[k]; + } + if (nf == 4) + { + SIP[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k] + 2.*CHP[k]; //(u+ub+d+db+s+sb+c+cb) -> all quarks + NS15P[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k] - 6.*CHP[k]; + NS24P[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k] + 2.*CHP[k]; + NS35P[k] = SIP[k]; + } + if (nf == 3) + { + SIP[k] = UVP[k] + DVP[k] + 2.*USP[k] + 2.*DSP[k] + 2.*SSP[k]; //(u+ub+d+db+s+sb) -> all quarks + NS15P[k] = SIP[k]; + NS24P[k] = SIP[k]; + NS35P[k] = SIP[k]; + } + } + } } // PROVIDES MOMENTS OF DENSITIES AT A GIVEN SCALE (INCLUDES EVOLUTION) // AND ANOMALOUS DIMENSIONS // EVERYTHING IN MELLIN SPACE // // input: I point in z grid, ISIGN (+ or -), IBEAM (1 or 2), SCALE2 (b), alps = alpqf (mass dependent) // output: FN, alpq = alps * alphasl(scale2) // dependence: FN(b, mass, I) // // *************************************** // Can completely cache the output FN in FN(I,IFIT,ISIG), and calculate in the init? No because of b // reno2 is cached and looped only on I, before entering the I1 I2 doube loop into cfx1 cfx2p cfx2m //Evolve the Mellin moment of PDF corresponding to the index i, sign sign, beam beam, from the scale q2 to the scale scale2(b) //Output in fx -void pdfevol::evolution(int i) //from reno2 +void pdfevol::evolution() //from reno2 { - // i is the index of the complex mellin moment in the z-space for the gaussian quadrature used for the mellin inversion //N flavour dependence int nf = resconst::NF; + //complex alpq, ALPr; + + //evolve only the moments of the positive branch, the negative branch is obtained by complex conjugation + int sign = mesq::positive; + + complex fx[11]; + //At LL there is no PDF evolution, PDFs are evaluated at the factorisation scale if (opts.order == 0) + for (int i = 0; i < mellinint::mdim; i++) + { + //XP[i] are moments of PDFs at the starting scale (factorisation scale) + fx[0+MAXNF] = GLP[i]; + fx[1+MAXNF] = UVP[i] + USP[i]; + fx[-1+MAXNF] = USP[i]; + fx[2+MAXNF] = DVP[i] + DSP[i]; + fx[-2+MAXNF] = DSP[i]; + fx[3+MAXNF] = SSP[i]; + fx[-3+MAXNF] = SSP[i]; + if (nf >= 4) + { + fx[4+MAXNF] = CHP[i]; + fx[-4+MAXNF] = CHP[i]; + } + else + { + fx[4+MAXNF] = 0.; + fx[-4+MAXNF] = 0.; + } + if (nf >= 5) + { + fx[5+MAXNF] = BOP[i]; + fx[-5+MAXNF] = BOP[i]; + } + else + { + fx[5+MAXNF] = 0.; + fx[-5+MAXNF] = 0.; + } + + storemoments(i, fx); + continue; + } + + //i is the index of the complex mellin moment in the z-space for the gaussian quadrature used for the mellin inversion + for (int i = 0; i < mellinint::mdim; i++) { - //XP[i] are moments of PDFs at the starting scale (factorisation scale) - complex fx[11]; - fx[0+MAXNF] = GLP[i]; - fx[1+MAXNF] = UVP[i] + USP[i]; - fx[-1+MAXNF] = USP[i]; - fx[2+MAXNF] = DVP[i] + DSP[i]; - fx[-2+MAXNF] = DSP[i]; - fx[3+MAXNF] = SSP[i]; - fx[-3+MAXNF] = SSP[i]; + /* + //Moments at the factorisation scale + //complex UVI, DVI, USI, DSI, SSI, GLI, CHI, BOI; + + //Moments of PDFs at the starting scale (factorisation scale) + UVI = UVP[i]; + DVI = DVP[i]; + USI = USP[i]; + DSI = DSP[i]; + SSI = SSP[i]; + GLI = GLP[i]; + CHI = CHP[i]; + BOI = BOP[i]; + + // *********************** + // this part can be precomputed + complex UVN = UVI; + complex DVN = DVI; + complex NS3N = UVI + 2.*USI - DVI - 2.*DSI; //(u+ub-d-db) u-d + complex NS8N = UVI + 2.*USI + DVI + 2.*DSI - 4.*SSI; //(u+ub+d+db-s-sb) u+d-s + complex GLN = GLI; + + complex SIN, NS15N, NS24N, NS35N; + if (nf == 5) + { + SIN = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI + 2.*BOI; //(u+ub+d+db+s+sb+c+cb+b+bb) -> all quarks u+d+s+c+b + NS15N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI - 6.*CHI; //(u+ub+d+db+s+sb-3(c+cb)) u+d+s-3c + NS24N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI - 8.*BOI; //(u+ub+d+db+s+sb+c+cb-4(b+bb)) u+d+s+c-4b + NS35N = SIN; + } + if (nf == 4) + { + SIN = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI; //(u+ub+d+db+s+sb+c+cb) -> all quarks + NS15N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI - 6.*CHI; + NS24N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI; + NS35N = SIN; + } + if (nf == 3) + { + SIN = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI; //(u+ub+d+db+s+sb) -> all quarks + NS15N = SIN; + NS24N = SIN; + NS35N = SIN; + } + complex SG = SIN; + complex GL = GLN; + */ + + // Singlet/non-singlet decomposition at the starting scale (factorisation scale) + complex UVN = UVP[i]; + complex DVN = DVP[i]; + complex GLN = GLP[i]; + complex SIN = SIP[i]; + complex NS3N = NS3P[i]; + complex NS8N = NS8P[i]; + complex NS15N = NS15P[i]; + complex NS24N = NS24P[i]; + complex NS35N = NS35P[i]; + + complex SG = SIN; + complex GL = GLN; + + + // retrieved values cached in anomalous.C + complex ANS = anomalous::ans[anomalous::index(i,sign)]; + complex AM = anomalous::am[anomalous::index(i,sign)]; + complex AP = anomalous::ap[anomalous::index(i,sign)]; + complex AL = anomalous::al[anomalous::index(i,sign)]; + complex BE = anomalous::be[anomalous::index(i,sign)]; + complex AB = anomalous::ab[anomalous::index(i,sign)]; + complex AC = 1. -AL; + complex RMIN = anomalous::rmin[anomalous::index(i,sign)]; + complex RPLUS = anomalous::rplus[anomalous::index(i,sign)]; + complex RMMQQ = anomalous::RMMQQ[anomalous::index(i,sign)]; + complex RMMQG = anomalous::RMMQG[anomalous::index(i,sign)]; + complex RMMGQ = anomalous::RMMGQ[anomalous::index(i,sign)]; + complex RMMGG = anomalous::RMMGG[anomalous::index(i,sign)]; + complex RMPQQ = anomalous::RMPQQ[anomalous::index(i,sign)]; + complex RMPQG = anomalous::RMPQG[anomalous::index(i,sign)]; + complex RMPGQ = anomalous::RMPGQ[anomalous::index(i,sign)]; + complex RMPGG = anomalous::RMPGG[anomalous::index(i,sign)]; + complex RPMQQ = anomalous::RPMQQ[anomalous::index(i,sign)]; + complex RPMQG = anomalous::RPMQG[anomalous::index(i,sign)]; + complex RPMGQ = anomalous::RPMGQ[anomalous::index(i,sign)]; + complex RPMGG = anomalous::RPMGG[anomalous::index(i,sign)]; + complex RPPQQ = anomalous::RPPQQ[anomalous::index(i,sign)]; + complex RPPQG = anomalous::RPPQG[anomalous::index(i,sign)]; + complex RPPGQ = anomalous::RPPGQ[anomalous::index(i,sign)]; + complex RPPGG = anomalous::RPPGG[anomalous::index(i,sign)]; + // ************************************** + + // ************************************** + // b-dependence + //resummation scale + // complex XL = 1./cx(alphasl_(fcx(scale2))); + // complex XL1 = 1.- XL; + // complex SALP = log(XL); + //--> SALP ~ log[alphas(Q)/alphas(b0/b)] + + complex S = SALP; + // cout << S << " " << < ENS = exp(-ANS*S); + complex EM = exp(-AM*S); + complex EP = exp(-AP*S); + complex EMP = EM/EP; + complex EPM = EP/EM; + + //...EVOLUTION OF LIGHT PARTON DENSITIES + //double q2s = q2/pow(resint::a,2); //resummation scale + //double alpqf = dyalphas_lhapdf_(sqrt(q2s))/4./M_PI; //alphas at the resummation scale + //complex alpq = alpqf * alphasl(scale2); //alphas at the resummation scale times alphas at 1/b + //complex alpr= alpq * 1 *(opts.order-1); + //--> alpr = 0 at NLL; alpr = alphas(Q) * alphasl ~ alphas(b0/b) at NNLL + + if (opts.order == 1) + { + UVN = UVN * ENS; + DVN = DVN * ENS; + NS3N = NS3N * ENS; + NS8N = NS8N * ENS; + + SIN = EM * (AL*SG + BE*GL) + EP * (AC*SG - BE*GL); + GLN = EM * (AB*SG + AC*GL) + EP * (-AB*SG + AL*GL); + + NS15N = NS15N * ENS; + NS24N = NS24N * ENS; + } + else if (opts.order == 2) + { + UVN = UVN * ENS * (1.+ alpr * XL1 * RMIN); + DVN = DVN * ENS * (1.+ alpr * XL1 * RMIN); + NS3N = NS3N * ENS * (1.+ alpr * XL1 * RPLUS); + NS8N = NS8N * ENS * (1.+ alpr * XL1 * RPLUS); + + SIN = EM * ((AL + alpr * (RMMQQ*XL1 + RMPQQ*(EPM-XL)))* SG + (BE + alpr * (RMMQG*XL1 + RMPQG*(EPM-XL))) * GL) + + EP * ((AC + alpr * (RPPQQ*XL1 + RPMQQ*(EMP-XL)))* SG + (-BE + alpr * (RPPQG*XL1 + RPMQG*(EMP-XL))) * GL); + GLN = EM * ((AB + alpr * (RMMGQ*XL1 + RMPGQ*(EPM-XL)))* SG + (AC + alpr * (RMMGG*XL1 + RMPGG*(EPM-XL))) * GL) + + EP *((-AB + alpr * (RPPGQ*XL1 + RPMGQ*(EMP-XL)))* SG + (AL + alpr * (RPPGG*XL1 + RPMGG*(EMP-XL))) * GL); + + NS15N = NS15N * ENS * (1.+ alpr * XL1 * RPLUS); + NS24N = NS24N * ENS * (1.+ alpr * XL1 * RPLUS); + } + + + NS35N = SIN; + + //... FLAVOUR DECOMPOSITION OF THE QUARK SEA : + complex SSN, DSN, USN, CHN, BON; + SSN = (10.* SIN + 2.* NS35N + 3.* NS24N + 5.* NS15N - 20.* NS8N) / 120.; + DSN = (10.* SIN + 2.* NS35N + 3.* NS24N + 5.* NS15N + 10.* NS8N - 30.* NS3N - 60.* DVN) / 120.; + USN = (10.* SIN + 2.* NS35N + 3.* NS24N + 5.* NS15N + 10.* NS8N + 30.* NS3N - 60.* UVN) / 120.; + CHN = (UVN + DVN + 2.*USN + 2.*DSN + 2.*SSN - NS15N)/6.; + BON = (UVN + DVN + 2.*USN + 2.*DSN + 2.*SSN + 2.*CHN - NS24N)/8.; + //equivalent to: + //CHN = (10.* SIN + 2. *NS35N + 3.* NS24N - 15.* NS15N) / 120.; + //BON = (10.* SIN + 2. *NS35N - 12.* NS24N) / 120.; + + if (nf == 3) //GRV + { + SSN= (20.* SIN - 20.* NS8N)/120.; + DSN = (20.* SIN + 10.* NS8N - 30.* NS3N - 60.* DVN) / 120.; + USN = (20.* SIN + 10.* NS8N + 30.* NS3N - 60.* UVN) / 120.; + CHN=0.; + BON=0.; + } + + //if (fabs(bstarscale) < LHAPDF::getThreshold(4)) + // CHN *= exp(-pow((LHAPDF::getThreshold(4)-fabs(bstarscale)),2)/pow((LHAPDF::getThreshold(4)/10.) ,2)); //=0 + // double delta = 1./5.; + // if (fabs(bstarscale) < LHAPDF::getThreshold(5)*(1+delta/2.)) + // BON *= exp(-pow((LHAPDF::getThreshold(5)*(1+delta/2.)-fabs(bstarscale)),2)/pow((LHAPDF::getThreshold(5)*delta),2)); //=0 + + // ************************************** + + // output: + // bbar cbar sbar dbar ubar gluon u d s c b + // -5 -4 -3 -2 -1 0 1 2 3 4 5 + + fx[0+MAXNF] = GLN; + fx[1+MAXNF] = UVN + USN; + fx[-1+MAXNF] = USN; + fx[2+MAXNF] = DVN + DSN; + fx[-2+MAXNF] = DSN; + fx[3+MAXNF] = SSN; + fx[-3+MAXNF] = SSN; if (nf >= 4) { - fx[4+MAXNF] = CHP[i]; - fx[-4+MAXNF] = CHP[i]; + fx[4+MAXNF] = CHN; + fx[-4+MAXNF] = CHN; } else { fx[4+MAXNF] = 0.; fx[-4+MAXNF] = 0.; } + if (nf >= 5) { - fx[5+MAXNF] = BOP[i]; - fx[-5+MAXNF] = BOP[i]; + fx[5+MAXNF] = BON; + fx[-5+MAXNF] = BON; } else { fx[5+MAXNF] = 0.; fx[-5+MAXNF] = 0.; } storemoments(i, fx); - return; - } - - complex alpq, ALPr; - - //Moments at the factorisation scale - complex UVI, DVI, USI, DSI, SSI, GLI, CHI, BOI; - - //evolve only the moments of the positive branch, the negative branch is obtained by complex conjugation - int sign = mesq::positive; - - //Moments of PDFs at the starting scale (factorisation scale) - UVI = UVP[i]; - DVI = DVP[i]; - USI = USP[i]; - DSI = DSP[i]; - SSI = SSP[i]; - GLI = GLP[i]; - CHI = CHP[i]; - BOI = BOP[i]; - - // *********************** - // this part can be precomputed - complex UVN = UVI; - complex DVN = DVI; - complex NS3N = UVI + 2.*USI - DVI - 2.*DSI; //(u+ub-d-db) u-d - complex NS8N = UVI + 2.*USI + DVI + 2.*DSI - 4.*SSI; //(u+ub+d+db-s-sb) u+d-s - complex GLN = GLI; - - complex SIN, NS15N, NS24N, NS35N; - if (nf == 5) - { - SIN = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI + 2.*BOI; //(u+ub+d+db+s+sb+c+cb+b+bb) -> all quarks u+d+s+c+b - NS15N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI - 6.*CHI; //(u+ub+d+db+s+sb-3(c+cb)) u+d+s-3c - NS24N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI - 8.*BOI; //(u+ub+d+db+s+sb+c+cb-4(b+bb)) u+d+s+c-4b - NS35N = SIN; - } - if (nf == 4) - { - SIN = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI; //(u+ub+d+db+s+sb+c+cb) -> all quarks - NS15N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI - 6.*CHI; - NS24N = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI + 2.*CHI; - NS35N = SIN; - } - if (nf == 3) - { - SIN = UVI + DVI + 2.*USI + 2.*DSI + 2.*SSI; //(u+ub+d+db+s+sb) -> all quarks - NS15N = SIN; - NS24N = SIN; - NS35N = SIN; - } - complex SG = SIN; - complex GL = GLN; - - // ************************************** - // retrieved values cached in anomalous.C - complex ANS = anomalous::ans[anomalous::index(i,sign)]; - complex AM = anomalous::am[anomalous::index(i,sign)]; - complex AP = anomalous::ap[anomalous::index(i,sign)]; - complex AL = anomalous::al[anomalous::index(i,sign)]; - complex BE = anomalous::be[anomalous::index(i,sign)]; - complex AB = anomalous::ab[anomalous::index(i,sign)]; - complex AC = 1. -AL; - complex RMIN = anomalous::rmin[anomalous::index(i,sign)]; - complex RPLUS = anomalous::rplus[anomalous::index(i,sign)]; - complex RMMQQ = anomalous::RMMQQ[anomalous::index(i,sign)]; - complex RMMQG = anomalous::RMMQG[anomalous::index(i,sign)]; - complex RMMGQ = anomalous::RMMGQ[anomalous::index(i,sign)]; - complex RMMGG = anomalous::RMMGG[anomalous::index(i,sign)]; - complex RMPQQ = anomalous::RMPQQ[anomalous::index(i,sign)]; - complex RMPQG = anomalous::RMPQG[anomalous::index(i,sign)]; - complex RMPGQ = anomalous::RMPGQ[anomalous::index(i,sign)]; - complex RMPGG = anomalous::RMPGG[anomalous::index(i,sign)]; - complex RPMQQ = anomalous::RPMQQ[anomalous::index(i,sign)]; - complex RPMQG = anomalous::RPMQG[anomalous::index(i,sign)]; - complex RPMGQ = anomalous::RPMGQ[anomalous::index(i,sign)]; - complex RPMGG = anomalous::RPMGG[anomalous::index(i,sign)]; - complex RPPQQ = anomalous::RPPQQ[anomalous::index(i,sign)]; - complex RPPQG = anomalous::RPPQG[anomalous::index(i,sign)]; - complex RPPGQ = anomalous::RPPGQ[anomalous::index(i,sign)]; - complex RPPGG = anomalous::RPPGG[anomalous::index(i,sign)]; - // ************************************** - - // ************************************** - // b-dependence - //resummation scale - // complex XL = 1./cx(alphasl_(fcx(scale2))); - // complex XL1 = 1.- XL; - // complex SALP = log(XL); - //--> SALP ~ log[alphas(Q)/alphas(b0/b)] - - complex S = SALP; - // cout << S << " " << < ENS = exp(-ANS*S); - complex EM = exp(-AM*S); - complex EP = exp(-AP*S); - complex EMP = EM/EP; - complex EPM = EP/EM; - - //...EVOLUTION OF LIGHT PARTON DENSITIES - //double q2s = q2/pow(resint::a,2); //resummation scale - //double alpqf = dyalphas_lhapdf_(sqrt(q2s))/4./M_PI; //alphas at the resummation scale - //complex alpq = alpqf * alphasl(scale2); //alphas at the resummation scale times alphas at 1/b - //complex alpr= alpq * 1 *(opts.order-1); - //--> alpr = 0 at NLL; alpr = alphas(Q) * alphasl ~ alphas(b0/b) at NNLL - - UVN = UVN * ENS * (1.+ alpr * XL1 * RMIN); - DVN = DVN * ENS * (1.+ alpr * XL1 * RMIN); - NS3N = NS3N * ENS * (1.+ alpr * XL1 * RPLUS); - NS8N = NS8N * ENS * (1.+ alpr * XL1 * RPLUS); - - SIN = EM * ((AL + alpr * (RMMQQ*XL1 + RMPQQ*(EPM-XL)))* SG + (BE + alpr * (RMMQG*XL1 + RMPQG*(EPM-XL))) * GL) - + EP * ((AC + alpr * (RPPQQ*XL1 + RPMQQ*(EMP-XL)))* SG + (-BE + alpr * (RPPQG*XL1 + RPMQG*(EMP-XL))) * GL); - GLN = EM * ((AB + alpr * (RMMGQ*XL1 + RMPGQ*(EPM-XL)))* SG + (AC + alpr * (RMMGG*XL1 + RMPGG*(EPM-XL))) * GL) - + EP *((-AB + alpr * (RPPGQ*XL1 + RPMGQ*(EMP-XL)))* SG + (AL + alpr * (RPPGG*XL1 + RPMGG*(EMP-XL))) * GL); - - NS15N = NS15N * ENS * (1.+ alpr * XL1 * RPLUS); - NS24N = NS24N * ENS * (1.+ alpr * XL1 * RPLUS); - NS35N = SIN; - - //... FLAVOUR DECOMPOSITION OF THE QUARK SEA : - complex SSN, DSN, USN, CHN, BON; - SSN = (10.* SIN + 2.* NS35N + 3.* NS24N + 5.* NS15N - 20.* NS8N) / 120.; - DSN = (10.* SIN + 2.* NS35N + 3.* NS24N + 5.* NS15N + 10.* NS8N - 30.* NS3N - 60.* DVN) / 120.; - USN = (10.* SIN + 2.* NS35N + 3.* NS24N + 5.* NS15N + 10.* NS8N + 30.* NS3N - 60.* UVN) / 120.; - CHN = (UVN + DVN + 2.*USN + 2.*DSN + 2.*SSN - NS15N)/6.; - BON = (UVN + DVN + 2.*USN + 2.*DSN + 2.*SSN + 2.*CHN - NS24N)/8.; - //equivalent to: - //CHN = (10.* SIN + 2. *NS35N + 3.* NS24N - 15.* NS15N) / 120.; - //BON = (10.* SIN + 2. *NS35N - 12.* NS24N) / 120.; - - if (nf == 3) //GRV - { - SSN= (20.* SIN - 20.* NS8N)/120.; - DSN = (20.* SIN + 10.* NS8N - 30.* NS3N - 60.* DVN) / 120.; - USN = (20.* SIN + 10.* NS8N + 30.* NS3N - 60.* UVN) / 120.; - CHN=0.; - BON=0.; - } - - //if (fabs(bstarscale) < LHAPDF::getThreshold(4)) - // CHN *= exp(-pow((LHAPDF::getThreshold(4)-fabs(bstarscale)),2)/pow((LHAPDF::getThreshold(4)/10.) ,2)); //=0 - // double delta = 1./5.; - // if (fabs(bstarscale) < LHAPDF::getThreshold(5)*(1+delta/2.)) - // BON *= exp(-pow((LHAPDF::getThreshold(5)*(1+delta/2.)-fabs(bstarscale)),2)/pow((LHAPDF::getThreshold(5)*delta),2)); //=0 - - // ************************************** - - // output: - // bbar cbar sbar dbar ubar gluon u d s c b - // -5 -4 -3 -2 -1 0 1 2 3 4 5 - - complex fx[11]; - fx[0+MAXNF] = GLN; - fx[1+MAXNF] = UVN + USN; - fx[-1+MAXNF] = USN; - fx[2+MAXNF] = DVN + DSN; - fx[-2+MAXNF] = DSN; - fx[3+MAXNF] = SSN; - fx[-3+MAXNF] = SSN; - if (nf >= 4) - { - fx[4+MAXNF] = CHN; - fx[-4+MAXNF] = CHN; } - else - { - fx[4+MAXNF] = 0.; - fx[-4+MAXNF] = 0.; - } - - if (nf >= 5) - { - fx[5+MAXNF] = BON; - fx[-5+MAXNF] = BON; - } - else - { - fx[5+MAXNF] = 0.; - fx[-5+MAXNF] = 0.; - } - - storemoments(i, fx); } //Old fortran code void pdfevol::calculate(int i) { //N flavour dependence int nf = resconst::NF; int hadron = 1; //double facscale = fabs(bscale); //double facscale = fabs(bstarscale); //double facscale = fabs(opts.muf); double facscale = fabs(pdfevol::bstartilde); //bstartilde = qbstar * resint::mures / sqrt(pow(qbstar,2) + resint::mures2); fcomplex XN = fcx(mellinint::Np[i]); double xmin = 1e-8; complex fx[11]; fcomplex uval,dval,usea,dsea,s,sbar,glu,charm,bot; pdfmoments_(hadron,facscale,XN,uval,dval,usea,dsea,s,sbar,glu,charm,bot,xmin); fx[0+MAXNF] = cx(glu); fx[1+MAXNF] = cx(uval) + cx(usea); fx[-1+MAXNF] = cx(usea); fx[2+MAXNF] = cx(dval) + cx(dsea); fx[-2+MAXNF] = cx(dsea); fx[3+MAXNF] = cx(s); fx[-3+MAXNF] = cx(sbar); if (nf >= 4) { fx[4+MAXNF] = cx(charm); fx[-4+MAXNF] = cx(charm); } else { fx[4+MAXNF] = 0.; fx[-4+MAXNF] = 0.; } if (nf >= 5) { fx[5+MAXNF] = cx(bot); fx[-5+MAXNF] = cx(bot); } else { fx[5+MAXNF] = 0.; fx[-5+MAXNF] = 0.; } storemoments(i, fx); } void pdfevol::calculate() { clock_t begin_time, end_time; //double facscale = fabs(bscale); //double facscale = fabs(bstarscale); //double facscale = fabs(opts.muf); double facscale = fabs(pdfevol::bstartilde); //bstartilde = qbstar * resint::mures / sqrt(pow(qbstar,2) + resint::mures2); //cout << facscale << endl; complex fx[11]; mellinpdf::allocate(); begin_time = clock(); mellinpdf::evalpdfs(facscale); if (opts.mellininv == 1 || opts.phi > 0.5) mellinpdf::laguerre_ipol(); else mellinpdf::gauss_quad(); end_time = clock(); //cout << "facscale is " << facscale << " x to N done in " << float(end_time - begin_time)/CLOCKS_PER_SEC*1000. << "ms" << endl; for (int k = 0; k < mellinint::mdim; k++) { fx[g] = mellinpdf::GL[k]; fx[u] = mellinpdf::UV[k]+mellinpdf::US[k]; fx[ub] = mellinpdf::US[k]; fx[d] = mellinpdf::DV[k]+mellinpdf::DS[k]; fx[db] = mellinpdf::DS[k]; fx[s] = mellinpdf::SP[k]; fx[sb] = mellinpdf::SM[k]; if (resconst::NF >= 4) { fx[c] = mellinpdf::CP[k]; fx[cb] = mellinpdf::CM[k]; } else { fx[c] = 0.; fx[cb] = 0.; } if (resconst::NF >= 5) { fx[b] = mellinpdf::BP[k]; fx[bb] = mellinpdf::BM[k]; } else { fx[b] = 0.; fx[bb] = 0.; } //cout << k << " " << fx[g] << endl; storemoments(k, fx); } mellinpdf::free(); } void pdfevol::storemoments(int i, complex fx[11]) { //Save the evolved PDFs into the fx1 and fx2 arrays int negidx = mellinint::mdim*11; int nf = 2*MAXNF+1; - //beam 1 positive - fx1[i*nf+bb] = fx[bb]; - fx1[i*nf+cb] = fx[cb]; - fx1[i*nf+sb] = fx[sb]; - fx1[i*nf+g ] = fx[g ]; - fx1[i*nf+s ] = fx[s ]; - fx1[i*nf+c ] = fx[c ]; - fx1[i*nf+b ] = fx[b ]; - if (opts.ih1 == 1) - { - fx1[i*nf+db] = fx[db]; - fx1[i*nf+ub] = fx[ub]; - fx1[i*nf+u ] = fx[u ]; - fx1[i*nf+d ] = fx[d ]; - } - else if (opts.ih1 == -1) + copy(fx, fx+nf, fx1+i*nf); + copy(fx, fx+nf, fx2+i*nf); + + if (opts.ih1 == -1) { fx1[i*nf+db] = fx[d ]; fx1[i*nf+ub] = fx[u ]; fx1[i*nf+u ] = fx[ub]; fx1[i*nf+d ] = fx[db]; } - - //beam 1 negative (never used) - fx1[negidx+i*nf+bb] = conj(fx[bb]); - fx1[negidx+i*nf+cb] = conj(fx[cb]); - fx1[negidx+i*nf+sb] = conj(fx[sb]); - fx1[negidx+i*nf+g ] = conj(fx[g ]); - fx1[negidx+i*nf+s ] = conj(fx[s ]); - fx1[negidx+i*nf+c ] = conj(fx[c ]); - fx1[negidx+i*nf+b ] = conj(fx[b ]); - if (opts.ih1 == 1) - { - fx1[negidx+i*nf+db] = conj(fx[db]); - fx1[negidx+i*nf+ub] = conj(fx[ub]); - fx1[negidx+i*nf+u ] = conj(fx[u ]); - fx1[negidx+i*nf+d ] = conj(fx[d ]); - } - else if (opts.ih1 == -1) - { - fx1[negidx+i*nf+db] = conj(fx[d ]); - fx1[negidx+i*nf+ub] = conj(fx[u ]); - fx1[negidx+i*nf+u ] = conj(fx[ub]); - fx1[negidx+i*nf+d ] = conj(fx[db]); - } - - //beam 2 positive - fx2[i*nf+bb] = fx[bb]; - fx2[i*nf+cb] = fx[cb]; - fx2[i*nf+sb] = fx[sb]; - fx2[i*nf+g ] = fx[g ]; - fx2[i*nf+s ] = fx[s ]; - fx2[i*nf+c ] = fx[c ]; - fx2[i*nf+b ] = fx[b ]; - if (opts.ih2 == 1) - { - fx2[i*nf+db] = fx[db]; - fx2[i*nf+ub] = fx[ub]; - fx2[i*nf+u ] = fx[u ]; - fx2[i*nf+d ] = fx[d ]; - } - else if (opts.ih2 == -1) + if (opts.ih2 == -1) { fx2[i*nf+db] = fx[d ]; fx2[i*nf+ub] = fx[u ]; fx2[i*nf+u ] = fx[ub]; fx2[i*nf+d ] = fx[db]; } +// //beam 1 positive +// fx1[i*nf+bb] = fx[bb]; +// fx1[i*nf+cb] = fx[cb]; +// fx1[i*nf+sb] = fx[sb]; +// fx1[i*nf+g ] = fx[g ]; +// fx1[i*nf+s ] = fx[s ]; +// fx1[i*nf+c ] = fx[c ]; +// fx1[i*nf+b ] = fx[b ]; +// if (opts.ih1 == 1) +// { +// fx1[i*nf+db] = fx[db]; +// fx1[i*nf+ub] = fx[ub]; +// fx1[i*nf+u ] = fx[u ]; +// fx1[i*nf+d ] = fx[d ]; +// } +// else if (opts.ih1 == -1) +// { +// fx1[i*nf+db] = fx[d ]; +// fx1[i*nf+ub] = fx[u ]; +// fx1[i*nf+u ] = fx[ub]; +// fx1[i*nf+d ] = fx[db]; +// } +// +// //beam 1 negative (never used) +// fx1[negidx+i*nf+bb] = conj(fx[bb]); +// fx1[negidx+i*nf+cb] = conj(fx[cb]); +// fx1[negidx+i*nf+sb] = conj(fx[sb]); +// fx1[negidx+i*nf+g ] = conj(fx[g ]); +// fx1[negidx+i*nf+s ] = conj(fx[s ]); +// fx1[negidx+i*nf+c ] = conj(fx[c ]); +// fx1[negidx+i*nf+b ] = conj(fx[b ]); +// if (opts.ih1 == 1) +// { +// fx1[negidx+i*nf+db] = conj(fx[db]); +// fx1[negidx+i*nf+ub] = conj(fx[ub]); +// fx1[negidx+i*nf+u ] = conj(fx[u ]); +// fx1[negidx+i*nf+d ] = conj(fx[d ]); +// } +// else if (opts.ih1 == -1) +// { +// fx1[negidx+i*nf+db] = conj(fx[d ]); +// fx1[negidx+i*nf+ub] = conj(fx[u ]); +// fx1[negidx+i*nf+u ] = conj(fx[ub]); +// fx1[negidx+i*nf+d ] = conj(fx[db]); +// } +// +// //beam 2 positive +// fx2[i*nf+bb] = fx[bb]; +// fx2[i*nf+cb] = fx[cb]; +// fx2[i*nf+sb] = fx[sb]; +// fx2[i*nf+g ] = fx[g ]; +// fx2[i*nf+s ] = fx[s ]; +// fx2[i*nf+c ] = fx[c ]; +// fx2[i*nf+b ] = fx[b ]; +// if (opts.ih2 == 1) +// { +// fx2[i*nf+db] = fx[db]; +// fx2[i*nf+ub] = fx[ub]; +// fx2[i*nf+u ] = fx[u ]; +// fx2[i*nf+d ] = fx[d ]; +// } +// else if (opts.ih2 == -1) +// { +// fx2[i*nf+db] = fx[d ]; +// fx2[i*nf+ub] = fx[u ]; +// fx2[i*nf+u ] = fx[ub]; +// fx2[i*nf+d ] = fx[db]; +// } + +// if (opts.mellin1d && opts.bprescription != 2) return; + //beam 2 negative fx2[negidx+i*nf+bb] = conj(fx[bb]); fx2[negidx+i*nf+cb] = conj(fx[cb]); fx2[negidx+i*nf+sb] = conj(fx[sb]); fx2[negidx+i*nf+g ] = conj(fx[g ]); fx2[negidx+i*nf+s ] = conj(fx[s ]); fx2[negidx+i*nf+c ] = conj(fx[c ]); fx2[negidx+i*nf+b ] = conj(fx[b ]); if (opts.ih2 == 1) { fx2[negidx+i*nf+db] = conj(fx[db]); fx2[negidx+i*nf+ub] = conj(fx[ub]); fx2[negidx+i*nf+u ] = conj(fx[u ]); fx2[negidx+i*nf+d ] = conj(fx[d ]); } else if (opts.ih2 == -1) { fx2[negidx+i*nf+db] = conj(fx[d ]); fx2[negidx+i*nf+ub] = conj(fx[u ]); fx2[negidx+i*nf+u ] = conj(fx[ub]); fx2[negidx+i*nf+d ] = conj(fx[db]); } } void pdfevol::storemoments_fortran(int i, complex fx[11]) { //Save the evolved PDFs into the fortran common block //beam 1 creno_.cfx1_[i][-5+MAXNF] = fcx(fx[-5+MAXNF]); creno_.cfx1_[i][-4+MAXNF] = fcx(fx[-4+MAXNF]); creno_.cfx1_[i][-3+MAXNF] = fcx(fx[-3+MAXNF]); creno_.cfx1_[i][ 0+MAXNF] = fcx(fx[ 0+MAXNF]); creno_.cfx1_[i][ 3+MAXNF] = fcx(fx[ 3+MAXNF]); creno_.cfx1_[i][ 4+MAXNF] = fcx(fx[ 4+MAXNF]); creno_.cfx1_[i][ 5+MAXNF] = fcx(fx[ 5+MAXNF]); if (opts.ih1 == 1) { creno_.cfx1_[i][-2+MAXNF] = fcx(fx[-2+MAXNF]); creno_.cfx1_[i][-1+MAXNF] = fcx(fx[-1+MAXNF]); creno_.cfx1_[i][ 1+MAXNF] = fcx(fx[ 1+MAXNF]); creno_.cfx1_[i][ 2+MAXNF] = fcx(fx[ 2+MAXNF]); } else if (opts.ih1 == -1) { creno_.cfx1_[i][-2+MAXNF] = fcx(fx[ 2+MAXNF]); creno_.cfx1_[i][-1+MAXNF] = fcx(fx[ 1+MAXNF]); creno_.cfx1_[i][ 1+MAXNF] = fcx(fx[-1+MAXNF]); creno_.cfx1_[i][ 2+MAXNF] = fcx(fx[-2+MAXNF]); } //beam 2 positive creno_.cfx2p_[i][-5+MAXNF] = fcx(fx[-5+MAXNF]); creno_.cfx2p_[i][-4+MAXNF] = fcx(fx[-4+MAXNF]); creno_.cfx2p_[i][-3+MAXNF] = fcx(fx[-3+MAXNF]); creno_.cfx2p_[i][ 0+MAXNF] = fcx(fx[ 0+MAXNF]); creno_.cfx2p_[i][ 3+MAXNF] = fcx(fx[ 3+MAXNF]); creno_.cfx2p_[i][ 4+MAXNF] = fcx(fx[ 4+MAXNF]); creno_.cfx2p_[i][ 5+MAXNF] = fcx(fx[ 5+MAXNF]); if (opts.ih2 == 1) { creno_.cfx2p_[i][-2+MAXNF] = fcx(fx[-2+MAXNF]); creno_.cfx2p_[i][-1+MAXNF] = fcx(fx[-1+MAXNF]); creno_.cfx2p_[i][ 1+MAXNF] = fcx(fx[ 1+MAXNF]); creno_.cfx2p_[i][ 2+MAXNF] = fcx(fx[ 2+MAXNF]); } else if (opts.ih2 == -1) { creno_.cfx2p_[i][-2+MAXNF] = fcx(fx[ 2+MAXNF]); creno_.cfx2p_[i][-1+MAXNF] = fcx(fx[ 1+MAXNF]); creno_.cfx2p_[i][ 1+MAXNF] = fcx(fx[-1+MAXNF]); creno_.cfx2p_[i][ 2+MAXNF] = fcx(fx[-2+MAXNF]); } //beam 2 negative creno_.cfx2m_[i][-5+MAXNF] = fcx(conj(fx[-5+MAXNF])); creno_.cfx2m_[i][-4+MAXNF] = fcx(conj(fx[-4+MAXNF])); creno_.cfx2m_[i][-3+MAXNF] = fcx(conj(fx[-3+MAXNF])); creno_.cfx2m_[i][ 0+MAXNF] = fcx(conj(fx[ 0+MAXNF])); creno_.cfx2m_[i][ 3+MAXNF] = fcx(conj(fx[ 3+MAXNF])); creno_.cfx2m_[i][ 4+MAXNF] = fcx(conj(fx[ 4+MAXNF])); creno_.cfx2m_[i][ 5+MAXNF] = fcx(conj(fx[ 5+MAXNF])); if (opts.ih2 == 1) { creno_.cfx2m_[i][-2+MAXNF] = fcx(conj(fx[-2+MAXNF])); creno_.cfx2m_[i][-1+MAXNF] = fcx(conj(fx[-1+MAXNF])); creno_.cfx2m_[i][ 1+MAXNF] = fcx(conj(fx[ 1+MAXNF])); creno_.cfx2m_[i][ 2+MAXNF] = fcx(conj(fx[ 2+MAXNF])); } else if (opts.ih2 == -1) { creno_.cfx2m_[i][-2+MAXNF] = fcx(conj(fx[ 2+MAXNF])); creno_.cfx2m_[i][-1+MAXNF] = fcx(conj(fx[ 1+MAXNF])); creno_.cfx2m_[i][ 1+MAXNF] = fcx(conj(fx[-1+MAXNF])); creno_.cfx2m_[i][ 2+MAXNF] = fcx(conj(fx[-2+MAXNF])); } } //Apply the flavour dependent form factors void pdfevol::flavour_kt() { - int negidx = mellinint::mdim*11; int nf = 2*MAXNF+1; + int negidx = mellinint::mdim*nf; complex uv, dv, us, ds; for (int i = 0; i < mellinint::mdim; i++) { //beam 1 positive fx1[i*nf+bb] *= npff::boff; fx1[i*nf+cb] *= npff::chff; fx1[i*nf+sb] *= npff::ssff; fx1[i*nf+g ] *= npff::glff; fx1[i*nf+s ] *= npff::ssff; fx1[i*nf+c ] *= npff::chff; fx1[i*nf+b ] *= npff::boff; if (opts.ih1 == 1) { uv = (fx1[i*nf+u] - fx1[i*nf+ub]); dv = (fx1[i*nf+d] - fx1[i*nf+db]); us = fx1[i*nf+ub]; ds = fx1[i*nf+db]; fx1[i*nf+u] = uv * npff::uvff + us * npff::usff; fx1[i*nf+d] = dv * npff::dvff + ds * npff::dsff; fx1[i*nf+ub] = us * npff::usff; fx1[i*nf+db] = ds * npff::dsff; } else if (opts.ih1 == -1) { uv = (fx1[i*nf+ub] - fx1[i*nf+u]); dv = (fx1[i*nf+db] - fx1[i*nf+d]); us = fx1[i*nf+u]; ds = fx1[i*nf+d]; fx1[i*nf+ub] = uv * npff::uvff + us * npff::usff; fx1[i*nf+db] = dv * npff::dvff + ds * npff::dsff; fx1[i*nf+u] = us * npff::usff; fx1[i*nf+d] = ds * npff::dsff; } //beam 2 positive fx2[i*nf+bb] *= npff::boff; fx2[i*nf+cb] *= npff::chff; fx2[i*nf+sb] *= npff::ssff; fx2[i*nf+g ] *= npff::glff; fx2[i*nf+s ] *= npff::ssff; fx2[i*nf+c ] *= npff::chff; fx2[i*nf+b ] *= npff::boff; if (opts.ih2 == 1) { uv = (fx2[i*nf+u] - fx2[i*nf+ub]); dv = (fx2[i*nf+d] - fx2[i*nf+db]); us = fx2[i*nf+ub]; ds = fx2[i*nf+db]; fx2[i*nf+u] = uv * npff::uvff + us * npff::usff; fx2[i*nf+d] = dv * npff::dvff + ds * npff::dsff; fx2[i*nf+ub] = us * npff::usff; fx2[i*nf+db] = ds * npff::dsff; } else if (opts.ih2 == -1) { uv = (fx2[i*nf+ub] - fx2[i*nf+u]); dv = (fx2[i*nf+db] - fx2[i*nf+d]); us = fx2[i*nf+u]; ds = fx2[i*nf+d]; fx2[i*nf+ub] = uv * npff::uvff + us * npff::usff; fx2[i*nf+db] = dv * npff::dvff + ds * npff::dsff; fx2[i*nf+u] = us * npff::usff; fx2[i*nf+d] = ds * npff::dsff; } + + //beam 2 negative + fx2[negidx+i*nf+bb] *= npff::boff; + fx2[negidx+i*nf+cb] *= npff::chff; + fx2[negidx+i*nf+sb] *= npff::ssff; + fx2[negidx+i*nf+g ] *= npff::glff; + fx2[negidx+i*nf+s ] *= npff::ssff; + fx2[negidx+i*nf+c ] *= npff::chff; + fx2[negidx+i*nf+b ] *= npff::boff; + if (opts.ih2 == 1) + { + uv = (fx2[negidx+i*nf+u] - fx2[negidx+i*nf+ub]); + dv = (fx2[negidx+i*nf+d] - fx2[negidx+i*nf+db]); + us = fx2[negidx+i*nf+ub]; + ds = fx2[negidx+i*nf+db]; + fx2[negidx+i*nf+u] = uv * npff::uvff + us * npff::usff; + fx2[negidx+i*nf+d] = dv * npff::dvff + ds * npff::dsff; + fx2[negidx+i*nf+ub] = us * npff::usff; + fx2[negidx+i*nf+db] = ds * npff::dsff; + } + else if (opts.ih2 == -1) + { + uv = (fx2[negidx+i*nf+ub] - fx2[negidx+i*nf+u]); + dv = (fx2[negidx+i*nf+db] - fx2[negidx+i*nf+d]); + us = fx2[negidx+i*nf+u]; + ds = fx2[negidx+i*nf+d]; + fx2[negidx+i*nf+ub] = uv * npff::uvff + us * npff::usff; + fx2[negidx+i*nf+db] = dv * npff::dvff + ds * npff::dsff; + fx2[negidx+i*nf+u] = us * npff::usff; + fx2[negidx+i*nf+d] = ds * npff::dsff; + } } } void pdfevol::retrieve(int i1, int i2, int sign) { int nf = 2*MAXNF+1; int negidx = mellinint::mdim*nf; fn1[bb] = fx1[i1*nf+bb]; fn1[cb] = fx1[i1*nf+cb]; fn1[sb] = fx1[i1*nf+sb]; fn1[db] = fx1[i1*nf+db]; fn1[ub] = fx1[i1*nf+ub]; fn1[g ] = fx1[i1*nf+g ]; fn1[u ] = fx1[i1*nf+u ]; fn1[d ] = fx1[i1*nf+d ]; fn1[s ] = fx1[i1*nf+s ]; fn1[c ] = fx1[i1*nf+c ]; fn1[b ] = fx1[i1*nf+b ]; if (sign == mesq::positive) { fn2[bb] = fx2[i2*nf+bb]; fn2[cb] = fx2[i2*nf+cb]; fn2[sb] = fx2[i2*nf+sb]; fn2[db] = fx2[i2*nf+db]; fn2[ub] = fx2[i2*nf+ub]; fn2[g ] = fx2[i2*nf+g ]; fn2[u ] = fx2[i2*nf+u ]; fn2[d ] = fx2[i2*nf+d ]; fn2[s ] = fx2[i2*nf+s ]; fn2[c ] = fx2[i2*nf+c ]; fn2[b ] = fx2[i2*nf+b ]; } else if (sign == mesq::negative) { fn2[bb] = fx2[negidx+i2*nf+bb]; fn2[cb] = fx2[negidx+i2*nf+cb]; fn2[sb] = fx2[negidx+i2*nf+sb]; fn2[db] = fx2[negidx+i2*nf+db]; fn2[ub] = fx2[negidx+i2*nf+ub]; fn2[g ] = fx2[negidx+i2*nf+g ]; fn2[u ] = fx2[negidx+i2*nf+u ]; fn2[d ] = fx2[negidx+i2*nf+d ]; fn2[s ] = fx2[negidx+i2*nf+s ]; fn2[c ] = fx2[negidx+i2*nf+c ]; fn2[b ] = fx2[negidx+i2*nf+b ]; } } void pdfevol::retrieve_fortran(int i1, int i2, int sign) { // cout << i1 << endl; // cout << creno_.cfx1_[i1][5].real << " " << creno_.cfx1_[i1][5].imag << endl; fn1[-5+MAXNF] = cx(creno_.cfx1_[i1][-5+MAXNF]); fn1[-4+MAXNF] = cx(creno_.cfx1_[i1][-4+MAXNF]); fn1[-3+MAXNF] = cx(creno_.cfx1_[i1][-3+MAXNF]); fn1[-2+MAXNF] = cx(creno_.cfx1_[i1][-2+MAXNF]); fn1[-1+MAXNF] = cx(creno_.cfx1_[i1][-1+MAXNF]); fn1[ 0+MAXNF] = cx(creno_.cfx1_[i1][ 0+MAXNF]); fn1[ 1+MAXNF] = cx(creno_.cfx1_[i1][ 1+MAXNF]); fn1[ 2+MAXNF] = cx(creno_.cfx1_[i1][ 2+MAXNF]); fn1[ 3+MAXNF] = cx(creno_.cfx1_[i1][ 3+MAXNF]); fn1[ 4+MAXNF] = cx(creno_.cfx1_[i1][ 4+MAXNF]); fn1[ 5+MAXNF] = cx(creno_.cfx1_[i1][ 5+MAXNF]); if (sign == mesq::positive) { fn2[-5+MAXNF] = cx(creno_.cfx2p_[i2][-5+MAXNF]); fn2[-4+MAXNF] = cx(creno_.cfx2p_[i2][-4+MAXNF]); fn2[-3+MAXNF] = cx(creno_.cfx2p_[i2][-3+MAXNF]); fn2[-2+MAXNF] = cx(creno_.cfx2p_[i2][-2+MAXNF]); fn2[-1+MAXNF] = cx(creno_.cfx2p_[i2][-1+MAXNF]); fn2[ 0+MAXNF] = cx(creno_.cfx2p_[i2][ 0+MAXNF]); fn2[ 1+MAXNF] = cx(creno_.cfx2p_[i2][ 1+MAXNF]); fn2[ 2+MAXNF] = cx(creno_.cfx2p_[i2][ 2+MAXNF]); fn2[ 3+MAXNF] = cx(creno_.cfx2p_[i2][ 3+MAXNF]); fn2[ 4+MAXNF] = cx(creno_.cfx2p_[i2][ 4+MAXNF]); fn2[ 5+MAXNF] = cx(creno_.cfx2p_[i2][ 5+MAXNF]); } else if (sign == mesq::negative) { fn2[-5+MAXNF] = cx(creno_.cfx2m_[i2][-5+MAXNF]); fn2[-4+MAXNF] = cx(creno_.cfx2m_[i2][-4+MAXNF]); fn2[-3+MAXNF] = cx(creno_.cfx2m_[i2][-3+MAXNF]); fn2[-2+MAXNF] = cx(creno_.cfx2m_[i2][-2+MAXNF]); fn2[-1+MAXNF] = cx(creno_.cfx2m_[i2][-1+MAXNF]); fn2[ 0+MAXNF] = cx(creno_.cfx2m_[i2][ 0+MAXNF]); fn2[ 1+MAXNF] = cx(creno_.cfx2m_[i2][ 1+MAXNF]); fn2[ 2+MAXNF] = cx(creno_.cfx2m_[i2][ 2+MAXNF]); fn2[ 3+MAXNF] = cx(creno_.cfx2m_[i2][ 3+MAXNF]); fn2[ 4+MAXNF] = cx(creno_.cfx2m_[i2][ 4+MAXNF]); fn2[ 5+MAXNF] = cx(creno_.cfx2m_[i2][ 5+MAXNF]); } } void pdfevol::retrieve_beam1(int i1) { int nf = 2*MAXNF+1; - memcpy(fn1, &(fx1[i1*nf]), nf*sizeof(complex)); + //memcpy(fn1, &(fx1[i1*nf]), nf*sizeof(complex)); + copy(fx1+i1*nf, fx1+(i1+1)*nf, fn1); } void pdfevol::retrieve_beam2_pos(int i2) { int nf = 2*MAXNF+1; - memcpy(fn2, &(fx2[i2*nf]), nf*sizeof(complex)); + //memcpy(fn2, &(fx2[i2*nf]), nf*sizeof(complex)); + copy(fx2+i2*nf, fx2+(i2+1)*nf, fn2); } void pdfevol::retrieve_beam2_neg() { - fn2[bb] = conj(fn2[bb]); fn2[cb] = conj(fn2[cb]); fn2[sb] = conj(fn2[sb]); fn2[db] = conj(fn2[db]); fn2[ub] = conj(fn2[ub]); fn2[g ] = conj(fn2[g ]); fn2[u ] = conj(fn2[u ]); fn2[d ] = conj(fn2[d ]); fn2[s ] = conj(fn2[s ]); fn2[c ] = conj(fn2[c ]); fn2[b ] = conj(fn2[b ]); } -void pdfevol::retrieve1d(int i, int sign) +void pdfevol::retrieve1d_pos(int i) { int nf = 2*MAXNF+1; - int negidx = mellinint::mdim*nf; - if (sign == mesq::positive) - { - fn1[bb] = fx1[i*nf+bb]; - fn1[cb] = fx1[i*nf+cb]; - fn1[sb] = fx1[i*nf+sb]; - fn1[db] = fx1[i*nf+db]; - fn1[ub] = fx1[i*nf+ub]; - fn1[g ] = fx1[i*nf+g ]; - fn1[u ] = fx1[i*nf+u ]; - fn1[d ] = fx1[i*nf+d ]; - fn1[s ] = fx1[i*nf+s ]; - fn1[c ] = fx1[i*nf+c ]; - fn1[b ] = fx1[i*nf+b ]; - } - else if (sign == mesq::negative) - { - fn1[bb] = fx1[negidx+i*nf+bb]; - fn1[cb] = fx1[negidx+i*nf+cb]; - fn1[sb] = fx1[negidx+i*nf+sb]; - fn1[db] = fx1[negidx+i*nf+db]; - fn1[ub] = fx1[negidx+i*nf+ub]; - fn1[g ] = fx1[negidx+i*nf+g ]; - fn1[u ] = fx1[negidx+i*nf+u ]; - fn1[d ] = fx1[negidx+i*nf+d ]; - fn1[s ] = fx1[negidx+i*nf+s ]; - fn1[c ] = fx1[negidx+i*nf+c ]; - fn1[b ] = fx1[negidx+i*nf+b ]; - } - if (sign == mesq::positive) - { - fn2[bb] = fx2[i*nf+bb]; - fn2[cb] = fx2[i*nf+cb]; - fn2[sb] = fx2[i*nf+sb]; - fn2[db] = fx2[i*nf+db]; - fn2[ub] = fx2[i*nf+ub]; - fn2[g ] = fx2[i*nf+g ]; - fn2[u ] = fx2[i*nf+u ]; - fn2[d ] = fx2[i*nf+d ]; - fn2[s ] = fx2[i*nf+s ]; - fn2[c ] = fx2[i*nf+c ]; - fn2[b ] = fx2[i*nf+b ]; - } - else if (sign == mesq::negative) - { - fn2[bb] = fx2[negidx+i*nf+bb]; - fn2[cb] = fx2[negidx+i*nf+cb]; - fn2[sb] = fx2[negidx+i*nf+sb]; - fn2[db] = fx2[negidx+i*nf+db]; - fn2[ub] = fx2[negidx+i*nf+ub]; - fn2[g ] = fx2[negidx+i*nf+g ]; - fn2[u ] = fx2[negidx+i*nf+u ]; - fn2[d ] = fx2[negidx+i*nf+d ]; - fn2[s ] = fx2[negidx+i*nf+s ]; - fn2[c ] = fx2[negidx+i*nf+c ]; - fn2[b ] = fx2[negidx+i*nf+b ]; - } + //memcpy(fn1, &(fx1[i*nf]), nf*sizeof(complex)); + //memcpy(fn2, &(fx2[i*nf]), nf*sizeof(complex)); + + copy(fx1+i*nf, fx1+(i+1)*nf, fn1); + copy(fx2+i*nf, fx2+(i+1)*nf, fn2); } - + +void pdfevol::retrieve1d_neg() +{ + //the negative branch is needed in mellin1d only for bprescription = 2 + fn1[bb] = conj(fn1[bb]); + fn1[cb] = conj(fn1[cb]); + fn1[sb] = conj(fn1[sb]); + fn1[db] = conj(fn1[db]); + fn1[ub] = conj(fn1[ub]); + fn1[g ] = conj(fn1[g ]); + fn1[u ] = conj(fn1[u ]); + fn1[d ] = conj(fn1[d ]); + fn1[s ] = conj(fn1[s ]); + fn1[c ] = conj(fn1[c ]); + fn1[b ] = conj(fn1[b ]); + + fn2[bb] = conj(fn2[bb]); + fn2[cb] = conj(fn2[cb]); + fn2[sb] = conj(fn2[sb]); + fn2[db] = conj(fn2[db]); + fn2[ub] = conj(fn2[ub]); + fn2[g ] = conj(fn2[g ]); + fn2[u ] = conj(fn2[u ]); + fn2[d ] = conj(fn2[d ]); + fn2[s ] = conj(fn2[s ]); + fn2[c ] = conj(fn2[c ]); + fn2[b ] = conj(fn2[b ]); + + // int negidx = mellinint::mdim*nf; + // copy(fx1+negidx+i*nf, fx1+negidx+(i+1)*nf, fn1); + // copy(fx2+negidx+i*nf, fx2+negidx+(i+1)*nf, fn2); +} + void pdfevol::retrieve1d_fortran(int i, int sign) { //cout << i << endl; //cout << creno_.cfx1_[i][5].real << " " << creno_.cfx1_[i][5].imag << endl; if (sign == mesq::positive) { fn1[-5+MAXNF] = cx(creno_.cfx1_[i][-5+MAXNF]); fn1[-4+MAXNF] = cx(creno_.cfx1_[i][-4+MAXNF]); fn1[-3+MAXNF] = cx(creno_.cfx1_[i][-3+MAXNF]); fn1[-2+MAXNF] = cx(creno_.cfx1_[i][-2+MAXNF]); fn1[-1+MAXNF] = cx(creno_.cfx1_[i][-1+MAXNF]); fn1[ 0+MAXNF] = cx(creno_.cfx1_[i][ 0+MAXNF]); fn1[ 1+MAXNF] = cx(creno_.cfx1_[i][ 1+MAXNF]); fn1[ 2+MAXNF] = cx(creno_.cfx1_[i][ 2+MAXNF]); fn1[ 3+MAXNF] = cx(creno_.cfx1_[i][ 3+MAXNF]); fn1[ 4+MAXNF] = cx(creno_.cfx1_[i][ 4+MAXNF]); fn1[ 5+MAXNF] = cx(creno_.cfx1_[i][ 5+MAXNF]); } else if (sign == mesq::negative) { fn1[-5+MAXNF] = conj(cx(creno_.cfx1_[i][-5+MAXNF])); fn1[-4+MAXNF] = conj(cx(creno_.cfx1_[i][-4+MAXNF])); fn1[-3+MAXNF] = conj(cx(creno_.cfx1_[i][-3+MAXNF])); fn1[-2+MAXNF] = conj(cx(creno_.cfx1_[i][-2+MAXNF])); fn1[-1+MAXNF] = conj(cx(creno_.cfx1_[i][-1+MAXNF])); fn1[ 0+MAXNF] = conj(cx(creno_.cfx1_[i][ 0+MAXNF])); fn1[ 1+MAXNF] = conj(cx(creno_.cfx1_[i][ 1+MAXNF])); fn1[ 2+MAXNF] = conj(cx(creno_.cfx1_[i][ 2+MAXNF])); fn1[ 3+MAXNF] = conj(cx(creno_.cfx1_[i][ 3+MAXNF])); fn1[ 4+MAXNF] = conj(cx(creno_.cfx1_[i][ 4+MAXNF])); fn1[ 5+MAXNF] = conj(cx(creno_.cfx1_[i][ 5+MAXNF])); } if (sign == mesq::positive) { fn2[-5+MAXNF] = cx(creno_.cfx2p_[i][-5+MAXNF]); fn2[-4+MAXNF] = cx(creno_.cfx2p_[i][-4+MAXNF]); fn2[-3+MAXNF] = cx(creno_.cfx2p_[i][-3+MAXNF]); fn2[-2+MAXNF] = cx(creno_.cfx2p_[i][-2+MAXNF]); fn2[-1+MAXNF] = cx(creno_.cfx2p_[i][-1+MAXNF]); fn2[ 0+MAXNF] = cx(creno_.cfx2p_[i][ 0+MAXNF]); fn2[ 1+MAXNF] = cx(creno_.cfx2p_[i][ 1+MAXNF]); fn2[ 2+MAXNF] = cx(creno_.cfx2p_[i][ 2+MAXNF]); fn2[ 3+MAXNF] = cx(creno_.cfx2p_[i][ 3+MAXNF]); fn2[ 4+MAXNF] = cx(creno_.cfx2p_[i][ 4+MAXNF]); fn2[ 5+MAXNF] = cx(creno_.cfx2p_[i][ 5+MAXNF]); } else if (sign == mesq::negative) { fn2[-5+MAXNF] = cx(creno_.cfx2m_[i][-5+MAXNF]); fn2[-4+MAXNF] = cx(creno_.cfx2m_[i][-4+MAXNF]); fn2[-3+MAXNF] = cx(creno_.cfx2m_[i][-3+MAXNF]); fn2[-2+MAXNF] = cx(creno_.cfx2m_[i][-2+MAXNF]); fn2[-1+MAXNF] = cx(creno_.cfx2m_[i][-1+MAXNF]); fn2[ 0+MAXNF] = cx(creno_.cfx2m_[i][ 0+MAXNF]); fn2[ 1+MAXNF] = cx(creno_.cfx2m_[i][ 1+MAXNF]); fn2[ 2+MAXNF] = cx(creno_.cfx2m_[i][ 2+MAXNF]); fn2[ 3+MAXNF] = cx(creno_.cfx2m_[i][ 3+MAXNF]); fn2[ 4+MAXNF] = cx(creno_.cfx2m_[i][ 4+MAXNF]); fn2[ 5+MAXNF] = cx(creno_.cfx2m_[i][ 5+MAXNF]); } } -//Retrieve PDFs at the starting scale (muf) -void pdfevol::retrievemuf(int i, int sign) +//Retrieve PDFs at the starting scale (muf), only positive branch +void pdfevol::retrievemuf(int i) { // i is the index of the complex mellin moment in the z-space for the gaussian quadrature used for the mellin inversion //N flavour dependence int nf = resconst::NF; //XP[i] are moments of PDFs at the starting scale (factorisation scale) complex fx[11]; fx[0+MAXNF] = GLP[i]; fx[1+MAXNF] = UVP[i] + USP[i]; fx[-1+MAXNF] = USP[i]; fx[2+MAXNF] = DVP[i] + DSP[i]; fx[-2+MAXNF] = DSP[i]; fx[3+MAXNF] = SSP[i]; fx[-3+MAXNF] = SSP[i]; if (nf >= 4) { fx[4+MAXNF] = CHP[i]; fx[-4+MAXNF] = CHP[i]; } else { fx[4+MAXNF] = 0.; fx[-4+MAXNF] = 0.; } if (nf >= 5) { fx[5+MAXNF] = BOP[i]; fx[-5+MAXNF] = BOP[i]; } else { fx[5+MAXNF] = 0.; fx[-5+MAXNF] = 0.; } storemoments(i, fx); - retrieve1d(i, sign); + retrieve1d_pos(i); // cout << i << " " << GLP[i] << " " << fx[0+MAXNF] << " " << fn1[MAXNF] << " " << fn2[MAXNF] << endl; return; } void pdfevol::truncate() { //Calculate truncated moments double x1 = phasespace::m/opts.sroot*exp(phasespace::ymin); double x2 = phasespace::m/opts.sroot*exp(-phasespace::ymax); //double x1 = 1e-8;//pow(phasespace::m/opts.sroot,2); //double x2 = 1e-8;//pow(phasespace::m/opts.sroot,2); double lx1 = log(x1); double lx2 = log(x2); //truncated moments (output) complex fx1[mellinint::mdim][2*MAXNF+1] = {0.}; complex fx2[mellinint::mdim][2*MAXNF+1] = {0.}; //cache x^(N) values complex x1n[mellinint::mdim]; complex x2n[mellinint::mdim]; for (int n = 0; n < mellinint::mdim; n++) { x1n[n] = pow(x1,mellinint::Np[n]); x2n[n] = pow(x2,mellinint::Np[n]); } //Normalisation times Jacobian complex facp = mellinint::CCp/2./M_PI/complex (0.,1); complex facm = mellinint::CCm/2./M_PI/complex (0.,1); //original moments times prefactor and weight complex fm1p[mellinint::mdim][2*MAXNF+1]; complex fm2p[mellinint::mdim][2*MAXNF+1]; complex fm1m[mellinint::mdim][2*MAXNF+1]; complex fm2m[mellinint::mdim][2*MAXNF+1]; for (int m = 0; m < mellinint::mdim; m++) for (int f = 0; f < 2*MAXNF+1; f++) { fm1p[m][f] = facp * cx(creno_.cfx1_[m][f] ) * mellinint::wn[m]; fm2p[m][f] = facp * cx(creno_.cfx2p_[m][f]) * mellinint::wn[m]; fm1m[m][f] = facm * conj(cx(creno_.cfx1_[m][f]) ) * mellinint::wn[m]; fm2m[m][f] = facm * conj(cx(creno_.cfx2p_[m][f])) * mellinint::wn[m]; } //cache factor (1-x^(N-M))/(N-M) which limit is ln(x) when N-M -> 0 complex llx1p[mellinint::mdim][mellinint::mdim]; complex llx2p[mellinint::mdim][mellinint::mdim]; complex llx1m[mellinint::mdim][mellinint::mdim]; complex llx2m[mellinint::mdim][mellinint::mdim]; for (int n = 0; n < mellinint::mdim; n++) for (int m = 0; m < mellinint::mdim; m++) { llx1p[n][m] = (1.-x1n[n]/x1n[m])/(mellinint::Np[n]-mellinint::Np[m]); llx2p[n][m] = (1.-x2n[n]/x2n[m])/(mellinint::Np[n]-mellinint::Np[m]); llx1m[n][m] = (1.-x1n[n]/conj(x1n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); llx2m[n][m] = (1.-x2n[n]/conj(x2n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); } //overwrite divergent diagonal part for (int n = 0; n < mellinint::mdim; n++) { llx1p[n][n] = -lx1; llx2p[n][n] = -lx2; } for (int n = 0; n < mellinint::mdim; n++) { //positive branch for (int m = 0; m < mellinint::mdim; m++) for (int f = 0; f < 2*MAXNF+1; f++) { /* if (m == n) { fx1[n][f] += fm1p[m][f] * (-lx1); fx2[n][f] += fm2p[m][f] * (-lx2); } else { fx1[n][f] += fm1p[m][f] * (1.-x1n[n]/x1n[m])/(mellinint::Np[n]-mellinint::Np[m]); fx2[n][f] += fm2p[m][f] * (1.-x2n[n]/x2n[m])/(mellinint::Np[n]-mellinint::Np[m]); } */ // fx1[n][f] += fm1p[m][f]*llx1p[n][m]; // fx2[n][f] += fm2p[m][f]*llx2p[n][m]; fx1[n][f] += fm1p[m][f]*llx1p[n][m] - fm1m[m][f]*llx1m[n][m]; fx2[n][f] += fm2p[m][f]*llx2p[n][m] - fm2m[m][f]*llx2m[n][m]; } /* //negative branch for (int m = 0; m < mellinint::mdim; m++) for (int f = 0; f < 2*MAXNF+1; f++) { // fx1[n][f] -= fm1m[m][f]*(1.-x1n[n]/conj(x1n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); // fx2[n][f] -= fm2m[m][f]*(1.-x2n[n]/conj(x2n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); fx1[n][f] -= fm1m[m][f]*llx1m[n][m]; fx2[n][f] -= fm2m[m][f]*llx2m[n][m]; } */ //cout << "truncated " << n << fx1[5] << endl; } //replace moments for (int n = 0; n < mellinint::mdim; n++) storemoments(n, fx1[n]); //storemoments(n, fx1[n], fx2[n]); } void pdfevol::uppertruncate() { //Calculate truncated moments double xmin1 = pow(phasespace::m/opts.sroot,2); double xmin2 = pow(phasespace::m/opts.sroot,2); double xmax1 = phasespace::m/opts.sroot*exp(phasespace::ymax); double xmax2 = phasespace::m/opts.sroot*exp(-phasespace::ymin); double lx1 = log(xmin1/xmax1); double lx2 = log(xmin2/xmax2); //truncated moments (output) complex fx1[mellinint::mdim][2*MAXNF+1] = {0.}; complex fx2[mellinint::mdim][2*MAXNF+1] = {0.}; //cache x^(N) values complex xmin1n[mellinint::mdim]; complex xmin2n[mellinint::mdim]; complex xmax1n[mellinint::mdim]; complex xmax2n[mellinint::mdim]; for (int n = 0; n < mellinint::mdim; n++) { xmin1n[n] = pow(xmin1,mellinint::Np[n]); xmin2n[n] = pow(xmin2,mellinint::Np[n]); xmax1n[n] = pow(xmax1,mellinint::Np[n]); xmax2n[n] = pow(xmax2,mellinint::Np[n]); } //Normalisation times Jacobian complex facp = mellinint::CCp/2./M_PI/complex (0.,1); complex facm = mellinint::CCm/2./M_PI/complex (0.,1); //original moments times prefactor and weight complex fm1p[mellinint::mdim][2*MAXNF+1]; complex fm2p[mellinint::mdim][2*MAXNF+1]; complex fm1m[mellinint::mdim][2*MAXNF+1]; complex fm2m[mellinint::mdim][2*MAXNF+1]; for (int m = 0; m < mellinint::mdim; m++) for (int f = 0; f < 2*MAXNF+1; f++) { fm1p[m][f] = facp * cx(creno_.cfx1_[m][f] ) * mellinint::wn[m]; fm2p[m][f] = facp * cx(creno_.cfx2p_[m][f]) * mellinint::wn[m]; fm1m[m][f] = facm * conj(cx(creno_.cfx1_[m][f]) ) * mellinint::wn[m]; fm2m[m][f] = facm * conj(cx(creno_.cfx2p_[m][f])) * mellinint::wn[m]; } //cache factor (xmax^(N-M)-xmin^(N-M))/(N-M) which limit is ln(xmin/xmax) when N-M -> 0 complex llx1p[mellinint::mdim][mellinint::mdim]; complex llx2p[mellinint::mdim][mellinint::mdim]; complex llx1m[mellinint::mdim][mellinint::mdim]; complex llx2m[mellinint::mdim][mellinint::mdim]; for (int n = 0; n < mellinint::mdim; n++) for (int m = 0; m < mellinint::mdim; m++) { llx1p[n][m] = (xmax1n[n]/xmax1n[m]-xmin1n[n]/xmin1n[m])/(mellinint::Np[n]-mellinint::Np[m]); llx2p[n][m] = (xmax2n[n]/xmax2n[m]-xmin2n[n]/xmin2n[m])/(mellinint::Np[n]-mellinint::Np[m]); llx1m[n][m] = (xmax1n[n]/conj(xmax1n[m])-xmin1n[n]/conj(xmin1n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); llx2m[n][m] = (xmax2n[n]/conj(xmax2n[m])-xmin2n[n]/conj(xmin2n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); } //overwrite divergent diagonal part for (int n = 0; n < mellinint::mdim; n++) { llx1p[n][n] = -lx1; llx2p[n][n] = -lx2; } for (int n = 0; n < mellinint::mdim; n++) { //positive branch for (int m = 0; m < mellinint::mdim; m++) for (int f = 0; f < 2*MAXNF+1; f++) { /* if (m == n) { fx1[n][f] += fm1p[m][f] * (-lx1); fx2[n][f] += fm2p[m][f] * (-lx2); } else { fx1[n][f] += fm1p[m][f] * (1.-x1n[n]/x1n[m])/(mellinint::Np[n]-mellinint::Np[m]); fx2[n][f] += fm2p[m][f] * (1.-x2n[n]/x2n[m])/(mellinint::Np[n]-mellinint::Np[m]); } */ // fx1[n][f] += fm1p[m][f]*llx1p[n][m]; // fx2[n][f] += fm2p[m][f]*llx2p[n][m]; fx1[n][f] += fm1p[m][f]*llx1p[n][m] - fm1m[m][f]*llx1m[n][m]; fx2[n][f] += fm2p[m][f]*llx2p[n][m] - fm2m[m][f]*llx2m[n][m]; } /* //negative branch for (int m = 0; m < mellinint::mdim; m++) for (int f = 0; f < 2*MAXNF+1; f++) { // fx1[n][f] -= fm1m[m][f]*(1.-x1n[n]/conj(x1n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); // fx2[n][f] -= fm2m[m][f]*(1.-x2n[n]/conj(x2n[m]))/(mellinint::Np[n]-mellinint::Nm[m]); fx1[n][f] -= fm1m[m][f]*llx1m[n][m]; fx2[n][f] -= fm2m[m][f]*llx2m[n][m]; } */ //cout << "truncated " << n << fx1[5] << endl; } //replace moments for (int n = 0; n < mellinint::mdim; n++) storemoments(n, fx1[n]); //storemoments(n, fx1[n], fx2[n]); } /* void pdfevol::invert(double x, double fx[2*MAXNF+1]) { fx = {0.}; complex fm[2*MAXNF+1]; complex fm2[2*MAXNF+1]; double lx = log(x); //positive branch for (int m = 0; m < mellinint::mdim; m++) { complex cexp = mellinint::CCp/2./M_PI/complex (0.,1) * exp(-mellinint::Np[m] * lx); for (int f = 0; f < 2*MAXNF+1; f++) { fm[f] = cx(creno_.cfx1_[m][f]); fx[f] += cexp * fm[f] * mellinint::wn[m]; } //negative branch for (int m = 0; m < mellinint::mdim; m++) { complex cexm = mellinint::CCm/2./M_PI/complex (0.,1) * exp(-mellinint::Nm[m] * lx); for (int f = 0; f < 2*MAXNF+1; f++) { fm[f] = cx(creno_.cfx2m_[m][f]); fx1[f] -= mellinint::CCm/2./M_PI/complex (0.,1) * conj(fm1[f])*(1.-pow(x1,mellinint::Np[n]-mellinint::Nm[m]))/(mellinint::Np[n]-mellinint::Nm[m]) * mellinint::wn[m]; fx2[f] -= mellinint::CCm/2./M_PI/complex (0.,1) * conj(fm2[f])*(1.-pow(x2,mellinint::Np[n]-mellinint::Nm[m]))/(mellinint::Np[n]-mellinint::Nm[m]) * mellinint::wn[m]; } cout << "truncated " << n << fx1[5] << endl; } */ diff --git a/resum/pdfevol.h b/resum/pdfevol.h index 6b8aa8b..d7a863d 100644 --- a/resum/pdfevol.h +++ b/resum/pdfevol.h @@ -1,108 +1,119 @@ #ifndef pdfevol_h #define pdfevol_h #include "interface.h" #include using namespace std; extern "C" { void pdfmoments_(int &beam, double &scale, fcomplex &N, fcomplex &UV, fcomplex &DV, fcomplex &US, fcomplex &DS, fcomplex &SP, fcomplex &SM, fcomplex &GL, fcomplex &CH, fcomplex &BO, double &xmin); //access dyres PDF in N-space extern struct { // fcomplex cfx1_[136][11]; // fcomplex cfx2p_[136][11]; // fcomplex cfx2m_[136][11]; fcomplex cfx1_[512][11]; fcomplex cfx2p_[512][11]; fcomplex cfx2m_[512][11]; } creno_; } #pragma omp threadprivate(creno_) namespace pdfevol { //PDFs mellin moments at the factorisation scale extern complex *UVP; extern complex *DVP; extern complex *USP; extern complex *DSP; extern complex *SSP; extern complex *GLP; extern complex *CHP; extern complex *BOP; #pragma omp threadprivate(UVP,DVP,USP,DSP,SSP,GLP,CHP,BOP) + //Singlet/non-singlet decomposition + extern complex *SIP; + extern complex *NS3P; + extern complex *NS8P; + extern complex *NS15P; + extern complex *NS24P; + extern complex *NS35P; +#pragma omp threadprivate(SIP,NS3P,NS8P,NS15P,NS24P,NS35P) + //evolved PDFs mellin moments extern complex *fx1; extern complex *fx2; #pragma omp threadprivate(fx1,fx2) //moments for the PDFs convolution extern complex fn1[2*MAXNF+1]; extern complex fn2[2*MAXNF+1]; #pragma omp threadprivate(fn1,fn2) //scales extern complex bscale; extern complex bstarscale; extern complex bstartilde; extern complex qbstar; extern complex bcomplex; extern complex XL; extern complex XL1; extern complex SALP; extern complex alpr; #pragma omp threadprivate(bscale,bstarscale,bstartilde,qbstar,bcomplex,XL,XL1,SALP,alpr) extern void allocate(); //allocate dynamic memory extern void free(); //free dynamic memory extern void allocate_fx(); //allocate dynamic memory extern void free_fx(); //free dynamic memory //initialise and compute Mellin moments of PDFs at the starting scale (factorisation scale) extern void init(); //release memory if fmufac = 0 extern void release(); //update Mellin moments of PDFs at the starting scale with a dynamic factorisation scale extern void update(); //evolve Mellin moments of PDFs from the factorisation scale to the scale ~1/b, set in bscale - extern void evolution(int i); + extern void evolution(); //calculate Mellin moments of PDFs at all scales by direct Mellin transform extern void calculate(int i); extern void calculate(); //store the moments in the fx1 and fx2 arrays extern void storemoments(int i, complex fx[11]); //store the moments in the Fortran common block extern void storemoments_fortran(int i, complex fx[11]); //retrieve moments corresponding to i1, i2, and sign from the fx1 and fx2 arrays and store them in fn1 and fn2 extern void retrieve(int i1, int i2, int sign); //retrieve moments corresponding to i1, i2, and sign from the Fortran common block and store them in fn1 and fn2 extern void retrieve_fortran(int i1, int i2, int sign); //retrieve only fn1 extern void retrieve_beam1(int i1); //retrieve only fn2 positive extern void retrieve_beam2_pos(int i2); //retrieve only fn2 negative extern void retrieve_beam2_neg(); //flavour dependent form factors extern void flavour_kt(); - //retrieve moments corresponding to i, and sign from fx1 and fx2 and store them in fn1 and fn2, to be used with the mellin1d option - extern void retrieve1d(int i, int sign); + //retrieve moments corresponding to i, and positive branch from fx1 and fx2 and store them in fn1 and fn2, to be used with the mellin1d option + extern void retrieve1d_pos(int i); + //retrieve moments corresponding to i, and negative branch by complex conjugation of the positive branch moments into fn1 fn2, to be used with the mellin1d option + extern void retrieve1d_neg(); //retrieve moments corresponding to i, and sign from the Fortran common block and store them in fn1 and fn2, to be used with the mellin1d option extern void retrieve1d_fortran(int i, int sign); - extern void retrievemuf(int i, int sign); + extern void retrievemuf(int i); extern void truncate(); extern void uppertruncate(); } #endif