From 910549800b427134e555ebea61759db8a571e522 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Fri, 8 Jul 2022 15:28:33 +0200 Subject: [PATCH] Backup more debug stuff --- debug.py | 17 ++++++++++++ sv.py | 28 ++++++++++++------- yadism_debug.py | 72 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 108 insertions(+), 9 deletions(-) create mode 100644 debug.py create mode 100644 yadism_debug.py diff --git a/debug.py b/debug.py new file mode 100644 index 00000000..9e506f05 --- /dev/null +++ b/debug.py @@ -0,0 +1,17 @@ +# -*- coding: utf-8 -*- +import eko +import lhapdf +import numpy as np +from ekobox import apply + +gonly = lhapdf.mkPDF("gonly", 0) + +ho = eko.output.Output.load_tar("data/ekos/3201/test.tar") +hpdf = apply.apply_pdf(ho, gonly) + +bo = eko.output.Output.load_tar("data/ekos/2201/test.tar") +bpdf = apply.apply_pdf(bo, gonly) + +plot(np.log(a[:, 2]), a[:, 4]) +plot(np.log(a[:, 2]), a[:, 6]) +plot(np.log(a[:, 2]), 4 * 5 / 9 * a[:, 2] * bpdf[300.0]["pdfs"][2]) diff --git a/sv.py b/sv.py index 8c101e15..77662cd9 100644 --- a/sv.py +++ b/sv.py @@ -3,28 +3,36 @@ import click import lhapdf +import numpy as np import pandas as pd import pineappl # global settings here = pathlib.Path(__file__).parent fk_path = here / "data" / "fktables" -fk_name = "HERA_CC_318GEV_EM_SIGMARED" +fk_name = "test" +# fk_name = "HERA_CC_318GEV_EM_SIGMARED" # fk_name = "HERA_NC_225GEV_EP_SIGMARED" is_dis = True @click.command() @click.argument("order", type=int) +@click.argument("value", type=float) @click.option("--abs", is_flag=True, help="keep absolute values") -def run(order, abs): +def run(order, value, abs): # derive order specific settings if order == 2: pdf_name = "NNPDF40_nnlo_as_01180" - theories = {"central": 200, "C": 301, "B": 302} + central = 200 + theories = {"A": 10200, "B": 20200, "C": 30200} else: - pdf_name = "NNPDF40_nlo_as_01180" - theories = {"central": 208, "C": 310, "B": 311} + pdf_name = "gonly" + central = 4001 + if np.isclose(value, np.sqrt(2), atol=0.05): + theories = {"A": 1140, "B": 2140, "C": 3140} + else: + theories = {"B": 2201, "C": 3201} # derive common settings pdf = lhapdf.mkPDF(pdf_name, 0) @@ -33,7 +41,7 @@ def run(order, abs): # collect central df = pd.DataFrame() fk = pineappl.fk_table.FkTable.read( - fk_path / str(theories["central"]) / f"{fk_name}.pineappl.lz4" + fk_path / str(central) / f"{fk_name}.pineappl.lz4" ) for bd in range(fk.bin_dimensions()): df[f"b{bd} l"] = fk.bin_left(bd) @@ -42,20 +50,22 @@ def run(order, abs): df["central"] = fk.convolute_with_one(pdgid, pdf.xfxQ2) # collect SVs - for tn in ("B", "C"): - ti = theories[tn] + for tn, ti in theories.items(): fk_sv = pineappl.fk_table.FkTable.read( fk_path / str(ti) / f"{fk_name}.pineappl.lz4" ) df[tn] = fk_sv.convolute_with_one(pdgid, pdf.xfxQ2) df[f"{tn}2central [%]"] = (df[tn] / df["central"] - 1.0) * 100.0 + # df["ΔA2ΔC [%]"] = ( + # (df["A"] - df["central"]) / (df["C"] - df["central"]) - 1.0 + # ) * 100.0 df["ΔB2ΔC [%]"] = ( (df["B"] - df["central"]) / (df["C"] - df["central"]) - 1.0 ) * 100.0 # remove predictions if not abs: - for tn in ("B", "C"): + for tn in theories.keys(): del df[tn] print(df.to_string()) diff --git a/yadism_debug.py b/yadism_debug.py new file mode 100644 index 00000000..8404e5dd --- /dev/null +++ b/yadism_debug.py @@ -0,0 +1,72 @@ +# -*- coding: utf-8 -*- +import numpy as np +import yadism +import yaml +from eko.anomalous_dimensions import as1 +from eko.harmonics import S1 +from eko.interpolation import InterpolatorDispatcher +from scipy.integrate import quad +from scipy.interpolate import interp1d + +# def mel(y, n): +# xs = out["interpolation_xgrid"] +# f = interp1d(xs, y) +# return quad(lambda x: x ** (n - 1.0) * f(x), min(xs), 1.0) + + +def mel2(out, ys, n): + interp = InterpolatorDispatcher.from_dict(out) + lnxmin = np.log(out["interpolation_xgrid"][0]) + res = 0.0 + for y, bf in zip(ys, interp): + pj = bf(n, lnxmin) * np.exp(n * lnxmin) + res += y * pj + return res + + +with open("210.yaml", "r") as f: + t = yaml.safe_load(f) +with open("test.yaml", "r") as f: + o = yaml.safe_load(f) +# out = yadism.output.Output.load_tar("test.tar") +out = yadism.run_yadism(t, o) +out.dump_tar("test.tar") + +locbf10 = [ + out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3][10] + for j in range(len(out["interpolation_xgrid"])) +] +pqqcbf10 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][-3][10] + for j in range(len(out["interpolation_xgrid"])) +] +pgqgbf10 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][7][10] + for j in range(len(out["interpolation_xgrid"])) +] + +mel2bf10 = np.array([mel2(out, [0] * 10 + [1] + [0] * 39, n) for n in range(1, 10)]) +mel2lobf10 = np.array([mel2(out, locbf10, n) for n in range(1, 10)]) +mel2pqqcbf10 = np.array([mel2(out, pqqcbf10, n) for n in range(1, 10)]) +mel2pgqgbf10 = np.array([mel2(out, pgqgbf10, n) for n in range(1, 10)]) + +locbf20 = [ + out["F2_light"][j].orders[(0, 0, 0, 0)][0][-3][20] + for j in range(len(out["interpolation_xgrid"])) +] +pqqcbf20 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][-3][20] + for j in range(len(out["interpolation_xgrid"])) +] +pgqgbf20 = [ + out["F2_light"][j].orders[(1, 0, 0, 1)][0][7][20] + for j in range(len(out["interpolation_xgrid"])) +] + +mel2bf20 = np.array([mel2(out, [0] * 20 + [1] + [0] * 29, n) for n in range(1, 10)]) +mel2lobf20 = np.array([mel2(out, locbf20, n) for n in range(1, 10)]) +mel2pqqcbf20 = np.array([mel2(out, pqqcbf20, n) for n in range(1, 10)]) +mel2pgqgbf20 = np.array([mel2(out, pgqgbf20, n) for n in range(1, 10)]) + +gns = np.array([as1.gamma_ns(n, S1(n)) for n in range(1, 10)]) +ggq = np.array([np.nan if n == 1 else as1.gamma_gq(n) for n in range(1, 10)])