diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index cdd87ed11..242fee66e 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -312,6 +312,7 @@ def quad_ker( sv_mode, is_threshold, n3lo_ad_variation, + use_fhmruvv, ) # recombine everything @@ -461,6 +462,7 @@ def quad_ker_qed( sv_mode, is_threshold, n3lo_ad_variation, + use_fhmruvv, ): """Raw evolution kernel inside quad. @@ -501,7 +503,9 @@ def quad_ker_qed( is_threshold : boolean is this an itermediate threshold operator? n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv : bool + if True use the |FHMRUVV| |N3LO| anomalous dimensions Returns ------- @@ -510,7 +514,9 @@ def quad_ker_qed( """ # compute the actual evolution kernel for QEDxQCD if ker_base.is_QEDsinglet: - gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation) + gamma_s = ad_us.gamma_singlet_qed( + order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_s = sv.exponentiated.gamma_variation_qed( @@ -538,7 +544,9 @@ def quad_ker_qed( ) @ np.ascontiguousarray(ker) ker = select_QEDsinglet_element(ker, mode0, mode1) elif ker_base.is_QEDvalence: - gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf) + gamma_v = ad_us.gamma_valence_qed( + order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_v = sv.exponentiated.gamma_variation_qed( @@ -563,7 +571,9 @@ def quad_ker_qed( ) @ np.ascontiguousarray(ker) ker = select_QEDvalence_element(ker, mode0, mode1) else: - gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf) + gamma_ns = ad_us.gamma_ns_qed( + order, mode0, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv + ) # scale var exponentiated is directly applied on gamma if sv_mode == sv.Modes.exponentiated: gamma_ns = sv.exponentiated.gamma_variation_qed( diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index 09f78e861..32687e489 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -137,20 +137,24 @@ def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=False): @nb.njit(cache=True) -def gamma_ns_qed(order, mode, n, nf): +def gamma_ns_qed(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False): r""" Compute the grid of the QED non-singlet anomalous dimensions. Parameters ---------- - order : tuple(int,int) - perturbative orders - mode : 10102 | 10103 | 10202 | 10203 - sector identifier - n : complex - Mellin variable - nf : int - Number of active flavors + order : tuple(int,int) + perturbative orders + mode : 10102 | 10103 | 10202 | 10203 + sector identifier + n : complex + Mellin variable + nf : int + Number of active flavors + n3lo_ad_variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions Returns ------- @@ -179,10 +183,20 @@ def gamma_ns_qed(order, mode, n, nf): elif mode in [10202, 10203]: gamma_ns[3, 0] = as3.gamma_nsm(n, nf, cache) if order[0] >= 4: - if mode in [10102, 10103]: - gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache) - elif mode in [10202, 10203]: - gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache) + if use_fhmruvv: + if mode in [10102, 10103]: + gamma_ns[4, 0] = as4.fhmruvv.gamma_nsp( + n, nf, cache, n3lo_ad_variation[4] + ) + elif mode in [10202, 10203]: + gamma_ns[4, 0] = as4.fhmruvv.gamma_nsm( + n, nf, cache, n3lo_ad_variation[5] + ) + else: + if mode in [10102, 10103]: + gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache) + elif mode in [10202, 10203]: + gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache) return gamma_ns @@ -278,7 +292,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache): @nb.njit(cache=True) -def gamma_singlet_qed(order, n, nf, n3lo_ad_variation): +def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False): r""" Compute the grid of the QED singlet anomalous dimensions matrices. @@ -291,7 +305,9 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation): nf : int Number of active flavors n3lo_ad_variation : tuple - |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)`` + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions Returns ------- @@ -311,12 +327,17 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation): if order[0] >= 3: gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, cache) if order[0] >= 4: - gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation) + if use_fhmruvv: + gamma_s[4, 0] = as4.fhmruvv.gamma_singlet_qed( + n, nf, cache, n3lo_ad_variation + ) + else: + gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation) return gamma_s @nb.njit(cache=True) -def gamma_valence_qed(order, n, nf): +def gamma_valence_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False): r""" Compute the grid of the QED valence anomalous dimensions matrices. @@ -328,6 +349,11 @@ def gamma_valence_qed(order, n, nf): Mellin variable nf : int Number of active flavors + n3lo_ad_variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)`` + use_fhmruvv: bool + if True use the |FHMRUVV| N3LO anomalous dimensions + Returns ------- @@ -347,5 +373,10 @@ def gamma_valence_qed(order, n, nf): if order[0] >= 3: gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, cache) if order[0] >= 4: - gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache) + if use_fhmruvv: + gamma_v[4, 0] = as4.fhmruvv.gamma_valence_qed( + n, nf, cache, n3lo_ad_variation + ) + else: + gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache) return gamma_v diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py index e36de43e8..e92a55ffa 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/__init__.py @@ -56,3 +56,88 @@ def gamma_singlet(N, nf, cache, variation): np.complex_, ) return gamma_S_0 + + +@nb.njit(cache=True) +def gamma_singlet_qed(N, nf, cache, variation): + r"""Compute the leading-order singlet anomalous dimension matrix for the unified evolution basis. + + .. math:: + \\gamma_S^{(3,0)} = \\left(\begin{array}{cccc} + \\gamma_{gg}^{(3,0)} & 0 & \\gamma_{gq}^{(3,0)} & 0\\ + 0 & 0 & 0 & 0 \\ + \\gamma_{qg}^{(3,0)} & 0 & \\gamma_{qq}^{(3,0)} & 0 \\ + 0 & 0 & 0 & \\gamma_{qq}^{(3,0)} \\ + \\end{array}\right) + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : tuple + |N3LO| anomalous dimension variation ``(gg, gq, qg, qq)`` + + Returns + ------- + numpy.ndarray + Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(3,0)}(N)` + + """ + gamma_np_p = gamma_nsp(N, nf, cache, variation[3]) + gamma_qq = gamma_np_p + gamma_ps(N, nf, cache, variation[3]) + gamma_S = np.array( + [ + [ + gamma_gg(N, nf, cache, variation[0]), + 0.0 + 0.0j, + gamma_gq(N, nf, cache, variation[1]), + 0.0 + 0.0j, + ], + [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j], + [gamma_qg(N, nf, cache, variation[2]), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j], + [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_np_p], + ], + np.complex_, + ) + return gamma_S + + +@nb.njit(cache=True) +def gamma_valence_qed(N, nf, cache, variation): + r"""Compute the leading-order valence anomalous dimension matrix for the unified evolution basis. + + .. math:: + \\gamma_V^{(3,0)} = \\left(\begin{array}{cc} + \\gamma_{nsV}^{(3,0)} & 0\\ + 0 & \\gamma_{ns-}^{(3,0)} + \\end{array}\right) + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + cache: numpy.ndarray + Harmonic sum cache + variation : tuple + |N3LO| anomalous dimension variation ``(nsm, nsv)`` + + Returns + ------- + numpy.ndarray + Leading-order singlet anomalous dimension matrix :math:`\\gamma_{V}^{(3,0)}(N)` + + """ + gamma_V = np.array( + [ + [gamma_nsv(N, nf, cache, variation[-1]), 0.0], + [0.0, gamma_nsm(N, nf, cache, variation[-2])], + ], + np.complex_, + ) + return gamma_V diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py index e6179ba99..bb18c3e17 100644 --- a/tests/eko/kernels/test_kernels_QEDns.py +++ b/tests/eko/kernels/test_kernels_QEDns.py @@ -77,7 +77,7 @@ def test_zero_true_gamma(): for qed in range(1, 2 + 1): order = (qcd, qed) n = np.random.rand() - gamma_ns = ad.gamma_ns_qed(order, mode, n, nf) + gamma_ns = ad.gamma_ns_qed(order, mode, n, nf, (0, 0, 0, 0, 0, 0, 0)) for method in methods: for running in running_alpha: np.testing.assert_allclose( @@ -248,7 +248,9 @@ def test_ode_true_gamma(): ) a1 = sc.a_s(mu2_to) n = 3 + np.random.rand() - gamma_ns = ad.gamma_ns_qed(order, mode, n, nf) + gamma_ns = ad.gamma_ns_qed( + order, mode, n, nf, (0, 0, 0, 0, 0, 0, 0) + ) gammatot = 0.0 for i in range(0, order[0] + 1): for j in range(0, order[1] + 1): diff --git a/tests/eko/kernels/test_kernels_QEDvalence.py b/tests/eko/kernels/test_kernels_QEDvalence.py index a27bc6ef6..29f0f3d7a 100644 --- a/tests/eko/kernels/test_kernels_QEDvalence.py +++ b/tests/eko/kernels/test_kernels_QEDvalence.py @@ -67,7 +67,7 @@ def test_zero_true_gamma(monkeypatch): order = (qcd, qed) n = np.random.rand() gamma_v = anomalous_dimensions.unpolarized.space_like.gamma_valence_qed( - order, n, nf + order, n, nf, (0, 0, 0, 0, 0, 0, 0) ) for method in methods: np.testing.assert_allclose( diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py index cb65f4ba6..5fdc1ee06 100644 --- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py +++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py @@ -147,66 +147,90 @@ def test_gamma_ns_qed(): nf = 3 # aem1 assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 2)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 2)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns+u"], 1, nf)[0, 1], + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns+u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + )[0, 1], 0, decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 1), br.non_singlet_pids_map["ns+d"], 1, nf)[0, 1], + ad_us.gamma_ns_qed( + (1, 1), br.non_singlet_pids_map["ns+d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + )[0, 1], 0, decimal=5, ) # as1aem1 assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) # aem2 assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((1, 2), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (1, 2), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((2, 3)), decimal=5, ) # ad_us.as2 assert_almost_equal( - ad_us.gamma_ns_qed((2, 1), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (2, 1), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((3, 2)), decimal=5, ) assert_almost_equal( - ad_us.gamma_ns_qed((2, 1), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (2, 1), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((3, 2)), decimal=5, ) # ad_us.as3 assert_almost_equal( - ad_us.gamma_ns_qed((3, 1), br.non_singlet_pids_map["ns-u"], 1, nf), + ad_us.gamma_ns_qed( + (3, 1), br.non_singlet_pids_map["ns-u"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((4, 2)), decimal=3, ) assert_almost_equal( - ad_us.gamma_ns_qed((3, 1), br.non_singlet_pids_map["ns-d"], 1, nf), + ad_us.gamma_ns_qed( + (3, 1), br.non_singlet_pids_map["ns-d"], 1, nf, (0, 0, 0, 0, 0, 0, 0) + ), np.zeros((4, 2)), decimal=3, ) @@ -238,7 +262,7 @@ def test_dim_valence(): nf = 3 N = 2 cache = h.cache.reset() - gamma_valence = ad_us.gamma_valence_qed((3, 2), N, nf) + gamma_valence = ad_us.gamma_valence_qed((3, 2), N, nf, (0, 0, 0, 0, 0, 0, 0)) assert gamma_valence.shape == (4, 3, 2, 2) gamma_valence_as1 = ad_us.as1.gamma_valence_qed(N, cache) assert gamma_valence_as1.shape == (2, 2) @@ -251,9 +275,9 @@ def test_dim_valence(): def test_dim_nsp(): nf = 3 N = 2 - gamma_nsup = ad_us.gamma_ns_qed((3, 2), 10102, N, nf) + gamma_nsup = ad_us.gamma_ns_qed((3, 2), 10102, N, nf, (0, 0, 0, 0, 0, 0, 0)) assert gamma_nsup.shape == (4, 3) - gamma_nsdp = ad_us.gamma_ns_qed((3, 2), 10103, N, nf) + gamma_nsdp = ad_us.gamma_ns_qed((3, 2), 10103, N, nf, (0, 0, 0, 0, 0, 0, 0)) assert gamma_nsdp.shape == (4, 3) with pytest.raises(NotImplementedError): - ad_us.gamma_ns_qed((2, 0), 10106, N, nf) + ad_us.gamma_ns_qed((2, 0), 10106, N, nf, (0, 0, 0, 0, 0, 0, 0))