From a8c1a0f140d4a32945e08beca340986d901c596e Mon Sep 17 00:00:00 2001 From: Giacomo Magni Date: Tue, 14 Feb 2023 13:54:38 +0100 Subject: [PATCH] implement approx tmc --- benchmarks/bench_yadism.py | 6 ++-- src/dis_tp/io.py | 26 +++++++++++---- src/dis_tp/runner.py | 49 ++++++++++++++++++++++----- src/dis_tp/tmc.py | 68 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 131 insertions(+), 18 deletions(-) create mode 100644 src/dis_tp/tmc.py diff --git a/benchmarks/bench_yadism.py b/benchmarks/bench_yadism.py index 8bee2be..08b7306 100644 --- a/benchmarks/bench_yadism.py +++ b/benchmarks/bench_yadism.py @@ -26,7 +26,6 @@ def __init__(self, pto): ) as file: th = yaml.safe_load(file) - th["TMC"] = 0 th["IC"] = 0 th["FactScaleVar"] = False th["RenScaleVar"] = False @@ -39,6 +38,7 @@ def yadism_like(self): def dis_tp_like(self): new_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["masses"] = [default_masses(4), default_masses(5), default_masses(6)] new_t_card["fns"] = "fonll" @@ -187,5 +187,5 @@ def benchmarkFONLL(pto, pdf_name, heavyness): # obj = benchmarkF_M_charm(pto=1, pdf_name=pdf_name) # obj = benchmarkFO_charm(pto=1, pdf_name=pdf_name) # benchmarkFONLL(pto=2, pdf_name=pdf_name, heavyness="charm") - # benchmarkFONLL(pto=3, pdf_name=pdf_name, heavyness="light") - benchmarkFONLL(pto=2, pdf_name=pdf_name, heavyness="total") + benchmarkFONLL(pto=3, pdf_name=pdf_name, heavyness="light") + # benchmarkFONLL(pto=2, pdf_name=pdf_name, heavyness="total") diff --git a/src/dis_tp/io.py b/src/dis_tp/io.py index cb4a7a4..52a7c07 100644 --- a/src/dis_tp/io.py +++ b/src/dis_tp/io.py @@ -10,13 +10,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 def yadism_like(self): return self._t_card @@ -33,11 +35,6 @@ def load_card(cls, configs, name): th = name # Disable some NNPDF features not included here - if "TMC" in th and th["TMC"] == 1: - 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 @@ -54,6 +51,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: @@ -84,7 +89,14 @@ def load_card(cls, configs, name): hqm_scheme=th.get("HQ", "POLE"), ) 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 b0fc2e5..e03b18c 100644 --- a/src/dis_tp/runner.py +++ b/src/dis_tp/runner.py @@ -9,6 +9,7 @@ from . import configs, io from .parameters import initialize_theory, number_active_flavors from .logging import console +from .tmc import TMC_StructureFunction mapfunc = { "F2": { @@ -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 @@ -104,14 +133,18 @@ 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