Backup more debug stuff

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])

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())

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)])