Add N3LO ren SV

src/pineko/scale_variations.py

@@ -51,16 +51,32 @@ def ren_sv_coeffs(m, max_as, logpart, which_part, nf):
     float
         renormalization scale variation contribution
     """
-    bcoeff = beta.beta_qcd((max_as - logpart - which_part + 2, 0), nf)
-    as_normalization = AS_NORM ** (max_as - which_part)
+    # nothing to do?
     if max_as == 0:
         return 0.0
-    if max_as == 2:
-        if which_part > 0:
-            m += 1
-        elif logpart > 1:
-            m = 0.5 * m * (m + 1)
-    return m * as_normalization * bcoeff
+    # eko uses as = alpha_s/(4pi), but pineappl just alpha_s
+    as_normalization = AS_NORM ** (max_as - which_part)
+    # the coefficients can be found in the NNLO MHOU paper
+    # (which also contains a generating MMa script)
+    beta0 = beta.beta_qcd_as2(nf)
+    beta1 = beta.beta_qcd_as3(nf)
+    beta2 = beta.beta_qcd_as4(nf)
+    ren_coeffs = {
+        # NLO
+        (1, 1, 0): m * beta0,
+        # NNLO
+        (2, 1, 1): (m + 1) * beta0,
+        (2, 1, 0): (m + 0) * beta1,
+        (2, 2, 0): m * (m + 1) / 2.0 * beta0**2,
+        # N3LO
+        (3, 1, 2): (m + 2) * beta0,
+        (3, 1, 1): (m + 1) * beta1,
+        (3, 1, 0): (m + 0) * beta2,
+        (3, 2, 1): (m + 1) * (m + 2) / 2.0 * beta0**2,
+        (3, 2, 0): m * (2 * m + 3) / 2.0 * beta0 * beta1,
+        (3, 3, 0): m * (m + 1) * (m + 2) / 6 * beta0**3,
+    }
+    return as_normalization * ren_coeffs[(max_as, logpart, which_part)]
 
 
 def requirements(m: int, max_as: int, al: int) -> Dict[OrderTuple, List[OrderTuple]]:

tests/test_scale_variations.py

@@ -29,6 +29,42 @@ def test_ren_sv_coeffs():
         scale_variations.ren_sv_coeffs(m=0, max_as=2, logpart=2, which_part=0, nf=5),
         0.0,
     )
+    # check all coeff are set
+    for m in (0, 1, 10):
+        for max_as in range(1, 3 + 1):
+            for logpart in range(1, max_as + 1):
+                for which_part in range(0, max_as - logpart + 1):
+                    for nf in (4, 5):
+                        c = scale_variations.ren_sv_coeffs(
+                            m=m,
+                            max_as=max_as,
+                            logpart=logpart,
+                            which_part=which_part,
+                            nf=nf,
+                        )
+                        assert np.isfinite(c)
+                        if which_part > 0 and not m == 0:
+                            assert c > 0.0
+                        else:
+                            assert c >= 0.0
+    # due to the exponential structure we can predict some things:
+    for nf in (4, 5):
+        # the highest log is always proportional beta0^k
+        for k in range(1, 3 + 1):
+            c = scale_variations.ren_sv_coeffs(
+                m=1, max_as=k, logpart=k, which_part=0, nf=nf
+            )
+            bare_c = c * (4.0 * np.pi / beta_qcd((2, 0), nf)) ** k
+            int_c = bare_c * np.math.factorial(k)
+            np.testing.assert_allclose(int_c, int(int_c))
+        # and even the second highest for the highest coeff
+        for k in range(2, 3 + 1):
+            c = scale_variations.ren_sv_coeffs(
+                m=1, max_as=k, logpart=k - 1, which_part=1, nf=nf
+            )
+            bare_c = c * (4.0 * np.pi / beta_qcd((2, 0), nf)) ** (k - 1)
+            int_c = bare_c * np.math.factorial(k)
+            np.testing.assert_allclose(int_c, int(int_c))
 
 
 def test_requirements():