Beam functions

src/eko/basis_rotation.py

@@ -149,6 +149,9 @@
     (non_singlet_pids_map["ns+u"], 0),
     (non_singlet_pids_map["ns-u"], 0),
 )
+
+scet_labels = ((1, 1), (1, 21), (21, 1), (21, 21), (1,2), (1,-1), (1,-2))
+
 full_labels = (*singlet_labels, *non_singlet_labels)
 full_unified_labels = (
     *singlet_unified_labels,

src/eko/evolution_operator/beam_function.py

@@ -0,0 +1,160 @@
+"""SCET I kernels"""
+
+import copy
+import functools
+import logging
+
+import numba as nb
+import numpy as np
+
+import ekore.scet_I as scet_I
+
+from .. import basis_rotation as br
+from .. import scale_variations as sv
+from ..matchings import Segment
+from . import Operator, QuadKerBase
+
+logger = logging.getLogger(__name__)
+
+@nb.njit(cache=True)
+def quad_ker(
+    u,
+    order,
+    space,
+    nf,
+    mode0,
+    mode1,
+    is_log,
+    logx,
+    areas,
+):
+    r"""Raw kernel inside quad.
+
+    Parameters
+    ----------
+    u : float
+        quad argument
+    order : tuple(int,int)
+        perturbation matching order
+    mode0 : int
+        pid for first element in the singlet sector
+    mode1 : int
+        pid for second element in the singlet sector
+    is_log : boolean
+        logarithmic interpolation
+    logx : float
+        Mellin inversion point
+    areas : tuple
+        basis function configuration
+
+    Returns
+    -------
+    ker : float
+        evaluated integration kernel
+
+    """
+    ker_base = QuadKerBase(u, is_log, logx, mode0)
+    integrand = ker_base.integrand(areas)
+    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, nf, ker_base.n)
+    # select the needed matrix element
+    ker = A[indices[mode0], indices[mode1]]
+
+    # recombine everything
+    return np.real(ker * integrand)
+
+
+class SCET_I(Operator):
+    r"""
+    Internal representation of a single |SCET1| mathcing kernel.
+
+    The actual matrices are computed upon calling :meth:`compute`.
+
+    Parameters
+    ----------
+    config : dict
+        configuration
+    managers : dict
+        managers
+    order: tuple (int, int)
+        order in alpha_s and L
+    """
+
+    log_label = "Scet_I"
+    full_labels = br.scet_labels
+
+    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):
+        """Necessary sector labels to compute.
+
+        Returns
+        -------
+        list(str)
+            sector labels
+        """
+        labels = []
+
+
+        labels.extend(
+            [
+                *br.scet_labels,
+            ]
+        )
+
+        return labels
+
+    def quad_ker(self, label, logx, areas):
+        """Return partially initialized integrand function.
+
+        Parameters
+        ----------
+        label: tuple
+            operator element pids
+        logx: float
+            Mellin inversion point
+        areas : tuple
+            basis function configuration
+
+        Returns
+        -------
+        functools.partial
+            partially initialized integration kernel
+        """
+        return functools.partial(
+            quad_ker,
+            order=self.order_scet,
+            space=self.space,
+            nf=self.nf,
+            mode0=label[0],
+            mode1=label[1],
+            is_log=self.int_disp.log,
+            logx=logx,
+            areas=areas,
+        )
+
+    @property
+    def a_s(self):
+        """Return the computed values for :math:`a_s`.
+
+        Note that here you need to use :math:`a_s^{n_f+1}`
+        """
+        sc = self.managers["couplings"]
+        return sc.a_s(
+            self.q2_from
+            * (self.xif2 if self.sv_mode == sv.Modes.exponentiated else 1.0),
+            nf_to=self.nf + 1,
+        )
+
+    def compute(self):
+        """Compute the actual operators (i.e. run the integrations)."""
+        self.initialize_op_members()
+        self.integrate()

src/eko/evolution_operator/scet_kernel.py

@@ -0,0 +1,60 @@
+"""Defines the SCET kernel."""
+
+from .. import basis_rotation as br
+from .. import member
+
+
+class ScetKernel(member.OperatorBase):
+    """
+    Scet kernel for |PDF|.
+
+    """
+
+    @classmethod
+    def split_ad_to_evol_map(
+        cls,
+        op_members,
+        q2_thr,
+    ):
+        """
+        Create the instance from the |OME|.
+
+        Parameters
+        ----------
+            op_members : eko.beam_fuctions.SCET_I.op_members
+                Attribute of :class:`~eko.beam_fuctions.SCET_I` containing the scet kernels
+            nf : int
+                number of active flavors *below* the threshold
+            q2_thr: float
+                dummy value
+        """
+
+        quark_names = br.flavor_basis_names
+        quarks = quark_names[8:]
+
+        m = {}
+
+        for q1 in quarks:
+            for q2 in quarks:
+                if q1 == q2:
+                    m[f"{q1}.{q1}"] = op_members[(1, 1)]
+                    m[f"{q1}bar.{q1}bar"] = op_members[(1, 1)]
+                    m[f"{q1}.{q1}bar"] = op_members[(1, -1)]
+                    m[f"{q1}bar.{q1}"] = op_members[(1, -1)]         
+                else:
+                    m[f"{q1}.{q2}"] = op_members[(1, 2)]
+                    m[f"{q2}.{q1}"] = op_members[(1, 2)]
+                    m[f"{q1}bar.{q2}bar"] = op_members[(1, 2)]
+                    m[f"{q2}bar.{q1}bar"] = op_members[(1, 2)]
+                    m[f"{q1}.{q2}bar"] = op_members[(1, -2)]
+                    m[f"{q2}bar.{q1}"] = op_members[(1, -2)]
+
+        for q in quarks:
+            m[f"{q}.g"] = op_members[(1, 21)]
+            m[f"{q}bar.g"] = op_members[(1, 21)]  
+            m[f"g.{q}"] = op_members[(21, 1)]
+            m[f"g.{q}bar"] = op_members[(21, 1)]
+
+        m[f"g.g"] = op_members[(21,21)]
+
+        return cls.promote_names(m, q2_thr)

src/eko/io/items.py

@@ -11,7 +11,7 @@
 from .. import matchings
 from . import exceptions
 from .types import EvolutionPoint as EPoint
-from .types import FlavorIndex, FlavorsNumber, SquaredScale
+from .types import FlavorIndex, FlavorsNumber, SquaredScale, Order, Space
 
 
 @dataclass(frozen=True)
@@ -98,6 +98,26 @@ def as_atlas(self) -> matchings.Matching:
         """The associated segment."""
         return matchings.Matching(self.scale, self.hq, self.inverse)
 
+@dataclass(frozen=True)
+class ScetKernel(Header):
+    """Information to compute a SCET matching operator.
+
+    """
+
+    order: Order
+    space: Space
+    nf: FlavorsNumber
+
+    @classmethod
+    def from_atlas(cls, scet_kernel: matchings.ScetKernel):
+        """Create instance from analogous :class:`eko.matchings.Atlas` object."""
+        return cls(**asdict(scet_kernel))
+
+    @property
+    def as_atlas(self) -> matchings.ScetKernel:
+        """The associated segment."""
+        return matchings.ScetKernel(self.order, self.space, self.nf)
+
 
 Recipe = Union[Evolution, Matching]
 

src/eko/io/paths.py

@@ -11,6 +11,7 @@
 PARTSDIR = "parts"
 MATCHINGDIR = "matching"
 OPERATORSDIR = "operators"
+SCETDIR = "scet"
 
 
 @dataclass
@@ -45,6 +46,11 @@ def recipes_matching(self):
         """Matching recipes folder."""
         return self.root / RECIPESDIR / MATCHINGDIR
 
+    @property
+    def recipes_scet(self):
+        """Scet recipes folder."""
+        return self.root / RECIPESDIR / SCETDIR
+
     @property
     def parts(self):
         """Parts folder."""
@@ -55,6 +61,11 @@ def parts_matching(self):
         """Matching parts folder."""
         return self.root / PARTSDIR / MATCHINGDIR
 
+    @property
+    def parts_scet(self):
+        """Scet kernel folder."""
+        return self.root / PARTSDIR / SCETDIR
+
     @property
     def operators(self):
         """Operators folder.
@@ -85,3 +96,5 @@ def bootstrap(self, theory: dict, operator: dict, metadata: dict):
         self.parts.mkdir()
         self.parts_matching.mkdir()
         self.operators.mkdir()
+        self.recipes_scet.mkdir()
+        self.parts_scet.mkdir()

src/eko/io/runcards.py

@@ -31,6 +31,8 @@
     ScaleVariationsMethod,
     SquaredScale,
     T,
+    Space,
+    FlavorsNumber,
 )
 
 
@@ -51,6 +53,12 @@ class TheoryCard(DictLike):
     """|N3LO| anomalous dimension variation: ``(gg_var, gq_var, qg_var, qq_var)``."""
     matching_order: Optional[Order] = None
     """Matching conditions perturbative order tuple, ``(QCD, QED)``."""
+    space: Optional[Space] = None 
+    """Space to compute beam function kernels"""
+    orders: Optional[List[Order]] = None
+    """Beam function order"""
+    nf: Optional[FlavorsNumber] = None
+    """Number of active flavors to use to compute beam function"""
 
 
 @dataclass

src/eko/io/struct.py

@@ -18,7 +18,7 @@
 from .access import AccessConfigs
 from .bases import Bases
 from .inventory import Inventory
-from .items import Evolution, Matching, Operator, Recipe, Target
+from .items import Evolution, Matching, Operator, Recipe, Target, ScetKernel
 from .metadata import Metadata
 from .paths import InternalPaths
 from .runcards import OperatorCard, TheoryCard
@@ -49,6 +49,14 @@ def inventories(path: pathlib.Path, access: AccessConfigs) -> dict:
             paths.parts_matching, access, Matching, name="matching-parts"
         ),
         operators=Inventory(paths.operators, access, Target, name="operators"),
+        recipes_scet=Inventory(
+            paths.recipes_scet,
+            access,
+            ScetKernel,
+            contentless=True,
+            name="scet-recipes",
+        ),
+        scet_kernels=Inventory(paths.parts_scet, access, ScetKernel, name="scet_kernels"),
     )
 
 
@@ -90,6 +98,8 @@ class EKO:
     parts: Inventory[Evolution]
     parts_matching: Inventory[Matching]
     operators: Inventory[Target]
+    recipes_scet: Inventory[ScetKernel]
+    scet_kernels: Inventory[ScetKernel]
 
     # public containers
     # -----------------
@@ -179,8 +189,11 @@ def load_recipes(self, recipes: List[Recipe]):
             # leverage auto-save
             if isinstance(recipe, Evolution):
                 self.recipes[recipe] = None
-            else:
+            elif isinstance(recipe, Matching):
                 self.recipes_matching[recipe] = None
+            else:
+                self.recipes_scet[recipe] = None
+
 
     # operator management
     # -------------------

src/eko/io/types.py

@@ -24,6 +24,7 @@
 FlavorIndex = int
 IntrinsicFlavors = typing.List[FlavorIndex]
 N3LOAdVariation = typing.Tuple[int, int, int, int]
+Space = str
 
 # Evolution coordinates
 # ---------------------