diff --git a/benchmarks/bench_yadism.py b/benchmarks/bench_yadism.py index 1610e2e..02d729c 100644 --- a/benchmarks/bench_yadism.py +++ b/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 diff --git a/src/dis_tp/io.py b/src/dis_tp/io.py index b123ac4..b087cf2 100644 --- a/src/dis_tp/io.py +++ b/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, ) diff --git a/src/dis_tp/runner.py b/src/dis_tp/runner.py index 2acd433..b683e36 100644 --- a/src/dis_tp/runner.py +++ b/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} ..." ) diff --git a/src/dis_tp/tmc.py b/src/dis_tp/tmc.py new file mode 100644 index 0000000..bd6c64b --- /dev/null +++ b/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