Decorrelate N3LO IHOU

src/yadism/coefficient_functions/asy/kernels.py

@@ -108,7 +108,17 @@ def generate_heavy_asy(esf, nf, pto_evol, ihq):
         )
         if not is_pv:
             n3lo_cf_variation = esf.info.theory["n3lo_cf_variation"]
+            if "F2" in kind:
+                n3lo_g_variation = n3lo_cf_variation[0]
+                n3lo_q_variation = n3lo_cf_variation[1]
+            elif "FL" in kind:
+                n3lo_g_variation = n3lo_cf_variation[2]
+                n3lo_q_variation = n3lo_cf_variation[3]
+            else:
+                n3lo_g_variation = 0
+                n3lo_q_variation = 0
             for c, channel in (("g", "Gluon"), ("s", "Singlet")):
+                variation = n3lo_g_variation if c == "g" else n3lo_q_variation
                 for res in range(pto_evol + 1):
                     name = "Asy" + ("N" * res) + "LL" + channel
                     for av in ("AA", "VV"):
@@ -119,7 +129,7 @@ def generate_heavy_asy(esf, nf, pto_evol, ihq):
                                     esf,
                                     nf,
                                     m2hq=m2hq,
-                                    n3lo_cf_variation=n3lo_cf_variation,
+                                    n3lo_cf_variation=variation,
                                 ),
                             )
                         )

src/yadism/coefficient_functions/heavy/kernels.py

@@ -70,15 +70,23 @@ def generate(esf, nf, ihq):
             ihq,
             is_pv,
         )
-        n3lo_cf_variation = esf.info.theory["n3lo_cf_variation"]
+        if "F2" in kind:
+            n3lo_g_variation = esf.info.theory["n3lo_cf_variation"][0]
+            n3lo_q_variation = esf.info.theory["n3lo_cf_variation"][1]
+        elif "FL" in kind:
+            n3lo_g_variation = esf.info.theory["n3lo_cf_variation"][2]
+            n3lo_q_variation = esf.info.theory["n3lo_cf_variation"][3]
+        else:
+            n3lo_g_variation = 0
+            n3lo_q_variation = 0
         gVV = kernels.Kernel(
             weights["gVV"],
-            pcs.GluonVV(esf, nf, m2hq=m2hq, n3lo_cf_variation=n3lo_cf_variation),
+            pcs.GluonVV(esf, nf, m2hq=m2hq, n3lo_cf_variation=n3lo_g_variation),
         )
         gAA = kernels.Kernel(weights["gAA"], pcs.GluonAA(esf, nf, m2hq=m2hq))
         sVV = kernels.Kernel(
             weights["sVV"],
-            pcs.SingletVV(esf, nf, m2hq=m2hq, n3lo_cf_variation=n3lo_cf_variation),
+            pcs.SingletVV(esf, nf, m2hq=m2hq, n3lo_cf_variation=n3lo_q_variation),
         )
         sAA = kernels.Kernel(weights["sAA"], pcs.SingletAA(esf, nf, m2hq=m2hq))
         return (gVV, gAA, sVV, sAA)

src/yadism/runner.py

@@ -121,6 +121,12 @@ def __init__(self, theory: dict, observables: dict):
             coupling_constants=coupling_constants,
             sv_manager=sv_manager,
         )
+
+        # The n3lo_cf_variations are ordered as (C2g, C2q, CLg, CLq)
+        # If only one number is passed, all are set to that
+        if isinstance(theory["n3lo_cf_variation"], int):
+            theory["n3lo_cf_variation"] = (theory["n3lo_cf_variation"], theory["n3lo_cf_variation"], theory["n3lo_cf_variation"], theory["n3lo_cf_variation"])
+
         # pass theory params
         theory_params = dict(
             pto=pto,
@@ -215,7 +221,6 @@ def replace_nans_with_0(self, out):
         for observable, points in out2.items():
             # Skip the keys that are not an observable
             if not observable_name.ObservableName.is_valid(observable):
-                logger.critical("Some NaNs are encountered and set to zero!")
                 continue
 
             # Loop over the kinematic points
@@ -229,6 +234,8 @@ def replace_nans_with_0(self, out):
                     for tup in range(2):
                         # Set any NaN or inf values in the array to 0
                         values[tup][~np.isfinite(values[tup])] = 0.0
+                        logger.critical(f"NaNs are encountered and set to zero at [x,Q2] = [{point.x}, {point.Q2}].!")
+
         return out2
 
     def get_result(self):

src/yadmark/data/observables.py

@@ -15,6 +15,7 @@
     TargetDIS="proton",
     NCPositivityCharge=None,
     observables={},
+    n3lo_cf_variation=(0, 0, 0, 0),
 )
 default_card = dict(sorted(default_card.items()))
 

tests/yadism/cf/test_nc_kernels.py

@@ -36,7 +36,7 @@ def __init__(self, n):
         self.coupling_constants = MockCouplingConstants()
         self.m2hq = [1.0, 2.0, 3.0]
         self.threshold = Atlas(matching_scales=self.m2hq, origin=(1.65**2, 4))
-        self.theory = {"n3lo_cf_variation": 0, "pto": 3}
+        self.theory = {"n3lo_cf_variation": [0, 0, 0, 0], "pto": 3}
 
 
 class MockESF:

tests/yadism/test_init.py

@@ -35,7 +35,7 @@
     "GF": 1.1663787e-05,
     "SIN2TW": 0.23126,
     "ModEv": "EXA",
-    "n3lo_cf_variation": 0,
+    "n3lo_cf_variation": (0, 0, 0, 0),
 }
 
 obs_dict = {

tests/yadism/test_runner.py

@@ -83,7 +83,7 @@ def fake_update(*args):
             MZ=100.0,
             GF=1.0,
             FONLLParts="full",
-            n3lo_cf_variation=0,
+            n3lo_cf_variation=(0, 0, 0, 0),
         )
         xgrid = np.geomspace(1e-5, 1.0, 6)
         observables = dict(