High-pT Ψ(2S) studies

.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

run/LHAPDFgrid.f

@@ -17,6 +17,7 @@ program Evolution
 *
       call ReadParameters(card)
 *
-      call LHAPDFgrid(2.8d0, "DeltaGluonFF")
+      call LHAPDFgrid(3.0d0, "DeltaGluonFF")
 *
+      stop
       end

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
-#########################################
+XMAXLHA  = 0.99999999999999D0           # x upper bound
+#########################################
+#########################################

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`")

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)