Fix scale variations scheme B

benchmarks/lha_paper_bench.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO)
+Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO).
 """
 import numpy as np
 from banana import register
@@ -35,9 +35,7 @@
 
 
 class LHABenchmark(Runner):
-    """
-    Globally set the external program to LHA
-    """
+    """Globally set the external program to LHA."""
 
     def __init__(self):
         super().__init__()
@@ -46,36 +44,34 @@ def __init__(self):
         self.rotate_to_evolution_basis = True
 
     def plain_theory(self, pto):
-        """
-        Plain theories at given PTO.
+        """Generate plain theories at given PTO.
 
         Parameters
         ----------
-            pto : int
-                perturbation order
+        pto : int
+            perturbation order
 
         Returns
         -------
-            list(dict)
-                theory updates
+        list(dict)
+            theory updates
         """
         th = self.theory.copy()
         th.update({"PTO": pto})
         return [th]
 
     def sv_theories(self, pto):
-        """
-        Scale variation theories.
+        """Generate scale variation theories.
 
         Parameters
         ----------
-            pto : int
-                perturbation order
+        pto : int
+            perturbation order
 
         Returns
         -------
-            list(dict)
-                theory updates
+        list(dict)
+            theory updates
         """
         low = self.theory.copy()
         low["PTO"] = pto
@@ -89,29 +85,27 @@ def sv_theories(self, pto):
 
     @staticmethod
     def skip_pdfs(_theory):
-        """
-        Adjust skip_pdf by the used theory
+        """Adjust skip_pdf by the used theory.
 
         Parameters
         ----------
-            theory_updates : list(dict)
-                theory updates
+        theory : dict
+            theory update
 
         Returns
         -------
-            list :
-                current skip_pdf
+        list
+            allowed PDFs in LHA
         """
         return default_skip_pdfs
 
     def run_lha(self, theory_updates):
-        """
-        Enforce operators and PDF
+        """Enforce operator grid and PDF.
 
         Parameters
         ----------
-            theory_updates : list(dict)
-                theory updates
+        theory_updates : list(dict)
+            theory updates
         """
         self.run(
             theory_updates,
@@ -126,16 +120,16 @@ def run_lha(self, theory_updates):
         )
 
     def benchmark_plain(self, pto):
-        """Plain configuration"""
+        """Run plain configuration."""
         self.run_lha(self.plain_theory(pto))
 
     def benchmark_sv(self, pto):
-        """Scale variations"""
+        """Run scale variations."""
         self.run_lha(self.sv_theories(pto))
 
 
 class BenchmarkVFNS(LHABenchmark):
-    """Variable Flavor Number Scheme"""
+    """Provide |VFNS| settings."""
 
     def __init__(self):
         super().__init__()
@@ -153,7 +147,7 @@ def __init__(self):
 
 
 class BenchmarkFFNS(LHABenchmark):
-    """Fixed Flavor Number Scheme"""
+    """Provide |FFNS| settings."""
 
     def __init__(self):
         super().__init__()
@@ -171,38 +165,42 @@ def __init__(self):
 
     @staticmethod
     def skip_pdfs(theory):
+        """Adjust skip_pdf by the used theory.
+
+        Parameters
+        ----------
+        theory : dict
+            theory update
+
+        Returns
+        -------
+        list
+            allowed PDFs in FFNS LHA
+        """
         ffns_skip_pdfs = default_skip_pdfs.copy()
         # remove bottom
         ffns_skip_pdfs.extend([-5, 5, "T24"])
-        # in NNLO V8 becomes available
+        # in NNLO also V8 gets removed
         if theory["PTO"] >= 2:
             ffns_skip_pdfs.remove("V8")
         return ffns_skip_pdfs
 
 
 class BenchmarkRunner(BenchmarkVFNS):
-    """
-    Generic benchmark runner using the LHA VFNS settings
-    """
+    """Generic benchmark runner using the LHA |VFNS| settings."""
 
     def __init__(self, external):
         super().__init__()
         self.external = external
         self.sandbox = True
 
     def benchmark_sv(self, pto):
-        """
-        Scale variations
+        """Run scale variations.
 
         Parameters
         ----------
-            pto : int
-                perturbation order
-
-        Returns
-        -------
-            list(dict)
-                theory updates
+        pto : int
+            perturbation order
         """
         high, low = self.sv_theories(pto)
 
@@ -223,15 +221,15 @@ def benchmark_sv(self, pto):
 if __name__ == "__main__":
 
     # Benchmark to LHA
-    # obj = BenchmarkVFNS()
-    obj = BenchmarkFFNS()
+    obj = BenchmarkVFNS()
+    # obj = BenchmarkFFNS()
 
     # obj.benchmark_plain(0)
-    # obj.benchmark_sv(2)
+    obj.benchmark_sv(2)
 
     # # VFNS benchmarks with LHA settings
-    programs = ["LHA", "pegasus", "apfel"]
-    for p in programs:
-        obj = BenchmarkRunner(p)
-        # obj.benchmark_plain(2)
-        obj.benchmark_sv(2)
+    # programs = ["LHA", "pegasus", "apfel"]
+    # for p in programs:
+    #     obj = BenchmarkRunner(p)
+    #     # obj.benchmark_plain(2)
+    #     obj.benchmark_sv(2)

doc/source/theory/pQCD.rst

@@ -156,11 +156,12 @@ This method provides many advantages:
     * it naturally incorporates renormalization group invariance,
       as the perturbative order increases, estimates of |MHOU| decrease;
     * the same procedure can be used for any perturbative process,
-      since the scale dependence of the strong coupling :math:`a_s(\mu^2)` and of PDFs is universal;
+      since the scale dependence of the strong coupling :math:`a_s(\mu^2)` and of |PDF| are universal;
 
 However, there is no unique prescription to determine the specific range of the scale variation,
-the most used prescription specify to vary the factor :math:`\mu_F/\mu_R` in the range:
-:math:`1/2 \le \mu_F/\mu_R \le 2`.
+the most common prescription specify to vary the factor :math:`\mu_F/\mu_R` in the range:
+:math:`1/2 \le \mu_F/\mu_R \le 2`. In the following we express this additional dependency as a function
+of :math:`k = \ln(\mu_F^2/\mu_R^2)`
 
 This variation can be performed at least at two different levels during the |PDF|
 evolution, always evaluating the strong coupling at :math:`\mu_R^2`.
@@ -169,43 +170,40 @@ evolution, always evaluating the strong coupling at :math:`\mu_R^2`.
       and the anomalous dimension are then modified using :cite:`Vogt:2004ns`:
 
         .. math ::
-            & \gamma^{(1)}(N) \to \gamma^{(1)}(N) - \beta_0 k \gamma^{(0)} \\
-            & \gamma^{(2)}(N) \to \gamma^{(2)}(N) - 2 \beta_0 k \gamma^{(1)} - ( \beta_1 k - \beta_0^2 k^2) \gamma^{(0)} \\
-            & \gamma^{(3)}(N) \to \gamma^{(3)}(N) - 3 \beta_0 k \gamma^{(2)} - ( 2 \beta_1 k - 3 \beta_0^2 k^2) \gamma^{(1)} - (\beta_2 k - \frac{5}{2} \beta_1 \beta_0 k^2 + \beta_0^3 k^3) \gamma^{(0)} \\
-            & k = \ln(\mu_F^2/\mu_R^2)
+            \gamma^{(1)}(N) &\to \gamma^{(1)}(N) - \beta_0 k \gamma^{(0)} \\
+            \gamma^{(2)}(N) &\to \gamma^{(2)}(N) - 2 \beta_0 k \gamma^{(1)} - ( \beta_1 k - \beta_0^2 k^2) \gamma^{(0)} \\
+            \gamma^{(3)}(N) &\to \gamma^{(3)}(N) - 3 \beta_0 k \gamma^{(2)} - ( 2 \beta_1 k - 3 \beta_0^2 k^2) \gamma^{(1)} - (\beta_2 k - \frac{5}{2} \beta_1 \beta_0 k^2 + \beta_0^3 k^3) \gamma^{(0)}
 
       This procedure corresponds to Eq. (3.32) of :cite:`AbdulKhalek:2019ihb`, and we recommend to use it along with
       ``ModEv='iterate-exact'`` in order to be in agreement with the treatment of the evolution integral expansion.
 
 
-    * In ``ModSV='expanded'`` the |EKO| is multiplied by an additional kernel, such that
-      the scale variation is applied to the whole evolution operator:
+    * In ``ModSV='expanded'`` the full |EKO| is multiplied by an additional kernel:
 
         .. math ::
             \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) & = \tilde{\mathbf{K}}(a_s) \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) \\
             \tilde{\mathbf{K}}(a_s) & = 1 - k \gamma + \frac{1}{2} k^2 \left ( \gamma^{2} - \beta \frac{\partial \gamma}{\partial a_s} \right ) \\
-            & + \frac{1}{6} k^3 \left [ - \beta \frac{\partial}{\partial a_s} \left( \beta \frac{\partial \gamma}{\partial a_s} \right) + 3 \beta \frac{\partial \gamma}{\partial a_s} \gamma - \gamma^3 \right ] + \mathcal{O}(k^4)
+            & \hspace{10pt} + \frac{1}{6} k^3 \left [ - \beta \frac{\partial}{\partial a_s} \left( \beta \frac{\partial \gamma}{\partial a_s} \right) + 3 \beta \frac{\partial \gamma}{\partial a_s} \gamma - \gamma^3 \right ] + \mathcal{O}(k^4)
 
-      where scale variation kernel is expanded consistently order by order in :math:`a_s`,
+      where the scale variation kernel :math:`\tilde{\mathbf{K}}` is expanded consistently order by order in :math:`a_s`,
       leading to:
 
         .. math ::
-            \tilde{\mathbf{K}}(a_s) \approx & 1 - a_s k \gamma^{(0)} + a_s^2 \left [ - k \gamma^{(1)} + \frac{1}{2} k^2 \gamma^{(0)} (\beta_0 + \gamma^{(0)}) \right ] \\
-            & + a_s^3 \left [ -k \gamma^{(2)} + \frac{1}{2} k^2 \left(\beta_1 \gamma^{(0)} + 2 \gamma^{(1)} (\beta_0 + \gamma^{(0)} ) \right) \right. \\
-            & \left. - \frac{1}{6} k^3 \gamma^{(0)} \left(2 \beta_0^2 + 3 \beta_0 \gamma^{(0)}+\gamma^{(0),2} \right) \right] + \mathcal{O}(a^4)
+            \tilde{\mathbf{K}}(a_s) &\approx 1 - a_s k \gamma^{(0)} + a_s^2 \left [ - k \gamma^{(1)} + \frac{1}{2} k^2 \gamma^{(0)} (\beta_0 + \gamma^{(0)}) \right ] \\
+            & \hspace{10pt} + a_s^3 \left [ -k \gamma^{(2)} + \frac{1}{2} k^2 \left(\beta_1 \gamma^{(0)} + 2 \beta_0\gamma^{(1)}  + \gamma^{(1)}\gamma^{(0)} + \gamma^{(0)}\gamma^{(1)} \right) \right. \\
+            & \hspace{35pt} \left. - \frac{1}{6} k^3 \gamma^{(0)} \left(2 \beta_0^2 + 3 \beta_0 \gamma^{(0)} + \left(\gamma^{(0)}\right)^2 \right) \right] + \mathcal{O}(a_s^4)
 
 
       In this way the dependence of the |EKO| on :math:`k` is factorized outside the unvaried evolution kernel.
       This procedure is repeated for each different flavor patch present in the evolution path.
       It corresponds to Eq. (3.35) of :cite:`AbdulKhalek:2019ihb`, and we recommend to use it along with
       ``ModEv='truncated'`` in order to keep consistency with the evolution integral expansion.
 
-
-By construction, the corrections of the order :math:`\mathcal{O}(k^n)` will appear
-at the order :math:`n` in the expansion :math:`a_s`.
-This happens because :math:`\beta \approx \mathcal{O}(a_s^2)`, :math:`\gamma \approx \mathcal{O}(a_s)`
-and the contribution proportional to :math:`\mathcal{O}(k^n)` is originated
-by the `n-th` derivative in :math:`\gamma` :cite:`AbdulKhalek:2019ihb`.
+      By construction, the corrections of the order :math:`\mathcal{O}(k^n)` will appear
+      at the order :math:`n` in the expansion :math:`a_s`.
+      This happens because :math:`\beta \approx \mathcal{O}(a_s^2)`, :math:`\gamma \approx \mathcal{O}(a_s)`
+      and the contribution proportional to :math:`\mathcal{O}(k^n)` is originated
+      by the `n-th` derivative in :math:`\gamma` :cite:`AbdulKhalek:2019ihb`.
 
 Furthermore the distance between the varied |EKO| and the unvaried one will decrease while
 keeping higher order terms in :math:`a_s`

src/eko/anomalous_dimensions/as4/ggg.py

@@ -137,17 +137,17 @@ def gamma_gg_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf2(n, sx):
+def gamma_gg_nf2(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf1(n, sx):
+def gamma_gg_nf1(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf0(n, sx):
+def gamma_gg_nf0(_n, _sx):
     return 0
 
 

src/eko/anomalous_dimensions/as4/gps.py

@@ -79,12 +79,12 @@ def gamma_ps_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_ps_nf2(n, sx):
+def gamma_ps_nf2(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_ps_nf1(n, sx):
+def gamma_ps_nf1(_n, _sx):
     return 0
 
 

src/eko/anomalous_dimensions/as4/gqg.py

@@ -329,17 +329,17 @@ def gamma_qg_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf2(n, sx):
+def gamma_qg_nf2(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf1(n, sx):
+def gamma_qg_nf1(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf0(n, sx):
+def gamma_qg_nf0(_n, _sx):
     return 0