diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/eko/benchmark_ad.py index b220cae5a..1087db0e0 100644 --- a/benchmarks/eko/benchmark_ad.py +++ b/benchmarks/eko/benchmark_ad.py @@ -3,8 +3,8 @@ import numpy as np import pytest +import eko.anomalous_dimensions.as2 as ad_as2 import eko.anomalous_dimensions.harmonics as h -import eko.anomalous_dimensions.nlo as ad_nlo from eko.constants import CA, CF, TR @@ -136,8 +136,8 @@ def check_gamma_1_pegasus(N, NF): P1NSP = CF * ((CF - CA / 2.0) * PNPA + CA * PNSB + TR * NF * PNSC) P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC) - np.testing.assert_allclose(ad_nlo.gamma_nsp_1(N, NF), -P1NSP) - np.testing.assert_allclose(ad_nlo.gamma_nsm_1(N, NF), -P1NSM) + np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF), -P1NSP) + np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF), -P1NSM) NS = N * N NT = NS * N @@ -254,7 +254,7 @@ def check_gamma_1_pegasus(N, NF): P1Sgq = (CF * CF * PGQA + CF * CA * PGQB + TR * NF * CF * PGQC) * 4.0 P1Sgg = (CA * CA * PGGA + TR * NF * (CA * PGGB + CF * PGGC)) * 4.0 - gS1 = ad_nlo.gamma_singlet_1(N, NF) + gS1 = ad_as2.gamma_singlet(N, NF) np.testing.assert_allclose(gS1[0, 0], -P1Sqq) np.testing.assert_allclose(gS1[0, 1], -P1Sqg) np.testing.assert_allclose(gS1[1, 0], -P1Sgq) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index df5540b34..79065061b 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/benchmarks/eko/benchmark_evol_to_unity.py @@ -2,6 +2,7 @@ import numpy as np import pytest +from eko import basis_rotation as br from eko.evolution_operator import Operator from eko.evolution_operator.grid import OperatorGrid from eko.interpolation import InterpolatorDispatcher @@ -76,38 +77,39 @@ def test_operator_grid( o.compute() o_back.compute() - dim = o_back.op_members["NS_v"].value.shape - for k in ["NS_v", "NS_m", "NS_p"]: + dim = o_back.op_members[(br.non_singlet_pids_map["nsV"], 0)].value.shape + for k in ["nsV", "ns-", "ns+"]: np.testing.assert_allclose( - o.op_members[k].value @ o_back.op_members[k].value, + o.op_members[(br.non_singlet_pids_map[k], 0)].value + @ o_back.op_members[(br.non_singlet_pids_map[k], 0)].value, np.eye(dim[0]), atol=7e-2, ) # qq np.testing.assert_allclose( - o_back.op_members["S_qq"].value @ o.op_members["S_qq"].value - + o_back.op_members["S_qg"].value @ o.op_members["S_gq"].value, + o_back.op_members[(100, 100)].value @ o.op_members[(100, 100)].value + + o_back.op_members[(100, 21)].value @ o.op_members[(21, 100)].value, np.eye(dim[0]), atol=7e-2, ) # qg np.testing.assert_allclose( - o_back.op_members["S_qq"].value @ o.op_members["S_qg"].value - + o_back.op_members["S_qg"].value @ o.op_members["S_gg"].value, + o_back.op_members[(100, 100)].value @ o.op_members[(100, 21)].value + + o_back.op_members[(100, 21)].value @ o.op_members[(21, 21)].value, np.zeros(dim), atol=7e-4, ) # gg np.testing.assert_allclose( - o_back.op_members["S_gg"].value @ o.op_members["S_gg"].value - + o_back.op_members["S_gq"].value @ o.op_members["S_qg"].value, + o_back.op_members[(21, 21)].value @ o.op_members[(21, 21)].value + + o_back.op_members[(21, 100)].value @ o.op_members[(100, 21)].value, np.eye(dim[0]), atol=9e-2, ) # gq np.testing.assert_allclose( - o_back.op_members["S_gg"].value @ o.op_members["S_gq"].value - + o_back.op_members["S_gq"].value @ o.op_members["S_qq"].value, + o_back.op_members[(21, 21)].value @ o.op_members[(21, 100)].value + + o_back.op_members[(21, 100)].value @ o.op_members[(100, 100)].value, np.zeros(dim), atol=2e-3, ) @@ -124,7 +126,7 @@ def test_operator_grid( # ome.compute( q2, L) # ome_back.compute(q2, L) - # dim = ome.ome_members["S_qq"].value.shape + # dim = ome.ome_members[(100, 100)].value.shape # ome_tensor = np.zeros((3,3,dim[0],dim[0])) # ome_tensor_back = ome_tensor # idx_dict = dict(zip(["g", "q", "H"],[0,1,2])) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index 624bcb5a6..8131027b0 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -420,6 +420,17 @@ @article{Bertone:2013vaa year = {2014} } +@phdthesis{Carrazza:2015dea, + author = "Carrazza, Stefano", + title = "{Parton distribution functions with QED corrections}", + eprint = "1509.00209", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + doi = "10.13130/carrazza-stefano_phd2015-07-06", + school = "Milan U.", + year = "2015" +} + @article{Vermaseren:1997fq, author = "Vermaseren, J. A. M. and Larin, S. A. and van Ritbergen, T.", title = "{The four loop quark mass anomalous dimension and the invariant quark mass}", diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py index 173fe5ead..bdf7d67ef 100644 --- a/src/eko/anomalous_dimensions/__init__.py +++ b/src/eko/anomalous_dimensions/__init__.py @@ -20,7 +20,8 @@ import numba as nb import numpy as np -from . import harmonics, lo, nlo, nnlo +from .. import basis_rotation as br +from . import as1, as2, as3, harmonics @nb.njit("Tuple((c16[:,:],c16,c16,c16[:,:],c16[:,:]))(c16[:,:])", cache=True) @@ -52,9 +53,9 @@ def exp_singlet(gamma_S): See Also -------- - eko.anomalous_dimensions.lo.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` - eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` - eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` + eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)` + eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)` """ # compute eigenvalues det = np.sqrt( @@ -71,7 +72,7 @@ def exp_singlet(gamma_S): return exp, lambda_p, lambda_m, e_p, e_m -@nb.njit("c16[:](u1,string,c16,u1)", cache=True) +@nb.njit("c16[:](u1,u2,c16,u1)", cache=True) def gamma_ns(order, mode, n, nf): r""" Computes the tower of the non-singlet anomalous dimensions @@ -80,7 +81,7 @@ def gamma_ns(order, mode, n, nf): ---------- order : int perturbative order - mode : "m" | "p" | "v" + mode : 10201 | 10101 | 10200 sector identifier n : complex Mellin variable @@ -94,37 +95,39 @@ def gamma_ns(order, mode, n, nf): See Also -------- - eko.anomalous_dimensions.lo.gamma_ns_0 : :math:`\gamma_{ns}^{(0)}(N)` - eko.anomalous_dimensions.nlo.gamma_nsp_1 : :math:`\gamma_{ns,+}^{(1)}(N)` - eko.anomalous_dimensions.nlo.gamma_nsm_1 : :math:`\gamma_{ns,-}^{(1)}(N)` - eko.anomalous_dimensions.nnlo.gamma_nsp_2 : :math:`\gamma_{ns,+}^{(2)}(N)` - eko.anomalous_dimensions.nnlo.gamma_nsm_2 : :math:`\gamma_{ns,-}^{(2)}(N)` - eko.anomalous_dimensions.nnlo.gamma_nsv_2 : :math:`\gamma_{ns,v}^{(2)}(N)` + eko.anomalous_dimensions.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)` + eko.anomalous_dimensions.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)` + eko.anomalous_dimensions.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)` + eko.anomalous_dimensions.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)` + eko.anomalous_dimensions.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)` """ # cache the s-es sx = np.full(1, harmonics.harmonic_S1(n)) # now combine gamma_ns = np.zeros(order + 1, np.complex_) - gamma_ns[0] = lo.gamma_ns_0(n, sx[0]) + gamma_ns[0] = as1.gamma_ns(n, sx[0]) # NLO and beyond if order >= 1: # TODO: pass the necessary harmonics to nlo gammas - if mode == "p": - gamma_ns_1 = nlo.gamma_nsp_1(n, nf) + if mode == 10101: + gamma_ns_1 = as2.gamma_nsp(n, nf) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here - elif mode in ["m", "v"]: - gamma_ns_1 = nlo.gamma_nsm_1(n, nf) + elif mode in [10201, 10200]: + gamma_ns_1 = as2.gamma_nsm(n, nf) + else: + raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 # NNLO and beyond if order >= 2: sx = np.append(sx, harmonics.harmonic_S2(n)) sx = np.append(sx, harmonics.harmonic_S3(n)) - if mode == "p": - gamma_ns_2 = -nnlo.gamma_nsp_2(n, nf, sx) - elif mode == "m": - gamma_ns_2 = -nnlo.gamma_nsm_2(n, nf, sx) - elif mode == "v": - gamma_ns_2 = -nnlo.gamma_nsv_2(n, nf, sx) + if mode == 10101: + gamma_ns_2 = -as3.gamma_nsp(n, nf, sx) + elif mode == 10201: + gamma_ns_2 = -as3.gamma_nsm(n, nf, sx) + elif mode == 10200: + gamma_ns_2 = -as3.gamma_nsv(n, nf, sx) gamma_ns[2] = gamma_ns_2 return gamma_ns @@ -150,9 +153,9 @@ def gamma_singlet(order, n, nf): See Also -------- - eko.anomalous_dimensions.lo.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)` - eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)` - eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)` + eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)` + eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)` + eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)` """ # cache the s-es sx = np.full(1, harmonics.harmonic_S1(n)) @@ -161,10 +164,10 @@ def gamma_singlet(order, n, nf): sx = np.append(sx, harmonics.harmonic_S3(n)) gamma_singlet = np.zeros((order + 1, 2, 2), np.complex_) - gamma_singlet[0] = lo.gamma_singlet_0(n, sx[0], nf) + gamma_singlet[0] = as1.gamma_singlet(n, sx[0], nf) if order >= 1: - gamma_singlet[1] = nlo.gamma_singlet_1(n, nf) + gamma_singlet[1] = as2.gamma_singlet(n, nf) if order == 2: sx = np.append(sx, harmonics.harmonic_S4(n)) - gamma_singlet[2] = -nnlo.gamma_singlet_2(n, nf, sx) + gamma_singlet[2] = -as3.gamma_singlet(n, nf, sx) return gamma_singlet diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py new file mode 100644 index 000000000..de620655a --- /dev/null +++ b/src/eko/anomalous_dimensions/aem1.py @@ -0,0 +1,93 @@ +# -*- coding: utf-8 -*- +import numba as nb + +from .. import constants +from . import as1 + + +@nb.njit("c16(c16)", cache=True) +def gamma_phq(N): + """ + Computes the leading-order photon-quark anomalous dimension + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + + Parameters + ---------- + N : complex + Mellin moment + + Returns + ------- + gamma_phq : complex + Leading-order photon-quark anomalous dimension :math:`\\gamma_{\\gamma q}^{(0)}(N)` + """ + + return as1.gamma_gq(N) / constants.CF + + +@nb.njit("c16(c16,u1)", cache=True) +def gamma_qph(N, nf): + """ + Computes the leading-order quark-photon anomalous dimension + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + But adding the :math:`N_C` and the :math:`2n_f` factors from :math:`\\theta` inside the + definition of :math:`\\gamma_{q \\gamma}^{(0)}(N)`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + Number of active flavors + + Returns + ------- + gamma_qph : complex + Leading-order quark-photon anomalous dimension :math:`\\gamma_{q \\gamma}^{(0)}(N)` + """ + return as1.gamma_qg(N, nf) / constants.TR * constants.NC + + +@nb.njit("c16(u1)", cache=True) +def gamma_phph(nf): + """ + Computes the leading-order photon-photon anomalous dimension + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + + Parameters + ---------- + nf : int + Number of active flavors + + Returns + ------- + gamma_phph : complex + Leading-order phton-photon anomalous dimension :math:`\\gamma_{\\gamma \\gamma}^{(0)}(N)` + """ + + return 2 / 3 * constants.NC * 2 * nf + + +@nb.njit("c16(c16,c16)", cache=True) +def gamma_ns(N, s1): + """ + Computes the leading-order non-singlet QED anomalous dimension. + + Implements Eq. (2.5) of :cite:`Carrazza:2015dea`. + + Parameters + ---------- + N : complex + Mellin moment + s1 : complex + S1(N) + + Returns + ------- + gamma_ns : complex + Leading-order non-singlet QED anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` + """ + return as1.gamma_ns(N, s1) / constants.CF diff --git a/src/eko/anomalous_dimensions/lo.py b/src/eko/anomalous_dimensions/as1.py similarity index 82% rename from src/eko/anomalous_dimensions/lo.py rename to src/eko/anomalous_dimensions/as1.py index 847fb4ba5..2efc0cae5 100644 --- a/src/eko/anomalous_dimensions/lo.py +++ b/src/eko/anomalous_dimensions/as1.py @@ -8,7 +8,7 @@ @nb.njit("c16(c16,c16)", cache=True) -def gamma_ns_0(N, s1): +def gamma_ns(N, s1): """ Computes the leading-order non-singlet anomalous dimension. @@ -23,7 +23,7 @@ def gamma_ns_0(N, s1): Returns ------- - gamma_ns_0 : complex + gamma_ns : complex Leading-order non-singlet anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` """ gamma = -(3.0 - 4.0 * s1 + 2.0 / N / (N + 1.0)) @@ -32,7 +32,7 @@ def gamma_ns_0(N, s1): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qg_0(N, nf): +def gamma_qg(N, nf): """ Computes the leading-order quark-gluon anomalous dimension @@ -47,7 +47,7 @@ def gamma_qg_0(N, nf): Returns ------- - gamma_qg_0 : complex + gamma_qg : complex Leading-order quark-gluon anomalous dimension :math:`\\gamma_{qg}^{(0)}(N)` """ gamma = -(N**2 + N + 2.0) / (N * (N + 1.0) * (N + 2.0)) @@ -56,7 +56,7 @@ def gamma_qg_0(N, nf): @nb.njit("c16(c16)", cache=True) -def gamma_gq_0(N): +def gamma_gq(N): """ Computes the leading-order gluon-quark anomalous dimension @@ -69,7 +69,7 @@ def gamma_gq_0(N): Returns ------- - gamma_gq_0 : complex + gamma_gq : complex Leading-order gluon-quark anomalous dimension :math:`\\gamma_{gq}^{(0)}(N)` """ gamma = -(N**2 + N + 2.0) / (N * (N + 1.0) * (N - 1.0)) @@ -78,7 +78,7 @@ def gamma_gq_0(N): @nb.njit("c16(c16,c16,u1)", cache=True) -def gamma_gg_0(N, s1, nf): +def gamma_gg(N, s1, nf): """ Computes the leading-order gluon-gluon anomalous dimension @@ -95,7 +95,7 @@ def gamma_gg_0(N, s1, nf): Returns ------- - gamma_gg_0 : complex + gamma_gg : complex Leading-order gluon-gluon anomalous dimension :math:`\\gamma_{gg}^{(0)}(N)` """ gamma = s1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0) @@ -104,7 +104,7 @@ def gamma_gg_0(N, s1, nf): @nb.njit("c16[:,:](c16,c16,u1)", cache=True) -def gamma_singlet_0(N, s1, nf): +def gamma_singlet(N, s1, nf): r""" Computes the leading-order singlet anomalous dimension matrix @@ -130,14 +130,13 @@ def gamma_singlet_0(N, s1, nf): See Also -------- - gamma_ns_0 : :math:`\gamma_{qq}^{(0)}` - gamma_qg_0 : :math:`\gamma_{qg}^{(0)}` - gamma_gq_0 : :math:`\gamma_{gq}^{(0)}` - gamma_gg_0 : :math:`\gamma_{gg}^{(0)}` + gamma_ns : :math:`\gamma_{qq}^{(0)}` + gamma_qg : :math:`\gamma_{qg}^{(0)}` + gamma_gq : :math:`\gamma_{gq}^{(0)}` + gamma_gg : :math:`\gamma_{gg}^{(0)}` """ - gamma_qq = gamma_ns_0(N, s1) - gamma_qg = gamma_qg_0(N, nf) - gamma_gq = gamma_gq_0(N) - gamma_gg = gamma_gg_0(N, s1, nf) - gamma_S_0 = np.array([[gamma_qq, gamma_qg], [gamma_gq, gamma_gg]], np.complex_) + gamma_qq = gamma_ns(N, s1) + gamma_S_0 = np.array( + [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, s1, nf)]], np.complex_ + ) return gamma_S_0 diff --git a/src/eko/anomalous_dimensions/nlo.py b/src/eko/anomalous_dimensions/as2.py similarity index 92% rename from src/eko/anomalous_dimensions/nlo.py rename to src/eko/anomalous_dimensions/as2.py index 3466be105..7c477c591 100644 --- a/src/eko/anomalous_dimensions/nlo.py +++ b/src/eko/anomalous_dimensions/as2.py @@ -15,7 +15,7 @@ @nb.njit("c16(c16,u1)", cache=True) -def gamma_nsm_1(n, nf): +def gamma_nsm(n, nf): """ Computes the |NLO| valence-like non-singlet anomalous dimension. @@ -30,7 +30,7 @@ def gamma_nsm_1(n, nf): Returns ------- - gamma_nsm_1 : complex + gamma_nsm : complex |NLO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(1)}(N)` """ @@ -58,7 +58,7 @@ def gamma_nsm_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_nsp_1(n, nf): +def gamma_nsp(n, nf): """ Computes the |NLO| singlet-like non-singlet anomalous dimension. @@ -73,7 +73,7 @@ def gamma_nsp_1(n, nf): Returns ------- - gamma_nsp_1 : complex + gamma_nsp : complex |NLO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(1)}(N)` """ @@ -99,7 +99,7 @@ def gamma_nsp_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_ps_1(n, nf): +def gamma_ps(n, nf): """ Computes the |NLO| pure-singlet quark-quark anomalous dimension. @@ -114,7 +114,7 @@ def gamma_ps_1(n, nf): Returns ------- - gamma_ps_1 : complex + gamma_ps : complex |NLO| pure-singlet quark-quark anomalous dimension :math:`\\gamma_{ps}^{(1)}(N)` """ @@ -126,7 +126,7 @@ def gamma_ps_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qg_1(n, nf): +def gamma_qg(n, nf): """ Computes the |NLO| quark-gluon singlet anomalous dimension. @@ -141,7 +141,7 @@ def gamma_qg_1(n, nf): Returns ------- - gamma_qg_1 : complex + gamma_qg : complex |NLO| quark-gluon singlet anomalous dimension :math:`\\gamma_{qg}^{(1)}(N)` """ @@ -159,7 +159,7 @@ def gamma_qg_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_gq_1(n, nf): +def gamma_gq(n, nf): """ Computes the |NLO| gluon-quark singlet anomalous dimension. @@ -174,7 +174,7 @@ def gamma_gq_1(n, nf): Returns ------- - gamma_gq_1 : complex + gamma_gq : complex |NLO| gluon-quark singlet anomalous dimension :math:`\\gamma_{gq}^{(1)}(N)` """ @@ -195,7 +195,7 @@ def gamma_gq_1(n, nf): @nb.njit("c16(c16,u1)", cache=True) -def gamma_gg_1(n, nf): +def gamma_gg(n, nf): """ Computes the |NLO| gluon-gluon singlet anomalous dimension. @@ -210,7 +210,7 @@ def gamma_gg_1(n, nf): Returns ------- - gamma_gg_1 : complex + gamma_gg : complex |NLO| gluon-gluon singlet anomalous dimension :math:`\\gamma_{gg}^{(1)}(N)` """ @@ -233,7 +233,7 @@ def gamma_gg_1(n, nf): @nb.njit("c16[:,:](c16,u1)", cache=True) -def gamma_singlet_1(N, nf): +def gamma_singlet(N, nf): r""" Computes the next-leading-order singlet anomalous dimension matrix @@ -257,15 +257,14 @@ def gamma_singlet_1(N, nf): See Also -------- - gamma_nsp_1 : :math:`\gamma_{qq}^{(1)}` - gamma_ps_1 : :math:`\gamma_{qq}^{(1)}` - gamma_qg_1 : :math:`\gamma_{qg}^{(1)}` - gamma_gq_1 : :math:`\gamma_{gq}^{(1)}` - gamma_gg_1 : :math:`\gamma_{gg}^{(1)}` + gamma_nsp : :math:`\gamma_{qq}^{(1)}` + gamma_ps : :math:`\gamma_{qq}^{(1)}` + gamma_qg : :math:`\gamma_{qg}^{(1)}` + gamma_gq : :math:`\gamma_{gq}^{(1)}` + gamma_gg : :math:`\gamma_{gg}^{(1)}` """ - gamma_qq = gamma_nsp_1(N, nf) + gamma_ps_1(N, nf) - gamma_qg = gamma_qg_1(N, nf) - gamma_gq = gamma_gq_1(N, nf) - gamma_gg = gamma_gg_1(N, nf) - gamma_S_0 = np.array([[gamma_qq, gamma_qg], [gamma_gq, gamma_gg]], np.complex_) + gamma_qq = gamma_nsp(N, nf) + gamma_ps(N, nf) + gamma_S_0 = np.array( + [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N, nf), gamma_gg(N, nf)]], np.complex_ + ) return gamma_S_0 diff --git a/src/eko/anomalous_dimensions/nnlo.py b/src/eko/anomalous_dimensions/as3.py similarity index 93% rename from src/eko/anomalous_dimensions/nnlo.py rename to src/eko/anomalous_dimensions/as3.py index 46a63eb5c..c0864d2ff 100644 --- a/src/eko/anomalous_dimensions/nnlo.py +++ b/src/eko/anomalous_dimensions/as3.py @@ -17,7 +17,7 @@ @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsm_2(n, nf, sx): +def gamma_nsm(n, nf, sx): """ Computes the |NNLO| valence-like non-singlet anomalous dimension. @@ -34,7 +34,7 @@ def gamma_nsm_2(n, nf, sx): Returns ------- - gamma_nsm_2 : complex + gamma_nsm : complex |NNLO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(2)}(N)` """ @@ -94,7 +94,7 @@ def gamma_nsm_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsp_2(n, nf, sx): +def gamma_nsp(n, nf, sx): """ Computes the |NNLO| singlet-like non-singlet anomalous dimension. @@ -111,7 +111,7 @@ def gamma_nsp_2(n, nf, sx): Returns ------- - gamma_nsp_2 : complex + gamma_nsp : complex |NNLO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(2)}(N)` """ @@ -171,7 +171,7 @@ def gamma_nsp_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsv_2(n, nf, sx): +def gamma_nsv(n, nf, sx): """ Computes the |NNLO| valence non-singlet anomalous dimension. @@ -188,7 +188,7 @@ def gamma_nsv_2(n, nf, sx): Returns ------- - gamma_nsv_2 : complex + gamma_nsv : complex |NNLO| valence non-singlet anomalous dimension :math:`\\gamma_{ns,v}^{(2)}(N)` """ @@ -221,12 +221,12 @@ def gamma_nsv_2(n, nf, sx): + 46.18 * E2 ) - result = gamma_nsm_2(n, nf, sx) + nf * ps2 + result = gamma_nsm(n, nf, sx) + nf * ps2 return result @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_ps_2(n, nf, sx): +def gamma_ps(n, nf, sx): """ Computes the |NNLO| pure-singlet quark-quark anomalous dimension. @@ -243,7 +243,7 @@ def gamma_ps_2(n, nf, sx): Returns ------- - gamma_ps_2 : complex + gamma_ps : complex |NNLO| pure-singlet quark-quark anomalous dimension :math:`\\gamma_{ps}^{(2)}(N)` """ @@ -298,7 +298,7 @@ def gamma_ps_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_qg_2(n, nf, sx): +def gamma_qg(n, nf, sx): """ Computes the |NNLO| quark-gluon singlet anomalous dimension. @@ -315,7 +315,7 @@ def gamma_qg_2(n, nf, sx): Returns ------- - gamma_qg_2 : complex + gamma_qg : complex |NNLO| quark-gluon singlet anomalous dimension :math:`\\gamma_{qg}^{(2)}(N)` """ @@ -372,7 +372,7 @@ def gamma_qg_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_gq_2(n, nf, sx): +def gamma_gq(n, nf, sx): """ Computes the |NNLO| gluon-quark singlet anomalous dimension. @@ -389,7 +389,7 @@ def gamma_gq_2(n, nf, sx): Returns ------- - gamma_gq_2 : complex + gamma_gq : complex |NNLO| gluon-quark singlet anomalous dimension :math:`\\gamma_{gq}^{(2)}(N)` """ @@ -462,7 +462,7 @@ def gamma_gq_2(n, nf, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_gg_2(n, nf, sx): +def gamma_gg(n, nf, sx): """ Computes the |NNLO| gluon-gluon singlet anomalous dimension. @@ -479,7 +479,7 @@ def gamma_gg_2(n, nf, sx): Returns ------- - gamma_gg_2 : complex + gamma_gg : complex |NNLO| gluon-gluon singlet anomalous dimension :math:`\\gamma_{gg}^{(2)}(N)` """ @@ -550,7 +550,7 @@ def gamma_gg_2(n, nf, sx): @nb.njit("c16[:,:](c16,u1,c16[:])", cache=True) -def gamma_singlet_2(N, nf, sx): +def gamma_singlet(N, nf, sx): r""" Computes the |NNLO| singlet anomalous dimension matrix @@ -578,15 +578,15 @@ def gamma_singlet_2(N, nf, sx): See Also -------- - gamma_nsp_2 : :math:`\gamma_{qq}^{(2)}` - gamma_ps_2 : :math:`\gamma_{qq}^{(2)}` - gamma_qg_2 : :math:`\gamma_{qg}^{(2)}` - gamma_gq_2 : :math:`\gamma_{gq}^{(2)}` - gamma_gg_2 : :math:`\gamma_{gg}^{(2)}` + gamma_nsp : :math:`\gamma_{qq}^{(2)}` + gamma_ps : :math:`\gamma_{qq}^{(2)}` + gamma_qg : :math:`\gamma_{qg}^{(2)}` + gamma_gq : :math:`\gamma_{gq}^{(2)}` + gamma_gg : :math:`\gamma_{gg}^{(2)}` """ - gamma_qq = gamma_nsp_2(N, nf, sx) + gamma_ps_2(N, nf, sx) - gamma_qg = gamma_qg_2(N, nf, sx) - gamma_gq = gamma_gq_2(N, nf, sx) - gamma_gg = gamma_gg_2(N, nf, sx) - gamma_S_0 = np.array([[gamma_qq, gamma_qg], [gamma_gq, gamma_gg]], np.complex_) + gamma_qq = gamma_nsp(N, nf, sx) + gamma_ps(N, nf, sx) + gamma_S_0 = np.array( + [[gamma_qq, gamma_qg(N, nf, sx)], [gamma_gq(N, nf, sx), gamma_gg(N, nf, sx)]], + np.complex_, + ) return gamma_S_0 diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py index 0aa256380..d3c89a769 100644 --- a/src/eko/basis_rotation.py +++ b/src/eko/basis_rotation.py @@ -69,13 +69,21 @@ ) """|pid| representation of :data:`evol_basis`.""" -singlet_labels = ("S_qq", "S_qg", "S_gq", "S_gg") -non_singlet_labels = ("NS_m", "NS_p", "NS_v") +non_singlet_pids_map = {"ns-": 10201, "ns+": 10101, "nsV": 10200} + +singlet_labels = ((100, 100), (100, 21), (21, 100), (21, 21)) +non_singlet_labels = ( + (non_singlet_pids_map["ns-"], 0), + (non_singlet_pids_map["ns+"], 0), + (non_singlet_pids_map["nsV"], 0), +) full_labels = (*singlet_labels, *non_singlet_labels) anomalous_dimensions_basis = full_labels r""" Sorted elements in Anomalous Dimensions Basis as :obj:`str`. """ +matching_hplus_pid = 90 +matching_hminus_pid = 91 # Tranformation from physical basis to QCD evolution basis rotate_flavor_to_evolution = np.array( @@ -101,13 +109,25 @@ """ map_ad_to_evolution = { - "S_qq": ["S.S"], - "S_qg": ["S.g"], - "S_gq": ["g.S"], - "S_gg": ["g.g"], - "NS_v": ["V.V"], - "NS_p": ["T3.T3", "T8.T8", "T15.T15", "T24.T24", "T35.T35"], - "NS_m": ["V3.V3", "V8.V8", "V15.V15", "V24.V24", "V35.V35"], + (100, 100): ["S.S"], + (100, 21): ["S.g"], + (21, 100): ["g.S"], + (21, 21): ["g.g"], + (non_singlet_pids_map["nsV"], 0): ["V.V"], + (non_singlet_pids_map["ns+"], 0): [ + "T3.T3", + "T8.T8", + "T15.T15", + "T24.T24", + "T35.T35", + ], + (non_singlet_pids_map["ns-"], 0): [ + "V3.V3", + "V8.V8", + "V15.V15", + "V24.V24", + "V35.V35", + ], } """ Map anomalous dimension sectors' names to their members diff --git a/src/eko/constants.py b/src/eko/constants.py index 06f96ce93..a9bb89027 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -2,25 +2,38 @@ r""" This files sets the physical constants. - -The constants are: - - - :math:`N_C` the number of colors - defaults to :math:`3` - - :math:`T_R` the normalization of fundamental generators - defaults to :math:`1/2` - - :math:`C_A` second Casimir constant in the adjoint representation - defaults to - :math:`N_C = 3` - - :math:`C_F` second Casimir constant in the fundamental representation - defaults to - :math:`\frac{N_C^2-1}{2N_C} = 4/3` """ NC = 3 """the number of colors""" TR = float(1.0 / 2.0) -"""the normalization of fundamental generators""" +"""the normalization of fundamental generators - defaults to :math:`T_R = 1/2`""" CA = float(NC) -"""second Casimir constant in the adjoint representation""" +"""second Casimir constant in the adjoint representation - defaults to :math:`N_C = 3`""" CF = float((NC * NC - 1.0) / (2.0 * NC)) -"""second Casimir constant in the fundamental representation""" +"""second Casimir constant in the fundamental representation - defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`""" + +eu2 = 4.0 / 9 +"""up quarks charge squared""" + +ed2 = 1.0 / 9 +"""down quarks charge squared""" + + +def update_colors(nc): + """ + Updates the number of colors to :math:`NC = nc` and the Casimirs for a generic value of :math:`NC` + + Parameters + ---------- + nc : int + Number of colors + """ + global NC, CA, CF # pylint: disable=global-statement + + NC = int(nc) + CA = float(NC) + CF = float((NC * NC - 1.0) / (2.0 * NC)) diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index ac18eb7b4..3ed1db0a5 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -15,9 +15,9 @@ from scipy import integrate from .. import anomalous_dimensions as ad -from .. import interpolation, mellin +from .. import basis_rotation as br +from .. import beta, interpolation, mellin from .. import scale_variations as sv -from ..basis_rotation import full_labels, singlet_labels from ..kernels import non_singlet as ns from ..kernels import singlet as s from ..member import OpMember @@ -32,25 +32,27 @@ ) -@nb.njit("c16(c16[:,:],string)") -def select_singlet_element(ker, mode): +@nb.njit("c16(c16[:,:],u2,u2)") +def select_singlet_element(ker, mode0, mode1): """ Select element of the singlet matrix Parameters ---------- - mode : str - sector element ker : numpy.ndarray singlet integration kernel - + mode0 : int + id for first sector element + mode1 : int + id for second sector element Returns ------- ker : complex singlet integration kernel element """ - k = 0 if mode[2] == "q" else 1 - l = 0 if mode[3] == "q" else 1 + + k = 0 if mode0 == 100 else 1 + l = 0 if mode1 == 100 else 1 return ker[k, l] @@ -79,8 +81,8 @@ class QuadKerBase: sector element """ - def __init__(self, u, is_log, logx, mode): - self.is_singlet = mode[0] == "S" + def __init__(self, u, is_log, logx, mode1): + self.is_singlet = mode1 != 0 self.is_log = is_log self.u = u self.logx = logx @@ -120,11 +122,12 @@ def integrand( return self.path.prefactor * pj * self.path.jac -@nb.njit("f8(f8,u1,string,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1,u1)", cache=True) +@nb.njit("f8(f8,u1,u2,u2,string,b1,f8,f8[:,:],f8,f8,f8,f8,u4,u1,u1)", cache=True) def quad_ker( u, order, - mode, + mode0, + mode1, method, is_log, logx, @@ -148,8 +151,10 @@ def quad_ker( perturbation order method : str method - mode : str - sector element + mode0: int + pid for first sector element + mode1 : int + pid for second sector element is_log : boolean is a logarithmic interpolation logx : float @@ -176,7 +181,7 @@ def quad_ker( ker : float evaluated integration kernel """ - ker_base = QuadKerBase(u, is_log, logx, mode) + ker_base = QuadKerBase(u, is_log, logx, mode1) integrand = ker_base.integrand(areas) if integrand == 0.0: return 0.0 @@ -197,9 +202,9 @@ def quad_ker( ker = np.ascontiguousarray(ker) @ np.ascontiguousarray( sv.expanded.singlet_variation(gamma_singlet, a1, order, nf, L) ) - ker = select_singlet_element(ker, mode) + ker = select_singlet_element(ker, mode0, mode1) else: - gamma_ns = ad.gamma_ns(order, mode[-1], ker_base.n, nf) + gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf) if sv_mode == sv.mode_exponentiated: gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L) ker = ns.dispatcher( @@ -297,16 +302,16 @@ def labels(self): logger.warning("Evolution: skipping non-singlet sector") else: # add + as default - labels.append("NS_p") + labels.append(br.non_singlet_labels[1]) if order >= 1: # - becomes different starting from NLO - labels.append("NS_m") + labels.append(br.non_singlet_labels[0]) if order >= 2: # v also becomes different starting from NNLO - labels.append("NS_v") + labels.append(br.non_singlet_labels[2]) # singlet sector is fixed if self.config["debug_skip_singlet"]: logger.warning("Evolution: skipping singlet sector") else: - labels.extend(singlet_labels) + labels.extend(br.singlet_labels) return labels @property @@ -320,10 +325,10 @@ def initialize_op_members(self): np.eye(self.grid_size), np.zeros((self.grid_size, self.grid_size)) ) zero = OpMember(*[np.zeros((self.grid_size, self.grid_size))] * 2) - for n in full_labels: + for n in br.full_labels: if n in self.labels: # off diag singlet are zero - if n in ["S_qg", "S_gq"]: + if n in br.singlet_labels and n[0] != n[1]: self.op_members[n] = zero.copy() else: self.op_members[n] = eye.copy() @@ -364,7 +369,8 @@ def run_op_integration( 1.0 - self._mellin_cut, args=( self.config["order"], - label, + label[0], + label[1], self.config["method"], self.int_disp.log, logx, @@ -448,9 +454,21 @@ def copy_ns_ops(self): """Copy non-singlet kernels, if necessary""" order = self.config["order"] if order == 0: # in LO +=-=v - for label in ["NS_v", "NS_m"]: - self.op_members[label].value = self.op_members["NS_p"].value.copy() - self.op_members[label].error = self.op_members["NS_p"].error.copy() + for label in ["nsV", "ns-"]: + self.op_members[ + (br.non_singlet_pids_map[label], 0) + ].value = self.op_members[ + (br.non_singlet_pids_map["ns+"], 0) + ].value.copy() + self.op_members[ + (br.non_singlet_pids_map[label], 0) + ].error = self.op_members[ + (br.non_singlet_pids_map["ns+"], 0) + ].error.copy() elif order == 1: # in NLO -=v - self.op_members["NS_v"].value = self.op_members["NS_m"].value.copy() - self.op_members["NS_v"].error = self.op_members["NS_m"].error.copy() + self.op_members[ + (br.non_singlet_pids_map["nsV"], 0) + ].value = self.op_members[(br.non_singlet_pids_map["ns-"], 0)].value.copy() + self.op_members[ + (br.non_singlet_pids_map["nsV"], 0) + ].error = self.op_members[(br.non_singlet_pids_map["ns-"], 0)].error.copy() diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py index 91548eac7..aed503520 100644 --- a/src/eko/evolution_operator/physical.py +++ b/src/eko/evolution_operator/physical.py @@ -48,21 +48,23 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range): """ # constant elements m = { - "S.S": op_members["S_qq"], - "S.g": op_members["S_qg"], - "g.g": op_members["S_gg"], - "g.S": op_members["S_gq"], - "V.V": op_members["NS_v"], + "S.S": op_members[(100, 100)], + "S.g": op_members[(100, 21)], + "g.g": op_members[(21, 21)], + "g.S": op_members[(21, 100)], + "V.V": op_members[(br.non_singlet_pids_map["nsV"], 0)], } # add elements which are already active for f in range(2, nf + 1): n = f**2 - 1 - m[f"V{n}.V{n}"] = op_members["NS_m"] - m[f"T{n}.T{n}"] = op_members["NS_p"] + m[f"V{n}.V{n}"] = op_members[(br.non_singlet_pids_map["ns-"], 0)] + m[f"T{n}.T{n}"] = op_members[(br.non_singlet_pids_map["ns+"], 0)] # deal with intrinsic heavy quark PDFs if intrinsic_range is not None: hqfl = "cbt" - op_id = member.OpMember.id_like(op_members["NS_v"]) + op_id = member.OpMember.id_like( + op_members[(br.non_singlet_pids_map["nsV"], 0)] + ) for intr_fl in intrinsic_range: if intr_fl <= nf: # light quarks are not intrinsic continue diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py index a0750f9e2..0229e75df 100644 --- a/src/eko/kernels/non_singlet.py +++ b/src/eko/kernels/non_singlet.py @@ -338,30 +338,30 @@ def dispatcher( if order == 0: return lo_exact(gamma_ns, a1, a0, nf) # NLO - elif order == 1: + if order == 1: if method in [ "iterate-expanded", "decompose-expanded", "perturbative-expanded", ]: return nlo_expanded(gamma_ns, a1, a0, nf) - elif method == "truncated": + if method == "truncated": return nlo_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) - elif method == "ordered-truncated": + if method == "ordered-truncated": return nlo_ordered_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) # if method in ["iterate-exact", "decompose-exact", "perturbative-exact"]: return nlo_exact(gamma_ns, a1, a0, nf) # NNLO - elif order == 2: + if order == 2: if method in [ "iterate-expanded", "decompose-expanded", "perturbative-expanded", ]: return nnlo_expanded(gamma_ns, a1, a0, nf) - elif method == "truncated": + if method == "truncated": return nnlo_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) - elif method == "ordered-truncated": + if method == "ordered-truncated": return nnlo_ordered_truncated(gamma_ns, a1, a0, nf, ev_op_iterations) # if method in ["iterate-exact", "decompose-exact", "perturbative-exact"]: return nnlo_exact(gamma_ns, a1, a0, nf) diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py index aa5d19f33..40ac12e0b 100644 --- a/src/eko/kernels/singlet.py +++ b/src/eko/kernels/singlet.py @@ -518,22 +518,22 @@ def dispatcher( # pylint: disable=too-many-return-statements # Common method for NLO and NNLO if method in ["iterate-exact", "iterate-expanded"]: return eko_iterate(gamma_singlet, a1, a0, nf, order, ev_op_iterations) - elif method == "perturbative-exact": + if method == "perturbative-exact": return eko_perturbative( gamma_singlet, a1, a0, nf, order, ev_op_iterations, ev_op_max_order, True ) - elif method == "perturbative-expanded": + if method == "perturbative-expanded": return eko_perturbative( gamma_singlet, a1, a0, nf, order, ev_op_iterations, ev_op_max_order, False ) - elif method in ["truncated", "ordered-truncated"]: + if method in ["truncated", "ordered-truncated"]: return eko_truncated(gamma_singlet, a1, a0, nf, order, ev_op_iterations) # These methods are scattered for nlo and nnlo - elif method == "decompose-exact": + if method == "decompose-exact": if order == 1: return nlo_decompose_exact(gamma_singlet, a1, a0, nf) return nnlo_decompose_exact(gamma_singlet, a1, a0, nf) - elif method == "decompose-expanded": + if method == "decompose-expanded": if order == 1: return nlo_decompose_expanded(gamma_singlet, a1, a0, nf) return nnlo_decompose_expanded(gamma_singlet, a1, a0, nf) diff --git a/src/eko/matching_conditions/__init__.py b/src/eko/matching_conditions/__init__.py index 96cf0356e..90324bb62 100644 --- a/src/eko/matching_conditions/__init__.py +++ b/src/eko/matching_conditions/__init__.py @@ -42,11 +42,11 @@ def split_ad_to_evol_map( """ m = { - "S.S": ome_members["S_qq"], - "S.g": ome_members["S_qg"], # This is always zero for the time being - "g.S": ome_members["S_gq"], - "g.g": ome_members["S_gg"], - "V.V": ome_members["NS_qq"], + "S.S": ome_members[(100, 100)], + "S.g": ome_members[(100, 21)], # This is always zero for the time being + "g.S": ome_members[(21, 100)], + "g.g": ome_members[(21, 21)], + "V.V": ome_members[(200, 200)], } # add elements which are already active @@ -59,15 +59,15 @@ def split_ad_to_evol_map( hq = br.quark_names[nf] m.update( { - f"{hq}-.V": ome_members["NS_Hq"], - f"{hq}+.S": ome_members["S_Hq"], - f"{hq}+.g": ome_members["S_Hg"], + f"{hq}-.V": ome_members[(br.matching_hminus_pid, 200)], + f"{hq}+.S": ome_members[(br.matching_hplus_pid, 100)], + f"{hq}+.g": ome_members[(br.matching_hplus_pid, 21)], } ) # intrinsic matching if len(intrinsic_range) != 0: - op_id = member.OpMember.id_like(ome_members["NS_qq"]) + op_id = member.OpMember.id_like(ome_members[(200, 200)]) for intr_fl in intrinsic_range: ihq = br.quark_names[intr_fl - 1] # find name if intr_fl > nf + 1: @@ -78,11 +78,15 @@ def split_ad_to_evol_map( # match the missing contribution from h+ and h- m.update( { - f"{ihq}+.{ihq}+": ome_members["S_HH"], - # f"S.{ihq}+": ome_members["S_qH"], - f"g.{ihq}+": ome_members["S_gH"], - f"{ihq}-.{ihq}-": ome_members["NS_HH"], - # f"V.{ihq}-": ome_members["NS_qH"], + f"{ihq}+.{ihq}+": ome_members[ + (br.matching_hplus_pid, br.matching_hplus_pid) + ], + # f"S.{ihq}+": ome_members[(100, br.matching_hplus_pid)], + f"g.{ihq}+": ome_members[(21, br.matching_hplus_pid)], + f"{ihq}-.{ihq}-": ome_members[ + (br.matching_hminus_pid, br.matching_hminus_pid) + ], + # f"V.{ihq}-": ome_members[(200, br.matching_hminus_pid)], } ) return cls.promote_names(m, q2_thr) diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index b610f10a0..978401234 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/src/eko/matching_conditions/operator_matrix_element.py @@ -11,9 +11,9 @@ import numpy as np from scipy import integrate +from .. import basis_rotation as br from .. import interpolation, mellin from ..anomalous_dimensions import harmonics -from ..basis_rotation import singlet_labels from ..member import OpMember from . import nlo, nnlo @@ -145,8 +145,10 @@ def build_ome(A, order, a_s, backward_method): return ome -@nb.njit("f8(f8,u1,string,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) -def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_msbar): +@nb.njit("f8(f8,u1,u2,u2,b1,f8,f8[:,:],f8,f8,string,b1)", cache=True) +def quad_ker( + u, order, mode0, mode1, is_log, logx, areas, a_s, L, backward_method, is_msbar +): """ Raw kernel inside quad @@ -156,8 +158,10 @@ def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_ms quad argument order : int perturbation order - mode : str - element in the singlet sector + mode0 : int + pid for first element in the singlet sector + mode1 : int + pid for second element in the singlet sector is_log : boolean logarithmic interpolation logx : float @@ -177,15 +181,15 @@ def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_ms ker : float evaluated integration kernel """ - is_singlet = mode[0] == "S" + is_singlet = mode0 in [100, 21, 90] # get transformation to N integral r = 0.4 * 16.0 / (1.0 - logx) if is_singlet: o = 1.0 - indeces = {"g": 0, "q": 1, "H": 2} + indices = {21: 0, 100: 1, 90: 2} else: o = 0.0 - indeces = {"q": 0, "H": 1} + indices = {200: 0, 91: 1} n = mellin.Talbot_path(u, r, o) jac = mellin.Talbot_jac(u, r, o) @@ -215,7 +219,7 @@ def quad_ker(u, order, mode, is_log, logx, areas, a_s, L, backward_method, is_ms ker = build_ome(A, order, a_s, backward_method) # select the need matrix element - ker = ker[indeces[mode[-2]], indeces[mode[-1]]] + ker = ker[indices[mode0], indices[mode1]] if ker == 0.0: return 0.0 @@ -267,10 +271,10 @@ def labels(self): if self.config["debug_skip_non_singlet"]: logger.warning("Matching: skipping non-singlet sector") else: - labels.extend(["NS_qq", "NS_Hq"]) + labels.extend([(200, 200), (br.matching_hminus_pid, 200)]) if self.is_intrinsic or self.backward_method != "": # intrisic labels, which are not zero at NLO - labels.append("NS_HH") + labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) # if self.backward_method == "exact": # # this contribution starts at NNLO, we don't have it for the moment # labels.append("NS_qH") @@ -279,9 +283,20 @@ def labels(self): if self.config["debug_skip_singlet"]: logger.warning("Matching: skipping singlet sector") else: - labels.extend([*singlet_labels, "S_Hg", "S_Hq"]) + labels.extend( + [ + *br.singlet_labels, + (br.matching_hplus_pid, 21), + (br.matching_hplus_pid, 100), + ] + ) if self.is_intrinsic or self.backward_method != "": - labels.extend(["S_gH", "S_HH"]) + labels.extend( + [ + (21, br.matching_hplus_pid), + (br.matching_hplus_pid, br.matching_hplus_pid), + ] + ) # if self.backward_method == "exact": # labels.extend(["S_qH"]) return labels @@ -306,7 +321,7 @@ def compute(self, q2, nf, L, is_msbar): grid_size = len(self.int_disp.xgrid) labels = self.labels() for n in labels: - if n[-1] == n[-2]: + if n[0] == n[1]: self.ome_members[n] = OpMember( np.eye(grid_size), np.zeros((grid_size, grid_size)) ) @@ -342,7 +357,8 @@ def compute(self, q2, nf, L, is_msbar): 1.0 - self._mellin_cut, args=( self.config["order"], - label, + label[0], + label[1], self.int_disp.log, logx, bf.areas_representation, @@ -375,7 +391,13 @@ def copy_ome(self): """Add the missing |OME|, if necessary""" grid_size = len(self.int_disp.xgrid) # basic labels skipped with skip debug - for label in ["NS_qq", "S_Hg", "S_Hq", "NS_Hq", *singlet_labels]: + for label in [ + (200, 200), + (br.matching_hplus_pid, 21), + (br.matching_hplus_pid, 100), + (br.matching_hminus_pid, 200), + *br.singlet_labels, + ]: if label not in self.ome_members: self.ome_members[label] = OpMember( np.zeros((grid_size, grid_size)), np.zeros((grid_size, grid_size)) @@ -383,7 +405,13 @@ def copy_ome(self): # intrinsic labels not computed yet if self.is_intrinsic: - for label in ["S_qH", "NS_qH", "NS_HH", "S_HH", "S_gH"]: + for label in [ + (100, br.matching_hplus_pid), + (200, br.matching_hminus_pid), + (br.matching_hminus_pid, br.matching_hminus_pid), + (br.matching_hplus_pid, br.matching_hplus_pid), + (21, br.matching_hplus_pid), + ]: if label not in self.ome_members: self.ome_members[label] = OpMember( np.zeros((grid_size, grid_size)), diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 0642ecdef..455a48b5d 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -23,7 +23,7 @@ def strong_coupling_mod_ev(mod_ev): """Map ModEv key to the available strong coupling evolution methods""" if mod_ev in ["EXA", "iterate-exact", "decompose-exact", "perturbative-exact"]: return "exact" - elif mod_ev in [ + if mod_ev in [ "TRN", "truncated", "ordered-truncated", diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index 40ccab8a7..1eabf6fe4 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -129,7 +129,7 @@ def run_external(self, theory, ocard, pdf): ocard, rotate_to_evolution_basis=self.rotate_to_evolution_basis, ) - elif self.external.lower() == "lhapdf": + if self.external.lower() == "lhapdf": from .external import lhapdf_utils # here theory card is not needed @@ -139,7 +139,7 @@ def run_external(self, theory, ocard, pdf): self.skip_pdfs(theory), rotate_to_evolution_basis=self.rotate_to_evolution_basis, ) - elif self.external.lower() == "pegasus": + if self.external.lower() == "pegasus": from .external import pegasus_utils return pegasus_utils.compute_pegasus_data( @@ -149,7 +149,7 @@ def run_external(self, theory, ocard, pdf): rotate_to_evolution_basis=self.rotate_to_evolution_basis, ) - elif self.external.lower() == "apfel": + if self.external.lower() == "apfel": from .external import apfel_utils return apfel_utils.compute_apfel_data( diff --git a/tests/eko/test_ad.py b/tests/eko/test_ad.py index 3363220c1..2627decd5 100644 --- a/tests/eko/test_ad.py +++ b/tests/eko/test_ad.py @@ -4,8 +4,9 @@ from numpy.testing import assert_allclose, assert_almost_equal from eko import anomalous_dimensions as ad +from eko import basis_rotation as br +from eko.anomalous_dimensions import as1 as ad_as1 from eko.anomalous_dimensions import harmonics -from eko.anomalous_dimensions import lo as ad_lo NF = 5 @@ -13,7 +14,7 @@ def test_eigensystem_gamma_singlet_0_values(): n = 3 s1 = harmonics.harmonic_S1(n) - gamma_S_0 = ad_lo.gamma_singlet_0(3, s1, NF) + gamma_S_0 = ad_as1.gamma_singlet(3, s1, NF) res = ad.exp_singlet(gamma_S_0) lambda_p = complex(12.273612971466964, 0) lambda_m = complex(5.015275917421917, 0) @@ -53,9 +54,15 @@ def test_eigensystem_gamma_singlet_projectors_EV(): def test_gamma_ns(): nf = 3 # LO - assert_almost_equal(ad.gamma_ns(2, "p", 1, nf)[0], 0.0) + assert_almost_equal(ad.gamma_ns(2, br.non_singlet_pids_map["ns+"], 1, nf)[0], 0.0) # NLO - assert_allclose(ad.gamma_ns(1, "m", 1, nf), np.zeros(2), atol=2e-6) + assert_allclose( + ad.gamma_ns(1, br.non_singlet_pids_map["ns-"], 1, nf), np.zeros(2), atol=2e-6 + ) # NNLO - assert_allclose(ad.gamma_ns(2, "m", 1, nf), np.zeros(3), atol=2e-4) - assert_allclose(ad.gamma_ns(2, "v", 1, nf), np.zeros(3), atol=8e-4) + assert_allclose( + ad.gamma_ns(2, br.non_singlet_pids_map["ns-"], 1, nf), np.zeros(3), atol=2e-4 + ) + assert_allclose( + ad.gamma_ns(2, br.non_singlet_pids_map["nsV"], 1, nf), np.zeros(3), atol=8e-4 + ) diff --git a/tests/eko/test_ad_lo.py b/tests/eko/test_ad_lo.py index 950116295..d08e79e32 100644 --- a/tests/eko/test_ad_lo.py +++ b/tests/eko/test_ad_lo.py @@ -2,7 +2,8 @@ # Test LO splitting functions import numpy as np -import eko.anomalous_dimensions.lo as ad_lo +import eko.anomalous_dimensions.aem1 as ad_aem1 +import eko.anomalous_dimensions.as1 as ad_as1 from eko.anomalous_dimensions import harmonics NF = 5 @@ -12,14 +13,21 @@ def test_number_conservation(): # number N = complex(1.0, 0.0) s1 = harmonics.harmonic_S1(N) - np.testing.assert_almost_equal(ad_lo.gamma_ns_0(N, s1), 0) + np.testing.assert_almost_equal(ad_as1.gamma_ns(N, s1), 0) def test_quark_momentum_conservation(): # quark momentum N = complex(2.0, 0.0) s1 = harmonics.harmonic_S1(N) - np.testing.assert_almost_equal(ad_lo.gamma_ns_0(N, s1) + ad_lo.gamma_gq_0(N), 0) + np.testing.assert_almost_equal( + ad_as1.gamma_ns(N, s1) + ad_as1.gamma_gq(N), + 0, + ) + np.testing.assert_almost_equal( + ad_aem1.gamma_ns(N, s1) + ad_aem1.gamma_phq(N), + 0, + ) def test_gluon_momentum_conservation(): @@ -27,24 +35,47 @@ def test_gluon_momentum_conservation(): N = complex(2.0, 0.0) s1 = harmonics.harmonic_S1(N) np.testing.assert_almost_equal( - ad_lo.gamma_qg_0(N, NF) + ad_lo.gamma_gg_0(N, s1, NF), 0 + ad_as1.gamma_qg(N, NF) + ad_as1.gamma_gg(N, s1, NF), 0 ) +def test_photon_momentum_conservation(): + # gluon momentum + N = complex(2.0, 0.0) + np.testing.assert_almost_equal(ad_aem1.gamma_qph(N, NF) + ad_aem1.gamma_phph(NF), 0) + + def test_gamma_qg_0(): N = complex(1.0, 0.0) res = complex(-20.0 / 3.0, 0.0) - np.testing.assert_almost_equal(ad_lo.gamma_qg_0(N, NF), res) + np.testing.assert_almost_equal(ad_as1.gamma_qg(N, NF), res) def test_gamma_gq_0(): N = complex(0.0, 1.0) res = complex(4.0, -4.0) / 3.0 - np.testing.assert_almost_equal(ad_lo.gamma_gq_0(N), res) + np.testing.assert_almost_equal(ad_as1.gamma_gq(N), res) def test_gamma_gg_0(): N = complex(0.0, 1.0) s1 = harmonics.harmonic_S1(N) res = complex(5.195725159621, 10.52008856962) - np.testing.assert_almost_equal(ad_lo.gamma_gg_0(N, s1, NF), res) + np.testing.assert_almost_equal(ad_as1.gamma_gg(N, s1, NF), res) + + +def test_gamma_phq_0(): + N = complex(0.0, 1.0) + res = complex(4.0, -4.0) / 3.0 / 4 * 3 + np.testing.assert_almost_equal(ad_aem1.gamma_phq(N), res) + + +def test_gamma_qph_0(): + N = complex(1.0, 0.0) + res = complex(-20.0 / 3.0, 0.0) * 3 / 0.5 + np.testing.assert_almost_equal(ad_aem1.gamma_qph(N, NF), res) + + +def test_gamma_phph_0(): + res = complex(2.0 / 3 * 3 * 2 * NF, 0.0) + np.testing.assert_almost_equal(ad_aem1.gamma_phph(NF), res) diff --git a/tests/eko/test_ad_nlo.py b/tests/eko/test_ad_nlo.py index 82033d3a7..23a541fae 100644 --- a/tests/eko/test_ad_nlo.py +++ b/tests/eko/test_ad_nlo.py @@ -2,18 +2,18 @@ # Test NLO anomalous dims import numpy as np +import eko.anomalous_dimensions.as2 as ad_as2 import eko.anomalous_dimensions.harmonics as h -import eko.anomalous_dimensions.nlo as ad_nlo -from eko.constants import CA, CF +from eko import constants as const NF = 5 def test_gamma_1(): # number conservation - np.testing.assert_allclose(ad_nlo.gamma_nsm_1(1, NF), 0.0, atol=2e-6) + np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF), 0.0, atol=2e-6) - gS1 = ad_nlo.gamma_singlet_1(2, NF) + gS1 = ad_as2.gamma_singlet(2, NF) # gluon momentum conservation # the CA*NF term seems to be tough to compute, so raise the constraint ... np.testing.assert_allclose(gS1[0, 1] + gS1[1, 1], 0.0, atol=4e-5) @@ -25,55 +25,88 @@ def test_gamma_1(): # reference values are obtained from MMa # Non singlet sector np.testing.assert_allclose( - ad_nlo.gamma_nsp_1(2, NF), (-112.0 * CF + 376.0 * CA - 64.0 * NF) * CF / 27.0 + ad_as2.gamma_nsp(2, NF), + (-112.0 * const.CF + 376.0 * const.CA - 64.0 * NF) * const.CF / 27.0, ) # singlet sector - np.testing.assert_allclose(ad_nlo.gamma_ps_1(2, NF), -40.0 * CF * NF / 27.0) - np.testing.assert_allclose(gS1[0, 1], (-74.0 * CF - 35.0 * CA) * NF / 27.0) # qg + np.testing.assert_allclose(ad_as2.gamma_ps(2, NF), -40.0 * const.CF * NF / 27.0) np.testing.assert_allclose( - gS1[1, 0], (112.0 * CF - 376.0 * CA + 104.0 * NF) * CF / 27.0 + gS1[0, 1], (-74.0 * const.CF - 35.0 * const.CA) * NF / 27.0 + ) # qg + np.testing.assert_allclose( + gS1[1, 0], (112.0 * const.CF - 376.0 * const.CA + 104.0 * NF) * const.CF / 27.0 ) # gq # add additional point at (analytical) continuation point np.testing.assert_allclose( - ad_nlo.gamma_nsm_1(2, NF), + ad_as2.gamma_nsm(2, NF), ( - (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.zeta3) * CF - + (373.0 / 9.0 - 34.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * CA + (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.zeta3) * const.CF + + (373.0 / 9.0 - 34.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * const.CA - 64.0 * NF / 27.0 ) - * CF, + * const.CF, ) np.testing.assert_allclose( - ad_nlo.gamma_nsp_1(3, NF), + ad_as2.gamma_nsp(3, NF), ( - (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.zeta3) * CF - + (459.0 / 8.0 - 43.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * CA + (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.zeta3) * const.CF + + (459.0 / 8.0 - 43.0 * np.pi**2 / 9.0 + 8.0 * h.zeta3) * const.CA - 415.0 * NF / 108.0 ) - * CF, + * const.CF, ) - np.testing.assert_allclose(ad_nlo.gamma_ps_1(3, NF), -1391.0 * CF * NF / 5400.0) - gS1 = ad_nlo.gamma_singlet_1(3, NF) + np.testing.assert_allclose(ad_as2.gamma_ps(3, NF), -1391.0 * const.CF * NF / 5400.0) + gS1 = ad_as2.gamma_singlet(3, NF) + np.testing.assert_allclose( + gS1[1, 0], + ( + 973.0 / 432.0 * const.CF + + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * const.CA + + 61.0 / 54.0 * NF + ) + * const.CF, + ) # gq + np.testing.assert_allclose( + gS1[1, 1], + ( + (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) + * const.CA**2 + - 967.0 / 270.0 * const.CA * NF + + 541.0 / 216.0 * const.CF * NF + ), + rtol=6e-7, + ) # gg + gS1 = ad_as2.gamma_singlet(4, NF) + np.testing.assert_allclose( + gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF + ) # qg + + const.update_colors(4) + + np.testing.assert_allclose(const.CA, 4.0) + gS1 = ad_as2.gamma_singlet(3, NF) np.testing.assert_allclose( gS1[1, 0], ( - 973.0 / 432.0 * CF - + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * CA + 973.0 / 432.0 * const.CF + + (2801.0 / 5400.0 - 7.0 * np.pi**2 / 9.0) * const.CA + 61.0 / 54.0 * NF ) - * CF, + * const.CF, ) # gq np.testing.assert_allclose( gS1[1, 1], ( - (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) * CA**2 - - 967.0 / 270.0 * CA * NF - + 541.0 / 216.0 * CF * NF + (-79909.0 / 3375.0 + 194.0 * np.pi**2 / 45.0 - 8.0 * h.zeta3) + * const.CA**2 + - 967.0 / 270.0 * const.CA * NF + + 541.0 / 216.0 * const.CF * NF ), rtol=6e-7, ) # gg - gS1 = ad_nlo.gamma_singlet_1(4, NF) + gS1 = ad_as2.gamma_singlet(4, NF) np.testing.assert_allclose( - gS1[0, 1], (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * NF + gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF ) # qg + const.update_colors(3) diff --git a/tests/eko/test_ad_nnlo.py b/tests/eko/test_ad_nnlo.py index 923d4a3c7..58381478a 100644 --- a/tests/eko/test_ad_nnlo.py +++ b/tests/eko/test_ad_nnlo.py @@ -2,7 +2,7 @@ # Test NNLO anomalous dims import numpy as np -import eko.anomalous_dimensions.nnlo as ad_nnlo +import eko.anomalous_dimensions.as3 as ad_as3 from eko.anomalous_dimensions import harmonics NF = 5 @@ -26,13 +26,13 @@ def test_gamma_2(): # number conservation - each is 0 on its own, see :cite:`Moch:2004pa` N = 1 sx = get_sx(N) - np.testing.assert_allclose(ad_nnlo.gamma_nsv_2(N, NF, sx), 0.000960586, rtol=3e-7) - np.testing.assert_allclose(ad_nnlo.gamma_nsm_2(N, NF, sx), -0.000594225, rtol=6e-7) + np.testing.assert_allclose(ad_as3.gamma_nsv(N, NF, sx), 0.000960586, rtol=3e-7) + np.testing.assert_allclose(ad_as3.gamma_nsm(N, NF, sx), -0.000594225, rtol=6e-7) # get singlet sector N = 2 sx = get_sx(N) - gS2 = ad_nnlo.gamma_singlet_2(N, NF, sx) + gS2 = ad_as3.gamma_singlet(N, NF, sx) # gluon momentum conservation np.testing.assert_allclose(gS2[0, 1] + gS2[1, 1], 0.00388726, rtol=2e-6) @@ -42,4 +42,4 @@ def test_gamma_2(): assert gS2.shape == (2, 2) # test nsv_2 equal to referece value - np.testing.assert_allclose(ad_nnlo.gamma_nsv_2(N, NF, sx), -188.325593, rtol=3e-7) + np.testing.assert_allclose(ad_as3.gamma_nsv(N, NF, sx), -188.325593, rtol=3e-7) diff --git a/tests/eko/test_basis_rotation.py b/tests/eko/test_basis_rotation.py index c55a49f78..faf320b73 100644 --- a/tests/eko/test_basis_rotation.py +++ b/tests/eko/test_basis_rotation.py @@ -10,19 +10,19 @@ def test_ad_projector(): g = br.rotate_flavor_to_evolution[2] v3 = br.rotate_flavor_to_evolution[br.evol_basis.index("V3")] - s_to_s = br.ad_projector("S_qq", nf=6) + s_to_s = br.ad_projector((100, 100), nf=6) np.testing.assert_allclose(s @ s_to_s, s) np.testing.assert_allclose(g @ s_to_s, 0.0) np.testing.assert_allclose(v3 @ s_to_s, 0.0) - g_to_s = br.ad_projector("S_gq", nf=6) + g_to_s = br.ad_projector((21, 100), nf=6) np.testing.assert_allclose(s @ g_to_s, 0.0) np.testing.assert_allclose(g @ g_to_s, s) np.testing.assert_allclose(v3 @ g_to_s, 0.0) - ns_m = br.ad_projector("NS_m", nf=6) + ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6) np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15) np.testing.assert_allclose(g @ ns_m, 0.0) diff --git a/tests/eko/test_ev_operator.py b/tests/eko/test_ev_operator.py index e242f5e17..556856a63 100644 --- a/tests/eko/test_ev_operator.py +++ b/tests/eko/test_ev_operator.py @@ -30,7 +30,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=0, - mode="NS_p", + mode0=br.non_singlet_pids_map["ns+"], + mode1=0, method="", is_log=is_log, logx=0.0, @@ -47,7 +48,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=0, - mode="S_qq", + mode0=100, + mode1=100, method="", is_log=is_log, logx=1.0, @@ -64,7 +66,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=0, - mode="S_qg", + mode0=100, + mode1=21, method="", is_log=is_log, logx=0.0, @@ -78,12 +81,13 @@ def test_quad_ker(monkeypatch): sv_mode=0, ) np.testing.assert_allclose(res_s, 0.0) - for mode in ["NS_p", "S_qq"]: + for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]: for sv in [1, 2]: res_sv = quad_ker( u=0, order=0, - mode=mode, + mode0=label[0], + mode1=label[1], method="", is_log=True, logx=1.0, @@ -102,7 +106,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=0, - mode="NS_p", + mode0=br.non_singlet_pids_map["ns+"], + mode1=0, method="", is_log=True, logx=0.0, @@ -127,9 +132,7 @@ def test_labels(self): 1, 2, ) - assert sorted(o.labels) == sorted( - ["NS_p", "NS_m", "NS_v", "S_qq", "S_qg", "S_gq", "S_gg"] - ) + assert sorted(o.labels) == sorted(br.full_labels) o = Operator( dict(order=1, debug_skip_non_singlet=True, debug_skip_singlet=True), {}, @@ -190,19 +193,25 @@ def test_compute(self, monkeypatch): ) # LO o.compute() - assert "NS_m" in o.op_members + assert (br.non_singlet_pids_map["ns-"], 0) in o.op_members np.testing.assert_allclose( - o.op_members["NS_m"].value, o.op_members["NS_p"].value + o.op_members[(br.non_singlet_pids_map["ns-"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, ) np.testing.assert_allclose( - o.op_members["NS_v"].value, o.op_members["NS_p"].value + o.op_members[(br.non_singlet_pids_map["nsV"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, ) # NLO o.config["order"] = 1 o.compute() - assert not np.allclose(o.op_members["NS_p"].value, o.op_members["NS_m"].value) + assert not np.allclose( + o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns-"], 0)].value, + ) np.testing.assert_allclose( - o.op_members["NS_v"].value, o.op_members["NS_m"].value + o.op_members[(br.non_singlet_pids_map["nsV"], 0)].value, + o.op_members[(br.non_singlet_pids_map["ns-"], 0)].value, ) # unity operators @@ -217,13 +226,13 @@ def test_compute(self, monkeypatch): def test_pegasus_path(): def quad_ker_pegasus( - u, order, mode, method, logx, areas, a1, a0, nf, ev_op_iterations + u, order, mode0, method, logx, areas, a1, a0, nf, ev_op_iterations ): # compute the mellin inversion as done in pegasus phi = 3 / 4 * np.pi c = 1.9 n = complex(c + u * np.exp(1j * phi)) - gamma_ns = ad.gamma_ns(order, mode[-1], n, nf) + gamma_ns = ad.gamma_ns(order, mode0, n, nf) ker = ns.dispatcher( order, method, @@ -241,7 +250,8 @@ def quad_ker_pegasus( xgrid = np.geomspace(1e-7, 1, 10) int_disp = InterpolatorDispatcher(xgrid, 1, True) order = 1 - mode = "NS_p" + mode0 = br.non_singlet_pids_map["ns+"] + mode1 = 0 method = "" logxs = np.log(int_disp.xgrid_raw) a1 = 1 @@ -257,7 +267,8 @@ def quad_ker_pegasus( 1.0, args=( order, - mode, + mode0, + mode1, method, int_disp.log, logx, @@ -282,7 +293,7 @@ def quad_ker_pegasus( np.inf, args=( order, - mode, + mode0, method, logx, bf.areas_representation, diff --git a/tests/eko/test_kernels_ns.py b/tests/eko/test_kernels_ns.py index 3f2a7dd9e..f97d8e7eb 100644 --- a/tests/eko/test_kernels_ns.py +++ b/tests/eko/test_kernels_ns.py @@ -2,6 +2,7 @@ import numpy as np import pytest +from eko import anomalous_dimensions as ad from eko import beta from eko.kernels import non_singlet as ns @@ -114,6 +115,8 @@ def test_ode_nlo(): def test_error(): with pytest.raises(NotImplementedError): ns.dispatcher(3, "iterate-exact", np.random.rand(3) + 0j, 0.2, 0.1, 3, 10) + with pytest.raises(NotImplementedError): + ad.gamma_ns(1, 10202, 1, 3) def test_gamma_usage(): diff --git a/tests/eko/test_matching.py b/tests/eko/test_matching.py index 02f1d0d04..6aecbc3e8 100644 --- a/tests/eko/test_matching.py +++ b/tests/eko/test_matching.py @@ -3,6 +3,7 @@ import numpy as np from numpy.testing import assert_almost_equal +from eko import basis_rotation as br from eko import member from eko.matching_conditions import MatchingCondition @@ -17,17 +18,25 @@ class TestMatchingCondition: def mkOME(self): ome = {} - for key in ["qq", "qg", "gq", "gg", "Hq", "Hg"]: - ome.update({f"S_{key}": mkOM(self.shape)}) - if "g" not in key: - ome.update({f"NS_{key}": mkOM(self.shape)}) + for key in [ + *br.singlet_labels, + (br.matching_hplus_pid, 100), + (br.matching_hplus_pid, 21), + (200, 200), + (br.matching_hminus_pid, 200), + ]: + ome.update({key: mkOM(self.shape)}) return ome def update_intrinsic_OME(self, ome): - for key in ["HH", "qH", "gH"]: - ome.update({f"S_{key}": mkOM(self.shape)}) - if "g" not in key: - ome.update({f"NS_{key}": mkOM(self.shape)}) + for key in [ + (br.matching_hplus_pid, br.matching_hplus_pid), + (br.matching_hminus_pid, br.matching_hminus_pid), + (200, br.matching_hminus_pid), + (100, br.matching_hplus_pid), + (21, br.matching_hplus_pid), + ]: + ome.update({key: mkOM(self.shape)}) def test_split_ad_to_evol_map(self): ome = self.mkOME() @@ -54,7 +63,7 @@ def test_split_ad_to_evol_map(self): ) assert_almost_equal( a.op_members[member.MemberName("V.V")].value, - ome["NS_qq"].value, + ome[(200, 200)].value, ) # # if alpha is zero, nothing non-trivial should happen b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, []) @@ -107,13 +116,13 @@ def test_split_ad_to_evol_map(self): ) assert_almost_equal( d.op_members[member.MemberName("b-.V")].value, - ome["NS_Hq"].value, + ome[(br.matching_hminus_pid, 200)].value, ) assert_almost_equal( d.op_members[member.MemberName("b+.S")].value, - ome["S_Hq"].value, + ome[(br.matching_hplus_pid, 100)].value, ) assert_almost_equal( d.op_members[member.MemberName("b+.g")].value, - ome["S_Hg"].value, + ome[(br.matching_hplus_pid, 21)].value, ) diff --git a/tests/eko/test_ome.py b/tests/eko/test_ome.py index dc8b6355c..1ded8b602 100644 --- a/tests/eko/test_ome.py +++ b/tests/eko/test_ome.py @@ -4,8 +4,8 @@ import numpy as np +from eko import basis_rotation as br from eko import interpolation, mellin -from eko.basis_rotation import singlet_labels from eko.evolution_operator.grid import OperatorGrid from eko.interpolation import InterpolatorDispatcher from eko.matching_conditions.operator_matrix_element import ( @@ -93,7 +93,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=2, - mode="NS_qq", + mode0=200, + mode1=200, is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -106,7 +107,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=2, - mode="S_qq", + mode0=100, + mode1=100, is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -119,7 +121,8 @@ def test_quad_ker(monkeypatch): res_s = quad_ker( u=0, order=2, - mode="S_qg", + mode0=100, + mode1=21, is_log=is_log, logx=0.0, areas=np.zeros(3), @@ -131,12 +134,13 @@ def test_quad_ker(monkeypatch): np.testing.assert_allclose(res_s, 0.0) # test expanded intrisic inverse kernels - labels = ["NS_qq", *singlet_labels] + labels = [(200, 200), *br.singlet_labels] for label in labels: res_ns = quad_ker( u=0, order=2, - mode=label, + mode0=label[0], + mode1=label[1], is_log=True, logx=0.0, areas=np.zeros(3), @@ -153,21 +157,22 @@ def test_quad_ker(monkeypatch): # test exact intrinsic inverse kernel labels.extend( [ - "S_Hq", - "S_Hg", - "S_HH", - "S_qH", - "S_gH", - "NS_qH", - "NS_HH", - "NS_Hq", + (br.matching_hplus_pid, 100), + (br.matching_hplus_pid, 21), + (br.matching_hplus_pid, br.matching_hplus_pid), + (100, br.matching_hplus_pid), + (21, br.matching_hplus_pid), + (200, br.matching_hminus_pid), + (br.matching_hminus_pid, br.matching_hminus_pid), + (br.matching_hminus_pid, 200), ] ) for label in labels: res_ns = quad_ker( u=0, order=2, - mode=label, + mode0=label[0], + mode1=label[1], is_log=True, logx=0.0, areas=np.zeros(3), @@ -185,7 +190,8 @@ def test_quad_ker(monkeypatch): res_ns = quad_ker( u=0, order=2, - mode="NS_qq", + mode0=200, + mode1=200, is_log=True, logx=0.0, areas=np.array([0.01, 0.1, 1.0]), @@ -246,13 +252,19 @@ def test_labels(self): ) o = OperatorMatrixElement(g.config, g.managers, is_backward=False) labels = o.labels() - test_labels = ["NS_qq", "NS_Hq"] + test_labels = [(200, 200), (br.matching_hminus_pid, 200)] for l in test_labels: if skip_ns: assert l not in labels else: assert l in labels - test_labels = ["S_qq", "S_Hq", "S_gg", "S_Hg", "S_gH"] + test_labels = [ + (100, 100), + (br.matching_hplus_pid, 100), + (21, 21), + (br.matching_hplus_pid, 21), + (21, br.matching_hplus_pid), + ] for l in test_labels: if skip_singlet: assert l not in labels @@ -281,24 +293,29 @@ def test_compute_lo(self): o = OperatorMatrixElement(g.config, g.managers, is_backward=False) o.compute(self.theory_card["mb"] ** 2, nf=4, L=0, is_msbar=False) - dim = o.ome_members["NS_qq"].value.shape - for idx in ["S", "NS"]: + dim = o.ome_members[(200, 200)].value.shape + for indices in [(100, br.matching_hplus_pid), (200, br.matching_hminus_pid)]: np.testing.assert_allclose( - o.ome_members[f"{idx}_qq"].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[(indices[0], indices[0])].value, np.eye(dim[0]), atol=1e-8 ) np.testing.assert_allclose( - o.ome_members[f"{idx}_HH"].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[(indices[1], indices[1])].value, np.eye(dim[0]), atol=1e-8 + ) + np.testing.assert_allclose( + o.ome_members[(indices[0], indices[1])].value, np.zeros(dim) + ) + np.testing.assert_allclose( + o.ome_members[(indices[1], indices[0])].value, np.zeros(dim) ) - np.testing.assert_allclose(o.ome_members[f"{idx}_qH"].value, np.zeros(dim)) - np.testing.assert_allclose(o.ome_members[f"{idx}_Hq"].value, np.zeros(dim)) np.testing.assert_allclose( - o.ome_members["S_gg"].value, np.eye(dim[0]), atol=1e-8 + o.ome_members[(21, 21)].value, np.eye(dim[0]), atol=1e-8 ) np.testing.assert_allclose( - o.ome_members["S_qg"].value, o.ome_members["S_gq"].value + o.ome_members[100, 21].value, o.ome_members[(21, 100)].value ) np.testing.assert_allclose( - o.ome_members["S_Hg"].value, o.ome_members["S_gH"].value + o.ome_members[(br.matching_hplus_pid, 21)].value, + o.ome_members[(21, br.matching_hplus_pid)].value, ) def test_compute_nlo(self): @@ -327,20 +344,20 @@ def test_compute_nlo(self): dim = len(operators_card["interpolation_xgrid"]) shape = (dim, dim) - for idx in ["S", "NS"]: - assert o.ome_members[f"{idx}_qq"].value.shape == shape - assert o.ome_members[f"{idx}_HH"].value.shape == shape - assert o.ome_members[f"{idx}_qH"].value.shape == shape - assert o.ome_members[f"{idx}_Hq"].value.shape == shape + for indices in [(100, br.matching_hplus_pid), (200, br.matching_hminus_pid)]: + assert o.ome_members[(indices[0], indices[0])].value.shape == shape + assert o.ome_members[(indices[1], indices[1])].value.shape == shape + assert o.ome_members[(indices[0], indices[1])].value.shape == shape + assert o.ome_members[(indices[1], indices[0])].value.shape == shape np.testing.assert_allclose( - o.ome_members[f"{idx}_qH"].value, np.zeros(shape) + o.ome_members[(indices[0], indices[1])].value, np.zeros(shape) ) np.testing.assert_allclose( - o.ome_members[f"{idx}_Hq"].value, np.zeros(shape) + o.ome_members[(indices[1], indices[0])].value, np.zeros(shape) ) - assert o.ome_members["S_gg"].value.shape == shape + assert o.ome_members[(21, 21)].value.shape == shape np.testing.assert_allclose( - o.ome_members["S_qg"].value, o.ome_members["S_gq"].value + o.ome_members[(100, 21)].value, o.ome_members[(21, 100)].value ) - assert o.ome_members["S_Hg"].value.shape == shape - assert o.ome_members["S_gH"].value.shape == shape + assert o.ome_members[(br.matching_hplus_pid, 21)].value.shape == shape + assert o.ome_members[(21, br.matching_hplus_pid)].value.shape == shape diff --git a/tests/eko/test_runner.py b/tests/eko/test_runner.py index 1523e4034..e74cae619 100644 --- a/tests/eko/test_runner.py +++ b/tests/eko/test_runner.py @@ -116,3 +116,22 @@ def check_shapes(o, txs, ixs, theory_card, operators_card): assert "operator_errors" in ops assert ops["operators"].shape == op_shape assert ops["operator_errors"].shape == op_shape + + +def test_vfns(): + # change targetbasis + tc = copy.deepcopy(theory_card) + oc = copy.deepcopy(operators_card) + tc["kcThr"] = 1.0 + tc["kbThr"] = 1.0 + tc["PTO"] = 2 + oc["debug_skip_non_singlet"] = False + r = eko.runner.Runner(tc, oc) + o = r.get_output() + check_shapes( + o, + o["interpolation_xgrid"], + o["interpolation_xgrid"], + tc, + oc, + ) diff --git a/tests/eko/test_sv_expanded.py b/tests/eko/test_sv_expanded.py index 076df98f0..df3cf6f2b 100644 --- a/tests/eko/test_sv_expanded.py +++ b/tests/eko/test_sv_expanded.py @@ -2,6 +2,7 @@ import numpy as np +from eko import basis_rotation as br from eko.anomalous_dimensions import gamma_ns, gamma_singlet from eko.beta import beta_0 from eko.kernels import non_singlet, singlet @@ -71,7 +72,7 @@ def scheme_diff(g, k, pto, is_singlet): for L in [np.log(0.5), np.log(2)]: for order in [1, 2]: # Non singlet kernels - gns = gamma_ns(order, "p", n, nf) + gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf) ker = non_singlet.dispatcher( order, method, gns, a1, a0, nf, ev_op_iterations=1 )