Interface QED and N3LO, n.2

src/eko/basis_rotation.py

@@ -275,7 +275,7 @@
 }
 
 
-def ad_projector(ad_lab, nf, qed=False):
+def ad_projector(ad_lab, nf, qed):
     """
     Build a projector (as a numpy array) for the given anomalous dimension sector.
 
@@ -355,7 +355,7 @@ def select_light_flavors_uni_ev(ad_lab, nf):
         return map_ad_to_evolution[ad_lab]
 
 
-def ad_projectors(nf, qed=False):
+def ad_projectors(nf, qed):
     """
     Build projectors tensor (as a numpy array), collecting all the individual sector projectors.
 

src/eko/evolution_operator/__init__.py

@@ -307,6 +307,7 @@ def quad_ker(
             ev_op_max_order,
             sv_mode,
             is_threshold,
+            n3lo_ad_variation,
         )
 
     # recombine everything
@@ -450,6 +451,7 @@ def quad_ker_qed(
     ev_op_max_order,
     sv_mode,
     is_threshold,
+    n3lo_ad_variation,
 ):
     """Raw evolution kernel inside quad.
 
@@ -489,6 +491,8 @@ def quad_ker_qed(
         scale variation mode, see `eko.scale_variations.Modes`
     is_threshold : boolean
         is this an itermediate threshold operator?
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -497,7 +501,7 @@ def quad_ker_qed(
     """
     # compute the actual evolution kernel for QEDxQCD
     if ker_base.is_QEDsinglet:
-        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf)
+        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
@@ -622,7 +626,7 @@ def __init__(
         self.alphaem_running = self.managers["couplings"].alphaem_running
         if self.log_label == "Evolution":
             self.a = self.compute_a()
-            self.compute_aem_list()
+            self.as_list, self.a_half_list = self.compute_aem_list()
 
     @property
     def n_pools(self):
@@ -697,8 +701,8 @@ def compute_aem_list(self):
         """
         ev_op_iterations = self.config["ev_op_iterations"]
         if self.order[1] == 0:
-            self.as_list = np.array([self.a_s[0], self.a_s[1]])
-            self.a_half_list = np.zeros((ev_op_iterations, 2))
+            as_list = np.array([self.a_s[0], self.a_s[1]])
+            a_half = np.zeros((ev_op_iterations, 2))
         else:
             as0 = self.a_s[0]
             as1 = self.a_s[1]
@@ -718,7 +722,7 @@ def compute_aem_list(self):
             couplings = self.managers["couplings"]
             mu2_steps = utils.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations)
             mu2_l = mu2_steps[0]
-            self.as_list = np.array(
+            as_list = np.array(
                 [
                     couplings.compute(
                         a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2
@@ -734,7 +738,7 @@ def compute_aem_list(self):
                 )
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
-            self.a_half_list = a_half
+        return as_list, a_half
 
     @property
     def labels(self):

src/eko/evolution_operator/flavors.py

@@ -52,7 +52,7 @@ def pids_from_intrinsic_evol(label, nf, normalize):
     return weights
 
 
-def get_range(evol_labels, qed=False):
+def get_range(evol_labels, qed):
     """Determine the number of light and heavy flavors participating in the input and output.
 
     Here, we assume that the T distributions (e.g. T15) appears *always*
@@ -73,7 +73,7 @@ def get_range(evol_labels, qed=False):
     nf_in = 3
     nf_out = 3
 
-    def update(label, qed=False):
+    def update(label, qed):
         nf = 3
         if label[0] == "T":
             if not qed:
@@ -129,7 +129,7 @@ def rotate_pm_to_flavor(label):
     return l
 
 
-def rotate_matching(nf, qed=False, inverse=False):
+def rotate_matching(nf, qed, inverse=False):
     """Rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
@@ -206,7 +206,7 @@ def rotate_matching(nf, qed=False, inverse=False):
     return l
 
 
-def rotate_matching_inverse(nf, qed=False):
+def rotate_matching_inverse(nf, qed):
     """Inverse rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters

src/eko/evolution_operator/matching_condition.py

@@ -21,7 +21,7 @@ def split_ad_to_evol_map(
         nf,
         q2_thr,
         intrinsic_range,
-        qed=False,
+        qed,
     ):
         """
         Create the instance from the |OME|.

src/eko/evolution_operator/physical.py

@@ -22,7 +22,7 @@ class PhysicalOperator(member.OperatorBase):
     """
 
     @classmethod
-    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False):
+    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed):
         """
         Obtain map between the 3-dimensional anomalous dimension basis and the 4-dimensional evolution basis.
 

src/eko/kernels/non_singlet.py

@@ -186,7 +186,7 @@ def nnlo_expanded(gamma_ns, a1, a0, beta):
 
 
 @nb.njit(cache=True)
-def n3lo_expanded(gamma_ns, a1, a0, nf):
+def n3lo_expanded(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet expanded EKO.
 
     Parameters
@@ -206,12 +206,8 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet expanded EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_expanded(a1, a0, beta0, b_list)
     j23 = as4_ei.j23_expanded(a1, a0, beta0, b_list)
@@ -225,7 +221,7 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
 
 
 @nb.njit(cache=True)
-def n3lo_exact(gamma_ns, a1, a0, nf):
+def n3lo_exact(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet exact EKO.
 
     Parameters
@@ -245,12 +241,8 @@ def n3lo_exact(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet exact EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     roots = as4_ei.roots(b_list)
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_exact(a1, a0, beta0, b_list, roots)
@@ -420,6 +412,6 @@ def dispatcher(
             "decompose-expanded",
             "perturbative-expanded",
         ]:
-            return n3lo_expanded(gamma_ns, a1, a0, nf)
-        return n3lo_exact(gamma_ns, a1, a0, nf)
+            return n3lo_expanded(gamma_ns, a1, a0, betalist)
+        return n3lo_exact(gamma_ns, a1, a0, betalist)
     raise NotImplementedError("Selected order is not implemented")

src/eko/kernels/non_singlet_qed.py

@@ -120,6 +120,29 @@ def as3_exact(gamma_ns, a1, a0, beta):
     return ns.nnlo_exact(gamma_ns, a1, a0, beta)
 
 
+@nb.njit(cache=True)
+def as4_exact(gamma_ns, a1, a0, beta):
+    """O(as3aem1) non-singlet exact EKO.
+
+    Parameters
+    ----------
+    gamma_ns : numpy.ndarray
+        non-singlet anomalous dimensions
+    a1 : float
+        target coupling value
+    a0 : float
+        initial coupling value
+    beta : list
+        list of the values of the beta functions
+
+    Returns
+    -------
+    e_ns^3 : complex
+        O(as4aem1) non-singlet exact EKO
+    """
+    return ns.n3lo_exact(gamma_ns, a1, a0, beta)
+
+
 @nb.njit(cache=True)
 def dispatcher(
     order,
@@ -216,6 +239,8 @@ def fixed_alphaem_exact(order, gamma_ns, a1, a0, aem, nf, mu2_from, mu2_to):
         qcd_only = as2_exact(gamma_ns_list[1:], a1, a0, betalist)
     elif order[0] == 3:
         qcd_only = as3_exact(gamma_ns_list[1:], a1, a0, betalist)
+    elif order[0] == 4:
+        qcd_only = as4_exact(gamma_ns_list[1:], a1, a0, betalist)
     else:
         raise NotImplementedError("Selected order is not implemented")
     return qcd_only * apply_qed(gamma_ns_list[0], mu2_from, mu2_to)

src/eko/member.py

@@ -266,7 +266,7 @@ def operator_multiply(left, right, operation):
                     new_oms[new_key] += operation(l_op, r_op)
         return new_oms
 
-    def to_flavor_basis_tensor(self, qed: bool = False):
+    def to_flavor_basis_tensor(self, qed: bool):
         """Convert the computations into an rank 4 tensor.
 
         A sparse tensor defined with dot-notation (e.g. ``S.g``) is converted

src/ekobox/apply.py

@@ -12,7 +12,6 @@ def apply_pdf(
     lhapdf_like,
     targetgrid=None,
     rotate_to_evolution_basis=False,
-    qed=False,
 ):
     """
     Apply all available operators to the input PDFs.
@@ -34,6 +33,7 @@ def apply_pdf(
         out_grid : dict
             output PDFs and their associated errors for the computed mu2grid
     """
+    qed = eko.theory_card.order[1] > 0
     if rotate_to_evolution_basis:
         if not qed:
             rotate_flavor_to_evolution = br.rotate_flavor_to_evolution
@@ -99,6 +99,7 @@ def apply_pdf_flavor(
         if error_final is not None:
             out_grid[ep]["errors"] = dict(zip(eko.bases.targetpids, error_final))
 
+    qed = eko.theory_card.order[1] > 0
     # rotate to evolution basis
     if flavor_rotation is not None:
         for q2, op in out_grid.items():

src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py

@@ -155,6 +155,11 @@ def gamma_ns_qed(order, mode, n, nf):
             gamma_ns[3, 0] = as3.gamma_nsp(n, nf, cache)
         elif mode in [10202, 10203]:
             gamma_ns[3, 0] = as3.gamma_nsm(n, nf, cache)
+    if order[0] >= 4:
+        if mode in [10102, 10103]:
+            gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache)
+        elif mode in [10202, 10203]:
+            gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache)
     return gamma_ns
 
 
@@ -250,7 +255,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache):
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(order, n, nf):
+def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
     r"""
     Compute the grid of the QED singlet anomalous dimensions matrices.
 
@@ -262,6 +267,8 @@ def gamma_singlet_qed(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -280,6 +287,8 @@ def gamma_singlet_qed(order, n, nf):
         gamma_s[0, 2] = aem2.gamma_singlet(n, nf, cache)
     if order[0] >= 3:
         gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, cache)
+    if order[0] >= 4:
+        gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation)
     return gamma_s
 
 
@@ -314,4 +323,6 @@ def gamma_valence_qed(order, n, nf):
         gamma_v[0, 2] = aem2.gamma_valence(n, nf, cache)
     if order[0] >= 3:
         gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, cache)
+    if order[0] >= 4:
+        gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache)
     return gamma_v