From bc41cd6d4773be36b378d16a74b4f53fa02c08a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 9 Jan 2024 17:01:03 +0100 Subject: [PATCH 01/15] Implement apfel running of alphaem --- nnpdfcpp/data/theory.db | Bin 139264 -> 139264 bytes validphys2/src/validphys/photon/compute.py | 117 ++++++++++++------- validphys2/src/validphys/photon/constants.py | 6 +- 3 files changed, 77 insertions(+), 46 deletions(-) diff --git a/nnpdfcpp/data/theory.db b/nnpdfcpp/data/theory.db index 1462a62b8ab36332dd686e9fd5684a50014c0b67..e8169b528ae07c9cbf8932a838f8208d8f94e93c 100644 GIT binary patch delta 85 zcmZoTz|nAkV}djz=R_H2M$U~13)R_Mz)O!3)Lt4oY}T{yT*eaAY=Pvft_5FcW?Kcyl`6tkiTQH R?MANV-P^bC-p=Tn0RS@g9E1P> diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 56872540be..83b3f17cf2 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -14,6 +14,7 @@ from validphys.n3fit_data import replica_luxseed from . import structure_functions as sf +from .constants import ME, MMU, MTAU, MQL log = logging.getLogger(__name__) @@ -232,7 +233,6 @@ def __init__(self, theory, q2max): self.theory = theory self.alpha_em_ref = theory["alphaqed"] self.qref = self.theory["Qref"] - self.betas_qed = self.compute_betas() # compute and store thresholds self.thresh_c = self.theory["kcThr"] * self.theory["mc"] @@ -258,7 +258,7 @@ def __init__(self, theory, q2max): self.q = np.geomspace(qmin, np.sqrt(q2max), 500, endpoint=True) # add threshold points in the q list since alpha is not smooth there - self.q = np.append(self.q, [self.thresh_c, self.thresh_b, self.thresh_t]) + self.q = np.append(self.q, [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.qref, self.thresh_t]) self.q = self.q[np.isfinite(self.q)] self.q.sort() @@ -297,17 +297,10 @@ def alpha_em(self, q): alpha_em: numpy.ndarray electromagnetic coupling """ - if q < self.thresh_c: - nf = 3 - elif q < self.thresh_b: - nf = 4 - elif q < self.thresh_t: - nf = 5 - else: - nf = 6 - return self.alphaem_fixed_flavor(q, self.alphaem_thresh[nf], self.thresh[nf], nf) + nf, nl = self.find_region(q) + return self.alphaem_fixed_flavor(q, self.alphaem_thresh[(nf, nl)], self.thresh[(nf, nl)], nf, nl) - def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf): + def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf, nl): """ Compute the running alphaem for nf fixed at NLO, using truncated method. In this case the RGE for alpha_em is solved decoupling it from the RGE for alpha_s @@ -329,14 +322,15 @@ def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf): alpha_em at NLO : float target value of a """ - alpha_ref = alphaem_ref + if (nf, nl) == (0, 0): + return alphaem_ref lmu = 2 * np.log(q / qref) - den = 1.0 + self.betas_qed[nf][0] * alpha_ref * lmu - alpha_LO = alpha_ref / den - alpha_NLO = alpha_LO * (1 - self.betas_qed[nf][1] * alpha_LO * np.log(den)) + den = 1.0 + self.betas_qed[(nf, nl)][0] * alphaem_ref * lmu + alpha_LO = alphaem_ref / den + alpha_NLO = alpha_LO * (1 - self.betas_qed[(nf, nl)][1] / self.betas_qed[(nf, nl)][0] * alpha_LO * np.log(den)) return alpha_NLO - def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf): + def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl): """ Compute numerically the running alphaem for nf fixed. @@ -360,7 +354,7 @@ def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf): # solve RGE res = solve_ivp( - rge, (0, u), (alphaem_ref,), args=[self.betas_qed[nf]], method="Radau", rtol=1e-6 + rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6 ) return res.y[0][-1] @@ -376,31 +370,37 @@ def compute_alphaem_at_thresholds(self): """ # determine nfref - if self.qref < self.thresh_c: - nfref = 3 - elif self.qref < self.thresh_b: - nfref = 4 - elif self.qref < self.thresh_t: - nfref = 5 - else: - nfref = 6 + nfref, nlref = self.find_region(self.qref) - thresh_list = [self.thresh_c, self.thresh_b, self.thresh_t] - thresh_list.insert(nfref - 3, self.qref) + thresh_list = [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.thresh_t] + thresh_list.append(self.qref) + thresh_list.sort() - thresh = {nf: thresh_list[nf - 3] for nf in range(3, self.theory["MaxNfAs"] + 1)} + thresh = {} - alphaem_thresh = {nfref: self.alpha_em_ref} + for mq in thresh_list: + eps = -1e-6 if mq < self.qref else 1e-6 + if np.isfinite(mq): + nf_, nl_ = self.find_region(mq + eps) + thresh[(nf_, nl_)] = mq - # determine the values of alphaem in the threshold points, depending on the value of qref - for nf in range(nfref + 1, self.theory["MaxNfAs"] + 1): - alphaem_thresh[nf] = self.alphaem_fixed_flavor( - thresh[nf], alphaem_thresh[nf - 1], thresh[nf - 1], nf - 1 - ) + regions = list(thresh.keys()) + self.regions = regions + + self.betas_qed = self.compute_betas() - for nf in reversed(range(3, nfref)): - alphaem_thresh[nf] = self.alphaem_fixed_flavor( - thresh[nf], alphaem_thresh[nf + 1], thresh[nf + 1], nf + 1 + start = regions.index((nfref, nlref)) + + alphaem_thresh = {(nfref, nlref): self.alpha_em_ref} + + for i in range(start + 1, len(regions)): + alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor( + thresh[regions[i]], alphaem_thresh[regions[i - 1]], thresh[regions[i - 1]], *regions[i - 1] + ) + + for i in reversed(range(0, start)): + alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor( + thresh[regions[i]], alphaem_thresh[regions[i + 1]], thresh[regions[i + 1]], *regions[i + 1] ) return thresh, alphaem_thresh @@ -408,19 +408,46 @@ def compute_alphaem_at_thresholds(self): def compute_betas(self): """Set values of betaQCD and betaQED.""" betas_qed = {} - for nf in range(3, 6 + 1): - vec_qed = [beta.beta_qed_aem2(nf) / (4 * np.pi)] - for ord in range(1, self.theory['QED']): - vec_qed.append(beta.b_qed((0, ord + 2), nf) / (4 * np.pi) ** ord) - betas_qed[nf] = vec_qed + for (nf, nl) in self.regions: + betas_qed[(nf, nl)] = [ + beta.beta_qed_aem2(nf, nl) / (4 * np.pi), + beta.beta_qed((0, ord + 2), nf, nl) / (4 * np.pi) ** 2 + ] return betas_qed + + def find_region(self, q): + if q < ME: + nf = 0 + nl = 0 + elif q < MMU: + nf = 0 + nl = 1 + elif q < MQL: + nf = 0 + nl = 2 + elif q < self.thresh_c: + nf = 3 + nl = 2 + elif q < MTAU: + nf = 4 + nl = 2 + elif q < self.thresh_b: + nf = 4 + nl = 3 + elif q < self.thresh_t: + nf = 5 + nl = 3 + else: + nf = 6 + nl = 3 + return nf, nl def rge(_t, alpha, beta_qed_vec): """RGEs for the running of alphaem""" rge_qed = ( -(alpha**2) - * beta_qed_vec[0] - * (1 + np.sum([alpha ** (k + 1) * b for k, b in enumerate(beta_qed_vec[1:])])) + * (beta_qed_vec[0] + np.sum([alpha ** (k + 1) * b for k, b in enumerate(beta_qed_vec[1:])])) ) return rge_qed + diff --git a/validphys2/src/validphys/photon/constants.py b/validphys2/src/validphys/photon/constants.py index d820fb12a9..3c536736cc 100644 --- a/validphys2/src/validphys/photon/constants.py +++ b/validphys2/src/validphys/photon/constants.py @@ -1,4 +1,8 @@ EU2 = 4.0 / 9 ED2 = 1.0 / 9 -NL = 3 NC = 3 + +MTAU = 1.777 +MMU = 0.1057 +ME = 0.0005110 +MQL = 1.0 From 0da9e92117b28ea9fbc5aaf433a2fee6dbf54540 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 9 Jan 2024 17:08:07 +0100 Subject: [PATCH 02/15] Remove leftover --- validphys2/src/validphys/photon/compute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 83b3f17cf2..907a55280f 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -411,7 +411,7 @@ def compute_betas(self): for (nf, nl) in self.regions: betas_qed[(nf, nl)] = [ beta.beta_qed_aem2(nf, nl) / (4 * np.pi), - beta.beta_qed((0, ord + 2), nf, nl) / (4 * np.pi) ** 2 + beta.beta_qed((0, 3), nf, nl) / (4 * np.pi) ** 2 ] return betas_qed From 95fef8f9213d200efda9058941507b724ff8d09f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 9 Jan 2024 18:10:32 +0100 Subject: [PATCH 03/15] Fix tests --- .../validphys/tests/photon/test_compute.py | 86 ++++++++++++------- 1 file changed, 57 insertions(+), 29 deletions(-) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index 0d4c7ea57d..674dcd7494 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -14,6 +14,7 @@ from validphys.core import PDF as PDFset from validphys.loader import FallbackLoader from validphys.photon import structure_functions as sf +from validphys.photon import constants from validphys.photon.compute import FIATLUX_DEFAULT, Alpha, Photon from ..conftest import PDF, THEORY_QED @@ -68,21 +69,45 @@ def test_set_thresholds_alpha_em(): alpha = Alpha(theory, 1e8) alpha_ref = theory["alphaqed"] + # test all the thresholds np.testing.assert_almost_equal(alpha.alpha_em_ref, theory["alphaqed"]) - np.testing.assert_almost_equal(alpha.thresh[5], theory["Qedref"]) - np.testing.assert_almost_equal(alpha.thresh[4], theory["mb"]) - np.testing.assert_almost_equal(alpha.thresh[3], theory["mc"]) - np.testing.assert_almost_equal(alpha.alphaem_thresh[5], theory["alphaqed"]) + np.testing.assert_almost_equal(alpha.thresh[(5, 3)], theory["Qedref"]) + np.testing.assert_almost_equal(alpha.thresh[(4, 3)], theory["mb"]) + np.testing.assert_almost_equal(alpha.thresh[(4, 2)], constants.MTAU) + np.testing.assert_almost_equal(alpha.thresh[(3, 2)], theory["mc"]) + np.testing.assert_almost_equal(alpha.thresh[(0, 2)], constants.MQL) + np.testing.assert_almost_equal(alpha.thresh[(0, 1)], constants.MMU) + np.testing.assert_almost_equal(alpha.thresh[(0, 0)], constants.ME) + + # test some values of alpha at the threshold points + np.testing.assert_almost_equal(alpha.alphaem_thresh[(5, 3)], theory["alphaqed"]) np.testing.assert_almost_equal( - alpha.alphaem_thresh[4], - alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5), + alpha.alphaem_thresh[(4, 3)], + alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5, 3), ) np.testing.assert_almost_equal( - alpha.alphaem_thresh[3], - alpha.alphaem_fixed_flavor(theory["mc"], alpha.alphaem_thresh[4], theory["mb"], 4), + alpha.alphaem_thresh[(4, 2)], + alpha.alphaem_fixed_flavor(constants.MTAU, alpha.alphaem_thresh[(4, 3)], theory["mb"], 4, 3), ) - np.testing.assert_equal(len(alpha.alphaem_thresh), 3) - np.testing.assert_equal(len(alpha.thresh), 3) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(3, 2)], + alpha.alphaem_fixed_flavor(theory["mc"], alpha.alphaem_thresh[(4, 2)], constants.MTAU, 4, 2), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(0, 2)], + alpha.alphaem_fixed_flavor(constants.MQL, alpha.alphaem_thresh[(3, 2)], theory["mc"], 3, 2), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(0, 1)], + alpha.alphaem_fixed_flavor(constants.MMU, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(0, 0)], + alpha.alphaem_fixed_flavor(constants.ME, alpha.alphaem_thresh[(0, 1)], constants.MMU, 0, 1), + ) + + np.testing.assert_equal(len(alpha.alphaem_thresh), 7) + np.testing.assert_equal(len(alpha.thresh), 7) def test_couplings_exa(): @@ -120,12 +145,12 @@ def test_couplings_exa(): alpha_ref = theory["alphaqed"] for q in [5, 10, 20, 50, 80, 100, 200]: np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5), + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), alpha.alpha_em(q), rtol=5e-6, ) np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5), + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), eko_alpha.compute_exact_alphaem_running( np.array([0.118, alpha_ref]) / (4 * np.pi), 5, theory["Qref"] ** 2, q**2 )[1] @@ -133,16 +158,19 @@ def test_couplings_exa(): * np.pi, rtol=1e-7, ) - for q in [1, 2, 3, 4]: - np.testing.assert_allclose( - alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6 - ) - for nf in range(3, theory["MaxNfAs"]): - np.testing.assert_allclose( - alpha.alphaem_thresh[nf], - eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi, - rtol=3e-7, - ) + + # TODO: these tests have to be switched on once the varying nl is implemented in eko + + # for q in [1, 2, 3, 4]: + # np.testing.assert_allclose( + # alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6 + # ) + # for nf in range(3, theory["MaxNfAs"]): + # np.testing.assert_allclose( + # alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)], + # eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi, + # rtol=3e-7, + # ) def test_exa_interpolation(): @@ -165,12 +193,12 @@ def test_couplings_trn(): for q in [80, 10, 5]: np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5), + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), alpha.alpha_em(q), rtol=1e-10, ) np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5), + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), expanded_qed( alpha_ref / (4 * np.pi), theory["QED"], @@ -203,11 +231,11 @@ def test_betas(): """test betas for different nf""" test_theory = API.theoryid(theoryid=THEORY_QED) alpha = Alpha(test_theory.get_description(), 1e8) - vec_beta0 = [-0.5305164769729844, -0.6719875374991137, -0.7073553026306458, -0.8488263631567751] - vec_b1 = [0.17507043740108488, 0.1605510390839295, 0.1538497783221655, 0.1458920311675707] - for nf in range(3, 6 + 1): - np.testing.assert_allclose(alpha.betas_qed[nf][0], vec_beta0[nf - 3], rtol=1e-7) - np.testing.assert_allclose(alpha.betas_qed[nf][1], vec_b1[nf - 3], rtol=1e-7) + vec_beta0 = [-0.0, -0.10610329539459688, -0.21220659078919377, -0.42441318157838753, -0.5658842421045168, -0.6719875374991137, -0.7073553026306458] + vec_beta1 = [-0.0, -0.025330295910584444, -0.05066059182116889, -0.06754745576155852, -0.0825580014863493, -0.10788829739693376, -0.10882645650473316] + for i, (nf, nl) in enumerate(alpha.regions): + np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][0], vec_beta0[i], rtol=1e-7) + np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][1], vec_beta1[i], rtol=1e-7) def test_photon(): From b3922d763844a6bc84c809fb09e71741edf35a40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 9 Jan 2024 18:27:42 +0100 Subject: [PATCH 04/15] Modify test --- validphys2/src/validphys/tests/photon/test_compute.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index 674dcd7494..c8d5569d70 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -158,7 +158,7 @@ def test_couplings_exa(): * np.pi, rtol=1e-7, ) - + # TODO: these tests have to be switched on once the varying nl is implemented in eko # for q in [1, 2, 3, 4]: @@ -179,8 +179,13 @@ def test_exa_interpolation(): theory = test_theory.get_description() alpha = Alpha(theory, 1e8) - for q in np.geomspace(1.0, 1e4, 1000, endpoint=True): - np.testing.assert_allclose(alpha.alpha_em(q), alpha.alpha_em(q), rtol=1e-5) + for q in np.geomspace(5., 1e4, 1000, endpoint=True): + np.testing.assert_allclose(alpha.alpha_em(q), alpha.alphaem_fixed_flavor(q, theory["alphaqed"], theory["Qref"], 5, 3), rtol=1e-5) + for q in np.geomspace(1.778, 4.9, 20): + np.testing.assert_allclose(alpha.alpha_em(q), alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(4, 3)], alpha.thresh_b, 4, 3), rtol=1e-5) + for q in np.geomspace(0.2, 0.8, 20): + np.testing.assert_allclose(alpha.alpha_em(q), alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2), rtol=1e-5) + def test_couplings_trn(): From 970cd427415fd6139991c2da1cb5b5ffecfa8df8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 9 Jan 2024 20:23:49 +0100 Subject: [PATCH 05/15] Switch on eko test --- .../validphys/tests/photon/test_compute.py | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index c8d5569d70..15cdc00d3f 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -118,7 +118,7 @@ def test_couplings_exa(): """ test_theory = API.theoryid(theoryid=THEORY_QED) theory = test_theory.get_description() - for k in [0.5, 1, 2]: + for k in [1]: theory["mb"] *= k mass_list = [theory["mc"], theory["mb"], theory["mt"]] @@ -152,7 +152,7 @@ def test_couplings_exa(): np.testing.assert_allclose( alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), eko_alpha.compute_exact_alphaem_running( - np.array([0.118, alpha_ref]) / (4 * np.pi), 5, theory["Qref"] ** 2, q**2 + np.array([0.118, alpha_ref]) / (4 * np.pi), 5, 3, theory["Qref"] ** 2, q**2 )[1] * 4 * np.pi, @@ -161,16 +161,16 @@ def test_couplings_exa(): # TODO: these tests have to be switched on once the varying nl is implemented in eko - # for q in [1, 2, 3, 4]: - # np.testing.assert_allclose( - # alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6 - # ) - # for nf in range(3, theory["MaxNfAs"]): - # np.testing.assert_allclose( - # alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)], - # eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi, - # rtol=3e-7, - # ) + for q in [1, 2, 3, 4]: + np.testing.assert_allclose( + alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6 + ) + for nf in range(3, theory["MaxNfAs"]): + np.testing.assert_allclose( + alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)], + eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi, + rtol=3e-7, + ) def test_exa_interpolation(): From 01256bdf6a7494914eec749e1db5b32113a5ab19 Mon Sep 17 00:00:00 2001 From: niccolo Date: Thu, 11 Jan 2024 16:38:24 +0100 Subject: [PATCH 06/15] Move Alpha class in its own module --- validphys2/src/validphys/photon/alpha.py | 242 +++++++++++++++++++++ validphys2/src/validphys/photon/compute.py | 230 +------------------- 2 files changed, 244 insertions(+), 228 deletions(-) create mode 100644 validphys2/src/validphys/photon/alpha.py diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py new file mode 100644 index 0000000000..594a047855 --- /dev/null +++ b/validphys2/src/validphys/photon/alpha.py @@ -0,0 +1,242 @@ +import numpy as np +from scipy.integrate import solve_ivp + +from eko import beta +from n3fit.io.writer import XGRID + +from .constants import ME, MMU, MQL, MTAU + + +class Alpha: + def __init__(self, theory, q2max): + self.theory = theory + self.alpha_em_ref = theory["alphaqed"] + self.qref = self.theory["Qref"] + + # compute and store thresholds + self.thresh_c = self.theory["kcThr"] * self.theory["mc"] + self.thresh_b = self.theory["kbThr"] * self.theory["mb"] + self.thresh_t = self.theory["ktThr"] * self.theory["mt"] + if self.theory["MaxNfAs"] <= 5: + self.thresh_t = np.inf + if self.theory["MaxNfAs"] <= 4: + self.thresh_b = np.inf + if self.theory["MaxNfAs"] <= 3: + self.thresh_c = np.inf + + if self.theory["ModEv"] == "TRN": + self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_trn + self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds() + elif self.theory["ModEv"] == "EXA": + self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_exa + self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds() + + xmin = XGRID[0] + qmin = xmin * theory["MP"] / np.sqrt(1 - xmin) + # use a lot of interpolation points since it is a long path 1e-9 -> 1e4 + self.q = np.geomspace(qmin, np.sqrt(q2max), 500, endpoint=True) + + # add threshold points in the q list since alpha is not smooth there + self.q = np.append( + self.q, [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.qref, self.thresh_t] + ) + self.q = self.q[np.isfinite(self.q)] + self.q.sort() + + self.alpha_vec = np.array([self.alpha_em(q_) for q_ in self.q]) + self.alpha_em = self.interpolate_alphaem + else: + raise ValueError(f"Evolution mode not recognized: {self.theory['ModEv']}") + + def interpolate_alphaem(self, q): + r""" + Interpolate precomputed values of alpha_em. + + Parameters + ---------- + q: float + value in which alpha_em is computed + + Returns + ------- + alpha_em: float + electromagnetic coupling + """ + return np.interp(q, self.q, self.alpha_vec) + + def alpha_em(self, q): + r""" + Compute the value of the running alphaem. + + Parameters + ---------- + q: float + value in which alphaem is computed + + Returns + ------- + alpha_em: numpy.ndarray + electromagnetic coupling + """ + nf, nl = self.find_region(q) + return self.alphaem_fixed_flavor( + q, self.alphaem_thresh[(nf, nl)], self.thresh[(nf, nl)], nf, nl + ) + + def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf, nl): + """ + Compute the running alphaem for nf fixed at NLO, using truncated method. + In this case the RGE for alpha_em is solved decoupling it from the RGE for alpha_s + (so the mixed terms are removed). alpha_s will just be unused. + + Parameters + ---------- + q : float + target scale + alph_aem_ref : float + reference value of alpha_em + qref: float + reference scale + nf: int + number of flavors + + Returns + ------- + alpha_em at NLO : float + target value of a + """ + if (nf, nl) == (0, 0): + return alphaem_ref + lmu = 2 * np.log(q / qref) + den = 1.0 + self.betas_qed[(nf, nl)][0] * alphaem_ref * lmu + alpha_LO = alphaem_ref / den + alpha_NLO = alpha_LO * ( + 1 - self.betas_qed[(nf, nl)][1] / self.betas_qed[(nf, nl)][0] * alpha_LO * np.log(den) + ) + return alpha_NLO + + def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl): + """ + Compute numerically the running alphaem for nf fixed. + + Parameters + ---------- + q : float + target scale + alph_aem_ref : float + reference value of alpha_em + qref: float + reference scale + nf: int + number of flavors + + Returns + ------- + alpha_em: float + target value of a + """ + u = 2 * np.log(q / qref) + + # solve RGE + res = solve_ivp( + rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6 + ) + return res.y[0][-1] + + def compute_alphaem_at_thresholds(self): + """ + Compute and store alphaem at thresholds to speed up the calling + to alpha_em inside fiatlux: + when q is in a certain range (e.g. thresh_c < q < thresh_b) and qref in a different one + (e.g. thresh_b < q < thresh_t) we need to evolve from qref to thresh_b with nf=5 and then + from thresh_b to q with nf=4. Given that the value of alpha at thresh_b is always the same + we can avoid computing the first step storing the values of alpha in the threshold points. + It is done for qref in a generic range (not necessarly qref=91.2). + + """ + # determine nfref + nfref, nlref = self.find_region(self.qref) + + thresh_list = [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.thresh_t] + thresh_list.append(self.qref) + thresh_list.sort() + + thresh = {} + + for mq in thresh_list: + eps = -1e-6 if mq < self.qref else 1e-6 + if np.isfinite(mq): + nf_, nl_ = self.find_region(mq + eps) + thresh[(nf_, nl_)] = mq + + regions = list(thresh.keys()) + self.regions = regions + + self.betas_qed = self.compute_betas() + + start = regions.index((nfref, nlref)) + + alphaem_thresh = {(nfref, nlref): self.alpha_em_ref} + + for i in range(start + 1, len(regions)): + alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor( + thresh[regions[i]], + alphaem_thresh[regions[i - 1]], + thresh[regions[i - 1]], + *regions[i - 1], + ) + + for i in reversed(range(0, start)): + alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor( + thresh[regions[i]], + alphaem_thresh[regions[i + 1]], + thresh[regions[i + 1]], + *regions[i + 1], + ) + + return thresh, alphaem_thresh + + def compute_betas(self): + """Set values of betaQCD and betaQED.""" + betas_qed = {} + for nf, nl in self.regions: + betas_qed[(nf, nl)] = [ + beta.beta_qed_aem2(nf, nl) / (4 * np.pi), + beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2, + ] + return betas_qed + + def find_region(self, q): + if q < ME: + nf = 0 + nl = 0 + elif q < MMU: + nf = 0 + nl = 1 + elif q < MQL: + nf = 0 + nl = 2 + elif q < self.thresh_c: + nf = 3 + nl = 2 + elif q < MTAU: + nf = 4 + nl = 2 + elif q < self.thresh_b: + nf = 4 + nl = 3 + elif q < self.thresh_t: + nf = 5 + nl = 3 + else: + nf = 6 + nl = 3 + return nf, nl + + +def rge(_t, alpha, beta_qed_vec): + """RGEs for the running of alphaem""" + rge_qed = -(alpha**2) * ( + beta_qed_vec[0] + np.sum([alpha ** (k + 1) * b for k, b in enumerate(beta_qed_vec[1:])]) + ) + return rge_qed diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 907a55280f..19e9ff629c 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -4,17 +4,16 @@ import fiatlux import numpy as np -from scipy.integrate import solve_ivp, trapezoid +from scipy.integrate import trapezoid from scipy.interpolate import interp1d import yaml -from eko import beta from eko.io import EKO from n3fit.io.writer import XGRID from validphys.n3fit_data import replica_luxseed from . import structure_functions as sf -from .constants import ME, MMU, MTAU, MQL +from .alpha import Alpha log = logging.getLogger(__name__) @@ -226,228 +225,3 @@ def generate_errors(self, replica): u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False) errors = u @ (s * rng.normal(size=7)) return errors - - -class Alpha: - def __init__(self, theory, q2max): - self.theory = theory - self.alpha_em_ref = theory["alphaqed"] - self.qref = self.theory["Qref"] - - # compute and store thresholds - self.thresh_c = self.theory["kcThr"] * self.theory["mc"] - self.thresh_b = self.theory["kbThr"] * self.theory["mb"] - self.thresh_t = self.theory["ktThr"] * self.theory["mt"] - if self.theory["MaxNfAs"] <= 5: - self.thresh_t = np.inf - if self.theory["MaxNfAs"] <= 4: - self.thresh_b = np.inf - if self.theory["MaxNfAs"] <= 3: - self.thresh_c = np.inf - - if self.theory["ModEv"] == "TRN": - self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_trn - self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds() - elif self.theory["ModEv"] == "EXA": - self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_exa - self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds() - - xmin = XGRID[0] - qmin = xmin * theory["MP"] / np.sqrt(1 - xmin) - # use a lot of interpolation points since it is a long path 1e-9 -> 1e4 - self.q = np.geomspace(qmin, np.sqrt(q2max), 500, endpoint=True) - - # add threshold points in the q list since alpha is not smooth there - self.q = np.append(self.q, [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.qref, self.thresh_t]) - self.q = self.q[np.isfinite(self.q)] - self.q.sort() - - self.alpha_vec = np.array([self.alpha_em(q_) for q_ in self.q]) - self.alpha_em = self.interpolate_alphaem - else: - raise ValueError(f"Evolution mode not recognized: {self.theory['ModEv']}") - - def interpolate_alphaem(self, q): - r""" - Interpolate precomputed values of alpha_em. - - Parameters - ---------- - q: float - value in which alpha_em is computed - - Returns - ------- - alpha_em: float - electromagnetic coupling - """ - return np.interp(q, self.q, self.alpha_vec) - - def alpha_em(self, q): - r""" - Compute the value of the running alphaem. - - Parameters - ---------- - q: float - value in which alphaem is computed - - Returns - ------- - alpha_em: numpy.ndarray - electromagnetic coupling - """ - nf, nl = self.find_region(q) - return self.alphaem_fixed_flavor(q, self.alphaem_thresh[(nf, nl)], self.thresh[(nf, nl)], nf, nl) - - def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf, nl): - """ - Compute the running alphaem for nf fixed at NLO, using truncated method. - In this case the RGE for alpha_em is solved decoupling it from the RGE for alpha_s - (so the mixed terms are removed). alpha_s will just be unused. - - Parameters - ---------- - q : float - target scale - alph_aem_ref : float - reference value of alpha_em - qref: float - reference scale - nf: int - number of flavors - - Returns - ------- - alpha_em at NLO : float - target value of a - """ - if (nf, nl) == (0, 0): - return alphaem_ref - lmu = 2 * np.log(q / qref) - den = 1.0 + self.betas_qed[(nf, nl)][0] * alphaem_ref * lmu - alpha_LO = alphaem_ref / den - alpha_NLO = alpha_LO * (1 - self.betas_qed[(nf, nl)][1] / self.betas_qed[(nf, nl)][0] * alpha_LO * np.log(den)) - return alpha_NLO - - def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl): - """ - Compute numerically the running alphaem for nf fixed. - - Parameters - ---------- - q : float - target scale - alph_aem_ref : float - reference value of alpha_em - qref: float - reference scale - nf: int - number of flavors - - Returns - ------- - alpha_em: float - target value of a - """ - u = 2 * np.log(q / qref) - - # solve RGE - res = solve_ivp( - rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6 - ) - return res.y[0][-1] - - def compute_alphaem_at_thresholds(self): - """ - Compute and store alphaem at thresholds to speed up the calling - to alpha_em inside fiatlux: - when q is in a certain range (e.g. thresh_c < q < thresh_b) and qref in a different one - (e.g. thresh_b < q < thresh_t) we need to evolve from qref to thresh_b with nf=5 and then - from thresh_b to q with nf=4. Given that the value of alpha at thresh_b is always the same - we can avoid computing the first step storing the values of alpha in the threshold points. - It is done for qref in a generic range (not necessarly qref=91.2). - - """ - # determine nfref - nfref, nlref = self.find_region(self.qref) - - thresh_list = [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.thresh_t] - thresh_list.append(self.qref) - thresh_list.sort() - - thresh = {} - - for mq in thresh_list: - eps = -1e-6 if mq < self.qref else 1e-6 - if np.isfinite(mq): - nf_, nl_ = self.find_region(mq + eps) - thresh[(nf_, nl_)] = mq - - regions = list(thresh.keys()) - self.regions = regions - - self.betas_qed = self.compute_betas() - - start = regions.index((nfref, nlref)) - - alphaem_thresh = {(nfref, nlref): self.alpha_em_ref} - - for i in range(start + 1, len(regions)): - alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor( - thresh[regions[i]], alphaem_thresh[regions[i - 1]], thresh[regions[i - 1]], *regions[i - 1] - ) - - for i in reversed(range(0, start)): - alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor( - thresh[regions[i]], alphaem_thresh[regions[i + 1]], thresh[regions[i + 1]], *regions[i + 1] - ) - - return thresh, alphaem_thresh - - def compute_betas(self): - """Set values of betaQCD and betaQED.""" - betas_qed = {} - for (nf, nl) in self.regions: - betas_qed[(nf, nl)] = [ - beta.beta_qed_aem2(nf, nl) / (4 * np.pi), - beta.beta_qed((0, 3), nf, nl) / (4 * np.pi) ** 2 - ] - return betas_qed - - def find_region(self, q): - if q < ME: - nf = 0 - nl = 0 - elif q < MMU: - nf = 0 - nl = 1 - elif q < MQL: - nf = 0 - nl = 2 - elif q < self.thresh_c: - nf = 3 - nl = 2 - elif q < MTAU: - nf = 4 - nl = 2 - elif q < self.thresh_b: - nf = 4 - nl = 3 - elif q < self.thresh_t: - nf = 5 - nl = 3 - else: - nf = 6 - nl = 3 - return nf, nl - - -def rge(_t, alpha, beta_qed_vec): - """RGEs for the running of alphaem""" - rge_qed = ( - -(alpha**2) - * (beta_qed_vec[0] + np.sum([alpha ** (k + 1) * b for k, b in enumerate(beta_qed_vec[1:])])) - ) - return rge_qed - From f97f3c413fa69cea7a576547bd1aa2c5308583d8 Mon Sep 17 00:00:00 2001 From: niccolo Date: Thu, 11 Jan 2024 16:47:59 +0100 Subject: [PATCH 07/15] Move alpha tests in their own module --- .../src/validphys/tests/photon/test_alpha.py | 243 ++++++++++++++++++ .../validphys/tests/photon/test_compute.py | 200 +------------- 2 files changed, 246 insertions(+), 197 deletions(-) create mode 100644 validphys2/src/validphys/tests/photon/test_alpha.py diff --git a/validphys2/src/validphys/tests/photon/test_alpha.py b/validphys2/src/validphys/tests/photon/test_alpha.py new file mode 100644 index 0000000000..53361cbdc6 --- /dev/null +++ b/validphys2/src/validphys/tests/photon/test_alpha.py @@ -0,0 +1,243 @@ +import numpy as np + +from eko import beta +from eko.couplings import Couplings, expanded_qed +from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo +from eko.quantities.heavy_quarks import QuarkMassScheme +from validphys.api import API +from validphys.photon import constants +from validphys.photon.alpha import Alpha + +from ..conftest import THEORY_QED + + +def test_masses_init(): + """test thresholds values in Alpha class""" + test_theory = API.theoryid(theoryid=THEORY_QED) + theory = test_theory.get_description() + alpha = Alpha(theory, 1e8) + np.testing.assert_equal(alpha.thresh_t, np.inf) + np.testing.assert_almost_equal(alpha.thresh_b, theory["mb"]) + np.testing.assert_almost_equal(alpha.thresh_c, theory["mc"]) + + +def test_set_thresholds_alpha_em(): + """test value of alpha_em at threshold values""" + test_theory = API.theoryid(theoryid=THEORY_QED) + theory = test_theory.get_description() + + for modev in ['EXA', 'TRN']: + theory['ModEv'] = modev + + alpha = Alpha(theory, 1e8) + alpha_ref = theory["alphaqed"] + + # test all the thresholds + np.testing.assert_almost_equal(alpha.alpha_em_ref, theory["alphaqed"]) + np.testing.assert_almost_equal(alpha.thresh[(5, 3)], theory["Qedref"]) + np.testing.assert_almost_equal(alpha.thresh[(4, 3)], theory["mb"]) + np.testing.assert_almost_equal(alpha.thresh[(4, 2)], constants.MTAU) + np.testing.assert_almost_equal(alpha.thresh[(3, 2)], theory["mc"]) + np.testing.assert_almost_equal(alpha.thresh[(0, 2)], constants.MQL) + np.testing.assert_almost_equal(alpha.thresh[(0, 1)], constants.MMU) + np.testing.assert_almost_equal(alpha.thresh[(0, 0)], constants.ME) + + # test some values of alpha at the threshold points + np.testing.assert_almost_equal(alpha.alphaem_thresh[(5, 3)], theory["alphaqed"]) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(4, 3)], + alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5, 3), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(4, 2)], + alpha.alphaem_fixed_flavor( + constants.MTAU, alpha.alphaem_thresh[(4, 3)], theory["mb"], 4, 3 + ), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(3, 2)], + alpha.alphaem_fixed_flavor( + theory["mc"], alpha.alphaem_thresh[(4, 2)], constants.MTAU, 4, 2 + ), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(0, 2)], + alpha.alphaem_fixed_flavor( + constants.MQL, alpha.alphaem_thresh[(3, 2)], theory["mc"], 3, 2 + ), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(0, 1)], + alpha.alphaem_fixed_flavor( + constants.MMU, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2 + ), + ) + np.testing.assert_almost_equal( + alpha.alphaem_thresh[(0, 0)], + alpha.alphaem_fixed_flavor( + constants.ME, alpha.alphaem_thresh[(0, 1)], constants.MMU, 0, 1 + ), + ) + + np.testing.assert_equal(len(alpha.alphaem_thresh), 7) + np.testing.assert_equal(len(alpha.thresh), 7) + + +def test_couplings_exa(): + """ + Benchmark the exact running of alpha: + Testing both the FFS and the VFS with the results from EKO. + Test also the values of the couplings at the threshold values. + """ + test_theory = API.theoryid(theoryid=THEORY_QED) + theory = test_theory.get_description() + for k in [1]: + theory["mb"] *= k + mass_list = [theory["mc"], theory["mb"], theory["mt"]] + + alpha = Alpha(theory, 1e8) + couplings = CouplingsInfo.from_dict( + dict( + alphas=theory["alphas"], + alphaem=theory["alphaqed"], + scale=theory["Qref"], + num_flavs_ref=None, + max_num_flavs=theory["MaxNfAs"], + em_running=True, + ) + ) + eko_alpha = Couplings( + couplings, + (theory["PTO"] + 1, theory["QED"]), + method=CouplingEvolutionMethod.EXACT, + masses=[m**2 for m in mass_list], + hqm_scheme=QuarkMassScheme.POLE, + thresholds_ratios=[1.0, 1.0, 1.0], + ) + eko_alpha.decoupled_running = True + alpha_ref = theory["alphaqed"] + for q in [5, 10, 20, 50, 80, 100, 200]: + np.testing.assert_allclose( + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), + alpha.alpha_em(q), + rtol=5e-6, + ) + np.testing.assert_allclose( + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), + eko_alpha.compute_exact_alphaem_running( + np.array([0.118, alpha_ref]) / (4 * np.pi), 5, 3, theory["Qref"] ** 2, q**2 + )[1] + * 4 + * np.pi, + rtol=1e-7, + ) + + # TODO: these tests have to be switched on once the varying nl is implemented in eko + + for q in [1, 2, 3, 4]: + np.testing.assert_allclose( + alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6 + ) + for nf in range(3, theory["MaxNfAs"]): + np.testing.assert_allclose( + alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)], + eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi, + rtol=3e-7, + ) + + +def test_exa_interpolation(): + """test the accuracy of the alphaem interpolation""" + test_theory = API.theoryid(theoryid=THEORY_QED) + theory = test_theory.get_description() + + alpha = Alpha(theory, 1e8) + for q in np.geomspace(5.0, 1e4, 1000, endpoint=True): + np.testing.assert_allclose( + alpha.alpha_em(q), + alpha.alphaem_fixed_flavor(q, theory["alphaqed"], theory["Qref"], 5, 3), + rtol=1e-5, + ) + for q in np.geomspace(1.778, 4.9, 20): + np.testing.assert_allclose( + alpha.alpha_em(q), + alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(4, 3)], alpha.thresh_b, 4, 3), + rtol=1e-5, + ) + for q in np.geomspace(0.2, 0.8, 20): + np.testing.assert_allclose( + alpha.alpha_em(q), + alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2), + rtol=1e-5, + ) + + +def test_couplings_trn(): + """benchmark the truncated running of alpha with EKO""" + test_theory = API.theoryid(theoryid=THEORY_QED) + theory = test_theory.get_description() + theory["ModEv"] = 'TRN' + alpha = Alpha(theory, 1e8) + alpha_ref = theory["alphaqed"] + + for q in [80, 10, 5]: + np.testing.assert_allclose( + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), + alpha.alpha_em(q), + rtol=1e-10, + ) + np.testing.assert_allclose( + alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), + expanded_qed( + alpha_ref / (4 * np.pi), + theory["QED"], + beta.beta_qed((0, 2), 5, 3), + [beta.b_qed((0, i + 2), 5, 3) for i in range(theory["QED"])], + 2 * np.log(q / theory["Qref"]), + ) + * 4 + * np.pi, + rtol=1e-7, + ) + + +def test_trn_vs_exa(): + """benchmark the exact running of alpha vs truncated""" + # does this test make sense? + test_theory = API.theoryid(theoryid=THEORY_QED) + theory_exa = test_theory.get_description() + theory_trn = theory_exa.copy() + theory_trn["ModEv"] = 'TRN' + + alpha_exa = Alpha(theory_exa, 1e8) + alpha_trn = Alpha(theory_trn, 1e8) + + for q in [1, 3, 10, 50, 80]: + np.testing.assert_allclose(alpha_exa.alpha_em(q), alpha_trn.alpha_em(q), rtol=2e-3) + + +def test_betas(): + """test betas for different nf""" + test_theory = API.theoryid(theoryid=THEORY_QED) + alpha = Alpha(test_theory.get_description(), 1e8) + vec_beta0 = [ + -0.0, + -0.10610329539459688, + -0.21220659078919377, + -0.42441318157838753, + -0.5658842421045168, + -0.6719875374991137, + -0.7073553026306458, + ] + vec_beta1 = [ + -0.0, + -0.025330295910584444, + -0.05066059182116889, + -0.06754745576155852, + -0.0825580014863493, + -0.10788829739693376, + -0.10882645650473316, + ] + for i, (nf, nl) in enumerate(alpha.regions): + np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][0], vec_beta0[i], rtol=1e-7) + np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][1], vec_beta1[i], rtol=1e-7) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index 15cdc00d3f..6c9bb11a95 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -13,9 +13,10 @@ from validphys.api import API from validphys.core import PDF as PDFset from validphys.loader import FallbackLoader -from validphys.photon import structure_functions as sf from validphys.photon import constants -from validphys.photon.compute import FIATLUX_DEFAULT, Alpha, Photon +from validphys.photon import structure_functions as sf +from validphys.photon.alpha import Alpha +from validphys.photon.compute import FIATLUX_DEFAULT, Photon from ..conftest import PDF, THEORY_QED @@ -48,201 +49,6 @@ def test_parameters_init(): np.testing.assert_equal(photon.q_in, 100.0) -def test_masses_init(): - """test thresholds values in Alpha class""" - test_theory = API.theoryid(theoryid=THEORY_QED) - theory = test_theory.get_description() - alpha = Alpha(theory, 1e8) - np.testing.assert_equal(alpha.thresh_t, np.inf) - np.testing.assert_almost_equal(alpha.thresh_b, theory["mb"]) - np.testing.assert_almost_equal(alpha.thresh_c, theory["mc"]) - - -def test_set_thresholds_alpha_em(): - """test value of alpha_em at threshold values""" - test_theory = API.theoryid(theoryid=THEORY_QED) - theory = test_theory.get_description() - - for modev in ['EXA', 'TRN']: - theory['ModEv'] = modev - - alpha = Alpha(theory, 1e8) - alpha_ref = theory["alphaqed"] - - # test all the thresholds - np.testing.assert_almost_equal(alpha.alpha_em_ref, theory["alphaqed"]) - np.testing.assert_almost_equal(alpha.thresh[(5, 3)], theory["Qedref"]) - np.testing.assert_almost_equal(alpha.thresh[(4, 3)], theory["mb"]) - np.testing.assert_almost_equal(alpha.thresh[(4, 2)], constants.MTAU) - np.testing.assert_almost_equal(alpha.thresh[(3, 2)], theory["mc"]) - np.testing.assert_almost_equal(alpha.thresh[(0, 2)], constants.MQL) - np.testing.assert_almost_equal(alpha.thresh[(0, 1)], constants.MMU) - np.testing.assert_almost_equal(alpha.thresh[(0, 0)], constants.ME) - - # test some values of alpha at the threshold points - np.testing.assert_almost_equal(alpha.alphaem_thresh[(5, 3)], theory["alphaqed"]) - np.testing.assert_almost_equal( - alpha.alphaem_thresh[(4, 3)], - alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5, 3), - ) - np.testing.assert_almost_equal( - alpha.alphaem_thresh[(4, 2)], - alpha.alphaem_fixed_flavor(constants.MTAU, alpha.alphaem_thresh[(4, 3)], theory["mb"], 4, 3), - ) - np.testing.assert_almost_equal( - alpha.alphaem_thresh[(3, 2)], - alpha.alphaem_fixed_flavor(theory["mc"], alpha.alphaem_thresh[(4, 2)], constants.MTAU, 4, 2), - ) - np.testing.assert_almost_equal( - alpha.alphaem_thresh[(0, 2)], - alpha.alphaem_fixed_flavor(constants.MQL, alpha.alphaem_thresh[(3, 2)], theory["mc"], 3, 2), - ) - np.testing.assert_almost_equal( - alpha.alphaem_thresh[(0, 1)], - alpha.alphaem_fixed_flavor(constants.MMU, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2), - ) - np.testing.assert_almost_equal( - alpha.alphaem_thresh[(0, 0)], - alpha.alphaem_fixed_flavor(constants.ME, alpha.alphaem_thresh[(0, 1)], constants.MMU, 0, 1), - ) - - np.testing.assert_equal(len(alpha.alphaem_thresh), 7) - np.testing.assert_equal(len(alpha.thresh), 7) - - -def test_couplings_exa(): - """ - Benchmark the exact running of alpha: - Testing both the FFS and the VFS with the results from EKO. - Test also the values of the couplings at the threshold values. - """ - test_theory = API.theoryid(theoryid=THEORY_QED) - theory = test_theory.get_description() - for k in [1]: - theory["mb"] *= k - mass_list = [theory["mc"], theory["mb"], theory["mt"]] - - alpha = Alpha(theory, 1e8) - couplings = CouplingsInfo.from_dict( - dict( - alphas=theory["alphas"], - alphaem=theory["alphaqed"], - scale=theory["Qref"], - num_flavs_ref=None, - max_num_flavs=theory["MaxNfAs"], - em_running=True, - ) - ) - eko_alpha = Couplings( - couplings, - (theory["PTO"] + 1, theory["QED"]), - method=CouplingEvolutionMethod.EXACT, - masses=[m**2 for m in mass_list], - hqm_scheme=QuarkMassScheme.POLE, - thresholds_ratios=[1.0, 1.0, 1.0], - ) - eko_alpha.decoupled_running = True - alpha_ref = theory["alphaqed"] - for q in [5, 10, 20, 50, 80, 100, 200]: - np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), - alpha.alpha_em(q), - rtol=5e-6, - ) - np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), - eko_alpha.compute_exact_alphaem_running( - np.array([0.118, alpha_ref]) / (4 * np.pi), 5, 3, theory["Qref"] ** 2, q**2 - )[1] - * 4 - * np.pi, - rtol=1e-7, - ) - - # TODO: these tests have to be switched on once the varying nl is implemented in eko - - for q in [1, 2, 3, 4]: - np.testing.assert_allclose( - alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6 - ) - for nf in range(3, theory["MaxNfAs"]): - np.testing.assert_allclose( - alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)], - eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi, - rtol=3e-7, - ) - - -def test_exa_interpolation(): - """test the accuracy of the alphaem interpolation""" - test_theory = API.theoryid(theoryid=THEORY_QED) - theory = test_theory.get_description() - - alpha = Alpha(theory, 1e8) - for q in np.geomspace(5., 1e4, 1000, endpoint=True): - np.testing.assert_allclose(alpha.alpha_em(q), alpha.alphaem_fixed_flavor(q, theory["alphaqed"], theory["Qref"], 5, 3), rtol=1e-5) - for q in np.geomspace(1.778, 4.9, 20): - np.testing.assert_allclose(alpha.alpha_em(q), alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(4, 3)], alpha.thresh_b, 4, 3), rtol=1e-5) - for q in np.geomspace(0.2, 0.8, 20): - np.testing.assert_allclose(alpha.alpha_em(q), alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2), rtol=1e-5) - - - -def test_couplings_trn(): - """benchmark the truncated running of alpha with EKO""" - test_theory = API.theoryid(theoryid=THEORY_QED) - theory = test_theory.get_description() - theory["ModEv"] = 'TRN' - alpha = Alpha(theory, 1e8) - alpha_ref = theory["alphaqed"] - - for q in [80, 10, 5]: - np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), - alpha.alpha_em(q), - rtol=1e-10, - ) - np.testing.assert_allclose( - alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3), - expanded_qed( - alpha_ref / (4 * np.pi), - theory["QED"], - beta.beta_qed((0, 2), 5), - [beta.b_qed((0, i + 2), 5) for i in range(theory["QED"])], - 2 * np.log(q / theory["Qref"]), - ) - * 4 - * np.pi, - rtol=1e-7, - ) - - -def test_trn_vs_exa(): - """benchmark the exact running of alpha vs truncated""" - # does this test make sense? - test_theory = API.theoryid(theoryid=THEORY_QED) - theory_exa = test_theory.get_description() - theory_trn = theory_exa.copy() - theory_trn["ModEv"] = 'TRN' - - alpha_exa = Alpha(theory_exa, 1e8) - alpha_trn = Alpha(theory_trn, 1e8) - - for q in [1, 3, 10, 50, 80]: - np.testing.assert_allclose(alpha_exa.alpha_em(q), alpha_trn.alpha_em(q), rtol=2e-3) - - -def test_betas(): - """test betas for different nf""" - test_theory = API.theoryid(theoryid=THEORY_QED) - alpha = Alpha(test_theory.get_description(), 1e8) - vec_beta0 = [-0.0, -0.10610329539459688, -0.21220659078919377, -0.42441318157838753, -0.5658842421045168, -0.6719875374991137, -0.7073553026306458] - vec_beta1 = [-0.0, -0.025330295910584444, -0.05066059182116889, -0.06754745576155852, -0.0825580014863493, -0.10788829739693376, -0.10882645650473316] - for i, (nf, nl) in enumerate(alpha.regions): - np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][0], vec_beta0[i], rtol=1e-7) - np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][1], vec_beta1[i], rtol=1e-7) - - def test_photon(): """ Test that photon coming out of Photon interpolator matches the photon array From 9a8b3647fa9c154c43f86caf63a2b6791a21dc77 Mon Sep 17 00:00:00 2001 From: niccolo Date: Thu, 11 Jan 2024 16:48:41 +0100 Subject: [PATCH 08/15] Call isort on test_compute --- validphys2/src/validphys/tests/photon/test_compute.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py index 6c9bb11a95..0cd51c561c 100644 --- a/validphys2/src/validphys/tests/photon/test_compute.py +++ b/validphys2/src/validphys/tests/photon/test_compute.py @@ -4,19 +4,14 @@ import numpy as np import yaml -from eko import beta -from eko.couplings import Couplings, expanded_qed from eko.io import EKO -from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo -from eko.quantities.heavy_quarks import QuarkMassScheme from n3fit.io.writer import XGRID from validphys.api import API from validphys.core import PDF as PDFset from validphys.loader import FallbackLoader -from validphys.photon import constants from validphys.photon import structure_functions as sf -from validphys.photon.alpha import Alpha from validphys.photon.compute import FIATLUX_DEFAULT, Photon +from validphys.photon.alpha import Alpha from ..conftest import PDF, THEORY_QED From b001ba369b0fee8716f869029e496bd473e5f305 Mon Sep 17 00:00:00 2001 From: niccolo Date: Fri, 12 Jan 2024 12:15:00 +0100 Subject: [PATCH 09/15] Update eko version --- conda-recipe/meta.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 23b2098e1e..60cd770671 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -52,7 +52,7 @@ requirements: - sphinxcontrib-bibtex - curio >=1.0 - pineappl >=0.6.2 - - eko >=0.14.0 + - eko >=0.14.1 - fiatlux test: From b0426a89356f5a0311ef84749b00958afa63adb8 Mon Sep 17 00:00:00 2001 From: niccolo Date: Fri, 12 Jan 2024 16:59:27 +0100 Subject: [PATCH 10/15] Small fix --- validphys2/src/validphys/photon/alpha.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py index 594a047855..0fda9c87b0 100644 --- a/validphys2/src/validphys/photon/alpha.py +++ b/validphys2/src/validphys/photon/alpha.py @@ -170,9 +170,8 @@ def compute_alphaem_at_thresholds(self): thresh[(nf_, nl_)] = mq regions = list(thresh.keys()) - self.regions = regions - self.betas_qed = self.compute_betas() + self.betas_qed = self.compute_betas(regions) start = regions.index((nfref, nlref)) @@ -196,10 +195,10 @@ def compute_alphaem_at_thresholds(self): return thresh, alphaem_thresh - def compute_betas(self): + def compute_betas(self, regions): """Set values of betaQCD and betaQED.""" betas_qed = {} - for nf, nl in self.regions: + for nf, nl in regions: betas_qed[(nf, nl)] = [ beta.beta_qed_aem2(nf, nl) / (4 * np.pi), beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2, From 82a4b70e04a0f18b11f68ad615188a5849692c77 Mon Sep 17 00:00:00 2001 From: niccolo Date: Fri, 12 Jan 2024 17:53:12 +0100 Subject: [PATCH 11/15] Fix tests broken with last commit --- validphys2/src/validphys/tests/photon/test_alpha.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validphys2/src/validphys/tests/photon/test_alpha.py b/validphys2/src/validphys/tests/photon/test_alpha.py index 53361cbdc6..7db54055b6 100644 --- a/validphys2/src/validphys/tests/photon/test_alpha.py +++ b/validphys2/src/validphys/tests/photon/test_alpha.py @@ -238,6 +238,6 @@ def test_betas(): -0.10788829739693376, -0.10882645650473316, ] - for i, (nf, nl) in enumerate(alpha.regions): + for i, (nf, nl) in enumerate([(0, 0), (0, 1), (0, 2), (3, 2), (4, 2), (4, 3), (5, 3)]): np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][0], vec_beta0[i], rtol=1e-7) np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][1], vec_beta1[i], rtol=1e-7) From 97490a8022cade6cfc575aefb4c9c97107f862be Mon Sep 17 00:00:00 2001 From: niccolo Date: Mon, 15 Jan 2024 16:17:30 +0100 Subject: [PATCH 12/15] Fix some stuff --- validphys2/src/validphys/photon/alpha.py | 23 +++++++++++++++++----- validphys2/src/validphys/photon/compute.py | 2 +- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py index 0fda9c87b0..94bb1cea7a 100644 --- a/validphys2/src/validphys/photon/alpha.py +++ b/validphys2/src/validphys/photon/alpha.py @@ -1,3 +1,4 @@ +"""Module that implements the alpha running""" import numpy as np from scipy.integrate import solve_ivp @@ -8,10 +9,12 @@ class Alpha: + """Class that solves the RGE for alpha_qed""" def __init__(self, theory, q2max): self.theory = theory self.alpha_em_ref = theory["alphaqed"] self.qref = self.theory["Qref"] + self.qed_order = theory['QED'] # compute and store thresholds self.thresh_c = self.theory["kcThr"] * self.theory["mc"] @@ -99,6 +102,8 @@ def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf, nl): reference scale nf: int number of flavors + nl: int + number of leptons Returns ------- @@ -129,18 +134,26 @@ def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl): reference scale nf: int number of flavors + nl: int + number of leptons Returns ------- alpha_em: float target value of a """ + # at LO in QED the TRN solution is exact + if self.qed_order == 1: + return self.alphaem_fixed_flavor_trn(q, alphaem_ref, qref, nf, nl) + u = 2 * np.log(q / qref) # solve RGE res = solve_ivp( - rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6 + _rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6 ) + # taking fist (and only) element of y since it is a 1-D differential equation. + # Then taking the last element since it is the requested value return res.y[0][-1] def compute_alphaem_at_thresholds(self): @@ -201,7 +214,7 @@ def compute_betas(self, regions): for nf, nl in regions: betas_qed[(nf, nl)] = [ beta.beta_qed_aem2(nf, nl) / (4 * np.pi), - beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2, + beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2 if self.theory['QED'] == 2 else 0., ] return betas_qed @@ -233,9 +246,9 @@ def find_region(self, q): return nf, nl -def rge(_t, alpha, beta_qed_vec): +def _rge(_t, alpha_t, beta_qed_vec): """RGEs for the running of alphaem""" - rge_qed = -(alpha**2) * ( - beta_qed_vec[0] + np.sum([alpha ** (k + 1) * b for k, b in enumerate(beta_qed_vec[1:])]) + rge_qed = -(alpha_t**2) * ( + beta_qed_vec[0] + np.sum([alpha_t ** (k + 1) * beta for k, beta in enumerate(beta_qed_vec[1:])]) ) return rge_qed diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py index 19e9ff629c..453f45c8ba 100644 --- a/validphys2/src/validphys/photon/compute.py +++ b/validphys2/src/validphys/photon/compute.py @@ -1,4 +1,4 @@ -"""Script that calls fiatlux to add the photon PDF.""" +"""Module that calls fiatlux to add the photon PDF.""" import logging import tempfile From ef042daf3ce75a6e083bf351db3e075d90aa4a7a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <57321974+niclaurenti@users.noreply.github.com> Date: Tue, 16 Jan 2024 14:16:12 +0100 Subject: [PATCH 13/15] Update validphys2/src/validphys/photon/alpha.py Co-authored-by: Juan M. Cruz-Martinez --- validphys2/src/validphys/photon/alpha.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py index 94bb1cea7a..964b38741b 100644 --- a/validphys2/src/validphys/photon/alpha.py +++ b/validphys2/src/validphys/photon/alpha.py @@ -153,7 +153,7 @@ def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl): _rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6 ) # taking fist (and only) element of y since it is a 1-D differential equation. - # Then taking the last element since it is the requested value + # then the last element of the array which corresponds to alpha(q) return res.y[0][-1] def compute_alphaem_at_thresholds(self): From 22715d80616953ee965fa988a855848cdbdb632b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <57321974+niclaurenti@users.noreply.github.com> Date: Tue, 16 Jan 2024 14:16:22 +0100 Subject: [PATCH 14/15] Update validphys2/src/validphys/photon/alpha.py Co-authored-by: Juan M. Cruz-Martinez --- validphys2/src/validphys/photon/alpha.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py index 964b38741b..dbede0b0b1 100644 --- a/validphys2/src/validphys/photon/alpha.py +++ b/validphys2/src/validphys/photon/alpha.py @@ -88,7 +88,7 @@ def alpha_em(self, q): def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf, nl): """ - Compute the running alphaem for nf fixed at NLO, using truncated method. + Compute the running alphaem for nf fixed at, at least, NLO, using truncated method. In this case the RGE for alpha_em is solved decoupling it from the RGE for alpha_s (so the mixed terms are removed). alpha_s will just be unused. From fe782343c6f3fd52d63d6fdf5baa52c262d770bb Mon Sep 17 00:00:00 2001 From: "Juan M. Cruz-Martinez" Date: Tue, 16 Jan 2024 15:46:18 +0100 Subject: [PATCH 15/15] Update validphys2/src/validphys/photon/alpha.py --- validphys2/src/validphys/photon/alpha.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py index dbede0b0b1..574f12d67c 100644 --- a/validphys2/src/validphys/photon/alpha.py +++ b/validphys2/src/validphys/photon/alpha.py @@ -13,7 +13,7 @@ class Alpha: def __init__(self, theory, q2max): self.theory = theory self.alpha_em_ref = theory["alphaqed"] - self.qref = self.theory["Qref"] + self.qref = theory["Qref"] self.qed_order = theory['QED'] # compute and store thresholds