diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index 40342db2..9c1b436f 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -357,6 +357,8 @@ pub struct GridAxes { pub x_grid: Vec, /// Parton IDs used in the grid. pub pids: Vec, + /// Interpolation grid for the renormalization scale of the `Grid`. + pub mur2_grid: Vec, /// Interpolation grid for the factorization scale of the `Grid`. pub muf2_grid: Vec, } @@ -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, /// 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, diff --git a/pineappl_py/pineappl/grid.py b/pineappl_py/pineappl/grid.py index 61887036..5c6a7ef2 100644 --- a/pineappl_py/pineappl/grid.py +++ b/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 ) ) diff --git a/pineappl_py/src/grid.rs b/pineappl_py/src/grid.rs index d377b2cb..dee63be7 100644 --- a/pineappl_py/src/grid.rs +++ b/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, &'py PyArray1, &'py PyArray1) { + ) -> ( + &'py PyArray1, + &'py PyArray1, + &'py PyArray1, + &'py PyArray1, + ) { 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, target_pids: PyReadonlyArray1, target_x_grid: PyReadonlyArray1, + mur2_grid: PyReadonlyArray1, muf2_grid: PyReadonlyArray1, operator: PyReadonlyArray5, lumi_id_types: String, order_mask: PyReadonlyArray1, + 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,