Merge pull request #135 from NNPDF/feature/qed-solution

Implement QED exact solution

benchmarks/CT14_bench.py

@@ -5,16 +5,18 @@
   whereas the lo family uses LO splitting functions with NLO alphas evolution
 """
 
+from math import nan
+
 from banana import register
 
-from eko import compatibility
 from ekomark.benchmark.runner import Runner
 
 register(__file__)
 
 
 base_theory = {
     "Qref": 91.1876,
+    "QrefQED": nan,
     "mc": 1.3,
     "mb": 4.75,
     "mt": 172,

benchmarks/CT18_bench.py

@@ -12,6 +12,7 @@
 
 base_theory = {
     "Qref": 91.1870,
+    "QrefQED": 91.1870,
     "mc": 1.3,
     "mb": 4.75,
     "mt": 172.0,
@@ -53,6 +54,28 @@ def benchmark_nnlo(self, Q0=1.295, Q2grid=(1e4,)):
         ]
         self.run([theory_card], [operator_card], ["CT18NNLO"])
 
+    def benchmark_nnlo_qed(self, Q0=1.295, Q2grid=(1e4,)):
+        theory_card = base_theory.copy()
+        theory_card.update(
+            {
+                "alphas": 0.118000,
+                "alphaqed": 0.007496,
+                "PTO": 2,
+                "QED": 1,
+                "Q0": Q0,
+                "MaxNfPdf": 5,
+                "MaxNfAs": 5,
+            }
+        )
+        operator_card = {"Q2grid": list(Q2grid)}
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run([theory_card], [operator_card], ["CT18qed"])
+
     def benchmark_znnlo(self, Q0=1.3, Q2grid=(1e4,)):
         theory_card = base_theory.copy()
         theory_card.update(

benchmarks/HERA20_bench.py

@@ -2,6 +2,8 @@
 Benchmark HERAPDF2.0 pdf family
 
 """
+from math import nan
+
 from banana import register
 
 from eko import interpolation
@@ -12,6 +14,7 @@
 
 base_theory = {
     "Qref": 91.1876,
+    "QrefQED": nan,
     "mc": 1.43,
     "mb": 4.5,
     "mt": 173.0,

benchmarks/NNPDF_bench.py

@@ -63,6 +63,33 @@ def benchmark_nlo(self, Q0=1.65, Q2grid=(100,)):
         self.run([theory_card], [operator_card], ["NNPDF31_nlo_as_0118"])
 
 
+class BenchmarkNNPDF31_luxqed(BenchmarkNNPDF):
+    """Benchmark NNPDF3.1_luxqed"""
+
+    def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)):
+        theory_card = {
+            **base_theory,
+            "PTO": 2,
+            "QED": 2,
+            "Q0": Q0,
+        }
+        theory_card.update({"ModEv": "iterate-exact", "FNS": "VFNS", "QrefQED": 91.2})
+
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+
+        operator_card = {
+            **base_operator,
+            "Q2grid": list(Q2grid),
+            "ev_op_iterations": 10,
+        }
+        self.run([theory_card], [operator_card], ["NNPDF31_nnlo_as_0118_luxqed"])
+
+
 class BenchmarkNNPDF40(BenchmarkNNPDF):
     """Benchmark NNPDF4.0"""
 
@@ -128,8 +155,10 @@ def benchmark(self, Q0=1.65, Q2grid=(100,)):
     # nn31.benchmark_nlo(Q0=np.sqrt(low2), Q2grid=[10])
     # # test backward
     # #nn31.benchmark_nlo(Q0=np.sqrt(high2), Q2grid=[low2])
-    # nn40 = BenchmarkNNPDF40()
-    # # nn40.benchmark_nnlo(Q2grid=[100])
+    nn40 = BenchmarkNNPDF40()
+    nn40.benchmark_nnlo(Q2grid=[100])
     # nn40.benchmark_nnlo(Q0=np.sqrt(high2), Q2grid=[low2])
-    nnpol = BenchmarkNNPDFpol11()
-    nnpol.benchmark(Q0=np.sqrt(low2), Q2grid=[high2])
+    # nnpol = BenchmarkNNPDFpol11()
+    # nnpol.benchmark(Q0=np.sqrt(low2), Q2grid=[high2])
+    # obj = BenchmarkNNPDF31_luxqed()
+    # obj.benchmark_nnlo(Q0=5.0)

benchmarks/apfel_bench.py

@@ -197,13 +197,137 @@ def benchmark_sv(self, pto, svmode):
         self.run(ts, operators.build(operators.apfel_config), ["ToyLH"])
 
 
-if __name__ == "__main__":
+class BenchmarkFFNS_qed(ApfelBenchmark):
+    """Benckmark FFNS"""
+
+    ffns_theory = {
+        "Qref": 91.1870,
+        "QrefQED": 91.1870,
+        "mc": 1.3,
+        "mb": 4.75,
+        "mt": 172.0,
+        "FNS": "FFNS",
+        "ModEv": [
+            "EXA",
+            # "EXP",
+            # "TRN",
+        ],
+        "NfFF": 5,
+        "kcThr": 0.0,
+        "kbThr": 0.0,
+        "ktThr": np.inf,
+        "Q0": 5.0,
+        "alphas": 0.118000,
+        "alphaqed": 0.007496,
+    }
+    ffns_theory = tolist(ffns_theory)
+
+    def benchmark_plain(self, pto, qed):
+        """Plain configuration"""
+
+        th = self.ffns_theory.copy()
+        th.update({"PTO": [pto], "QED": [qed]})
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run(
+            cartesian_product(th),
+            operators.build(operators.apfel_config),
+            ["NNPDF31_nnlo_as_0118_luxqed"],
+        )
+
+    def benchmark_sv(self, pto, qed, svmode):
+        """Scale Variation"""
 
+        ts = []
+        th = self.ffns_theory.copy()
+        th.update(
+            {
+                "PTO": [pto],
+                "QED": [qed],
+                "XIR": [np.sqrt(0.5)],
+                "fact_to_ren_scale_ratio": [np.sqrt(2.0)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
+            }
+        )
+        ts.extend(cartesian_product(th))
+        th = self.ffns_theory.copy()
+        th.update(
+            {
+                "PTO": [pto],
+                "QED": [qed],
+                "XIR": [np.sqrt(2.0)],
+                "fact_to_ren_scale_ratio": [np.sqrt(0.5)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
+            }
+        )
+        ts.extend(cartesian_product(th))
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run(ts, operators.build(operators.apfel_config), ["ToyLH"])
+
+
+class BenchmarkVFNS_qed(ApfelBenchmark):
+    """Benckmark FFNS"""
+
+    vfns_theory = {
+        "Qref": 91.1870,
+        "QrefQED": 91.1870,
+        "mc": 1.3,
+        "mb": 4.75,
+        "mt": 172.0,
+        "FNS": "VFNS",
+        "ModEv": [
+            "EXA",
+            # "EXP",
+            # "TRN",
+        ],
+        "kcThr": 1.0,
+        "kbThr": 1.0,
+        "ktThr": 1.0,
+        "Q0": 1.25,
+        "alphas": 0.118000,
+        "alphaqed": 0.007496,
+    }
+    vfns_theory = tolist(vfns_theory)
+
+    def benchmark_plain(self, pto, qed):
+        """Plain configuration"""
+
+        th = self.vfns_theory.copy()
+        th.update({"PTO": [pto], "QED": [qed]})
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run(
+            cartesian_product(th),
+            operators.build(operators.apfel_config),
+            ["NNPDF31_nnlo_as_0118_luxqed"],
+        )
+
+
+if __name__ == "__main__":
     # obj = BenchmarkVFNS()
-    obj = BenchmarkFFNS()
+    # obj = BenchmarkFFNS()
 
-    obj.benchmark_plain_pol(2)
+    # obj.benchmark_plain_pol(2)
     # obj.benchmark_plain(2)
     # obj.benchmark_sv(2, "exponentiated")
     # obj.benchmark_kthr(2)
     # obj.benchmark_msbar(2)
+
+    # obj = BenchmarkFFNS_qed()
+    obj = BenchmarkFFNS_qed()
+    obj.benchmark_plain(2, 2)

benchmarks/lha_paper_bench.py

@@ -1,6 +1,8 @@
 """
 Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO).
 """
+from math import nan
+
 import numpy as np
 from banana import register
 from banana.data import cartesian_product
@@ -27,6 +29,7 @@
     "alphas": 0.35,  # Eq. (4.55) :cite:`Dittmar:2005ed`
     "alphaqed": 0.007496,
     "QED": 0,
+    "QrefQED": nan,
 }
 """Global theory settings"""
 

benchmarks/pegasus_bench.py

@@ -1,6 +1,8 @@
 """
 Benchmark to Pegasus :cite:`Vogt:2004ns`
 """
+from math import nan
+
 import numpy as np
 from banana import register
 from banana.data import cartesian_product
@@ -54,6 +56,7 @@ class BenchmarkVFNS(PegasusBenchmark):
         "kbThr": 1.0,
         "ktThr": 1.0,
         "Qref": np.sqrt(2.0),
+        "QrefQED": nan,
         "alphas": 0.35,
         "alphaqed": 0.007496,
         "QED": 0,
@@ -107,7 +110,9 @@ class BenchmarkFFNS(PegasusBenchmark):
         "kbThr": np.inf,
         "ktThr": np.inf,
         "Qref": np.sqrt(2.0),
+        "QrefQED": nan,
         "alphas": 0.35,
+        "alphaqed": 0.007496,
         "Q0": np.sqrt(2.0),
     }
     ffns_theory = tolist(ffns_theory)
@@ -151,7 +156,6 @@ def benchmark_sv(self, pto, svmode):
 
 
 if __name__ == "__main__":
-
     # obj = BenchmarkVFNS()
     obj = BenchmarkFFNS()
     obj.benchmark_plain_pol(1)

benchmarks/sandbox.py

@@ -1,4 +1,6 @@
 # pylint: skip-file
+from math import nan
+
 import numpy as np
 from banana import register
 
@@ -25,6 +27,7 @@
 }
 nnpdf_base_theory = {
     "Qref": 91.2,
+    "QrefQED": nan,
     "mc": 1.51,
     "mb": 4.92,
     "mt": 172.5,

doc/source/code/IO.rst

@@ -69,7 +69,7 @@ and environment. The benchmark settings are available at :mod:`banana.data.theor
   * - ``alphaqed``
     - :py:obj:`float`
     - Reference value of the electromagnetic coupling :math:`\alpha_{em}`.
-  * - ``Qedref``
+  * - ``QrefQED``
     - :py:obj:`float`
     - Reference scale at which the ``alphaqed`` value is given (in GeV).
   * - ``HQ``

doc/source/theory/DGLAP.rst

@@ -172,7 +172,7 @@ Here the strategies are:
 - for ``method in ['iterate-exact', 'iterate-expanded']`` we use a discretized path-ordering :cite:`Bonvini:2012sh`:
 
 .. math::
-    \ESk{1}{a_s}{a_s^0} = \prod\limits_{k=n}^{0} \ESk{1}{a_s^{k+1}}{a_s^{k}}\quad \text{with} a_s^{n+1} = a_s
+    \ESk{1}{a_s}{a_s^0} = \prod\limits_{k=n}^{0} \ESk{1}{a_s^{k+1}}{a_s^{k}}\quad \text{with}\quad a_s^{n+1} = a_s
 
 where the order of the product is such that later |EKO| are to the left and
 
@@ -522,3 +522,36 @@ evolution is simply an identity operation: e.g. for an intrinsic charm distribut
 After :doc:`crossing the mass threshold </theory/Matching>` (charm in this example) the |PDF| can not be considered intrinsic
 any longer and hence, they have to be rejoined with their evolution basis elements and take then again
 part in the ordinary collinear evolution.
+
+Mixed |QCD| :math:`\otimes` |QED| evolution
+-----------------------------------------
+
+For the moment in this case only the `exact` evolution is implemented.
+
+Singlet
+^^^^^^^
+
+The evolution is obtained in the same way of the pure |QCD| case, with the only difference that
+now both :math:`\gamma` and :math:`\beta_{qcd}` contain the |QED| corrections.
+
+In the case in which :math:`\alpha_{em}` is running, at every step of the iteration the corresponding value
+of :math:`a_{em}(a_s)` is used.
+
+Non singlet
+^^^^^^^^^^^
+
+For the non singlet, being it diagonal, the solution is straightforward.
+When :math:`\alpha_{em}` is fixed, the terms proportional to it are just a constant in the splitting functions, and therefore
+they can be integrated directly. For example at ``order=(1,1)`` we have
+
+.. math::
+    \tilde E^{(1,1)}_{ns}(a_s \leftarrow a_s^0) &= \exp \Bigl( -\int_{\log \mu_0^2}^{\log \mu^2}d\log\mu^2 \gamma_{ns}^{(1,0)} a_s(\log\mu^2) + \gamma_{ns}^{(1,1)} a_s(\log\mu^2) a_{em} + \gamma_{ns}^{(0,1)} a_em \Bigr) \\
+    & = \exp \Bigl( \int_{a_s^0}^{a_s}da_s\frac{\gamma_{ns}^{(1,0)} a_s + \gamma_{ns}^{(1,1)} a_s a_{em} + \gamma_{ns}^{(0,1)} a_em}{a_s^2(\beta_0 + \beta_0^{mix} a_{em})}  -\int_{\log \mu_0^2}^{\log \mu^2}d\log\mu^2 \gamma_{ns}^{(0,1)} a_em\Bigr)
+
+In the last expression, the first term can be integrated with the :math:`j^{(n,m)` functions, while the second term is trivial.
+
+In the case of :math:`\alpha_{em}` running, the :math:`a_s` integration integral is divided in steps, such that in every step
+:math:`\alpha_{em}` is considered constant. In this way the solution will be the product of the solutions of every integration step:
+
+.. math::
+    \tilde E^{(1,1)}_{ns}(a_s \leftarrow a_s^0) = \prod\limits_{k=n}^{0} E^{(1,1)}_{ns}(a_s^{k+1} \leftarrow a_s^k, a_{em}^k)