diff --git a/debug.py b/debug.py
deleted file mode 100644
index 9e506f05..00000000
--- a/debug.py
+++ /dev/null
@@ -1,17 +0,0 @@
-# -*- 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/eko_debug.py b/eko_debug.py
deleted file mode 100644
index e21afe2f..00000000
--- a/eko_debug.py
+++ /dev/null
@@ -1,47 +0,0 @@
-# -*- coding: utf-8 -*-
-import eko
-import lhapdf
-import matplotlib.pyplot as plt
-import numpy as np
-from eko.anomalous_dimensions.as1 import gamma_ns
-from eko.harmonics import S1
-from ekomark import apply
-from scipy.integrate import quad
-from scipy.interpolate import interp1d
-
-
-def mel(y, n):
-    f = interp1d(o["targetgrid"], y)
-    return quad(lambda x: x ** (n - 1.0) * f(x), min(o["targetgrid"]), 1.0)
-
-
-o = eko.output.Output.load_tar("data/ekos/2201/noevol-withK-as0p35.tar")
-
-resu = apply.apply_pdf(o, lhapdf.mkPDF("uonly"))
-resd = apply.apply_pdf(o, lhapdf.mkPDF("donly"))
-rest3 = apply.apply_pdf(o, lhapdf.mkPDF("T3only"))
-restoyt3 = apply.apply_pdf(o, lhapdf.mkPDF("toyt3only"))
-resS = apply.apply_pdf(o, lhapdf.mkPDF("Sonly"))
-resNN40 = apply.apply_pdf(o, lhapdf.mkPDF("NNPDF40_nnlo_as_01180"))
-
-u_u = [lhapdf.mkPDF("uonly").xfxQ2(2, x, 10) / x for x in o["targetgrid"]]
-d_d = [lhapdf.mkPDF("donly").xfxQ2(1, x, 10) / x for x in o["targetgrid"]]
-t3_u = [lhapdf.mkPDF("T3only").xfxQ2(2, x, 10) / x for x in o["targetgrid"]]
-toyt3_u = [lhapdf.mkPDF("toyt3only").xfxQ2(2, x, 10) / x for x in o["targetgrid"]]
-toyt3_d = [lhapdf.mkPDF("toyt3only").xfxQ2(1, x, 10) / x for x in o["targetgrid"]]
-toyS_d = [lhapdf.mkPDF("Sonly").xfxQ2(1, x, 10) / x for x in o["targetgrid"]]
-nn40_u = [
-    lhapdf.mkPDF("NNPDF40_nnlo_as_01180").xfxQ2(2, x, 10) / x for x in o["targetgrid"]
-]
-
-g = np.array([gamma_ns(n, S1(n)) for n in range(1, 10)])
-
-K = 1.0 - 0.35 / (4.0 * np.pi) * g * np.log(2**2)
-
-melu_u = np.array([mel(u_u, n)[0] for n in range(1, 10)])
-meld_d = np.array([mel(d_d, n)[0] for n in range(1, 10)])
-melnn40_u = np.array([mel(nn40_u, n)[0] for n in range(1, 10)])
-
-melresu_u = np.array([mel(resu[300.0]["pdfs"][2], n)[0] for n in range(1, 10)])
-melresd_d = np.array([mel(resd[300.0]["pdfs"][1], n)[0] for n in range(1, 10)])
-melresnn40_u = np.array([mel(resNN40[300.0]["pdfs"][2], n)[0] for n in range(1, 10)])
diff --git a/sv.py b/sv.py
deleted file mode 100644
index 77662cd9..00000000
--- a/sv.py
+++ /dev/null
@@ -1,75 +0,0 @@
-# -*- coding: utf-8 -*-
-import pathlib
-
-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 = "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, value, abs):
-    # derive order specific settings
-    if order == 2:
-        pdf_name = "NNPDF40_nnlo_as_01180"
-        central = 200
-        theories = {"A": 10200, "B": 20200, "C": 30200}
-    else:
-        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)
-    pdgid = int(pdf.set().get_entry("Particle"))
-
-    # collect central
-    df = pd.DataFrame()
-    fk = pineappl.fk_table.FkTable.read(
-        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)
-        if not is_dis:
-            df[f"b{bd} r"] = fk.bin_right(bd)
-    df["central"] = fk.convolute_with_one(pdgid, pdf.xfxQ2)
-
-    # collect SVs
-    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 theories.keys():
-            del df[tn]
-
-    print(df.to_string())
-
-
-if __name__ == "__main__":
-    run()
diff --git a/yadism_debug.py b/yadism_debug.py
deleted file mode 100644
index 8404e5dd..00000000
--- a/yadism_debug.py
+++ /dev/null
@@ -1,72 +0,0 @@
-# -*- 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)])