Scale ratios in convolute_eko

pineappl/src/grid.rs

@@ -357,6 +357,8 @@ pub struct GridAxes {
     pub x_grid: Vec<f64>,
     /// Parton IDs used in the grid.
     pub pids: Vec<i32>,
+    /// Interpolation grid for the renormalization scale of the `Grid`.
+    pub mur2_grid: Vec<f64>,
     /// Interpolation grid for the factorization scale of the `Grid`.
     pub muf2_grid: Vec<f64>,
 }
@@ -366,7 +368,7 @@ pub struct GridAxes {
 pub struct EkoInfo {
     /// Scale of the FkTable.
     pub muf2_0: f64,
-    /// Strong coupling constants for the factorization scales in the same ordering as given in
+    /// Strong coupling constants for the renormalization scales in the same ordering as given in
     /// [`GridAxes`].
     pub alphas: Vec<f64>,
     /// Renormalization scale variation.
@@ -1302,6 +1304,7 @@ impl Grid {
             .iter()
             .all(|entry| entry.entry().iter().all(|&(_, b, _)| b == initial_state_2));
 
+        let mut mur2_grid = Vec::new();
         let mut muf2_grid = Vec::new();
         let mut x_grid = Vec::new();
         let pids = Vec::new();
@@ -1341,6 +1344,7 @@ impl Grid {
 
             // not all luminosities are equal (some appear only at higher orders)
             for subgrid in lane.iter() {
+                mur2_grid.append(&mut subgrid.mu2_grid().iter().map(|mu2| mu2.ren).collect());
                 muf2_grid.append(&mut subgrid.mu2_grid().iter().map(|mu2| mu2.fac).collect());
                 if has_pdf1 {
                     x_grid.extend_from_slice(&subgrid.x1_grid());
@@ -1351,14 +1355,18 @@ impl Grid {
             }
         }
 
-        muf2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!()));
-        muf2_grid.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = 64));
-        x_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!()));
+        // make grids unique
+        x_grid.sort_by(|a, b| a.partial_cmp(b).unwrap());
         x_grid.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = 64));
+        mur2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        mur2_grid.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = 64));
+        muf2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        muf2_grid.dedup_by(|a, b| approx_eq!(f64, *a, *b, ulps = 64));
 
         Some(GridAxes {
             x_grid,
             pids, // TODO: for the time being they are just empty, but we might use them for slicing the eko
+            mur2_grid,
             muf2_grid,
         })
     }
@@ -1381,6 +1389,7 @@ impl Grid {
         // Check operator layout
         let dim = operator.shape();
 
+        assert_eq!(dim[0], eko_info.grid_axes.mur2_grid.len());
         assert_eq!(dim[0], eko_info.grid_axes.muf2_grid.len());
         assert_eq!(dim[1], eko_info.target_pids.len());
         assert_eq!(dim[3], eko_info.grid_axes.pids.len());
@@ -1567,15 +1576,15 @@ impl Grid {
                             .unwrap();
                         let als_iq2 = eko_info
                             .grid_axes
-                            .muf2_grid
+                            .mur2_grid
                             .iter()
                             .position(|&q2| {
                                 approx_eq!(f64, q2, eko_info.xir * eko_info.xir * scale, ulps = 64)
                             })
                             .unwrap_or_else(|| {
                                 panic!(
-                                    "Couldn't find q2: {:?} with xir: {:?} and muf2_grid: {:?}",
-                                    scale, eko_info.xir, eko_info.grid_axes.muf2_grid
+                                    "Couldn't find mur2: {:?} with xir: {:?} and mur2_grid: {:?}",
+                                    scale, eko_info.xir, eko_info.grid_axes.mur2_grid
                                 )
                             });
 
@@ -1605,8 +1614,13 @@ impl Grid {
                             .grid_axes
                             .muf2_grid
                             .iter()
-                            .position(|&q2| approx_eq!(f64, q2, src_q2, ulps = 64))
-                            .unwrap()
+                            .position(|&q2| approx_eq!(f64, q2, eko_info.xif * eko_info.xif * src_q2, ulps = 64))
+                            .unwrap_or_else(|| {
+                                panic!(
+                                    "Couldn't find muf2: {:?} with xif: {:?} and muf2_grid: {:?}",
+                                    src_q2, eko_info.xif, eko_info.grid_axes.muf2_grid
+                                )
+                            })
                     })
                     .collect();
                 // Iterate target lumis
@@ -2283,6 +2297,7 @@ mod tests {
     // TODO: properly test axes returned
 
     fn simple_grid() -> (Grid, GridAxes) {
+        let mur2_grid = vec![20.];
         let muf2_grid = vec![20.];
         let x_grid = vec![0.1, 0.5, 1.];
 
@@ -2323,6 +2338,7 @@ mod tests {
             GridAxes {
                 x_grid,
                 pids,
+                mur2_grid,
                 muf2_grid,
             },
         )
@@ -2334,6 +2350,7 @@ mod tests {
 
         let ret_axes = grid.axes().unwrap();
         assert_eq!(ret_axes.x_grid, axes.x_grid);
+        assert_eq!(ret_axes.mur2_grid, axes.mur2_grid);
         assert_eq!(ret_axes.muf2_grid, axes.muf2_grid);
         assert_eq!(ret_axes.pids, vec![]);
     }
@@ -2356,6 +2373,7 @@ mod tests {
             grid_axes: GridAxes {
                 x_grid: axes.x_grid.clone(),
                 pids: axes.pids.clone(),
+                mur2_grid: axes.mur2_grid.clone(),
                 muf2_grid: axes.muf2_grid.clone(),
             },
             lumi_id_types,

pineappl_py/pineappl/grid.py

@@ -245,7 +245,7 @@ def convolute_with_one(
             xi,
         )
 
-    def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids", order_mask=()):
+    def convolute_eko(self, operators, mur2_grid, alphas_values, lumi_id_types="pdg_mc_ids", order_mask=(), xi=(1.0, 1.0)):
         """
         Create an FKTable with the EKO.
 
@@ -255,12 +255,22 @@ def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids", order_mask=()):
         ----------
             operators : dict
                 EKO Output
+            mur2_grid : list[float]
+                renormalization scales
+            alphas_values : list[float]
+                alpha_s values associated to the renormalization scales
             lumi_id_types : str
                 kind of lumi types (e.g. "pdg_mc_ids" for flavor basis, "evol"
                 for evolution basis)
             order_mask : list(bool)
                 Mask for selecting specific orders. The value `True` means the corresponding order
                 is included. An empty list corresponds to all orders being enabled.
+            xi : (float, float)
+                A tuple with the scale variation factors that should be used.
+                The first entry of a tuple corresponds to the variation of
+                the renormalization scale, the second entry to the variation of the factorization
+                scale. If only results for the central scale are need the tuple should be
+                `(1.0, 1.0)`.
 
         Returns
         ------
@@ -271,19 +281,20 @@ def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids", order_mask=()):
             [op["operators"] for op in operators["Q2grid"].values()]
         )
         q2grid = list(operators["Q2grid"].keys())
-        alphas_values = [op["alphas"] for op in operators["Q2grid"].values()]
         return FkTable(
             self.raw.convolute_eko(
                 operators["q2_ref"],
-                np.array(alphas_values),
+                np.array(alphas_values, dtype=np.float64),
                 np.array(operators["targetpids"], dtype=np.int32),
                 np.array(operators["targetgrid"]),
                 np.array(operators["inputpids"], dtype=np.int32),
                 np.array(operators["inputgrid"]),
+                np.array(mur2_grid, dtype=np.float64),
                 np.array(q2grid, dtype=np.float64),
                 np.array(operator_grid),
                 lumi_id_types,
                 np.array(order_mask, dtype=bool),
+                xi
             )
         )
 

pineappl_py/src/grid.rs

@@ -297,20 +297,29 @@ impl PyGrid {
     ///         interpolation grid
     ///     pids: numpy.ndarray(int)
     ///         particle ids
+    ///     mur2_grid : numpy.ndarray(float)
+    ///         factorization scale list
     ///     muf2_grid : numpy.ndarray(float)
     ///         factorization scale list
     pub fn axes<'py>(
         &self,
         py: Python<'py>,
-    ) -> (&'py PyArray1<f64>, &'py PyArray1<i32>, &'py PyArray1<f64>) {
+    ) -> (
+        &'py PyArray1<f64>,
+        &'py PyArray1<i32>,
+        &'py PyArray1<f64>,
+        &'py PyArray1<f64>,
+    ) {
         let GridAxes {
             x_grid,
             pids,
+            mur2_grid,
             muf2_grid,
         } = self.grid.axes().unwrap();
         (
             x_grid.into_pyarray(py),
             pids.into_pyarray(py),
+            mur2_grid.into_pyarray(py),
             muf2_grid.into_pyarray(py),
         )
     }
@@ -392,8 +401,10 @@ impl PyGrid {
     ///         sorting of the particles in the tensor for final FkTable
     ///     target_x_grid : numpy.ndarray(float)
     ///         final FKTable interpolation grid
+    ///     mur2_grid : numpy.ndarray(float)
+    ///         list of renormalization scales
     ///     muf2_grid : numpy.ndarray(float)
-    ///         list of process scales
+    ///         list of factorization scales
     ///     operator : numpy.ndarray(int, rank=5)
     ///         evolution tensor
     ///     orders_mask: numpy.ndarray(bool)
@@ -411,21 +422,24 @@ impl PyGrid {
         x_grid: PyReadonlyArray1<f64>,
         target_pids: PyReadonlyArray1<i32>,
         target_x_grid: PyReadonlyArray1<f64>,
+        mur2_grid: PyReadonlyArray1<f64>,
         muf2_grid: PyReadonlyArray1<f64>,
         operator: PyReadonlyArray5<f64>,
         lumi_id_types: String,
         order_mask: PyReadonlyArray1<bool>,
+        xi: (f64, f64),
     ) -> PyFkTable {
         let eko_info = EkoInfo {
             muf2_0,
             alphas: alphas.to_vec().unwrap(),
-            xir: 1.,
-            xif: 1.,
+            xir: xi.0,
+            xif: xi.1,
             target_pids: target_pids.to_vec().unwrap(),
             target_x_grid: target_x_grid.to_vec().unwrap(),
             grid_axes: GridAxes {
                 x_grid: x_grid.to_vec().unwrap(),
                 pids: pids.to_vec().unwrap(),
+                mur2_grid: mur2_grid.to_vec().unwrap(),
                 muf2_grid: muf2_grid.to_vec().unwrap(),
             },
             lumi_id_types,