refactoring and adding nf dependence

src/eko/evolution_operator/beam_function.py

@@ -21,6 +21,7 @@ def quad_ker(
     u,
     order,
     space,
+    nf,
     mode0,
     mode1,
     is_log,
@@ -57,7 +58,7 @@ def quad_ker(
     if integrand == 0.0:
         return 0.0
     indices = {21: 0, 1: 1, -1: 2, 2: 3, -2: 4}
-    A = scet_I.SCET_I_entry(order, space, ker_base.n)
+    A = scet_I.SCET_I_entry(order, space, nf, ker_base.n)
     # select the needed matrix element
     ker = A[indices[mode0], indices[mode1]]
 
@@ -84,11 +85,12 @@ class SCET_I(Operator):
     log_label = "Scet_I"
     full_labels = br.scet_labels
 
-    def __init__(self, config, managers, order, space):
-        super().__init__(config, managers, Segment(origin=1, target=1, nf=5))
+    def __init__(self, config, managers, order, space, nf):
+        super().__init__(config, managers, Segment(origin=1, target=1, nf=nf))
         # order (alpha_s, L) of the SCET kernel
         self.order_scet = order
         self.space = space
+        self.nf = nf
 
     @property
     def labels(self):
@@ -131,6 +133,7 @@ def quad_ker(self, label, logx, areas):
             quad_ker,
             order=self.order_scet,
             space=self.space,
+            nf=self.nf,
             mode0=label[0],
             mode1=label[1],
             is_log=self.int_disp.log,

src/eko/runner/parts.py

@@ -178,13 +178,14 @@ def scetI(eko: EKO, recipe: ScetKernel) -> Operator:
     """Compute scetI matching kernels"""
     order = recipe.order
     space = recipe.space
-    #nf = recipe.nf #starting from NNLO from here you can recover the flavor info
+    nf = recipe.nf # needed starting from NNLO 
 
     kernel = bf.SCET_I(
         scet_configs(eko),
         managers(eko),
         order,
-        space
+        space,
+        nf
     )
     kernel.compute()
 

src/ekore/scet_I/__init__.py

@@ -4,11 +4,59 @@
 import numpy as np
 
 from ..harmonics import cache as c
-from .k_space import as1 as k_as1
-from .tau_space import as1 as tau_as1
+from . import k_space as k
+from . import tau_space as tau
 
 @nb.njit(cache=True)
-def SCET_I_entry(order, space, n):
+def A_entries(n, nf, order, space, cache):
+    r"""Compute the beam function matching kernel at the given order.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    numpy.ndarray
+        `
+
+    """
+    if space=='k':
+        Agg = k.A_gg(n, nf, order, cache)
+        Agq = k.A_gq(n, nf, order, cache)
+        Aqg = k.A_qg(n, nf, order, cache)
+        Aqq = k.A_qq(n, nf, order, cache)
+        AqQ2 = k.A_qQ2(n, nf, order, cache)
+        Aqqbar = k.A_qqbar(n, nf, order, cache)
+        AqQ2bar = k.A_qQ2bar(n, nf, order, cache)
+
+    if space=='tau':
+        Agg = tau.A_gg(n, nf, order, cache)
+        Agq = tau.A_gq(n, nf, order, cache)
+        Aqg = tau.A_qg(n, nf, order, cache)
+        Aqq = tau.A_qq(n, nf, order, cache)
+        AqQ2 = tau.A_qQ2(n, nf, order, cache)
+        Aqqbar = tau.A_qqbar(n, nf, order, cache)
+        AqQ2bar = tau.A_qQ2bar(n, nf, order, cache)
+
+    A_S = np.array(
+        [
+            [Agg, Agq, Agq, Agq, Agq],
+            [Aqg, Aqq, Aqqbar, AqQ2, AqQ2bar],
+            [Aqg, Aqqbar, Aqq, AqQ2bar, AqQ2],
+            [Aqg, AqQ2, AqQ2bar, Aqq, Aqqbar],
+            [Aqg, AqQ2bar, AqQ2bar, Aqqbar, Aqq],
+        ],
+        np.complex_,
+    )
+    return A_S
+
+
+@nb.njit(cache=True)
+def SCET_I_entry(order, space, nf, n):
     r"""Compute the tower of the singlet |OME|.
 
     Parameters
@@ -30,9 +78,6 @@ def SCET_I_entry(order, space, n):
     """
     cache = c.reset()
     A = np.zeros((5, 5), np.complex_)
-    if space=='k':
-        A = k_as1.A_entries(n, order, cache)
-    if space=='tau':
-        A = tau_as1.A_entries(n, order, cache)
+    A = k.A_entries(n, nf, order, space, cache)
 
     return A

src/ekore/scet_I/k_space/__init__.py

@@ -0,0 +1,267 @@
+import numba as nb
+import numpy as np
+
+from eko.constants import CF, zeta2
+from . import as1, as2
+from ...harmonics import cache as c
+
+
+@nb.njit(cache=True)
+def A_gg(n, nf, order, cache):
+    r"""
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    complex
+        :math:`A_{gg}`
+
+    """
+
+    if order == (1, 0):
+        res = as1.Agg10(n, cache)
+
+    if order == (1, 1):
+        res = as1.Agg11(n, cache)
+
+    if order == (1, 2):
+        res = as1.Agg12(n, cache)
+
+    if order == (2, 0):
+        res = as2.Agg20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.Agg21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.Agg22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.Agg23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.Agg24(n, nf, cache)
+
+    return res
+
+@nb.njit(cache=True)
+def A_gq(n, nf, order, cache):
+    r"""
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    complex
+        :math:`A_{gq}`
+
+    """
+
+    if order == (1, 0):
+        res = as1.Agq10(n, cache)
+
+    if order == (1, 1):
+        res = as1.Agq11(n, cache)
+
+    if order == (1, 2):
+        res = as1.Agq12(n, cache)
+
+    if order == (2, 0):
+        res = as2.Agq20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.Agq21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.Agq22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.Agq23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.Agq24(n, nf, cache)
+
+    return res
+
+@nb.njit(cache=True)
+def A_qg(n, nf, order, cache):
+    r"""
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    complex
+        |NLO| :math:`A_{qg}`
+
+    """
+
+    if order == (1, 0):
+        res = as1.Aqg10(n, cache)
+
+    if order == (1, 1):
+        res = as1.Aqg11(n, cache)
+
+    if order == (1, 2):
+        res = as1.Aqg12(n, cache)
+
+    if order == (2, 0):
+        res = as2.Aqg20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.Aqg21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.Aqg22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.Aqg23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.Aqg24(n, nf, cache)
+
+    return res
+
+@nb.njit(cache=True)
+def A_qq(n, nf, order, cache):
+    r"""
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    complex
+        |NLO| :math:`A_{qq}`
+
+    """
+
+    if order == (1, 0):
+        res = as1.Aqq10(n, cache)
+
+    if order == (1, 1):
+        res = as1.Aqq11(n, cache)
+
+    if order == (1, 2):
+        res = as1.Aqq12(n, cache)
+
+    if order == (2, 0):
+        res = as2.Aqq20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.Aqq21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.Aqq22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.Aqq23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.Aqq24(n, nf, cache)
+
+    return res
+
+@nb.njit(cache=True)
+def A_qQ2(n, nf, order, cache):
+    if order == (1, 0):
+        res = 0.0
+
+    if order == (1, 1):
+        res = 0.0
+
+    if order == (1, 2):
+        res = 0.0
+
+    if order == (2, 0):
+        res = as2.AqQ20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.AqQ21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.AqQ22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.AqQ23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.AqQ24(n, nf, cache)
+
+    return res
+
+
+@nb.njit(cache=True)
+def A_qQ2bar(n, nf, order, cache):
+    if order == (1, 0):
+        res = 0.0
+
+    if order == (1, 1):
+        res = 0.0
+
+    if order == (1, 2):
+        res = 0.0
+
+    if order == (2, 0):
+        res = as2.AqQbar20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.AqQbar21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.AqQbar22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.AqQbar23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.AqQbar24(n, nf, cache)
+
+    return res
+
+@nb.njit(cache=True)
+def A_qqbar(n, nf, order, cache):
+    if order == (1, 0):
+        res = 0.0
+
+    if order == (1, 1):
+        res = 0.0
+
+    if order == (1, 2):
+        res = 0.0
+
+    if order == (2, 0):
+        res = as2.Aqqbar20(n, nf, cache)
+
+    if order == (2, 1):
+        res = as2.Aqqbar21(n, nf, cache)
+
+    if order == (2, 2):
+        res = as2.Aqqbar22(n, nf, cache)
+
+    if order == (2, 3):
+        res = as2.Aqqbar23(n, nf, cache)
+
+    if order == (2, 4):
+        res = as2.Aqqbar24(n, nf, cache)
+
+    return res
+
+
+