diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 0205e3ff8..361cef5eb 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -29,7 +29,7 @@
 from ..kernels import valence_qed as qed_v
 from ..matchings import Segment
 from ..member import OpMember
-from .quad_ker import QuadCargs, c_quad_ker_qcd, cb_quad_ker_qcd, cb_quad_ker_qcd_T
+from .quad_ker import QuadQCDargs, c_quad_ker_qcd, cb_quad_ker_qcd, cb_quad_ker_qcd_T
 
 logger = logging.getLogger(__name__)
 
@@ -306,50 +306,6 @@ def labels(self):
                 labels.extend(br.singlet_unified_labels)
         return labels
 
-    def quad_ker(self, label, logx, areas):
-        """Return partially initialized integrand function.
-
-        Parameters
-        ----------
-        label: tuple
-            operator element pids
-        logx: float
-            Mellin inversion point
-        areas : tuple
-            basis function configuration
-
-        Returns
-        -------
-        functools.partial
-            partially initialized integration kernel
-
-        """
-        return functools.partial(
-            quad_ker,
-            order=self.order,
-            mode0=label[0],
-            mode1=label[1],
-            method=self.config["method"],
-            is_log=self.int_disp.log,
-            logx=logx,
-            areas=areas,
-            as_list=self.as_list,
-            mu2_from=self.q2_from,
-            mu2_to=self.q2_to,
-            a_half=self.a_half_list,
-            alphaem_running=self.alphaem_running,
-            nf=self.nf,
-            L=np.log(self.xif2),
-            ev_op_iterations=self.config["ev_op_iterations"],
-            ev_op_max_order=tuple(self.config["ev_op_max_order"]),
-            sv_mode=self.sv_mode,
-            is_threshold=self.is_threshold,
-            n3lo_ad_variation=self.config["n3lo_ad_variation"],
-            is_polarized=self.config["polarized"],
-            is_time_like=self.config["time_like"],
-            quad_ker_qcd_ns_address=cb_quad_ker_qcd.address,
-        )
-
     def initialize_op_members(self):
         """Init all operators with the identity or zeros."""
         eye = OpMember(
@@ -372,10 +328,7 @@ def initialize_op_members(self):
             else:
                 self.op_members[n] = zero.copy()
 
-    def run_op_integration(
-        self,
-        log_grid,
-    ):
+    def run_op_integration(self, log_grid):
         """Run the integration for each grid point.
 
         Parameters
@@ -390,10 +343,11 @@ def run_op_integration(
         """
         column = []
         k, logx = log_grid
+        # call(!) self.labels only once
         labels = self.labels
         start_time = time.perf_counter()
-        # iterate basis functions
-        cfg = QuadCargs(
+        # start preparing C arguments
+        cfg = QuadQCDargs(
             order_qcd=self.order[0],
             is_polarized=self.config["polarized"],
             is_time_like=self.config["time_like"],
@@ -410,13 +364,16 @@ def run_op_integration(
             sv_mode_num=1,
             is_threshold=self.is_threshold,
         )
+        # iterate basis functions
         for l, bf in enumerate(self.int_disp):
             if k == l and l == self.grid_size - 1:
                 continue
+            # add emtpy labels with 0s
             if bf.is_below_x(np.exp(logx)):
                 column.append({label: (0.0, 0.0) for label in labels})
                 continue
             temp_dict = {}
+            # prepare areas for C
             curareas = bf.areas_representation
             cfg.areas = (ctypes.c_double * (curareas.shape[0] * curareas.shape[1]))(
                 *curareas.flatten().tolist()
@@ -427,6 +384,7 @@ def run_op_integration(
             for label in labels:
                 cfg.mode0 = label[0]
                 cfg.mode1 = label[1]
+                # construct the low level object
                 user_data = ctypes.cast(ctypes.pointer(cfg), ctypes.c_void_p)
                 func = LowLevelCallable(c_quad_ker_qcd, user_data)
                 res = integrate.quad(
diff --git a/src/eko/evolution_operator/c_quad_ker.cpp b/src/eko/evolution_operator/c_quad_ker.cpp
index fc3dc5aa1..de4cb47e7 100644
--- a/src/eko/evolution_operator/c_quad_ker.cpp
+++ b/src/eko/evolution_operator/c_quad_ker.cpp
@@ -1,4 +1,4 @@
-// g++ -shared -o c_quad_ker.so -fPIC src/eko/evolution_operator/c_quad_ker.cpp
+// g++ -shared -o c_quad_ker.so -fPIC src/eko/evolution_operator/c_quad_ker.cpp -O2
 
 #include <exception>
 #include <vector>
@@ -6,13 +6,30 @@
 #include "./../../ekorepp/harmonics/cache.hpp"
 #include "./../../ekorepp/anomalous_dimensions/unpolarized/space_like.hpp"
 
+/** @brief Talbot path in Mellin inversion */
 class TalbotPath {
+    /** @brief integration variable */
     double t;
+
+    /** @brief bending variable */
     double r;
+
+    /** @brief real offset */
     double o;
+
+    /** @brief auxilary angle */
+    double theta() const { return M_PI * (2.0 * t - 1.0); }
+
 public:
+    /**
+     * @brief Constructor from parameters
+     * @param t integration variable
+     * @param logx log of inversion point
+     * @param is_singlet add real offset?
+     */
     TalbotPath(const double t, const double logx, const bool is_singlet) : t(t), r(0.4 * 16.0 / (1.0 - logx)), o(is_singlet ? 1. : 0.) {}
-    double theta() const { return M_PI * (2.0 * t - 1.0); }
+
+    /** @brief Mellin-N */
     ekorepp::cmplx n() const {
         const double theta = this->theta();
         // treat singular point separately
@@ -20,6 +37,8 @@ class TalbotPath {
         const double im = theta;
         return o + r * ekorepp::cmplx(re, im);
     }
+
+    /** @brief transformation jacobian */
     ekorepp::cmplx jac() const {
         const double theta = this->theta();
         // treat singular point separately
@@ -27,11 +46,14 @@ class TalbotPath {
         const double im = 1.0;
         return r * M_PI * 2.0 * ekorepp::cmplx(re, im);
     }
+
+    /** @brief Mellin inversion prefactor */
     ekorepp::cmplx prefactor() const {
         return ekorepp::cmplx(0,-1./M_PI);
     }
 };
 
+/** @brief Python callback signature  */
 typedef double (*py_quad_ker_qcd)(
                         double* re_gamma, double* im_gamma,
                         double re_n, double im_n,
@@ -52,7 +74,8 @@ typedef double (*py_quad_ker_qcd)(
                         unsigned int ev_op_iterations, unsigned int ev_op_max_order_qcd,
                         unsigned int sv_mode_num, bool is_threshold);
 
-struct QuadCargs {
+/** @brief additional integration parameters */
+struct QuadQCDargs {
     unsigned int order_qcd;
     unsigned int mode0;
     unsigned int mode1;
@@ -75,8 +98,14 @@ struct QuadCargs {
     bool is_threshold;
 };
 
+/**
+ * @brief intergration kernel inside quad
+ * @param u interation variable
+ * @param rargs raw args, to be casted to QuadQCDargs
+ * @return intergration value
+ */
 extern "C" double c_quad_ker_qcd(const double u, void* rargs) {
-    QuadCargs args = *(QuadCargs* )rargs;
+    QuadQCDargs args = *(QuadQCDargs* )rargs;
     const bool is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0);
     // prepare gamma
     const TalbotPath path(u, args.logx, is_singlet);
diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py
index 599934faf..93011b131 100644
--- a/src/eko/evolution_operator/quad_ker.py
+++ b/src/eko/evolution_operator/quad_ker.py
@@ -123,7 +123,7 @@ def select_QEDvalence_element(ker, mode0, mode1):
     ctypes.c_uint,  # nf
     ctypes.c_bool,  # is_log
     ctypes.c_double,  # logx
-    ctypes.c_double * 12,  # areas_raw
+    ctypes.POINTER(ctypes.c_double),  # areas_raw
     ctypes.c_uint,  # polynomial_degreee
     ctypes.c_double,  # L
     ctypes.c_uint,  # method_num
@@ -139,7 +139,7 @@ def select_QEDvalence_element(ker, mode0, mode1):
 c_quad_ker_qcd = libc.c_quad_ker_qcd
 
 
-class QuadCargs(ctypes.Structure):
+class QuadQCDargs(ctypes.Structure):
     """Arguments to C call."""
 
     _fields_ = [
@@ -197,7 +197,7 @@ class QuadCargs(ctypes.Structure):
         nb.types.uint,  # sv_mode_num
         nb.types.bool_,  # is_threshold
     ),
-    cache=False,
+    cache=True,
 )
 def cb_quad_ker_qcd(
     re_gamma_raw,
@@ -226,15 +226,19 @@ def cb_quad_ker_qcd(
     is_threshold,
 ):
     """C Callback inside integration kernel."""
+    # recover complex variables
     n = re_n + im_n * 1j
     jac = re_jac + im_jac * 1j
+    # combute basis functions
     areas = nb.carray(areas_raw, (areas_x, areas_y))
     pj = interpolation.evaluate_grid(n, is_log, logx, areas)
+    # TODO recover parameters
     method = "iterate-expanded"
     sv_mode = sv.Modes.exponentiated
     order = (order_qcd, 0)
     ev_op_max_order = (ev_op_max_order_qcd, 0)
     if is_singlet:
+        # reconstruct singlet matrices
         re_gamma_singlet = nb.carray(re_gamma_raw, (order_qcd, 2, 2))
         im_gamma_singlet = nb.carray(im_gamma_raw, (order_qcd, 2, 2))
         gamma_singlet = re_gamma_singlet + im_gamma_singlet * 1j
@@ -242,6 +246,7 @@ def cb_quad_ker_qcd(
             gamma_singlet = sv.exponentiated.gamma_variation(
                 gamma_singlet, order, nf, L
             )
+        # construct eko
         ker = s.dispatcher(
             order,
             method,
@@ -259,11 +264,13 @@ def cb_quad_ker_qcd(
             ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
+        # construct non-singlet matrices
         re_gamma_ns = nb.carray(re_gamma_raw, order_qcd)
         im_gamma_ns = nb.carray(im_gamma_raw, order_qcd)
         gamma_ns = re_gamma_ns + im_gamma_ns * 1j
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L)
+        # construct eko
         ker = ns.dispatcher(
             order,
             method,