From 7569391cbe775957acc7f95131f0037c149d869e Mon Sep 17 00:00:00 2001 From: t7phy Date: Fri, 4 Oct 2024 16:26:24 +0200 Subject: [PATCH] Add support for a fragmentation scale --- pineappl/src/convolutions.rs | 74 +++++++++++++++++++++++++++++++----- pineappl/src/grid.rs | 25 ++++++------ 2 files changed, 78 insertions(+), 21 deletions(-) diff --git a/pineappl/src/convolutions.rs b/pineappl/src/convolutions.rs index f4a7f5b6..98a2b594 100644 --- a/pineappl/src/convolutions.rs +++ b/pineappl/src/convolutions.rs @@ -16,9 +16,11 @@ pub struct ConvolutionCache<'a> { alphas_cache: Vec, mur2_grid: Vec, muf2_grid: Vec, + mua2_grid: Vec, x_grid: Vec, imur2: Vec, imuf2: Vec, + imua2: Vec, ix: Vec>, pdg: Vec, perm: Vec>, @@ -38,16 +40,18 @@ impl<'a> ConvolutionCache<'a> { alphas_cache: Vec::new(), mur2_grid: Vec::new(), muf2_grid: Vec::new(), + mua2_grid: Vec::new(), x_grid: Vec::new(), imur2: Vec::new(), imuf2: Vec::new(), + imua2: Vec::new(), ix: Vec::new(), pdg, perm: Vec::new(), } } - pub(crate) fn setup(&mut self, grid: &Grid, xi: &[(f64, f64)]) -> Result<(), ()> { + pub(crate) fn setup(&mut self, grid: &Grid, xi: &[(f64, f64, f64)]) -> Result<(), ()> { let convolutions = grid.convolutions(); // TODO: the following code only works with exactly two convolutions @@ -129,7 +133,7 @@ impl<'a> ConvolutionCache<'a> { .flatten() .flat_map(|ren| { xi.iter() - .map(|(xir, _)| xir * xir * ren) + .map(|(xir, _, _)| xir * xir * ren) .collect::>() }) .collect(); @@ -161,17 +165,50 @@ impl<'a> ConvolutionCache<'a> { .flatten() .flat_map(|fac| { xi.iter() - .map(|(_, xif)| xif * xif * fac) + .map(|(_, xif, _)| xif * xif * fac) .collect::>() }) .collect(); muf2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!())); muf2_grid.dedup(); + let mut mua2_grid: Vec<_> = grid + .subgrids() + .iter() + .filter_map(|subgrid| { + if subgrid.is_empty() { + None + } else { + Some( + grid.kinematics() + .iter() + .zip(subgrid.node_values()) + .find_map(|(kin, node_values)| { + // TODO: generalize this for arbitrary scales + matches!(kin, &Kinematics::Scale(idx) if idx == 0) + .then_some(node_values) + }) + // TODO: convert this into an error + .unwrap() + .values(), + ) + } + }) + .flatten() + .flat_map(|frg| { + xi.iter() + .map(|(_, _, xia)| xia * xia * frg) + .collect::>() + }) + .collect(); + mua2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!())); + mua2_grid.dedup(); + self.alphas_cache = mur2_grid.iter().map(|&mur2| (self.alphas)(mur2)).collect(); self.mur2_grid = mur2_grid; self.muf2_grid = muf2_grid; + self.mua2_grid = mua2_grid; self.x_grid = x_grid; Ok(()) @@ -191,8 +228,7 @@ impl<'a> ConvolutionCache<'a> { .filter_map(|(index, (perm, &pdg_id))| { perm.map(|(idx, cc)| { let ix = self.ix[index][indices[index + 1]]; - let imuf2 = self.imuf2[indices[0]]; - let muf2 = self.muf2_grid[imuf2]; + let pid = if cc { pids::charge_conjugate_pdg_pid(pdg_id) } else { @@ -200,9 +236,20 @@ impl<'a> ConvolutionCache<'a> { }; let xfx = &mut self.xfx[idx]; let xfx_cache = &mut self.xfx_cache[idx]; - *xfx_cache.entry((pid, ix, imuf2)).or_insert_with(|| { + let (imu2, mu2) = match self.pdg[idx] { + Convolution::UnpolPDF(_) | Convolution::PolPDF(_) => { + let imuf2 = self.imuf2[indices[0]]; + (imuf2, self.muf2_grid[imuf2]) + } + Convolution::UnpolFF(_) | Convolution::PolFF(_) => { + let imua2 = self.imua2[indices[0]]; + (imua2, self.mua2_grid[imua2]) + } + Convolution::None => unreachable!(), + }; + *xfx_cache.entry((pid, ix, imu2)).or_insert_with(|| { let x = self.x_grid[ix]; - xfx(pid, x, muf2) / x + xfx(pid, x, mu2) / x }) }) }) @@ -224,6 +271,7 @@ impl<'a> ConvolutionCache<'a> { } self.mur2_grid.clear(); self.muf2_grid.clear(); + self.mua2_grid.clear(); self.x_grid.clear(); } @@ -237,9 +285,6 @@ impl<'a> ConvolutionCache<'a> { xif: f64, xia: f64, ) { - // TODO: generalize this for fragmentation functions - assert_eq!(xia, 1.0); - self.imur2 = mu2_grid .iter() .map(|ren| { @@ -258,6 +303,15 @@ impl<'a> ConvolutionCache<'a> { .unwrap_or_else(|| unreachable!()) }) .collect(); + self.imua2 = mu2_grid + .iter() + .map(|frg| { + self.mua2_grid + .iter() + .position(|&mua2| mua2 == xia * xia * frg) + .unwrap_or_else(|| unreachable!()) + }) + .collect(); // TODO: generalize this for arbitrary orderings of x self.ix = (0..grid.convolutions().len()) diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs index e4a0493c..2f176cc4 100644 --- a/pineappl/src/grid.rs +++ b/pineappl/src/grid.rs @@ -232,13 +232,6 @@ impl Grid { channel_mask: &[bool], xi: &[(f64, f64, f64)], ) -> Vec { - assert!(xi - .iter() - .all(|&(_, _, xia)| approx_eq!(f64, xia, 1.0, ulps = 2))); - let xi = xi - .iter() - .map(|&(xir, xif, _)| (xir, xif)) - .collect::>(); convolution_cache.setup(self, &xi).unwrap(); let bin_indices = if bin_indices.is_empty() { @@ -250,11 +243,14 @@ impl Grid { let normalizations = self.bin_info().normalizations(); let pdg_channels = self.channels_pdg(); - for (xi_index, &(xir, xif)) in xi.iter().enumerate() { + for (xi_index, &(xir, xif, xia)) in xi.iter().enumerate() { for ((ord, bin, chan), subgrid) in self.subgrids.indexed_iter() { let order = &self.orders[ord]; - if ((order.logxir > 0) && (xir == 1.0)) || ((order.logxif > 0) && (xif == 1.0)) { + if ((order.logxir > 0) && (xir == 1.0)) + || ((order.logxif > 0) && (xif == 1.0)) + || ((order.logxia > 0) && (xia == 1.0)) + { continue; } @@ -316,6 +312,10 @@ impl Grid { value *= (xif * xif).ln().powi(order.logxif.try_into().unwrap()); } + if order.logxia > 0 { + value *= (xia * xia).ln().powi(order.logxia.try_into().unwrap()); + } + bins[xi_index + xi.len() * bin_index] += value / normalizations[bin]; } } @@ -339,8 +339,7 @@ impl Grid { channel: usize, (xir, xif, xia): (f64, f64, f64), ) -> ArrayD { - assert_eq!(xia, 1.0); - convolution_cache.setup(self, &[(xir, xif)]).unwrap(); + convolution_cache.setup(self, &[(xir, xif, xia)]).unwrap(); let normalizations = self.bin_info().normalizations(); let pdg_channels = self.channels_pdg(); @@ -395,6 +394,10 @@ impl Grid { array *= (xif * xif).ln().powi(order.logxif.try_into().unwrap()); } + if order.logxia > 0 { + array *= (xia * xia).ln().powi(order.logxia.try_into().unwrap()); + } + array /= normalizations[bin]; array }