Implement QED exact solution

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(
@@ -81,4 +104,4 @@ def benchmark_znnlo(self, Q0=1.3, Q2grid=(1e4,)):
 
 if __name__ == "__main__":
     b = BenchmarkCT18()
-    b.benchmark_nnlo()
+    b.benchmark_nnlo_qed()

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"""
 
@@ -131,5 +158,7 @@ def benchmark(self, Q0=1.65, 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)

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)

doc/source/theory/Matching.rst

@@ -135,3 +135,34 @@ EKO implements two different strategies to perform this operation, that can be s
 We emphasize that in the backward evolution, below the threshold, the remaining high quark PDFs are always intrinsic and do not evolve anymore.
 In fact, if the initial PDFs (above threshold) do contain an intrinsic contribution, this has to be evolved below the threshold otherwise momentum sum rules
 can be violated.
+
+QED Matching
+------------
+
+In the QED case the matching is changed only because of the change of the evolution basis, therefore the only different part will be the basis rotation.
+In fact, the |OME| :math:`\mathbf{A}^{(n_f)}(\mu_{h}^2)` don't have |QED| corrections. The matching of the singlet sector is unchanged since it
+remains the same with respect to the |QCD| case. The same happens for the matching of the valence component. All the elements :math:`V_i` and :math:`T_i`
+are non-singlet components, therefore they are matched with :math:`A_{ns}`. In the end, the new components :math:`\Sigma_\Delta` and :math:`V_\Delta` are matched
+with :math:`A_{ns}` since they are both non-singlets.
+
+QED basis rotation
+------------------
+
+For the basis rotation we have to consider that we are using the intrinsic unified evolution basis. Here it will be discussed only the rotation to be applied
+to the sector :math:`(\Sigma, \Sigma_\Delta, T_i)`, being the rotation of the sector :math:`(V, V_\Delta, V_i)` completely equivalent.
+The rotation matrix is given by:
+
+.. math ::
+    \begin{pmatrix} \Sigma_{(n_f)} \\ \Sigma_{\Delta,(n_f)} \\ T_{i,(nf)} \end{pmatrix}^{(n_f+1)} =
+    \begin{pmatrix} 1 & 0 & 1 \\ a(n_f) & b(n_f) & c(n_f) \\ d(n_f) & e(n_f) & f(n_f) \end{pmatrix}
+    \begin{pmatrix} \Sigma_{(n_f)} \\ \Sigma_{\Delta,(n_f)} \\ h^+ \end{pmatrix}^{(n_f)}
+
+where
+
+.. math ::
+    a(n_f) & = \frac{1}{n_f}\Bigl(\frac{n_d(n_f+1)}{n_u(n_f+1)}n_u(n_f)-n_d(n_f)\Bigr) \\
+    b(n_f) & = \frac{n_f+1}{n_u(n_f+1)}\frac{n_u(n_f)}{n_f} \\
+    c(n_f) & = \begin{cases} \frac{n_d(n_f+1)}{n_u(n_f+1)} \quad \text{if $h$ is up-like}\\-1  \quad \text{if $h$ is down-like}\end{cases} \\
+    d(n_f) & = \begin{cases} &\frac{n_u(n_f)}{n_f} \quad \text{if $h$ is up-like ($n_f$=3,5)} \\ &\frac{n_d(n_f)}{n_f} \quad \text{if $h$ is down-like ($n_f$=2,4)} \end{cases} \\
+    e(n_f) & = \begin{cases} &\frac{n_u(n_f)}{n_f} \quad \text{if $h$ is up-like} \\ &-\frac{n_u(n_f)}{n_f} \quad \text{if $h$ is down-like} \end{cases} \\
+    f(n_f) & = \begin{cases} &-1\quad \text{if $h$ is $s$, $c$  ($n_f$=2,3)} \\ &-2 \quad \text{if $h$ is $b$, $t$  ($n_f$=4,5)} \end{cases}

doc/source/theory/pQCD.rst

@@ -129,10 +129,11 @@ Order specification
 In the code ``order=tuple(int,int)`` specifies the |QCD| and |QED| perturbative orders of the splitting functions in terms
 of :math:`a_s = \alpha_s/(4\pi)` and :math:`a_{em} = \alpha_{em}/(4\pi)`. The available perturbative expansions are the following:
 
-- ``order=(0,0)``: it is the non evolution case in which :math:`a_s` and :math:`a_{em}` are kept fixed and the splitting functions are null.
 - ``order=(n,0)``: with :math:`n=1,2,3,4` correspond to the pure |QCD| evolution at |LO|, |NLO|, |NNLO| and |N3LO| in which the |QCD| splitting functions are expanded up to :math:`\mathcal{O}(a_s^n)` and the strong coupling is evolved using the n-th coefficient of the beta function, i.e. :math:`\beta_{n-1}`.
 - ``order=(n,m)``; with :math:`n=1,2,3,4` and :math:`m=1,2` corresponds to the mixed |QED| :math:`\otimes` |QCD| evolution in which the splitting functions are expanded up to :math:`\mathcal{O}(a_s^na_{em}^m)`, the stromg coupling is evolved using up to the n-th coefficient of the beta function and the electromagnetic coupling is kept fixed.
 
+Observe that the case :math:`n=0` is not allowed, since it would correspond to the pure |QED| evolution or (if :math:`m > 0`) no evolution at all.
+
 Sum Rules
 ---------