diff --git a/src/yadism/coefficient_functions/coupling_constants.py b/src/yadism/coefficient_functions/coupling_constants.py
index 1db5e083..99c1023d 100644
--- a/src/yadism/coefficient_functions/coupling_constants.py
+++ b/src/yadism/coefficient_functions/coupling_constants.py
@@ -132,9 +132,7 @@ def nc_partonic_coupling(self, pid):
         qZ = {"V": self.vectorial_coupling(pid), "A": self.weak_isospin_3[pid]}
         return qph, qZ
 
-    def partonic_coupling_fl2(
-        self, mode, pid, quark_coupling_type, cc_mask=None, nf=None
-    ):
+    def partonic_coupling(self, mode, pid, quark_coupling_type, cc_mask=None):
         """Compute the coupling of the boson to the parton.
 
         Parameters
@@ -176,7 +174,13 @@ def partonic_coupling_fl11(self, mode, pid, nf, quark_coupling_type):
         """Compute the coupling of the boson to the parton for the flavor class :math:`fl_{11}`.
 
         This is a generalization of :cite:`Larin:1996wd` (Table 2) for NC.
-        See also pag. 27. This coupling is given by :math:`Tr(Q_f) Q_f`.
+        See also pag. 27. This coupling is given by :
+
+        .. :math:
+            W_{q,bb'} = \frac{Tr(Q_b)}{n_f} Q_b'
+
+        the gluon and pure singlet coupling are then build after summing on all the
+        different electroweak channels.
 
         Parameters
         ----------
@@ -184,6 +188,8 @@ def partonic_coupling_fl11(self, mode, pid, nf, quark_coupling_type):
             scattered bosons
         pid : int
             parton identifier
+        nf : int
+            number of active flavors
         quark_coupling_type : str
             flag to distinguish for heavy quarks between vectorial and axial-vectorial
             coupling
@@ -193,30 +199,28 @@ def partonic_coupling_fl11(self, mode, pid, nf, quark_coupling_type):
         float
             partonic coupling
         """
-        # select quark coupling
+        # In CC there is not such flavor class
+        if mode == "WW":
+            return 0
+
         first = quark_coupling_type[0]
         second = quark_coupling_type[1]
 
-        def switch_mode(quark):
+        def switch_mode(quark, coupling_type):
             qph, qZ = self.nc_partonic_coupling(quark)
-            if mode == "phph":
-                return qph[first], qph[second]
-            if mode == "phZ":
-                return qph[first], qZ[second]
-            if mode == "ZZ":
-                return qZ[first], qZ[second]
-            return ValueError(f"Unknown mode {mode}")
-
-        # load linear coupling
-        g1, _ = switch_mode(abs(pid))
-
-        # compute trace
+            if mode in ["phph", "Zph"]:
+                return qph[coupling_type]
+            if mode in ["phZ", "ZZ"]:
+                return qZ[coupling_type]
+            raise ValueError(f"Unknown mode: {mode}")
+
+        # first coupling contribute with a mean
         pids = range(1, nf + 1)
-        trace = 0
-        for quark in pids:
-            _, g2 = switch_mode(quark)
-            trace += g2
-        return g1 * trace
+        g1 = np.mean([switch_mode(quark, first) for quark in pids])
+        # the second coupling is just the linear coupling
+        g2 = switch_mode(abs(pid), second)
+
+        return g1 * g2
 
     def propagator_factor(self, mode, Q2):
         r"""Compute propagator correction to account for different bosons (:math:`\\eta` in PDG).
@@ -254,9 +258,12 @@ def propagator_factor(self, mode, Q2):
             return eta_W
         raise ValueError(f"Unknown mode: {mode}")
 
-    def get_weight(self, pid, Q2, quark_coupling_type, cc_mask=None, nf=None):
+    def get_weight(self, pid, Q2, quark_coupling_type, cc_mask=None):
         """Compute the weight for the pid contributions to the structure function.
 
+        Combine the charges, both on the leptonic side and the hadronic side, as well
+        as propagator changes and/or corrections.
+
         Parameters
         ----------
         pid : int
@@ -297,10 +304,33 @@ def get_weight(self, pid, Q2, quark_coupling_type, cc_mask=None, nf=None):
             if abs(pid) != pos_pid:
                 return 0.0
 
-        partonic_coupling = self.partonic_coupling_fl2
-        return self.get_nc_weight(partonic_coupling, pid, Q2, quark_coupling_type)
+        w_phph = (
+            self.leptonic_coupling("phph", quark_coupling_type)
+            * self.propagator_factor("phph", Q2)
+            * self.partonic_coupling("phph", pid, quark_coupling_type)
+        )
+        # pure photon exchange
+        if self.obs_config["process"] == "EM":
+            return w_phph
+        # allow Z to be mixed in
+        if self.obs_config["process"] == "NC":
+            # photon-Z interference
+            w_phZ = (
+                2
+                * self.leptonic_coupling("phZ", quark_coupling_type)
+                * self.propagator_factor("phZ", Q2)
+                * self.partonic_coupling("phZ", pid, quark_coupling_type)
+            )
+            # true Z contributions
+            w_ZZ = (
+                self.leptonic_coupling("ZZ", quark_coupling_type)
+                * self.propagator_factor("ZZ", Q2)
+                * self.partonic_coupling("ZZ", pid, quark_coupling_type)
+            )
+            return w_phph + w_phZ + w_ZZ
+        raise ValueError(f"Unknown process: {self.obs_config['process']}")
 
-    def get_fl11_weight(self, pid, Q2, quark_coupling_type, nf):
+    def get_fl11_weight(self, pid, Q2, nf, quark_coupling_type):
         """Same as :func:`get_weight`but now for the NC flavor class :math:`fl_{11}`.
 
         Combine the charges, both on the leptonic side and the hadronic side,
@@ -312,11 +342,11 @@ def get_fl11_weight(self, pid, Q2, quark_coupling_type, nf):
             particle identifier
         Q2 : float
             DIS virtuality
+        nf : int
+            number of active flavors
         quark_coupling_type : str
             flag to distinguish for heavy quarks between vectorial and axial-vectorial
             coupling
-        nf : float
-            number of active flavors
 
         Returns
         -------
@@ -327,43 +357,41 @@ def get_fl11_weight(self, pid, Q2, quark_coupling_type, nf):
         if self.obs_config["process"] == "CC":
             return 0.0
 
-        # if we are are we in positivity mode, we fallback to the standard coupling
+        # TODO: what does positivity mode means for this coupling ?
         if (
             self.obs_config["nc_pos_charge"] is not None
             and self.obs_config["nc_pos_charge"] != "all"
         ):
-            return self.get_weight(pid, Q2, quark_coupling_type)
-
-        partonic_coupling = self.partonic_coupling_fl11
-        return self.get_nc_weight(partonic_coupling, pid, Q2, quark_coupling_type, nf)
+            return 0.0
 
-    def get_nc_weight(self, partonic_coupling, pid, Q2, quark_coupling_type, nf=None):
-        """Combine the NC couplings."""
-        # NC or EM
         w_phph = (
             self.leptonic_coupling("phph", quark_coupling_type)
             * self.propagator_factor("phph", Q2)
-            * partonic_coupling("phph", pid, quark_coupling_type, nf=nf)
+            * self.partonic_coupling_fl11("phph", pid, nf, quark_coupling_type)
         )
         # pure photon exchange
         if self.obs_config["process"] == "EM":
             return w_phph
         # allow Z to be mixed in
         if self.obs_config["process"] == "NC":
-            # photon-Z interference
+            # photon-Z interference, this class is not symmetric
             w_phZ = (
-                2
-                * self.leptonic_coupling("phZ", quark_coupling_type)
+                self.leptonic_coupling("phZ", quark_coupling_type)
                 * self.propagator_factor("phZ", Q2)
-                * partonic_coupling("phZ", pid, quark_coupling_type, nf=nf)
+                * self.partonic_coupling_fl11("phZ", pid, nf, quark_coupling_type)
+            )
+            w_Zph = (
+                self.leptonic_coupling("phZ", quark_coupling_type)
+                * self.propagator_factor("phZ", Q2)
+                * self.partonic_coupling_fl11("Zph", pid, nf, quark_coupling_type)
             )
             # true Z contributions
             w_ZZ = (
                 self.leptonic_coupling("ZZ", quark_coupling_type)
                 * self.propagator_factor("ZZ", Q2)
-                * partonic_coupling("ZZ", pid, nf, quark_coupling_type, nf=nf)
+                * self.partonic_coupling_fl11("ZZ", pid, nf, quark_coupling_type)
             )
-            return w_phph + w_phZ + w_ZZ
+            return w_phph + w_phZ + w_Zph + w_ZZ
         raise ValueError(f"Unknown process: {self.obs_config['process']}")
 
     @classmethod
diff --git a/src/yadism/coefficient_functions/light/kernels.py b/src/yadism/coefficient_functions/light/kernels.py
index 8dfafd81..678cce5d 100644
--- a/src/yadism/coefficient_functions/light/kernels.py
+++ b/src/yadism/coefficient_functions/light/kernels.py
@@ -88,49 +88,46 @@ def generate(esf, nf):
     kernels_list = [ns, g, s]
 
     # at N3LO we need to add also the fl11 diagrams
-    if 3 in esf.orders:
-        weights_fl11 = nc_weights(
-            esf.info.coupling_constants, esf.Q2, nf, is_pv=False, is_fl11=True
-        )
+    if esf.info.configs.theory["pto"] == 3:
+
+        gluon_fl11 = pcs.GluonFL11(esf, nf)
+        quark_fl11 = pcs.QuarkFL11(esf, nf)
+        weights_fl11 = nc_fl11_weights(esf.info.coupling_constants, esf.Q2, nf)
         ns_fl11 = kernels.Kernel(
             weights_fl11["ns"],
-            pcs.NonSingletFL11(esf, nf),
-        )
-        g_fl11 = kernels.Kernel(
-            weights_fl11["g"],
-            pcs.GluonFL11(esf, nf),
+            quark_fl11,
         )
+        g_fl11 = kernels.Kernel(weights_fl11["g"], gluon_fl11)
         s_fl11 = kernels.Kernel(
             weights_fl11["s"],
-            pcs.SingletFL11(esf, nf),
+            quark_fl11,
         )
         kernels_list.extend([ns_fl11, g_fl11, s_fl11])
     return kernels_list
 
 
-def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False, is_fl11=False):
+def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False):
     """
     Compute light NC weights.
 
     Parameters
     ----------
-        coupling_constants : CouplingConstants
-            manager for coupling constants
-        Q2 : float
-            W virtuality
-        nf : int
-            number of light flavors
-        is_pv: bool
-            True if observable violates parity conservation
-        skip_heavylight : bool
-            prevent the last quark to couple to the boson
-        is_fl11 : bool
-            if True compute the coupling of the flavor class :math:`fl_{11}`
+    coupling_constants : CouplingConstants
+        manager for coupling constants
+    Q2 : float
+        W virtuality
+    nf : int
+        number of light flavors
+    is_pv: bool
+        True if observable violates parity conservation
+    skip_heavylight : bool
+        prevent the last quark to couple to the boson
+
 
     Returns
     -------
-        weights : dict
-            mapping pid -> weight for ns, g and s channel
+    weights : dict
+        mapping pid -> weight for ns, g and s channel
     """
     # quark couplings
     tot_ch_sq = 0
@@ -147,12 +144,9 @@ def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False, is_fl11
                 q, Q2, "VA"
             ) + coupling_constants.get_weight(q, Q2, "AV")
         else:
-            coupling = coupling_constants.get_weight
-
-            # TODO: do AA contribute to fl11 ??
-            if is_fl11:
-                coupling = coupling_constants.get_fl11_weight
-            w = coupling(q, Q2, "VV", nf=nf) + coupling(q, Q2, "AA", nf=nf)
+            w = coupling_constants.get_weight(
+                q, Q2, "VV", nf=nf
+            ) + coupling_constants.get_weight(q, Q2, "AA", nf=nf)
         ns_partons[q] = w
         ns_partons[-q] = w if not is_pv else -w
         tot_ch_sq += w
@@ -166,3 +160,52 @@ def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False, is_fl11
     # gluon and singlet coupling = charge average (omitting the *2/2)
     s_partons = {q: ch_av for q in [*pids, *(-q for q in pids)]}
     return {"ns": ns_partons, "g": {21: ch_av}, "s": s_partons}
+
+
+def nc_fl11_weights(coupling_constants, Q2, nf, skip_heavylight=False):
+    """Compute the NC weights for the flavor class :math:`fl_{11}`.
+
+    For the time being we don't have such diagrams for parity violating
+    structure functions.
+
+    Parameters
+    ----------
+    coupling_constants : CouplingConstants
+        manager for coupling constants
+    Q2 : float
+        W virtuality
+    nf : int
+        number of light flavors
+    skip_heavylight : bool
+        prevent the last quark to couple to the boson
+
+    Returns
+    -------
+    weights : dict
+        mapping pid -> weight for ns, g and s channel
+    """
+    # quark couplings
+    ns_partons = {}
+    tot_ch_sq = 0
+    pids = range(1, nf + 1)
+
+    for q in pids:
+        # we don't want this quark to couple to the photon (because it is treated separately),
+        # but still let it take part in the average
+        if skip_heavylight and q == nf:
+            continue
+        w = coupling_constants.get_fl11_weight(
+            q, Q2, nf, "VV"
+        ) + coupling_constants.get_fl11_weight(q, Q2, nf, "AA")
+        ns_partons[q] = w
+        ns_partons[-q] = w
+        tot_ch_sq += w
+
+    # compute gluon, singlet (S) (omitting the *2/2)
+    # NOTE: what we call singlet (S) is actually the pure singlet (PS) coefficient.
+    # Here the contribution of the singlet PDF to the NS coefficient
+    # and the contribution of the singlet PDF to the S coefficient
+    # do not cancel completely.
+    ch_av = tot_ch_sq / len(pids)
+    s_partons = {q: ch_av - ns_partons[q] for q in [*pids, *(-q for q in pids)]}
+    return {"ns": ns_partons, "g": {21: ch_av}, "s": s_partons}