TMC

benchmarks/bench_yadism.py

@@ -24,7 +24,6 @@ def __init__(self, pto):
         ) as file:
             th = yaml.safe_load(file)
 
-        th["TMC"] = 0
         th["IC"] = 0
         th["FactScaleVar"] = False
         th["RenScaleVar"] = False
@@ -37,6 +36,8 @@ def yadism_like(self):
     def dis_tp_like(self):
         new_t_card = self.t_card
         new_t_card["grids"] = True
+        new_t_card["TMC"] = self.t_card["TMC"]
+        # new_t_card["NfFF"] = self.t_card["NfFF"]
         new_t_card["fns"] = "fonll"
         return new_t_card
 

src/dis_tp/io.py

@@ -11,13 +11,15 @@
 class TheoryParameters:
     """Class containing all the theory parameters."""
 
-    def __init__(self, order, fns, masses, sc, grids, full_card=None):
+    def __init__(self, order, fns, masses, sc, grids, tmc, target_mass, full_card=None):
         self.order = order
         self.fns = fns
         self.masses = masses
         self.grids = grids
         self._t_card = full_card
         self.strong_coupling = sc
+        self.tmc = tmc
+        self.target_mass = target_mass
 
         if not self.grids:
             console.log(
@@ -39,11 +41,6 @@ def load_card(cls, configs, name):
             th = name
 
         # Disable some NNPDF features not included here
-        if "TMC" in th and th["TMC"] != 0:
-            console.log(
-                "[red underline]Warning, disable Target Mass Corrections:", "TMC=0"
-            )
-            th["TMC"] = 0
         if "IC" in th and th["IC"] == 1:
             console.log("[red underline]Warning, disable Intrinsic Charm:", "IC=0")
             th["IC"] = 0
@@ -60,6 +57,14 @@ def load_card(cls, configs, name):
             )
             th["RenScaleVar"] = False
 
+        tmc = th.get("TMC", 0)
+        if tmc in [1,3]:
+            tmc = 2
+            console.log(
+                f"[blue underline]Warning, TMC = 2 (approx) is the only one implemented."
+            )
+        target_mass = th.get("MP", 0.938)
+
         # compatibility layer
         # PTO
         if "order" in th:
@@ -94,7 +99,14 @@ def load_card(cls, configs, name):
             thresholds_ratios=[1.0, 1.0, 1.0],
         )
         return cls(
-            order=order, fns=fns, grids=grids, masses=masses, sc=sc, full_card=th
+            order=order,
+            fns=fns,
+            grids=grids,
+            masses=masses,
+            sc=sc,
+            tmc=tmc,
+            target_mass=target_mass,
+            full_card=th,
         )
 
 

src/dis_tp/runner.py

@@ -7,6 +7,7 @@
 
 from . import configs, io
 from .logging import console
+from .tmc import TMC_StructureFunction
 from .parameters import initialize_theory, number_active_flavors
 from .structure_functions import f2, fl
 
@@ -88,6 +89,34 @@ def compute_sf(self, kins):
         # console.log(f"x={x}, Q={q}")
         return float(self.partial_sf(x=x, Q=q))
 
+    def setup_tmc(self, func, nf, ob):
+        if "FL" in func.__name__:
+            partial_fl = functools.partial(
+                func,
+                order=self.t_par.order,
+                meth=self.t_par.fns,
+                pdf=ob.pdf,
+                h_id=nf,
+                target_dict=self.o_par.target_dict,
+            )
+            f2 = mapfunc["F2"][ob.restype][0]
+        else:
+            f2 = func
+            partial_fl = None
+        partial_f2 = functools.partial(
+            f2,
+            order=self.t_par.order,
+            meth=self.t_par.fns,
+            pdf=ob.pdf,
+            h_id=nf,
+            target_dict=self.o_par.target_dict,
+        )
+        return TMC_StructureFunction(
+            f2=partial_f2,
+            target_mass=self.t_par.target_mass,
+            fl=partial_fl,
+        )
+
     def compute(self, n_cores):
         # loop on observables
         for ob in self.o_par.obs:
@@ -102,14 +131,19 @@ def compute(self, n_cores):
 
             # loop on SF
             for func in func_to_call:
-                self.partial_sf = functools.partial(
-                    func,
-                    order=self.t_par.order,
-                    meth=self.t_par.fns,
-                    pdf=ob.pdf,
-                    h_id=nf,
-                    target_dict=self.o_par.target_dict,
-                )
+
+                if self.t_par.tmc == 2:
+                    self.partial_sf = self.setup_tmc(func, nf, ob)
+                else:
+                    self.partial_sf = functools.partial(
+                        func,
+                        order=self.t_par.order,
+                        meth=self.t_par.fns,
+                        pdf=ob.pdf,
+                        h_id=nf,
+                        target_dict=self.o_par.target_dict,
+                    )
+
                 console.log(
                     f"[green]Computing {func.__name__} @ order: {self.t_par.order} ..."
                 )

src/dis_tp/tmc.py

@@ -0,0 +1,68 @@
+import numpy as np
+from scipy import integrate
+
+
+class TMC_StructureFunction:
+    """Compute TMC structure function."""
+
+    def __init__(self, target_mass, f2, fl=None):
+        self.f2 = f2
+        self.fl = fl
+        self.target_mass = target_mass
+        self.kind = "F2" if self.fl is None else "FL"
+
+    def __call__(self, x, Q):
+        self.x = x
+        self.Q = Q
+
+        # compute variables
+        self.mu = self.target_mass**2 / self.Q**2
+        self.rho = np.sqrt(1 + 4 * self.x**2 * self.mu)  # = r = sqrt(tau)
+        self.xi = 2 * self.x / (1 + self.rho)
+
+        # run actual computaion
+        if self.fl is None:
+            return self.tmc_f2()
+        return self.tmc_fl()
+
+    def h2(self):
+        # TODO: here you need to use interpolation this is too slow
+        return integrate.quad(
+            lambda u: self.f2(x=u, Q=self.Q) / u,
+            self.xi,
+            1,
+        )
+
+    # def tmc_f2_apfel(self):
+    #     _factor_shifted = self.x**2 / (self.xi**2 * self.rho**3)
+    #     _factor_h2 = 6.0 * self.mu * self.x**3 / (self.rho**4)
+    #     F2out = self.f2(x=self.xi, Q=self.Q)
+    #     h2out = self.h2()
+    #     return _factor_shifted * F2out + _factor_h2 * h2out
+
+    # def tmc_fl_apfel(self):
+    #     _factor_shifted = self.x**2 / (self.xi**2 * self.rho)
+    #     _factor_h2 = 4.0 * self.mu * self.x**3 / (self.rho**2)
+    #     FLout = self.fl(x=self.xi, Q=self.Q)
+    #     h2out = self.h2()
+    #     return _factor_shifted * FLout + _factor_h2 * h2out
+
+    def tmc_f2(self):
+        _factor_shifted = self.x**2 / (self.xi**2 * self.rho**3)
+        approx_prefactor = _factor_shifted * (
+            1 + ((6 * self.mu * self.x * self.xi) / self.rho) * (1 - self.xi) ** 2
+        )
+        F2out = self.f2(x=self.xi, Q=self.Q)
+        return approx_prefactor * F2out
+
+    def tmc_fl(self):
+        _factor_shifted = self.x**2 / (self.xi**2 * self.rho)
+        approx_prefactor_F2 = _factor_shifted * (
+            (4 * self.mu * self.x * self.xi) / self.rho * (1 - self.xi)
+            + 8
+            * (self.mu * self.x * self.xi / self.rho) ** 2
+            * (-np.log(self.xi) - 1 + self.xi)
+        )
+        FLout = self.fl(x=self.xi, Q=self.Q)
+        F2out = self.f2(x=self.xi, Q=self.Q)
+        return _factor_shifted * FLout + approx_prefactor_F2 * F2out