From 6edfc484987098acf4b2b8388c78ebf8c30f4f81 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Sat, 18 Nov 2023 09:46:05 +0100 Subject: [PATCH 1/7] add FFNS to apfel++ --- .../benchmark/external/apfelpy_utils.py | 122 +++++++++++++----- 1 file changed, 89 insertions(+), 33 deletions(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index 72c48d58..85949332 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -17,6 +17,69 @@ } +def map_apfelpy_sf(init, observables, fns): + """Get required structure function. + + Parameters + ---------- + init: ap.initializers + Apfel++ initializers + observables: dict + observables runcard + fns: str + flavor number scheme + + Returns + ------- + Apfel++ structure function initalizer + """ + MAP_ZM_CC = { + "F2m": init.InitializeF2CCMinusObjectsZM, + "FLm": init.InitializeFLCCMinusObjectsZM, + "F3m": init.InitializeF3CCMinusObjectsZM, + "F2p": init.InitializeF2CCPlusObjectsZM, + "FLp": init.InitializeFLCCPlusObjectsZM, + "F3p": init.InitializeF3CCPlusObjectsZM, + } + + MAP_ZM_NC = { + "F2": init.InitializeF2NCObjectsZM, + "FL": init.InitializeFLNCObjectsZM, + "F3": init.InitializeF3NCObjectsZM, + "g1": init.Initializeg1NCObjectsZM, + "gL": init.InitializegLNCObjectsZM, + "g4": init.Initializeg4NCObjectsZM, + } + + MAP_FFNS_NC = { + "F2": init.InitializeF2NCObjectsMassive, + "FL": init.InitializeFLNCObjectsMassive, + } + MAP_FFNS0_NC = { + "F2": init.InitializeF2NCObjectsMassiveZero, + "FL": init.InitializeFLNCObjectsMassiveZero, + } + + if observables["prDIS"] == "CC": + projectile_pids = { + "electron": 11, + "positron": -11, + "neutrino": 12, + "antineutrino": -12, + } + if projectile_pids[observables["ProjectileDIS"]] > 0: + apfelpy_structure_functions = MAP_ZM_CC + else: + if fns == "ZM-VFNS": + apfelpy_structure_functions = MAP_ZM_NC + elif fns == "FFNS": + apfelpy_structure_functions = MAP_FFNS_NC + elif fns == "FFNS0": + apfelpy_structure_functions = MAP_FFNS0_NC + + return apfelpy_structure_functions + + def couplings(ap, pids, proc_type, obs_name): """Return the corresponding coupling given a process type and an observable. @@ -226,55 +289,42 @@ def compute_apfelpy_data(theory, observables, pdf): theory["mb"] * theory["kbThr"], theory["mt"] * theory["ktThr"], ] + masses = [ + 0, + 0, + 0, + theory["mc"], + theory["mb"], + theory["mt"], + ] # Setting the theory fns = theory["FNS"] - if fns != "ZM-VFNS": - raise ValueError("APFEL++ benchmark contains only ZM-VFNS.") + if fns not in ["ZM-VFNS", "FFNS", "FFNS0"]: + raise ValueError(f"APFEL++ does not contain {fns}.") # Perturbative Order pto = theory["PTO"] # Couplings - alphas = ap.AlphaQCD(theory["alphas"], theory["Qref"], thrs, pto) + alphas = ap.AlphaQCD(theory["alphas"], theory["Qref"], masses, thrs, pto) alphas = ap.TabulateObject(alphas, 2 * NQ, QMin, QMax, 3) # Map Yadism observables to Apfel++ Objects - init = ap.initializers - if observables["prDIS"] == "CC": - projectile_pids = { - "electron": 11, - "positron": -11, - "neutrino": 12, - "antineutrino": -12, - } - if projectile_pids[observables["ProjectileDIS"]] > 0: - apfelpy_structure_functions = { - "F2m": init.InitializeF2CCMinusObjectsZM, - "FLm": init.InitializeFLCCMinusObjectsZM, - "F3m": init.InitializeF3CCMinusObjectsZM, - "F2p": init.InitializeF2CCPlusObjectsZM, - "FLp": init.InitializeFLCCPlusObjectsZM, - "F3p": init.InitializeF3CCPlusObjectsZM, - } - else: - apfelpy_structure_functions = { - "F2": init.InitializeF2NCObjectsZM, - "FL": init.InitializeFLNCObjectsZM, - "F3": init.InitializeF3NCObjectsZM, - "g1": init.Initializeg1NCObjectsZM, - "gL": init.InitializegLNCObjectsZM, - "g4": init.Initializeg4NCObjectsZM, - } + apfelpy_structure_functions = map_apfelpy_sf(ap.initializers, observables, fns) # Initialize DGLAP Object - dglapobj = ap.initializers.InitializeDglapObjectsQCD(xgrid, thrs) + dglapobj = ap.initializers.InitializeDglapObjectsQCD(xgrid, masses, thrs) apf_tab = {} for obs_name, kinematics in observables["observables"].items(): apf_tab[obs_name] = [] sf_name, heaviness = obs_name.split("_") + if fns != "ZM-VFNS" and heaviness == "total": + raise ValueError( + "total is not provided in APFEL++ for massive calculations." + ) pids = MAP_HEAVINESS[heaviness] # Define the couplings @@ -303,11 +353,12 @@ def compute_apfelpy_data(theory, observables, pdf): tabulatedPDFs = ap.TabulateObjectSetD(evolvedPDFs, NQ, QMin, QMax, 3) # Initialize structure functions + eff_thrs = thrs if fns == "ZM_VFNS" else masses if observables["prDIS"] == "CC": sfobj = [] for sign in ["p", "m"]: tmp_sf = apfelpy_structure_functions[f"{sf_name}{sign}"]( - xgrid, thrs + xgrid, eff_thrs ) sfobj.append( ap.builders.BuildStructureFunctions( @@ -322,7 +373,7 @@ def compute_apfelpy_data(theory, observables, pdf): ) tab_sf = tabulate_cc(ap, obs_name, sfobj, NQ, QMin, QMax, thrs) else: - sfobj = apfelpy_structure_functions[sf_name](xgrid, thrs) + sfobj = apfelpy_structure_functions[sf_name](xgrid, eff_thrs) sfobj = ap.builders.BuildStructureFunctions( sfobj, tabulatedPDFs.EvaluateMapxQ, @@ -334,8 +385,13 @@ def compute_apfelpy_data(theory, observables, pdf): ) tab_sf = tabulate_nc(ap, obs_name, sfobj, NQ, QMin, QMax, thrs) + # shift convolution for massive + eta = 1.0 + if fns == "FFNS": + m_h = masses[3] if heaviness == "charm" else masses[4] + eta = Q2 / (Q2 + 4 * m_h**2) # compute the actual result - result = tab_sf.EvaluatexQ(x, np.sqrt(Q2)) + result = tab_sf.EvaluatexQ(x / eta, np.sqrt(Q2)) # take over the kinematics r = kin.copy() From 781ed78e1926c9d1592aa07e3831788351d526e5 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Mon, 20 Nov 2023 11:19:00 +0100 Subject: [PATCH 2/7] improve name --- src/yadmark/benchmark/external/apfelpy_utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index 85949332..9e16f126 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -3,9 +3,9 @@ from eko import basis_rotation as br # Q2 knots specs -NQ = 250 +NQ = 1200 QMin = 1 -QMax = 200 +QMax = 100 # Map observables names to APFEL++ methods MAP_HEAVINESS = { @@ -386,12 +386,12 @@ def compute_apfelpy_data(theory, observables, pdf): tab_sf = tabulate_nc(ap, obs_name, sfobj, NQ, QMin, QMax, thrs) # shift convolution for massive - eta = 1.0 + x_eval = 1.0 if fns == "FFNS": m_h = masses[3] if heaviness == "charm" else masses[4] - eta = Q2 / (Q2 + 4 * m_h**2) + x_eval = Q2 / (Q2 + 4 * m_h**2) # compute the actual result - result = tab_sf.EvaluatexQ(x / eta, np.sqrt(Q2)) + result = tab_sf.EvaluatexQ(x / x_eval, np.sqrt(Q2)) # take over the kinematics r = kin.copy() From 1dd9fcd3218f2f5cc129e8bb0296ddf7e8632241 Mon Sep 17 00:00:00 2001 From: Giacomo Magni <39065935+giacomomagni@users.noreply.github.com> Date: Mon, 20 Nov 2023 11:37:34 +0100 Subject: [PATCH 3/7] Update src/yadmark/benchmark/external/apfelpy_utils.py Co-authored-by: Alessandro Candido --- src/yadmark/benchmark/external/apfelpy_utils.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index 9e16f126..f8c99aac 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -68,16 +68,14 @@ def map_apfelpy_sf(init, observables, fns): "antineutrino": -12, } if projectile_pids[observables["ProjectileDIS"]] > 0: - apfelpy_structure_functions = MAP_ZM_CC + return MAP_ZM_CC else: if fns == "ZM-VFNS": - apfelpy_structure_functions = MAP_ZM_NC - elif fns == "FFNS": - apfelpy_structure_functions = MAP_FFNS_NC - elif fns == "FFNS0": - apfelpy_structure_functions = MAP_FFNS0_NC - - return apfelpy_structure_functions + return MAP_ZM_NC + if fns == "FFNS": + return MAP_FFNS_NC + if fns == "FFNS0": + return MAP_FFNS0_NC def couplings(ap, pids, proc_type, obs_name): From b853d88caf9d2880689233371b0fde09d12944cf Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Mon, 20 Nov 2023 11:54:05 +0100 Subject: [PATCH 4/7] address some style comments --- .../benchmark/external/apfelpy_utils.py | 58 ++++++++++--------- 1 file changed, 30 insertions(+), 28 deletions(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index f8c99aac..579e6a2c 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -3,9 +3,9 @@ from eko import basis_rotation as br # Q2 knots specs -NQ = 1200 -QMin = 1 -QMax = 100 +NQ = 250 +QMIN = 1 +QMAX = 200 # Map observables names to APFEL++ methods MAP_HEAVINESS = { @@ -16,6 +16,13 @@ "total": 0, } +PROJECTILE_PIDS = { + "electron": 11, + "positron": -11, + "neutrino": 12, + "antineutrino": -12, +} + def map_apfelpy_sf(init, observables, fns): """Get required structure function. @@ -61,21 +68,16 @@ def map_apfelpy_sf(init, observables, fns): } if observables["prDIS"] == "CC": - projectile_pids = { - "electron": 11, - "positron": -11, - "neutrino": 12, - "antineutrino": -12, - } - if projectile_pids[observables["ProjectileDIS"]] > 0: + if PROJECTILE_PIDS[observables["ProjectileDIS"]] > 0: return MAP_ZM_CC - else: - if fns == "ZM-VFNS": - return MAP_ZM_NC - if fns == "FFNS": - return MAP_FFNS_NC - if fns == "FFNS0": - return MAP_FFNS0_NC + return ValueError(f"{observables['ProjectileDIS']} not available in Apfel++") + + if fns == "ZM-VFNS": + return MAP_ZM_NC + if fns == "FFNS": + return MAP_FFNS_NC + if fns == "FFNS0": + return MAP_FFNS0_NC def couplings(ap, pids, proc_type, obs_name): @@ -137,9 +139,9 @@ def tabulate_nc(ap, obs_name, sfobj, nq, qmin, qmax, thrs): SF objects in Zero-Mass Flavour Number Scheme nq: int number of Q points - QMin: float + qmin: float minimal value of Q - QMax: float + qmax: float maximal value of Q thrs: list list of quark masses and thresholds @@ -180,9 +182,9 @@ def tabulate_cc(ap, obs_name, sfobj, nq, qmin, qmax, thrs): list of SF objects in Zero-Mass Flavour Number Scheme nq: int number of Q points - QMin: float + qmin: float minimal value of Q - QMax: float + qmax: float maximal value of Q thrs: list list of quark masses and thresholds @@ -306,7 +308,7 @@ def compute_apfelpy_data(theory, observables, pdf): # Couplings alphas = ap.AlphaQCD(theory["alphas"], theory["Qref"], masses, thrs, pto) - alphas = ap.TabulateObject(alphas, 2 * NQ, QMin, QMax, 3) + alphas = ap.TabulateObject(alphas, 2 * NQ, QMIN, QMAX, 3) # Map Yadism observables to Apfel++ Objects apfelpy_structure_functions = map_apfelpy_sf(ap.initializers, observables, fns) @@ -348,7 +350,7 @@ def compute_apfelpy_data(theory, observables, pdf): alphas.Evaluate, ) # Tabulate PDFs - tabulatedPDFs = ap.TabulateObjectSetD(evolvedPDFs, NQ, QMin, QMax, 3) + tabulatedPDFs = ap.TabulateObjectSetD(evolvedPDFs, NQ, QMIN, QMAX, 3) # Initialize structure functions eff_thrs = thrs if fns == "ZM_VFNS" else masses @@ -369,7 +371,7 @@ def compute_apfelpy_data(theory, observables, pdf): theory["XIF"], ) ) - tab_sf = tabulate_cc(ap, obs_name, sfobj, NQ, QMin, QMax, thrs) + tab_sf = tabulate_cc(ap, obs_name, sfobj, NQ, QMIN, QMAX, thrs) else: sfobj = apfelpy_structure_functions[sf_name](xgrid, eff_thrs) sfobj = ap.builders.BuildStructureFunctions( @@ -381,15 +383,15 @@ def compute_apfelpy_data(theory, observables, pdf): theory["XIR"], theory["XIF"], ) - tab_sf = tabulate_nc(ap, obs_name, sfobj, NQ, QMin, QMax, thrs) + tab_sf = tabulate_nc(ap, obs_name, sfobj, NQ, QMIN, QMAX, thrs) # shift convolution for massive - x_eval = 1.0 + x_eval = x if fns == "FFNS": m_h = masses[3] if heaviness == "charm" else masses[4] - x_eval = Q2 / (Q2 + 4 * m_h**2) + x_eval = x / (Q2 / (Q2 + 4 * m_h**2)) # compute the actual result - result = tab_sf.EvaluatexQ(x / x_eval, np.sqrt(Q2)) + result = tab_sf.EvaluatexQ(x_eval, np.sqrt(Q2)) # take over the kinematics r = kin.copy() From 823a0b1523b68736b4cf0496badf5406e6c339d1 Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Mon, 20 Nov 2023 13:09:48 +0200 Subject: [PATCH 5/7] Update src/yadmark/benchmark/external/apfelpy_utils.py --- src/yadmark/benchmark/external/apfelpy_utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index 579e6a2c..aa555d76 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -389,7 +389,8 @@ def compute_apfelpy_data(theory, observables, pdf): x_eval = x if fns == "FFNS": m_h = masses[3] if heaviness == "charm" else masses[4] - x_eval = x / (Q2 / (Q2 + 4 * m_h**2)) + eta = Q2 / (Q2 + 4 * m_h**2) + x_eval = x / eta # compute the actual result result = tab_sf.EvaluatexQ(x_eval, np.sqrt(Q2)) From 374d39cfa403e5fcadb0b74291eef69f9010f20d Mon Sep 17 00:00:00 2001 From: Giacomo Magni <39065935+giacomomagni@users.noreply.github.com> Date: Mon, 20 Nov 2023 12:54:12 +0100 Subject: [PATCH 6/7] Update src/yadmark/benchmark/external/apfelpy_utils.py Co-authored-by: Alessandro Candido --- src/yadmark/benchmark/external/apfelpy_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index aa555d76..c6602b8d 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -70,7 +70,7 @@ def map_apfelpy_sf(init, observables, fns): if observables["prDIS"] == "CC": if PROJECTILE_PIDS[observables["ProjectileDIS"]] > 0: return MAP_ZM_CC - return ValueError(f"{observables['ProjectileDIS']} not available in Apfel++") + raise ValueError(f"{observables['ProjectileDIS']} not available in Apfel++") if fns == "ZM-VFNS": return MAP_ZM_NC From b64de00c50377cf5c1c857bb3e5026cd7321b331 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Mon, 20 Nov 2023 22:02:25 +0100 Subject: [PATCH 7/7] small fixes --- src/yadmark/benchmark/external/apfelpy_utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/yadmark/benchmark/external/apfelpy_utils.py b/src/yadmark/benchmark/external/apfelpy_utils.py index c6602b8d..43200b93 100644 --- a/src/yadmark/benchmark/external/apfelpy_utils.py +++ b/src/yadmark/benchmark/external/apfelpy_utils.py @@ -387,8 +387,8 @@ def compute_apfelpy_data(theory, observables, pdf): # shift convolution for massive x_eval = x - if fns == "FFNS": - m_h = masses[3] if heaviness == "charm" else masses[4] + if fns == "FFNS" and heaviness != "light": + m_h = masses[MAP_HEAVINESS[heaviness] - 1] eta = Q2 / (Q2 + 4 * m_h**2) x_eval = x / eta # compute the actual result