diff --git a/.circleci/config.yml b/.circleci/config.yml deleted file mode 100644 index 36acdf0..0000000 --- a/.circleci/config.yml +++ /dev/null @@ -1,13 +0,0 @@ -version: 2 -jobs: - build: - machine: - image: ubuntu-2004:202010-01 - steps: - - checkout - - run: sudo apt-get update - - run: sudo apt-get install gfortran g++ gcc - - run: cmake . && make -j8 - - run: cd run/; echo "Reference.ini" | ./Evolution; cd ..; - - run: cd run/; echo "Reference.ini" | ./Timing; cd ..; - - run: cd run/; echo "Reference.ini" | ./StructureFunctions; cd ..; diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..441f1e0 --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,25 @@ +name: Build and test + +on: + push: + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Install dependencies + run: sudo apt-get update && sudo apt-get install -y gfortran g++ gcc + + - name: Build + run: cmake . && make -j$(nproc) + + - name: Run Evolution + run: cd run && echo "Reference.ini" | ./Evolution + + - name: Run Timing + run: cd run && echo "Reference.ini" | ./Timing + + - name: Run StructureFunctions + run: cd run && echo "Reference.ini" | ./StructureFunctions diff --git a/run/LHAPDFgrid.f b/run/LHAPDFgrid.f index 181271a..e8979cd 100644 --- a/run/LHAPDFgrid.f +++ b/run/LHAPDFgrid.f @@ -17,6 +17,7 @@ program Evolution * call ReadParameters(card) * - call LHAPDFgrid(2.8d0, "DeltaGluonFF") + call LHAPDFgrid(3.0d0, "DeltaGluonFF") * + stop end diff --git a/run/Reference.ini b/run/Reference.ini index d35370e..1b60afb 100644 --- a/run/Reference.ini +++ b/run/Reference.ini @@ -7,24 +7,32 @@ ######################################### # Evolution parameters # ######################################### -IPT = 2 # Perturbative order (0,1,2) +IPT = 1 # Perturbative order (0,1,2) MODEV = "ITE" # Solution of the DGLAP equation ("ITE", "TRN", "PTH", "GFN") NS = "VFNS" # Mass scheme used for the LHA evoultion ("VFNS","FFNS") -NFMAX = 6 # Maximum number of flavours for the LHA evolution (3,4,5,6) -NFFN = 3 # Number of active flavours in the FFNS (3,4,5,6) -HQMASS = "POLE" # Heavy quark mass scheme ("POLE","MSBAR") +NFMAX = 5 # Maximum number of flavours for the LHA evolution (3,4,5,6) +NFFN = 5 # Number of active flavours in the FFNS (3,4,5,6) +HQMASS = "MSBAR" # Heavy quark mass scheme ("POLE","MSBAR") KRF = 1d0 # muR / muF -EVOL = "SPACE" # Evolution ("SPACE" = space-like (PDFs),"TIME" = time-like (FFs)) -POL = "OFF" # Polarization of the evolution ("OFF" = unpolarized, "ON" = polarized) +EVOL = "TIME" # Evolution ("SPACE" = space-like (PDFs),"TIME" = time-like (FFs)) +POL = "OFF" # Polarization of the evolution ("OFF" = unpolarized, "ON" = polarized) +DISTF = "FF3S11v4" # Distribution function +# options "deltagNLO","deltag","pwave0","pwave1","pwave2","c3s11LO","c3s11NLO", "constant", "Gc3S11","Fc3S11""g3s18MA" +#"d0g3S11" ; "c3S11LO"; g3S11c0 +# FENG FENG : FENGg3S11 +# full FF3S11 +# full FF3S11v4 +# full FF3S11c c +# full FF3S11g g ######################################### # Coupling # ######################################### -QREF = 1.414213562363095d0, 0d0 # Reference scale for alpha_s [GeV] -ASREF = 0.35d0, 0d0 # alpha_s at the refence scale +QREF = 91.1876d0, 0d0 # Reference scale for alpha_s [GeV] +ASREF = 0.118d0, 0d0 # alpha_s at the refence scale ######################################### # Heavy quark masses # ######################################### -MC = 1.414213562373095d0 # Charm threshold [GeV] +MC = 1.5d0 # Charm threshold [GeV] MB = 4.5d0 # Bottom threshold [GeV] MT = 175d0 # Top threshold [GeV] ######################################### @@ -36,12 +44,13 @@ KFACQ = 1.0d0 # muF / Q ######################################### # LHAPDF grid parameters # ######################################### -NXLHA = 100 # Number of x-space grid nodes -NXMLHA = 50 # Number of log x-spage nodes -NQ2LHA = 50 # Number of Q-space grid nodes -Q2MINLHA = 1D0 # Q2 lower bound +NXLHA = 300 # Number of x-space grid nodes +NXMLHA = 150 # Number of log x-spage nodes +NQ2LHA = 200 # Number of Q-space grid nodes +Q2MINLHA = 9D0 # Q2 lower bound Q2MAXLHA = 1D8 # Q2 upper bound XMINLHA = 1D-5 # x lower bound XMLHA = 1D-1 # x transition point -XMAXLHA = 1D0 # x upper bound -######################################### \ No newline at end of file +XMAXLHA = 0.99999999999999D0 # x upper bound +######################################### +######################################### diff --git a/scripts/compare_Nscan.py b/scripts/compare_Nscan.py new file mode 100644 index 0000000..66bee39 --- /dev/null +++ b/scripts/compare_Nscan.py @@ -0,0 +1,114 @@ +""" +Compare MELA N-space results with numerical Mellin transforms +of the LHAPDF x-space grid. + +Reads `mela_Nscan.dat` (produced by LHAPDFgrid) and computes the +corresponding Mellin moments from the LHAPDF grid for comparison. +""" +import numpy as np +import lhapdf +import mpmath as mp +import matplotlib.pyplot as plt + +lhapdf.setVerbosity(0) +PDF = lhapdf.mkPDF("DeltaGluonFF") + +data = [] +with open("../build/run/mela_Nscan.dat") as f: + for line in f: + if line.startswith("#"): + continue + parts = line.split() + N = float(parts[0]) + Q = float(parts[1]) + gluon = float(parts[2]) + charm = float(parts[3]) + if not (np.isnan(gluon) or np.isnan(charm)): + data.append((N, Q, gluon, charm)) + +data = np.array(data) +Ns = data[:, 0] +Qs = data[:, 1] +mela_gluon = data[:, 2] +mela_charm = data[:, 3] + +# Compute LHAPDF Mellin transforms at each (N, Q) point +lhapdf_gluon = np.zeros(len(data)) +lhapdf_charm = np.zeros(len(data)) + +print(f"Computing {len(data)} Mellin integrals...") +for i, (N, Q, mg, mc) in enumerate(data): + def g_integrand(z, nn=N, qq=Q): + return z**(nn - 1) * PDF.xfxQ(0, float(z), qq) / float(z) + + def c_integrand(z, nn=N, qq=Q): + return z**(nn - 1) * PDF.xfxQ(4, float(z), qq) / float(z) + + lhapdf_gluon[i] = float(mp.quad(g_integrand, [PDF.xMin, PDF.xMax])) + lhapdf_charm[i] = float(mp.quad(c_integrand, [PDF.xMin, PDF.xMax])) +print("Done.") + +ratio_gluon = lhapdf_gluon / mela_gluon +ratio_charm = lhapdf_charm / mela_charm + +diff_gluon_percent = (lhapdf_gluon - mela_gluon) / mela_gluon * 100 +diff_charm_percent = (lhapdf_charm - mela_charm) / mela_charm * 100 + +print(f"\n{'Q':>10s} {'MELA_g':>14s} {'LHAPDF_g':>14s} {'diff_g %':>10s}" + f" {'ratio_g':>10s} {'MELA_c':>14s} {'LHAPDF_c':>14s}" + f" {'diff_c %':>10s} {'ratio_c':>10s}") +print("-" * 120) + +for i in range(len(data)): + print(f"{Qs[i]:10.2f} {mela_gluon[i]:14.7e} {lhapdf_gluon[i]:14.7e}" + f" {diff_gluon_percent[i]:10.2f} {ratio_gluon[i]:10.4f} {mela_charm[i]:14.7e}" + f" {lhapdf_charm[i]:14.7e} {diff_charm_percent[i]:10.2f} {ratio_charm[i]:10.4f}") + +fig, axes = plt.subplots(2, 2, figsize=(12, 8)) + +# Top left: gluon absolute values +ax = axes[0, 0] +ax.plot(Qs, mela_gluon, "b-", label="MELA (N-space)") +ax.plot(Qs, lhapdf_gluon, "r--", label="LHAPDF (Mellin)") +ax.set_xlabel("Q [GeV]") +ax.set_ylabel("g(N=6.2, Q)") +ax.set_title("Gluon Mellin moment") +ax.legend() +ax.set_xscale("log") + +# Top right: charm absolute values +ax = axes[0, 1] +ax.plot(Qs, mela_charm, "b-", label="MELA (N-space)") +ax.plot(Qs, lhapdf_charm, "r--", label="LHAPDF (Mellin)") +ax.set_xlabel("Q [GeV]") +ax.set_ylabel("c(N=6.2, Q)") +ax.set_title("Charm Mellin moment") +ax.legend() +ax.set_xscale("log") + +# Bottom left: gluon ratio +ax = axes[1, 0] +ax.axhline(y=1, color="gray", linestyle=":", alpha=0.5) +ax.axhspan(0.99, 1.01, color="green", alpha=0.1, label="1% band") +ax.plot(Qs, ratio_gluon, "b.-") +ax.set_xlabel("Q [GeV]") +ax.set_ylabel("LHAPDF / MELA") +ax.set_title("Gluon ratio") +ax.set_xscale("log") +ax.legend() + +# Bottom right: charm ratio +ax = axes[1, 1] +ax.axhline(y=1, color="gray", linestyle=":", alpha=0.5) +ax.axhspan(0.99, 1.01, color="green", alpha=0.1, label="1% band") +ax.plot(Qs, ratio_charm, "r.-") +ax.set_xlabel("Q [GeV]") +ax.set_ylabel("LHAPDF / MELA") +ax.set_title("Charm ratio") +ax.set_xscale("log") +ax.legend() + +plt.suptitle("MELA N-space vs LHAPDF Mellin transform (N=6.2)", fontsize=14) +plt.tight_layout() +plt.savefig("compare_Nscan.pdf", dpi=350) +print("\nPlots saved to `compare_Nscan.pdf`") diff --git a/src/commons/colfact.h b/src/commons/colfact.h index aa3e4b7..207fbaf 100644 --- a/src/commons/colfact.h +++ b/src/commons/colfact.h @@ -2,8 +2,9 @@ * * QCD colour factors * - double precision CA,CF,TR + double precision CA,CF,TR,NUMFL * parameter(CA = 3d0) parameter(CF = 4d0/3d0) parameter(TR = 1d0/2d0) + parameter(NUMFL = 3d0) diff --git a/src/core/hypgeo.f b/src/core/hypgeo.f index 4827641..46d53a1 100644 --- a/src/core/hypgeo.f +++ b/src/core/hypgeo.f @@ -1,5 +1,321 @@ ************************************************************************ * +* KL series implementation of MT of Log(2-z) +* z^m Log(2-z) = Log(2) - sum_1^inf (1/n)(z/2)^n +* log[2]/(N+m) - sum_(k=1)^40 (1/k)(1/2)^k 1/(N+k+m) +* +************************************************************************ + FUNCTION LOG2KL(N,m) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX N + DOUBLE PRECISION m +* +* Output Variable + DOUBLE COMPLEX LOG2KL +* +* Internal Variables + INTEGER k,kcheck + DOUBLE COMPLEX TERM,SUM + DOUBLE PRECISION THRESH +* +* Initialization + SUM = (0.0D0, 0.0D0) + THRESH = 1.0D-6 +* +* Series Summation + DO k = 1, 40 + TERM = 1d0/k*0.5d0**k/(N+m+k) + SUM = SUM + TERM + END DO + LOG2KL=0.693147d0/(N+m)-SUM +* +* Check the N=41 term + kcheck=41 + TERM = 1/kcheck*0.5d0**kcheck/(N+m+kcheck) +C + IF (ABS(TERM) .GT. THRESH) THEN + WRITE(6,*) 'Warning: Term 101 =', TERM, ' exceeds threshold.' + END IF +C + RETURN + END +************************************************************************ +* +* KL series implementation +* sum_n=0^infty (9+k)!/k! (1/2)^k/(N+m+k) * 1/2^10 * 1/9! +* Mellin Transform of (2-z)^-10 +* +************************************************************************ + FUNCTION SMKL10(N,m) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX N + DOUBLE PRECISION m +* +* Output Variable + DOUBLE COMPLEX SMKL10 +* +* Internal Variables + INTEGER k,kcheck + DOUBLE COMPLEX TERM + DOUBLE PRECISION THRESH +* +* Initialization + SMKL10 = (0.0D0, 0.0D0) + THRESH = 1.0D-6 +* +* Series Summation + DO k = 0, 100 + TERM = (9d0+k)*(8d0+k)*(7d0+k)*(6d0+k)*(5d0+k) + & *(4d0+k)*(3d0+k)*(2d0+k) + & *(1d0+k)*0.5d0**k/(N+m+k) + SMKL10 = SMKL10 + TERM + END DO + SMKL10=SMKL10*1d0/362880d0/1024d0 +* +* Check the N=101 term + kcheck=101 + TERM = (9d0+kcheck)*(8d0+kcheck) + & *(7d0+kcheck)*(6d0+kcheck) + & *(5d0+kcheck)*(4d0+kcheck)*(3d0+kcheck) + & *(2d0+kcheck)*(1d0+kcheck)*0.5d0**kcheck/(N+m+kcheck) +C + IF (ABS(TERM) .GT. THRESH) THEN + WRITE(6,*) 'Warning: Term 101 =', TERM, ' exceeds threshold.' + END IF +C + RETURN + END +************************************************************************ +* +* KL series implementation +* sum_n=0^infty (7+k)!/k! (1/2)^k/(N+m+k) * 1/2^8 * 1/7! +* Mellin Transform of (2-z)^-8 +* +************************************************************************ + FUNCTION SMKL8(N,m) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX N + DOUBLE PRECISION m +* +* Output Variable + DOUBLE COMPLEX SMKL8 +* +* Internal Variables + INTEGER k,kcheck + DOUBLE COMPLEX TERM + DOUBLE PRECISION THRESH +* +* Initialization + SMKL8 = (0.0D0, 0.0D0) + THRESH = 1.0D-6 +* +* Series Summation + DO k = 0, 100 + TERM = (7d0+k)*(6d0+k)*(5d0+k) + & *(4d0+k)*(3d0+k)*(2d0+k) + & *(1d0+k)*0.5d0**k/(N+m+k) + SMKL8 = SMKL8 + TERM + END DO + SMKL8=SMKL8*1d0/5040d0/256d0 +* +* Check the N=101 term + kcheck=101 + TERM = (7d0+kcheck)*(6d0+kcheck) + & *(5d0+kcheck)*(4d0+kcheck)*(3d0+kcheck) + & *(2d0+kcheck)*(1d0+kcheck)*0.5d0**kcheck/(N+m+kcheck) +C + IF (ABS(TERM) .GT. THRESH) THEN + WRITE(6,*) 'Warning: Term 101 =', TERM, ' exceeds threshold.' + END IF +C + RETURN + END +************************************************************************ +* +* KL series implementation +* sum_n=0^infty (4+k)!/k! (1/2)^k/(N+m+k) * 1/2^5 * 1/4! +* Mellin Transform of (2-z)^-5 +* +************************************************************************ + FUNCTION SMKL5(N,m) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX N + DOUBLE PRECISION m +* +* Output Variable + DOUBLE COMPLEX SMKL5 +* +* Internal Variables + INTEGER k,kcheck + DOUBLE COMPLEX TERM + DOUBLE PRECISION THRESH +* +* Initialization + SMKL5 = (0.0D0, 0.0D0) + THRESH = 1.0D-6 +* +* Series Summation + DO k = 0, 100 + TERM = (4d0+k)*(3d0+k)*(2d0+k) + & *(1d0+k)*0.5d0**k/(N+m+k) + SMKL5 = SMKL5 + TERM + END DO + SMKL5=SMKL5*1d0/24d0/32d0 +* +* Check the N=101 term + kcheck=101 + TERM = (4d0+kcheck)*(3d0+kcheck) + & *(2d0+kcheck)*(1d0+kcheck)*0.5d0**kcheck/(N+m+kcheck) +C + IF (ABS(TERM) .GT. THRESH) THEN + WRITE(6,*) 'Warning: Term 101 =', TERM, ' exceeds threshold.' + END IF +C + RETURN + END +************************************************************************ +* +* KL series implementation +* sum_n=0^infty (5+k)!/k! (1/2)^k/(N+m+k) * 1/2^6 * 1/5! +* Mellin Transform of (2-z)^-6 +* +************************************************************************ + FUNCTION SMKL6(N,m) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX N + DOUBLE PRECISION m +* +* Output Variable + DOUBLE COMPLEX SMKL6 +* +* Internal Variables + INTEGER k,kcheck + DOUBLE COMPLEX TERM + DOUBLE PRECISION THRESH +* +* Initialization + SMKL6 = (0.0D0, 0.0D0) + THRESH = 1.0D-6 +* +* Series Summation + DO k = 0, 100 + TERM = (5d0+k)*(4d0+k)*(3d0+k)*(2d0+k) + & *(1d0+k)*0.5d0**k/(N+m+k) + SMKL6 = SMKL6 + TERM + END DO + SMKL6=SMKL6*1d0/120d0/64d0 +* +* Check the N=101 term + kcheck=101 + TERM = (5d0+kcheck)*(4d0+kcheck)*(3d0+kcheck) + & *(2d0+kcheck)*(1d0+kcheck)*0.5d0**kcheck/(N+m+kcheck) +* Part that should be here to check poles (N+m+k)=0 +* IF (A.EQ.INT(REAL(A)).AND.REAL(A).LE.0D0.AND.AIMAG(A).EQ.0D0)THEN +* WRITE(6,*) 'Warning: A =', A, ' A!=0,-1,-2,...' +* END IF +C + IF (ABS(TERM) .GT. THRESH) THEN + WRITE(6,*) 'Warning: Term 101 =', TERM, ' exceeds threshold.' + END IF +C + RETURN + END +************************************************************************ +* +* LerchTranscendent.f: +* +* It returns the Lerch Transcendent function. +* +* The function is written in terms of an infinate series +* we truncate at N=100 +* +************************************************************************ + FUNCTION LERCH(Z,S,A) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX Z, S, A +* +* Output Variable + DOUBLE COMPLEX LERCH +* +* Internal Variables + INTEGER N + DOUBLE COMPLEX TERM + DOUBLE PRECISION THRESH +* +* Initialization + LERCH = (0.0D0, 0.0D0) + THRESH = 1.0D-6 +* +* Series Summation + DO N = 0, 100 + TERM = Z**N / (DBLE(N) + A)**S + LERCH = LERCH + TERM + END DO +* +* Check the N=101 term + TERM = Z**101 / (101.0D0 + A)**S + IF (A.EQ.INT(REAL(A)).AND.REAL(A).LE.0D0.AND.AIMAG(A).EQ.0D0)THEN + WRITE(6,*) 'Warning: A =', A, ' A!=0,-1,-2,...' + END IF +C + IF (ABS(TERM) .GT. THRESH) THEN + WRITE(6,*) 'Warning: Term 101 =', TERM, ' exceeds threshold.' + END IF +C + RETURN + END +* =================================================================av== +* ===================================================================== +* +* KL implementation S_3 = zeta(3) - (-1)^l/(l-1)! Gamma^(l-1)(N-1) +* +* ===================================================================== +* + FUNCTION STHREE (Z) +* + IMPLICIT NONE +* +* Input Variables + DOUBLE COMPLEX Z +* +* Output Variable + DOUBLE COMPLEX STHREE +* +* Internal Variables + DOUBLE COMPLEX LERCH,L1,L2,L3 + DOUBLE PRECISION THRESH + + +* ..Shift of the argument using the functional equation +* +* + L1 =(1d0,0d0) + L2 =(3d0,0d0) + L3= Z-L1 + STHREE = 1.2020569031d0 - LERCH(L1,L2,L3) +* + RETURN + END +************************************************************************ +* * hypgeo.f: * * It returns the hypergeometric function/series. diff --git a/src/evolution/LHAPDFgrid.f b/src/evolution/LHAPDFgrid.f index 181da5b..60ac181 100644 --- a/src/evolution/LHAPDFgrid.f +++ b/src/evolution/LHAPDFgrid.f @@ -15,11 +15,20 @@ subroutine LHAPDFgrid(Qin,fname) include "../commons/ns.h" include "../commons/nfmax.h" include "../commons/nffn.h" +* KL + include "../commons/distf.h" +** +* KL parameters +* + double complex QKL(2) + integer iKL,jKL + double complex xfevNfinalKL(-6:6),xfevNKL(13) + double complex NKL ** * Input Variables * - double precision Qin - character*50 fname + double precision Qin,QKLi + character*(*) fname ** * Intenal Variables * @@ -53,17 +62,18 @@ subroutine LHAPDFgrid(Qin,fname) 1 " points in [",q2minLHA,":",q2maxLHA,"] GeV^2" * * Compute x-space grid -* +* + do ix=1,nxLHA - if(ix.le.nxmLHA)then - xbLHA(ix) = xminLHA * ( xmLHA / xminLHA ) - 1 **( dble(ix-1) / dble(nxmLHA-1) ) - else - xbLHA(ix) = xmLHA + ( xmaxLHA - xmLHA ) - 1 * ( dble(ix - nxmLHA) - 2 / dble(nxLHA - nxmLHA) ) - endif - enddo + if (ix.le.nxmLHA) then + xbLHA(ix) = xminLHA * ( xmLHA / xminLHA ) + & ** ( dble(ix-1) / dble(nxmLHA-1) ) + else + xbLHA(ix) = xmLHA + ( xmaxLHA - xmLHA ) + & * ( dble(ix - nxmLHA) + & / dble(nxLHA - nxmLHA) ) + endif + enddo * * Compute Q2 grid. * Use a distribution of the Q2 nodes uniform in ln(ln(Q2/Lambda2)) @@ -143,7 +153,7 @@ subroutine LHAPDFgrid(Qin,fname) 3 * dlog( lnQmax / lnQmin ) ) ) enddo lnQmin = dlog( ( q2LHA(iq2c) ) / Lambda2 ) - if(isg.eq.nffi-1)then + if(isg.eq.nffi)then lnQmax = dlog( q2maxLHA / Lambda2 ) else lnQmax = dlog( ( dble(q2th(isg+2)) ) / Lambda2 ) @@ -209,7 +219,8 @@ subroutine LHAPDFgrid(Qin,fname) * LHAPDF6 output * ln = index(fname,char(0)) - 1 - if(ln.eq.-1) ln = index(fname,char(32)) - 1 + if(ln.le.0) ln = index(fname,char(32)) - 1 + if(ln.le.0) ln = len_trim(fname) * creating main folder call mkdir(fname(1:ln)) * creating info file @@ -293,5 +304,44 @@ subroutine LHAPDFgrid(Qin,fname) * 40 format(13(es14.7,1x)) * + + write(6,*) "------------------KL edit-------------------" + jKL=2 + distf = "FF3S11v4" + QKLi = 3d0 + QKL(1) = dcmplx(QKLi, 0d0) +* +* Write MELA N-space scan to file for comparison +* + open(unit=14,status="unknown",file="mela_Nscan.dat") + write(14,"(a)") "# N Q gluon charm" +* +* Scan 50 Q points from 3.1 to 99 (within LHAPDF grid) +* + do iKL=1,50 + QKL(2) = dcmplx(3d0 + (99d0 - 3d0) + & * dble(iKL) / 50d0, 0d0) +* +* N = 1 + 1e-10 +* + NKL = dcmplx(1d0+1d-10,0d0) + call NDistributions(NKL,jKL,QKL,xfevNKL) + call evln2lhac(xfevNKL,xfevNfinalKL) + write(14,"(f6.2,1x,f10.4,1x,es22.14,1x,es22.14)") + & REAL(NKL),REAL(QKL(2)), + & REAL(xfevNfinalKL(0)),REAL(xfevNfinalKL(4)) +* +* N = 6.2 +* + NKL = dcmplx(6.2d0,0d0) + call NDistributions(NKL,jKL,QKL,xfevNKL) + call evln2lhac(xfevNKL,xfevNfinalKL) + write(14,"(f6.2,1x,f10.4,1x,es22.14,1x,es22.14)") + & REAL(NKL),REAL(QKL(2)), + & REAL(xfevNfinalKL(0)),REAL(xfevNfinalKL(4)) + enddo + close(14) + write(6,*) "MELA N-space scan written to mela_Nscan.dat" + return end diff --git a/src/evolution/NDistributions.f b/src/evolution/NDistributions.f index 264e10d..e065bca 100644 --- a/src/evolution/NDistributions.f +++ b/src/evolution/NDistributions.f @@ -26,6 +26,8 @@ subroutine NDistributions(N,nQ,Q,xfevN) * integer iQ integer i,j +* KL mod + double complex aQCD double complex zfuncns(3),zfuncsg(2,2) double complex zfuncns15(2),zfuncns24(2),zfuncns35(2) double complex zfuncnsv15,zfuncnsv24,zfuncnsv35 @@ -57,8 +59,76 @@ subroutine NDistributions(N,nQ,Q,xfevN) elseif(evol(1:4).eq."TIME")then call HKNSFFsn(N,xfph0N) ! HKNS for pi+ at NLO endif - elseif(distf.eq."deltag")then - call GluonDelta(xfph0N) +* elseif(distf.eq."g3S11v4")then +* call g3S11v4(N,Q(1),xfph0N) +* write(6,*) "v4; alphas(Q)=",12.56637d0*asEvolIni(1) + elseif(distf.eq."deltagNLO")then + call GluonDeltaNLO(N,Q(1),xfph0N) +* write(6,*) "NLO function! Q=", Q(1) + write(6,*) "NLO; alphas(Q)=",12.56637d0*asEvolIni(1) + elseif(distf(1:6).eq."deltag")then + call GluonDelta(N,Q,xfph0N) +* write(6,*) "LO function! alphas(Q)=",12.56637d0*asEvolIni(1) +* write(6,*) "LO function! alphas(Q)=",12.56637d0*asEvolFi(1) + elseif(distf.eq."constant")then + call Constant(N,Q(1),xfph0N) +* write(6,*) "Constant function! Q=", Q(1) + elseif(distf.eq."pwave0")then + call P0wave(N,Q(1),xfph0N) +* write(6,*) "J=0! Q=", Q(1) + elseif(distf.eq."pwave1")then + call P1wave(N,Q(1),xfph0N) +* write(6,*) "J=1! Q=", Q(1) + elseif(distf.eq."pwave2")then + call P2wave(N,Q(1),xfph0N) +* write(6,*) "J=2! Q=", Q(1) + elseif(distf.eq."g3s18")then + call g3s18NLO(N,Q(1),xfph0N) + write(6,*) "alphas(Q)=",12.56637d0*asEvolIni(1) + write(6,*) "Ma function! Q=", Q(1) +* write(6,*) "LO function! alphas(Q)=",12.56637d0*asEvolIni(1) + elseif(distf.eq."g3s18MA")then + call Ma3S18(N,Q(1),xfph0N) + write(6,*) "Ma function! Q=", Q(1) + write(6,*) "alphas(Q)=",12.56637d0*asEvolIni(1) + +*GLUON FF + elseif(distf.eq."d0g3S11")then + call d0g3S11(N,Q(1),xfph0N) + elseif(distf.eq."d2g3S11")then + call d2g3S11(N,Q(1),xfph0N) + elseif(distf.eq."d4g3S11")then + call d4g3S11(N,Q(1),xfph0N) + elseif(distf.eq."g3S11c0")then + call g3S11c0(N,Q(1),xfph0N) + elseif(distf.eq."g3S11c1")then + call g3S11c1(N,Q(1),xfph0N) + elseif(distf.eq."FENGg3S11")then + call FENGg3S11(N,Q(1),xfph0N) +*CHARM FF + elseif(distf.eq."c3S11LO")then + call c3S11LO(N,Q(1),xfph0N) + elseif(distf.eq."c3S11REL")then + call c3S11REL(N,Q(1),xfph0N) + elseif(distf.eq."c3S11PQQ")then + call c3S11PQQ(N,Q(1),xfph0N) + elseif(distf.eq."c3S11NLO")then + call c3S11NLO(N,Q(1),xfph0N) +*ALL together + elseif(distf.eq."FF3S11")then + call FF3S11(N,Q(1),xfph0N) + elseif(distf.eq."FF3S11v4")then + call FF3S11v4(N,Q(1),xfph0N) +* JUST V4 charm and gluon + elseif(distf.eq."d4gc3S11")then + call d4gc3S11(N,Q(1),xfph0N) +*MISC + elseif(distf.eq."TEST")then + call TEST(N,Q(1),xfph0N) + elseif(distf.eq."TESTF")then + call TESTF(N,Q(1),xfph0N) + elseif(distf.eq."TESTS")then + call TESTS(N,Q(1),xfph0N) elseif(distf.eq."XFitter")then call XFitterParametrization(N-1d0,xfph0N) ! XFitter PDFs elseif(distf(1:9).eq."ZeroScale")then @@ -68,6 +138,7 @@ subroutine NDistributions(N,nQ,Q,xfevN) else write(6,*) "Unknown input distributions, distf = ",distf endif +* write(6,*) "alphas(Q=", Q(1), "= ", 12.56637d0*asEvolIni(1) * * Rotate distribution into the evolution basis * @@ -84,8 +155,14 @@ subroutine NDistributions(N,nQ,Q,xfevN) * * Values of alphas at the initial and final scales * +* from XDist KL + asEvolIni(1) = aQCD(Q(1)**2d0) + asEvolFin(1) = aQCD(Q(2)**2d0) + as0 = asEvolIni(iQ) asq = asEvolFin(iQ) +* write(6,*) as0*12.56637d0 +* write(6,*) asq*12.56637d0 * * Call evolution kernels * @@ -163,10 +240,10 @@ subroutine NDistributions(N,nQ,Q,xfevN) * if(iQ.lt.nQ-1)then do i=1,13 - xfev0N(i) = xfevN(i) + xfev0N(i) = xfevN(i) enddo endif enddo -* + return end diff --git a/src/evolution/toyLHPDFs.f b/src/evolution/toyLHPDFs.f index 2e8986b..b4e8608 100644 --- a/src/evolution/toyLHPDFs.f +++ b/src/evolution/toyLHPDFs.f @@ -438,17 +438,22 @@ subroutine SetXFitterParametersMELA(ubar,dbar, ************************************************************************ * * Gluon delta function multiplied by alpha_s +* KL 5/11/24 * ************************************************************************ - subroutine GluonDelta(nff) + subroutine GluonDelta(N,Q,nff) * implicit none * include "../commons/alphas.h" + include "../commons/consts.h" ** * Internal Variables * + double complex N,Q integer iff + double complex alphas + double precision mc,LDME ** * Output Variables * @@ -456,10 +461,1741 @@ subroutine GluonDelta(nff) * * Initialize PDFs to zero * + LDME=1d0 + mc=1.5d0 + alphas=PI*4d0*asEvolIni(1) do iff=-6,6 nff(iff) = (0d0,0d0) enddo - nff(0) = 12.566370614359173d0 * asEvolIni(1) * (1d0, 0d0) + nff(0) = alphas*PI*LDME/(mc**3*24d0) +* write(*,*) "KL inside function !nff(0)", nff(0) * return end +************************************************************************ +* +* constant +* KL 5/11/24 +* +************************************************************************ + subroutine Constant(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" +** +* Internal Variables +* + double complex N,Q + integer iff +** +* Output Variables +* + double complex nff(-6:6) +* +* Initialize PDFs to zero +* + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + nff(0) = 1/N +* write(*,*) "CONSTANT KL inside function !nff(0)", nff(0) +* + return + end + +************************************************************************ +* +* Gluon delta function at NLO KL 08/10/24 +* +************************************************************************ + subroutine GluonDeltaNLO(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex DPSI,PSI,S1,S2,dL,dLD,dT,dTLO,dTNLO,FdT + double complex SHarm,S1N,S1N1,S1N2,S1N22,S2N2,PGG,FRAC +** +* Internal Variables +* + integer iff + double precision mc,A,beta + double complex N,Q,logscale,alphas +** +* Output Variables +* + double complex nff(-6:6) + alphas=PI*4d0*asEvolIni(1) + mc = 1.5d0 + beta = (11d0*CA - 4d0*TR*NUMFL)/3d0 +* +* In the file "andim_nlo.f" the harmonic functions are defined +* S[L=1,N] S1 = EMC + PSI(N+1d0) +* S[L=2,N] S2 = ZETA2 - DPSI(N+1d0,1) +* + S1 = EMC + PSI(N+1d0) + S1N = (EMC + PSI(N+1d0))/N + S1N1 = (EMC + PSI(N+2d0))/(N+1d0) + S1N2 = (EMC + PSI(N+3d0))/(N+2d0) + S1N22 = (EMC + PSI(N-1d0))**2 + S2N2 = ZETA2 - DPSI(N-1d0,1) + +* Longitudinal Component + dLD = 8d0 * mc**3 * N * (N-1d0) + dL = (alphas)**2 * dLD**(-1) +* +* Transverse component + FdT = alphas * PI / (48d0 * mc**3) + dTLO = FdT * (1d0, 0d0) + logscale = log(Q/(2*mc)) + A = beta * (logscale + 13d0/6d0) + 4.15975994110758d0 + SHarm = -6d0*(2*S1N - S1N1 + S1N2) - 3d0*(S1N22 + S2N2 ) + FRAC = 2d0 * (N**2 + N + 1d0)/(N * (N**2 - 1d0) * (N+2d0) ) + PGG = (logscale -1d0/2d0)*6d0*(-S1 +beta/6d0 + FRAC) + dTNLO = SHarm + A + PGG + dT = dTLO + dTNLO*FdT*asEvolIni(1) * 12.566370614359173d0/PI +* Initialize PDFs to zero +* + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* + nff(0) = 2*dT + dL + + return + end + +************************************************************************ +* +* P-wave KL 14/03/25 +* J=0 +* +************************************************************************ + subroutine P0wave(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex DPSI,PSI,S1 + double complex T1,T2,T3 +** +* Internal Variables +* + integer iff + double precision Nc,QJ + double precision mc,mulambda + double complex N,Q,logscale,alphas +** +* Output Variables +* + double complex nff(-6:6) + alphas=PI*4d0*asEvolIni(1) + Nc=3d0 + QJ=1/4d0 + mc = 1.5d0 + S1 = EMC + PSI(N+1d0) + mulambda=mc + T1=QJ-0.5d0*log(mulambda**2/(4*mc**2)) + T2=-S1 + T3=85/8*1/(N+1)-26/8*1/(N+2)+45/4*1/(N)-27/4*1/(N+1) +* + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* + nff(0) = (4d0,0d0)/(9d0*Nc)*(T1+T2+T3)*alphas**2/mc**5 + + return + end +************************************************************************ +* +* P-wave KL 14/03/25 +* J=1 +* +************************************************************************ + subroutine P1wave(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex DPSI,PSI,S1 + double complex T1,T2,T3 +** +* Internal Variables +* + integer iff + double precision Nc,QJ + double precision mc,mulambda + double complex N,Q,logscale,alphas +** +* Output Variables +* + double complex nff(-6:6) + alphas=PI*4d0*asEvolIni(1) + Nc=3d0 + QJ=3/8d0/(3d0) + mc = 1.5d0 + S1 = EMC + PSI(N+1d0) + mulambda=mc + T1=QJ-0.5d0*log(mulambda**2/(4*mc**2)) + T2=-S1 + T3=-3/4*(1/(N+1)+4/(N+2))/3d0 +* + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* + nff(0) = (4d0,0d0)/(9d0*Nc)*(T1+T2+T3)*alphas**2/mc**5 + + return + end +************************************************************************ +* +* P-wave KL 14/03/25 +* J=2 +* +************************************************************************ + subroutine P2wave(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex DPSI,PSI,S1,S2 + double complex T1,T2,T3 +** +* Internal Variables +* + integer iff + double precision Nc,QJ + double precision mc,mulambda + double complex N,Q,logscale,alphas +** +* Output Variables +* + double complex nff(-6:6) + alphas=PI*4d0*asEvolIni(1) + Nc=3d0 + QJ=7/8d0/(5d0) + mc = 1.5d0 + S1 = EMC + PSI(N+1d0) + S1 = EMC + PSI(N+2d0) + mulambda=mc + T1=QJ-0.5d0*log(mulambda**2/(4*mc**2)) + T2=-S1 + T3=5/4*(11/(N+1)-4/(N+2)) +9*(-2*S1/N + S2/(N+1)) +* + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* + nff(0) = (4d0,0d0)/(9d0*Nc)*(T1+T2+T3/5d0)*alphas**2/mc**5 + + return + end +************************************************************************ +* +* Gluon delta function at NLO KL 05/06/25 +* +************************************************************************ + subroutine g3s18NLO(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex DPSI,PSI,S1,S2,S3,S4,S5 + double complex PGG,S,LOGR + double complex T1,T2,T3,T4,T5 +** +* Internal Variables +* + integer iff + double precision mc,beta,NC,alphas + double complex N,Q +** +* Output Variables +* + double complex nff(-6:6) +* +* In the file "andim_nlo.f" the harmonic functions are defined +* S[L=1,N] S1 = EMC + PSI(N+1d0) +* S[L=2,N] S2 = ZETA2 - DPSI(N+1d0,1) +* + NC=3d0 + alphas=PI*4d0*asEvolIni(1) + mc = 1.5d0 + beta = (11d0*CA - 4d0*TR*3d0)/3d0 + LOGR=log(Q**2/(4d0*mc**2)) + + T1=PI*alphas/24d0 + + T2=beta/NC*(LOGR+13d0/3d0)+4d0/NC**2-PI**2/3d0+16d0/3d0*log(2d0) + + S=EMC + PSI(N+1d0) + PGG=2*NC*(-S+1/(N-1d0)-1/(N)+1/(N+1d0)-1/(N+2d0)+beta/(2*NC)) + T3=1/NC*PGG*(LOGR-1d0) + + T4=2d0/(N-1d0)-2d0/(N) + + S1 = (EMC + PSI(N-1d0))**2 + ZETA2 - DPSI(N-1d0,1) + S2 = (EMC + PSI(N))**2 + ZETA2 - DPSI(N,1) + S3 = (EMC + PSI(N+1d0))**2 + ZETA2 - DPSI(N+1d0,1) + S4 = (EMC + PSI(N+2d0))**2 + ZETA2 - DPSI(N+2d0,1) + S5 = (EMC + PSI(N+3d0))**2 + ZETA2 - DPSI(N+3d0,1) + T5= -2d0*(S1-2d0*S2+3d0*S3-2d0*S4+S5) +* Initialize PDFs to zero +* + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* + nff(0)=(T1+alphas**2/12d0/CF*(T2+T3+T4+T5))/mc**3 + + return + end + +************************************************************************ +* +* TEST FUNC +* +* +************************************************************************ + subroutine TEST(N,Q,nff) +* +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q + double complex PSI,DPSI + double complex S1N,S1Np1,S1Np2,S1Np3 + double complex Talpha + double complex Ta0,Ta1,Ta2,Ta3 + double precision alpha1,alpha2,alpha3 +** +* Output Variables +* + double complex nff(-6:6) + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* +* define params from table 3 +* + alpha1=-4.9866d-3 + alpha2=9.8448d-3 + alpha3=1.9512d-2 + +* In the file "andim_nlo.f" the harmonic functions are defined +* S[L=1,N] S1 = EMC + PSI(N+1d0) +* S[L=2,N] S2 = ZETA2 - DPSI(N+1d0,1) +* +*--- Harmonic sums +* + S1N = EMC + PSI(N+1d0) + S1Np1 = EMC + PSI(N+2d0) + S1Np2 = EMC + PSI(N+3d0) + S1Np3 = EMC + PSI(N+4d0) +* +* co-ef of alpha part + Ta0=-S1N/N + Ta1=-S1Np1/(N+1d0) + Ta2=-S1Np2/(N+2d0) + Ta3=-S1Np3/(N+3d0) + + +* making alpha, beta, etc. + Talpha=(alpha1+alpha2+alpha3)*Ta0 + & +(-alpha1-2d0*alpha2-3d0*alpha3)*Ta1 + & +(alpha2+3d0*alpha3)*Ta2 + & -alpha3*Ta3 + + + nff(0)=Talpha + + return + end +************************************************************************ +* +* TEST FUNC 2 +* +* +************************************************************************ + subroutine TESTS(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +** +* Internal Variables +* + integer iff + double complex SMKL8 + double complex N,Q + double precision m +** +** +* Output Variables +* + double complex nff(-6:6) + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + m=0d0 + nff(0)=SMKL8(N,m) + + return + end +************************************************************************ +* +* g->3S18 NLO from MA: Phys. Rev. D, 89(9):094029, 2014 +* +* +************************************************************************ + subroutine Ma3S18(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +** +* Internal Variables +* + integer iff + double complex N,Q,LL + double complex mc,alphas,beta0 + double complex PSI,S,PGG + double complex D1,D2,D21,D22,D23,D24 +** +** +* Output Variables +* + double complex nff(-6:6) + double complex GH + external GH + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + mc=1.5d0 + LL=Q**2/(4d0*mc**2) +* beta def with nf=3 + beta0=(11d0*3d0-2d0*3d0)/6d0 + alphas=PI*4d0*asEvolIni(1) +* S[L=1,N] S1 = EMC + PSI(N+1d0) + S=EMC + PSI(N+1d0) + PGG=6d0*(-S+1d0/(N-1d0)-1d0/N+1d0/(N+1d0) + & -1d0/(N+2d0)+beta0/6d0) + D1=(1d0,0d0)/(24d0) + D21=beta0/3d0*(log(LL)+13/3) + & +4d0/3d0**2-PI**2/3d0+16d0/3d0*log(2d0) + D22=1d0/3d0*PGG*(log(LL)-1d0) + D23=2d0/(N-1d0)-2d0/N + D24=-2*(GH(N+3d0)-2*GH(N+2d0)+3*GH(N+1d0)-2*GH(N)+GH(N-1d0)) + D2=(1d0,0d0)*(D21+D22+D23+D24)/(12d0*CF) +* D2=(1d0,0d0)/(12d0*CF)*(D22) + + nff(0)=((PI*alphas*D1) + (alphas**2*D2)) / (mc**3) +* nff(0)=(PI*alphas*D1)/mc**3 + +* nff(0)=(alphas**2*D2)/mc**3 + return + end +C===================================================================== +C Function GHARM(N): computes S[L=1,N]**2 + S[L=2,N] +C===================================================================== +* * +* In the file "andim_nlo.f" the harmonic functions are defined +* S[L=1,N] S1 = EMC + PSI(N+1d0) +* S[L=2,N] S2 = ZETA2 - DPSI(N+1d0,1) +* S[L=1,N-1] S1 = EMC + PSI(N) +* S[L=2,N-1] S2 = ZETA2 - DPSI(N,1) +* + double complex function GH(X) + implicit none + double complex X + double complex PSI, DPSI + include "../commons/consts.h" + include "../commons/alphas.h" + include "../commons/colfact.h" + + GH = (EMC + PSI(X))**2 + (ZETA2 - DPSI(X,1)) + + return + end +************************************************************************ +* +* v4 gluon divergent part +* JHEP11 2012 020, eq. 7.7 +* Does not include finite piece +* +************************************************************************ + double complex function G4G3S11DIST(N,muL) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* + double complex DPSI,PSI + double complex S1N,S1Nm1,S1Np1,S1Np2,S1Nm2 + double complex S2N,S2Nm1,S2Np1 + double complex PolyN,PolyNp1 + double complex bracket + double complex N + double complex T1,T2,T3,T4,T5,T6,T7,T8,T9 + double complex T91,T92,T93,T93i,T93ii + double precision mc,NC,dimd + double complex muL,Lmu + double precision prefacN,prefacD,prefac +* +*--- parameters +* + NC = 3d0 + mc = 1.5d0 + dimd = 4d0 +* +* + Lmu = log(muL/(2d0*mc)) + + +* In the file "andim_nlo.f" the harmonic functions are defined +* S[L=1,N] S1 = EMC + PSI(N+1d0) +* S[L=2,N] S2 = ZETA2 - DPSI(N+1d0,1) +* + S1N = EMC + PSI(N+1d0) + S1Nm1 = EMC + PSI(N ) + S1Np1 = EMC + PSI(N+2d0) + S1Np2 = EMC + PSI(N+3d0) + S1Nm2 = EMC + PSI(N-1d0) +* + S2N = ZETA2 - DPSI(N+1d0,1) + S2Np1 = ZETA2 - DPSI(N+2d0,1) + S2Nm1 = ZETA2 - DPSI(N ,1) +* +*--- Polygamma combinations +* + + PolyN=(S1N**2-S2N+2d0*ZETA2)/N-2d0*S1N/N**2 + PolyNp1=(S1Np1**2-S2Np1+2d0*ZETA2)/(N+1d0)-2d0*S1Np1/(N+1d0)**2 + +* +*--- build bracket +* + T1=1d0/24d0 - ZETA2 - Lmu/3d0 + Lmu*Lmu + + T2= - S1Nm1*(1d0/3d0 - 2d0*Lmu) + T3= S1Nm1*S1Nm1 + S2Nm1 + + T4=(-1d0/24d0)*(104d0/N - 29d0/(N+1d0) - 10d0/(N+2d0)) + + T5=(7d0/2d0)*(1d0/(N-1d0) - S1Nm2/(N-2d0) - S1Nm1/(N-1d0)) + + T6=(8d0/N - 21d0/(N+1d0) + 10d0/(N+2d0))*Lmu/4d0 + + T7=1/4d0*(5d0*S1N/(N)+36d0*S1Nm1/(N-1d0)-25d0*S1Np1/(N+1d0) + * + 6d0*S1Np2/(N+2d0)) + + T8=39d0/(4d0*(N+1d0)**2) - 3d0/(2d0*(N+2d0)**2) + & - 2d0*(S2N - ZETA2) + + T91= (13d0/2d0)*PolyN - (7d0/2d0)*PolyNp1 + T92= Lmu*((13d0/2d0)*S1N/N - (7d0/2d0)*S1Np1/(N+1d0)) + T93i=ZETA2/N -S1N/(N**2) + T93ii=ZETA2/(N+1d0) -S1Np1/(N+1d0)**2 + T93=-13d0/2d0*T93i+7d0/2d0*T93ii +* +*--- prefactor +* + prefacN = 2d0*(NC*NC - 4d0) + prefacD = (3d0*PI*(dimd-1d0)**3*NC**3) + prefac = prefacN/prefacD + T9 = T91+T92+T93 + + bracket=T1+T2+T3+T4+T5+T6+T7+T8+T9 + + G4G3S11DIST=prefac*bracket +* + return + end +************************************************************************ +* +* function v0 gluon paramterised FF kernel +* JHEP11 2012 020, eq. B2 +* +************************************************************************ + double complex function G0G3S11(N) +* + implicit none +* + double complex N + double complex PSI,DPSI + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables +* + double complex S1N,S1Np1,S1Np2,S1Np3 + double complex S2N,S2Np1,S2Np2 + double complex Talpha,Tbeta,Tmu,Tnu,Tomega + double complex Ta0,Ta1,Ta2,Ta3 + double complex Tb0,Tb1,Tb2 + double complex Tnu1,Tnu2,Tnu3,Tnu4 + double complex To11,To12 + double precision alpha1,alpha2,alpha3,beta1,beta2 + double precision mu1,nu1,nu2,nu3,nu4,omega11,omega12 +* +* Parameters (Table 3) +* + alpha1=-4.9866d-3 + alpha2= 9.8448d-3 + alpha3= 1.9512d-2 + beta1 =-5.1697d-6 + beta2 = 1.0462d-2 + mu1 =-1.8921d-3 + nu1 = 1.2154d-3 + nu2 = 1.3039d-3 + nu3 =-2.7246d-3 + nu4 =-1.4814d-3 + omega11=-1.6910d-2 + omega12= 3.8110d-2 +* +*--- Harmonic sums +* + S1N = EMC + PSI(N+1d0) + S1Np1 = EMC + PSI(N+2d0) + S1Np2 = EMC + PSI(N+3d0) + S1Np3 = EMC + PSI(N+4d0) + + S2N = ZETA2 - DPSI(N+1d0,1) + S2Np1 = ZETA2 - DPSI(N+2d0,1) + S2Np2 = ZETA2 - DPSI(N+3d0,1) +* +*--- Coefficients +* + Ta0=-S1N/N + Ta1=-S1Np1/(N+1d0) + Ta2=-S1Np2/(N+2d0) + Ta3=-S1Np3/(N+3d0) + + Tb0=(S1N*S1N+S2N)/N + Tb1=(S1Np1*S1Np1+S2Np1)/(N+1d0) + Tb2=(S1Np2*S1Np2+S2Np2)/(N+2d0) + + Tnu1=2d0/(N+1d0)**3 + Tnu2=2d0/(N+2d0)**3 + Tnu3=2d0/(N+3d0)**3 + Tnu4=2d0/(N+4d0)**3 + + To11=1d0/(N+1d0)-1d0/(N+2d0) + To12=1d0/(N+1d0)-2d0/(N+2d0)+1d0/(N+3d0) +* +*--- Assemble pieces +* + Talpha=(alpha1+alpha2+alpha3)*Ta0 + & +(-alpha1-2d0*alpha2-3d0*alpha3)*Ta1 + & +(alpha2+3d0*alpha3)*Ta2 + & -alpha3*Ta3 + + Tbeta=(beta1+beta2)*Tb0 + & +(-beta1-2d0*beta2)*Tb1 + & + beta2*Tb2 + + Tmu=-mu1/(N+1d0)**2 + Tnu=nu1*Tnu1+nu2*Tnu2+nu3*Tnu3+nu4*Tnu4 + Tomega=omega11*To11+omega12*To12 +* +*--- Final result +* + G0G3S11 = Talpha + Tbeta + Tmu + Tnu + Tomega + + return + end +************************************************************************ +* +* Subroutine calling gluon FF kernel +* +************************************************************************ + subroutine d0g3S11(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Function declaration +* + double complex G0G3S11 + external G0G3S11 +* +* Internal variables +* + integer iff + double complex N,Q +* +* Output +* + double complex nff(-6:6) + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* +* Call function +* + nff(0) = G0G3S11(N) + + return + end +************************************************************************ +* +* function v2 gluon paramterised FF kernel +* JHEP11 2012 020, eq. B2 +* +************************************************************************ + double complex function G2G3S11(N) +* + implicit none +* + double complex N + double complex PSI,DPSI + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables +* + double complex S1N,S1Np1,S1Np2,S1Np3 + double complex S2N,S2Np1,S2Np2 + double complex Tb,Talpha,Tbeta,Tmu,Tnu,Tomega + double complex Ta0,Ta1 + double complex Tb0,Tb1 + double complex Tnu1 + double complex To11,To12,To13,To14 + double precision b,alpha1,beta1 + double precision mu1,mu2,nu1,omega11,omega12,omega13,omega14 +* +* Parameters (Table 3) +* + b=1.4710d-2 + alpha1=-1.8127d-2 + beta1 =-1.2282d-2 + mu1 =4.1069d-3 + mu2 =1.8341d-2 + nu1 =-1.3653d-3 + omega11=-2.1832d-2 + omega12=4.1531d-3 + omega13=8.4949d-4 + omega14=-3.8207d-3 +* +*--- Harmonic sums +* + S1N = EMC + PSI(N+1d0) + S1Np1 = EMC + PSI(N+2d0) + S1Np2 = EMC + PSI(N+3d0) + S1Np3 = EMC + PSI(N+4d0) + + S2N = ZETA2 - DPSI(N+1d0,1) + S2Np1 = ZETA2 - DPSI(N+2d0,1) + S2Np2 = ZETA2 - DPSI(N+3d0,1) +* +*--- Coefficients +* + Ta0=-S1N/N + Ta1=-S1Np1/(N+1d0) + + Tb0=(S1N*S1N+S2N)/N + Tb1=(S1Np1*S1Np1+S2Np1)/(N+1d0) + + Tnu1=2d0/(N+1d0)**3 + + To11=1d0/(N+1d0)-1d0/(N+2d0) + To12=1d0/(N+1d0)-2d0/(N+2d0)+1d0/(N+3d0) + To13 = 1d0/(N+1d0) - 3d0/(N+2d0) + 3d0/(N+3d0) - 1d0/(N+4d0) + To14=1d0/(N+1d0)-4d0/(N+2d0)+6d0/(N+3d0)-4d0/(N+4d0)+1d0/(N+5d0) + +* +*--- Assemble pieces +* + Tb=b/(N+1d0) + + Talpha=alpha1*Ta0-alpha1*Ta1 + + Tbeta=beta1*Tb0-beta1*Tb1 + + Tmu=-mu1/(N+1d0)**2-mu2/(N+2d0)**2 + + Tnu=nu1*Tnu1 + Tomega=omega11*To11+omega12*To12+omega13*To13+omega14*To14 +* +*--- Final result +* + G2G3S11 = Tb + Talpha + Tbeta + Tmu + Tnu + Tomega + + return + end +************************************************************************ +* +* Subroutine calling gluon FF kernel +* +************************************************************************ + subroutine d2g3S11(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Function declaration +* + double complex G2G3S11 + external G2G3S11 +* +* Internal variables +* + integer iff + double complex N,Q +* +* Output +* + double complex nff(-6:6) + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* +* Call function +* + nff(0) = G2G3S11(N) + + return + end +************************************************************************ +* +* function v4 gluon paramterised FF kernel [Finite part] +* JHEP11 2012 020, eq. B2 +* +************************************************************************ + double complex function G4G3S11(N) +* + implicit none +* + double complex N + double complex PSI,DPSI + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables +* + double complex S1N,S1Np1,S1Np2,S1Np3 + double complex S2N,S2Np1,S2Np2 + double complex Tb,Talpha,Tbeta,Tmu,Tnu,Tomega + double complex Ta0,Ta1,Ta2 + double complex Tb0,Tb1 + double complex Tnu1 + double complex To11,To12,To13,To14 + double precision b,alpha1,alpha2,beta1 + double precision mu1,mu2,mu3,mu4,nu1,omega11 +* +* Parameters (Table 3) +* + b=-3.6531d-2 + alpha1=5.2157d-2 + alpha2=2.4588d-3 + beta1 =3.6020d-2 + mu1 =-7.3565d-3 + mu2 =-2.5387d-2 + mu3 =-9.3477d-3 + mu4 =-3.6926d-3 + nu1 =1.3839d-3 + omega11=7.7264d-2 +* +*--- Harmonic sums +* + S1N = EMC + PSI(N+1d0) + S1Np1 = EMC + PSI(N+2d0) + S1Np2 = EMC + PSI(N+3d0) + S1Np3 = EMC + PSI(N+4d0) + + S2N = ZETA2 - DPSI(N+1d0,1) + S2Np1 = ZETA2 - DPSI(N+2d0,1) + S2Np2 = ZETA2 - DPSI(N+3d0,1) +* +*--- Coefficients +* + Ta0=-S1N/N + Ta1=-S1Np1/(N+1d0) + Ta2=-S1Np2/(N+2d0) + + Tb0=(S1N*S1N+S2N)/N + Tb1=(S1Np1*S1Np1+S2Np1)/(N+1d0) + + Tnu1=2d0/(N+1d0)**3 + + To11=1d0/(N+1d0)-1d0/(N+2d0) + +* +*--- Assemble pieces +* + Tb=b/(N+1d0) + Talpha=(alpha1+alpha2)*Ta0 + & +(-alpha1-2d0*alpha2)*Ta1 + & +alpha2*Ta2 + + + Tbeta=beta1*Tb0-beta1*Tb1 + + Tmu=-mu1/(N+1d0)**2-mu2/(N+2d0)**2 + & -mu3/(N+3d0)**2 + & -mu4/(N+4d0)**2 + + Tnu=nu1*Tnu1 + Tomega=omega11*To11 +* +*--- Final result +* + G4G3S11 = Tb + Talpha + Tbeta + Tmu + Tnu + Tomega + + return + end +************************************************************************ +* +* c->3S11 at v^4 KL (Function) 04/02/26 +* +************************************************************************ + double complex function Fv4c3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex SMKL10 + double complex N + double complex T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11 + double precision m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11 + double precision a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11 + +* powers (m) + m1=1d0 + m2=2d0 + m3=3d0 + m4=4d0 + m5=5d0 + m6=6d0 + m7=7d0 + m8=8d0 + m9=9d0 + m10=10d0 + m11=11d0 + +* coefs of powers (a) + a1=1639680d0 + a2=-12648960d0 + a3=53962240d0 + a4=-141141120d0 + a5=225693536d0 + a6=-222599040d0 + a7=136706160d0 + a8=-52504200d0 + a9=12461885d0 + a10=-1668066d0 + a11=97885d0 +* + T1=a1*SMKL10(N,m1) + T2=a2*SMKL10(N,m2) + T3=a3*SMKL10(N,m3) + T4=a4*SMKL10(N,m4) + T5=a5*SMKL10(N,m5) + T6=a6*SMKL10(N,m6) + T7=a7*SMKL10(N,m7) + T8=a8*SMKL10(N,m8) + T9=a9*SMKL10(N,m9) + T10=a10*SMKL10(N,m10) + T11=a11*SMKL10(N,m11) + Fv4c3S11=2d0/164025d0*(T1+T2+T3+T4+T5+T6+T7+T8+T9+T10+T11) +* + return + end +************************************************************************ +* +* Subroutine calling gluon FF kernel +* +************************************************************************ + subroutine d4g3S11(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Function declaration +* + double complex G4G3S11,G4G3S11DIST + external G4G3S11,G4G3S11DIST +* +* Internal variables +* + integer iff + double complex N,Q,muL + double precision mc +* +* Output +* + double complex nff(-6:6) + + mc=1.5d0 + muL=(1.5d0,0d0) + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* +* Call function +* + nff(0) = G4G3S11(N)+G4G3S11DIST(N,muL) + + return + end + + +************************************************************************ +* +* Subroutine calling v4 corr to charm AND gluon +* +************************************************************************ + subroutine d4gc3S11(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Function declaration +* + double complex G4G3S11,G4G3S11DIST,Fv4c3S11 + external G4G3S11,G4G3S11DIST,Fv4c3S11 +* +* Internal variables +* + integer iff + double complex N,Q,muL + double precision mc,alphas +* +* Output +* + double complex nff(-6:6) + alphas=0.31494871292d0 + mc=1.5d0 + muL=(1.5d0,0d0) + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo +* +* Call function +* + nff(0) = alphas**3*(G4G3S11(N)+G4G3S11DIST(N,muL)) + nff(4) = alphas**2*(Fv4c3S11(N)) + nff(-4) = nff(4) + return + end +************************************************************************ +* +* testing function which works +* +************************************************************************ + double complex function TESTFUNC(X) + implicit none + double complex X + double complex PSI, DPSI + include "../commons/consts.h" + include "../commons/alphas.h" + include "../commons/colfact.h" + + TESTFUNC = 1/X + + return + end +* +* subroutine calling above function +* + subroutine TESTF(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q + double complex TESTFUNC + external TESTFUNC +** +** +* Output Variables +* + double complex nff(-6:6) + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + + nff(0)=TESTFUNC(N) + + return + end +************************************************************************ +* +* c->3S11 at LO KL (Function F) 12/05/25 +* +************************************************************************ + double complex function Fc3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex SMKL6 + double complex N,Q + double complex T1,T2,T3,T4,T5,T6,T7 + double precision m1,m2,m3,m4,m5,m6,m7 + double precision a1,a2,a3,a4,a5,a6,a7 + +* powers (m) + m1=1d0 + m2=2d0 + m3=3d0 + m4=4d0 + m5=5d0 + m6=6d0 + m7=7d0 + +* coefs of powers (a) + a1=16d0 + a2=-64d0 + a3=152d0 + a4=-208d0 + a5=141d0 + a6=-42d0 + a7=5d0 +* + T1=a1*SMKL6(N,m1) + T2=a2*SMKL6(N,m2) + T3=a3*SMKL6(N,m3) + T4=a4*SMKL6(N,m4) + T5=a5*SMKL6(N,m5) + T6=a6*SMKL6(N,m6) + T7=a7*SMKL6(N,m7) + Fc3S11=16d0/243d0*(T1+T2+T3+T4+T5+T6+T7) +* + return + end +* +************************************************************************ +* +* subroutine calling above function +* +************************************************************************ + subroutine c3S11LO(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables +* + integer iff + double complex N, Q + double complex Fc3S11 + external Fc3S11 +* +* Output +* + double complex nff(-6:6) +* + do iff = -6, 6 + nff(iff) = (0d0,0d0) + end do +* + nff(4) = Fc3S11(N) + nff(-4) = nff(4) +* + return + end +************************************************************************ +* +* c->3S11 at v^2 KL (Function G) 12/05/25 +* +************************************************************************ + double complex function Gc3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex SMKL8 + double complex N + double complex T1,T2,T3,T4,T5,T6,T7,T8,T9 + double precision m1,m2,m3,m4,m5,m6,m7,m8,m9 + double precision a1,a2,a3,a4,a5,a6,a7,a8,a9 + +* powers (m) + m1=1d0 + m2=2d0 + m3=3d0 + m4=4d0 + m5=5d0 + m6=6d0 + m7=7d0 + m8=8d0 + m9=9d0 + +* coefs of powers (a) + a1=-10752d0 + a2=62976d0 + a3=-209024d0 + a4=417536d0 + a5=-477536d0 + a6=304928d0 + a7=-106328d0 + a8=19664d0 + a9=-1464d0 +* + T1=a1*SMKL8(N,m1) + T2=a2*SMKL8(N,m2) + T3=a3*SMKL8(N,m3) + T4=a4*SMKL8(N,m4) + T5=a5*SMKL8(N,m5) + T6=a6*SMKL8(N,m6) + T7=a7*SMKL8(N,m7) + T8=a8*SMKL8(N,m8) + T9=a9*SMKL8(N,m9) + Gc3S11=1d0/2187d0*(T1+T2+T3+T4+T5+T6+T7+T8+T9) +* + return + end +************************************************************************ +* +* subroutine calling above function +* +************************************************************************ + subroutine c3S11REL(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables +* + integer iff + double complex N, Q + double complex Gc3S11 + external Gc3S11 +* +* Output +* + double complex nff(-6:6) +* + do iff = -6, 6 + nff(iff) = (0d0,0d0) + end do +* + nff(4) = Gc3S11(N) + nff(-4) = nff(4) +* + return + end +************************************************************************ +* +* c->3S11 NLO (PQQ) KL 12/05/25 +* +************************************************************************ + double complex function PQQc3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex PSI + double complex N + double complex S1Np1,S1Nm1,DELTA + +* S[L=1,N] S1 = EMC + PSI(N+1d0) + + S1Nm1 = EMC + PSI(N) + S1Np1 = EMC + PSI(N+2d0) + DELTA = (1.5d0,0d0) + PQQc3S11=2d0*CF*(-S1Nm1-S1Np1+DELTA) +* + return + end +************************************************************************ +* +* subroutine calling above function +* +************************************************************************ + subroutine c3S11PQQ(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables +* + integer iff + double complex N, Q + double complex PQQc3S11,Fc3S11,P,F + external PQQc3S11,Fc3S11 +* +* Output +* + double complex nff(-6:6) +* + do iff = -6, 6 + nff(iff) = (0d0,0d0) + end do +* + P=PQQc3S11(N) + F=Fc3S11(N) + nff(4) = P*F + nff(-4) = nff(4) +* + return + end +************************************************************************ +* +* c->3S11 NLO (POLY) KL 12/05/25 +* +************************************************************************ + double complex function POLYc3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex N + double complex FACT,F + + FACT=0.6981317008d0 + F=-9.01726/(N+10d0)+18.22777/(N+9d0) + F=F+16.11858/(N+8d0)-82.54936/(N+7d0) + F=F+106.57565/(N+6d0)-72.30107/(N+5d0) + F=F+28.85798/(N+4d0)-6.70607/(N+3d0) + F=F+0.84950/(N+2d0)-0.05376/(N+1d0)-0.00205/N + POLYc3S11=F*FACT +* + return + end +************************************************************************ +* +* subroutine calling above function +* +************************************************************************ + subroutine c3S11NLO(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" +* +* Internal variables + integer iff + double complex N, Q + double complex muR0C,mc,beta0 + double complex F,cAS3 + double complex POLYc3S11,Fc3S11,PQQc3S11 + external POLYc3S11,Fc3S11,PQQc3S11 +* +* Output +* + double complex nff(-6:6) +* + do iff = -6, 6 + nff(iff) = (0d0,0d0) + end do +* + + beta0=(11d0*CA - 4d0*TR*3d0)/3d0 + mc = 1.5d0 + muR0C = 2*mc + F = Fc3S11(N) + cAS3 = ( + * +F/(2d0*PI)*beta0*log(muR0C**2/(2d0*mc)**2) + * +PQQc3S11(N)*F/(4d0*PI)*log(Q**2/(3d0*mc)**2) + * +POLYc3S11(N) + * ) + nff(4) = cAS3 + nff(-4) = nff(4) +* + return + end +************************************************************************ +* +* g->3S11 FENG 12/05/25 +* +************************************************************************ + double complex function c0g3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex N + double precision FACT + double precision a0,a1,a2,a3,a4,a5,a6,a7 + double complex T20,T21,T22,T23,T24,T25,T26,T27 + double complex T1,T2,T3 + double complex SMKL5,LOG2KL + external SMKL5,LOG2KL + + + FACT = 1/(PI*4d0*3645d0) +* a's for (z-2)^-5 term (T2) + a0=-6112d0 + a1=22632d0 + a2=-35176d0 + a3=30160d0 + a4=-15685d0 + a5=5032d0 + a6=-927d0 + a7=76d0 + +* Terms for T2 + T20=a0*SMKL5(N,0d0) + T21=a1*SMKL5(N,1d0) + T22=a2*SMKL5(N,2d0) + T23=a3*SMKL5(N,3d0) + T24=a4*SMKL5(N,4d0) + T25=a5*SMKL5(N,5d0) + T26=a6*SMKL5(N,6d0) + T27=a7*SMKL5(N,7d0) + +* Terms + T1=-480d0*(LOG2KL(N,2d0)-19d0*LOG2KL(N,1d0)+36d0*LOG2KL(N,0d0)) + T2=-64d0*(T20+T21+T22+T23+T24+T25+T26+T27) + T3=480d0*(-1d0/(N+2d0)**2-1d0/(N+1d0)**2) + + c0g3S11=FACT*(T1+T2+T3) +* + return + end + + double complex function c1g3S11(N) + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + + double complex N + double complex F + double precision a0,a1,a2,a3,a4,a5,a6,a7,a8 + + a0=0.068789941142313182 + a1=-0.28294758545234799 + a2=0.43755027838051319 + a3=-0.24403801281005144 + a4=-0.12154015467967838 + a5=0.27499430533498526 + a6=-0.20288020418956876 + a7=0.092142697889357805 + a8=-0.022071291547035798 + + F=a8/(N+8d0)+a7/(N+7d0)+a6/(N+6d0)+a5/(N+5d0) + F=F+a4/(N+4d0)+a3/(N+3d0)+a2/(N+2d0)+a1/(N+1d0)+a0/N + c1g3S11=F +* + return + end +************************************************************* +* c1,c0 subroutines +************************************************************ + subroutine g3S11c0(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q + double complex c0g3S11 + external c0g3S11 +** +** +* Output Variables +* + double complex nff(-6:6) + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + + nff(0)=c0g3S11(N) + + return + end + + subroutine g3S11c1(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q + double complex c1g3S11 + external c1g3S11 +** +** +* Output Variables +* + double complex nff(-6:6) + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + + nff(0)=c1g3S11(N) + + return + end + + subroutine FENGg3S11(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q + double complex c1g3S11,c0g3S11 + external c1g3S11, c0g3S11 +** +** +* Output Variables +* + double complex nff(-6:6) + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + + nff(0)=c0g3S11(N)*log(Q/1.5d0) + c1g3S11(N) + + return + end + +************************************************************* +* FULL sigma up to v^2 +************************************************************ + subroutine FF3S11(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q + double complex gFENG,gLO,g + double complex CAS2,CAS3,CAS2v2,c + + double precision mc,v2,beta0 + double precision alphasG,alphasCLO,alphasCNLO,muR0G,muR0C + + double complex c0g3S11,c1g3S11,G0G3S11,G2G3S11 + double complex F,Fc3S11,POLYc3S11,Gc3S11,PQQc3S11 + external c0g3S11,c1g3S11,G0G3S11,G2G3S11 + external Fc3S11,POLYc3S11,Gc3S11,PQQc3S11 +** +** +* Output Variables +* + double complex nff(-6:6) + + +* set scales and alphas + + beta0=(11d0*CA - 4d0*TR*3d0)/3d0 + + mc = 1.5d0 + muR0G=2d0*mc + muR0C=2d0*mc + + v2=0.5d0 + + alphasG=0.31494871292 + alphasCLO=0.31494871292 + alphasCNLO=0.251822 + +* gluon induced part + + gFENG = alphasG**3*(c0g3S11(N)*log(Q/mc) + c1g3S11(N)) + gLO = alphasG**3*(G0G3S11(N) + v2*G2G3S11(N)) + g = (gLO + gFENG)/mc**3 + +* charm induced part + + F = Fc3S11(N) + cAS2 = alphasCNLO**2*F + cAS3 = alphasCNLO**3*( + * +F/(2d0*PI)*beta0*log(muR0C**2/(2d0*mc)**2) + * +PQQc3S11(N)*F/(4d0*PI)*log(Q**2/(3d0*mc)**2) + * +POLYc3S11(N) + * ) + CAS2v2 =Gc3S11(N)*v2*alphasCLO**2 + c = (CAS2+CAS3+CAS2v2)/mc**3 + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + + nff(0)=g + nff(4)=c + nff(-4)=c + + return + end +************************************************************* +* FULL sigma up to v4 +************************************************************ + subroutine FF3S11v4(N,Q,nff) +* + implicit none +* + include "../commons/alphas.h" + include "../commons/consts.h" + include "../commons/colfact.h" + +** +* Internal Variables +* + integer iff + double complex N,Q,muL + double complex gFENG,gLO,g + double complex CAS2,CAS3,CAS2v2,c + + double precision mc,v2,beta0 + double precision alphasG,alphasCLO,alphasCNLO,muR0G,muR0C + + double complex c0g3S11,c1g3S11,G0G3S11,G2G3S11 + double complex F,Fc3S11,POLYc3S11,Gc3S11,PQQc3S11 + double complex G4G3S11,G4G3S11DIST,gv4 + external c0g3S11,c1g3S11,G0G3S11,G2G3S11 + external Fc3S11,POLYc3S11,Gc3S11,PQQc3S11 + external G4G3S11,G4G3S11DIST +** +** +* Output Variables +* + double complex nff(-6:6) + + +* set scales and alphas + + beta0=(11d0*CA - 4d0*TR*3d0)/3d0 + + mc = 1.5d0 + muL = (1.5d0,0d0) + muR0G=2d0*mc + muR0C=2d0*mc + + v2=0.5d0 + + alphasG=0.31494871292 + alphasCLO=0.31494871292 + alphasCNLO=0.251822 + +* gluon induced part + + gFENG = alphasG**3*(c0g3S11(N)*log(Q/mc) + c1g3S11(N)) + gLO = alphasG**3*(G0G3S11(N) + v2*G2G3S11(N)) + gv4 = alphasG**3*v2**2*(G4G3S11(N)+G4G3S11DIST(N,muL)) + g = (gLO + gFENG + gv4)/mc**3 + +* charm induced part + + F = Fc3S11(N) + cAS2 = alphasCNLO**2*F + cAS3 = alphasCNLO**3*( + * +F/(2d0*PI)*beta0*log(muR0C**2/(2d0*mc)**2) + * +PQQc3S11(N)*F/(4d0*PI)*log(Q**2/(3d0*mc)**2) + * +POLYc3S11(N) + * ) + CAS2v2 =Gc3S11(N)*v2*alphasCLO**2 + c = (CAS2+CAS3+CAS2v2)/mc**3 + + do iff=-6,6 + nff(iff) = (0d0,0d0) + enddo + + nff(0)=g + nff(4)=c + nff(-4)=c + + return + end diff --git a/src/evolution/xDistributions.f b/src/evolution/xDistributions.f index faf0caa..4dea6fe 100644 --- a/src/evolution/xDistributions.f +++ b/src/evolution/xDistributions.f @@ -43,10 +43,12 @@ subroutine xDistributions(x,nQ,Q,xfph) enddo * t = - dlog(x) - m = 33 ! Must be odd + m = 65 ! Must be odd (increased from 33) r = 2d0 * m / 5d0 / t * - rmax = 10d0 +* Adaptive rmax: larger for small t (large x), bounded + rmax = 18d0 / t - 1d0 + rmax = max(min(rmax, 25d0), 1d0) * if(r.gt.rmax) r = rmax *