|
1 | 1 | //! Module for everything related to luminosity functions. |
2 | 2 |
|
3 | | -use super::grid::Grid; |
4 | 3 | use super::pids; |
5 | | -use super::subgrid::{Mu2, Subgrid}; |
6 | | -use rustc_hash::FxHashMap; |
7 | | - |
8 | | -enum Pdfs<'a> { |
9 | | - Two { |
10 | | - xfx1: &'a mut dyn FnMut(i32, f64, f64) -> f64, |
11 | | - xfx1_cache: FxHashMap<(i32, usize, usize), f64>, |
12 | | - xfx2: &'a mut dyn FnMut(i32, f64, f64) -> f64, |
13 | | - xfx2_cache: FxHashMap<(i32, usize, usize), f64>, |
14 | | - }, |
15 | | - One { |
16 | | - xfx: &'a mut dyn FnMut(i32, f64, f64) -> f64, |
17 | | - xfx_cache: FxHashMap<(i32, usize, usize), f64>, |
18 | | - }, |
19 | | -} |
20 | | - |
21 | | -impl<'a> Pdfs<'a> { |
22 | | - pub fn clear(&mut self) { |
23 | | - match self { |
24 | | - Self::One { xfx_cache, .. } => xfx_cache.clear(), |
25 | | - Self::Two { |
26 | | - xfx1_cache, |
27 | | - xfx2_cache, |
28 | | - .. |
29 | | - } => { |
30 | | - xfx1_cache.clear(); |
31 | | - xfx2_cache.clear(); |
32 | | - } |
33 | | - } |
34 | | - } |
35 | | -} |
36 | | - |
37 | | -/// A cache for evaluating PDFs. Methods like [`Grid::convolve`] accept instances of this `struct` |
38 | | -/// instead of the PDFs themselves. |
39 | | -pub struct LumiCache<'a> { |
40 | | - pdfs: Pdfs<'a>, |
41 | | - alphas: &'a mut dyn FnMut(f64) -> f64, |
42 | | - alphas_cache: Vec<f64>, |
43 | | - mur2_grid: Vec<f64>, |
44 | | - muf2_grid: Vec<f64>, |
45 | | - x_grid: Vec<f64>, |
46 | | - imur2: Vec<usize>, |
47 | | - imuf2: Vec<usize>, |
48 | | - ix1: Vec<usize>, |
49 | | - ix2: Vec<usize>, |
50 | | - pdg1: i32, |
51 | | - pdg2: i32, |
52 | | - cc1: i32, |
53 | | - cc2: i32, |
54 | | -} |
55 | | - |
56 | | -impl<'a> LumiCache<'a> { |
57 | | - /// Construct a luminosity cache with two PDFs, `xfx1` and `xfx2`. The types of hadrons the |
58 | | - /// PDFs correspond to must be given as `pdg1` and `pdg2`. The function to evaluate the |
59 | | - /// strong coupling must be given as `alphas`. The grid that the cache will be used with must |
60 | | - /// be given as `grid`; this parameter determines which of the initial states are hadronic, and |
61 | | - /// if an initial states is not hadronic the corresponding 'PDF' is set to `xfx = x`. If some |
62 | | - /// of the PDFs must be charge-conjugated, this is automatically done in this function. |
63 | | - pub fn with_two( |
64 | | - pdg1: i32, |
65 | | - xfx1: &'a mut dyn FnMut(i32, f64, f64) -> f64, |
66 | | - pdg2: i32, |
67 | | - xfx2: &'a mut dyn FnMut(i32, f64, f64) -> f64, |
68 | | - alphas: &'a mut dyn FnMut(f64) -> f64, |
69 | | - ) -> Self { |
70 | | - Self { |
71 | | - pdfs: Pdfs::Two { |
72 | | - xfx1, |
73 | | - xfx1_cache: FxHashMap::default(), |
74 | | - xfx2, |
75 | | - xfx2_cache: FxHashMap::default(), |
76 | | - }, |
77 | | - alphas, |
78 | | - alphas_cache: vec![], |
79 | | - mur2_grid: vec![], |
80 | | - muf2_grid: vec![], |
81 | | - x_grid: vec![], |
82 | | - imur2: Vec::new(), |
83 | | - imuf2: Vec::new(), |
84 | | - ix1: Vec::new(), |
85 | | - ix2: Vec::new(), |
86 | | - pdg1, |
87 | | - pdg2, |
88 | | - cc1: 0, |
89 | | - cc2: 0, |
90 | | - } |
91 | | - } |
92 | | - |
93 | | - /// Construct a luminosity cache with a single PDF `xfx`. The type of hadron the PDF |
94 | | - /// corresponds to must be given as `pdg`. The function to evaluate the strong coupling must be |
95 | | - /// given as `alphas`. The grid that the cache should be used with must be given as `grid`; |
96 | | - /// this parameter determines which of the initial states are hadronic, and if an initial |
97 | | - /// states is not hadronic the corresponding 'PDF' is set to `xfx = x`. If some of the PDFs |
98 | | - /// must be charge-conjugated, this is automatically done in this function. |
99 | | - pub fn with_one( |
100 | | - pdg: i32, |
101 | | - xfx: &'a mut dyn FnMut(i32, f64, f64) -> f64, |
102 | | - alphas: &'a mut dyn FnMut(f64) -> f64, |
103 | | - ) -> Self { |
104 | | - Self { |
105 | | - pdfs: Pdfs::One { |
106 | | - xfx, |
107 | | - xfx_cache: FxHashMap::default(), |
108 | | - }, |
109 | | - alphas, |
110 | | - alphas_cache: vec![], |
111 | | - mur2_grid: vec![], |
112 | | - muf2_grid: vec![], |
113 | | - x_grid: vec![], |
114 | | - imur2: Vec::new(), |
115 | | - imuf2: Vec::new(), |
116 | | - ix1: Vec::new(), |
117 | | - ix2: Vec::new(), |
118 | | - pdg1: pdg, |
119 | | - pdg2: pdg, |
120 | | - cc1: 0, |
121 | | - cc2: 0, |
122 | | - } |
123 | | - } |
124 | | - |
125 | | - pub(crate) fn setup(&mut self, grid: &Grid, xi: &[(f64, f64)]) -> Result<(), ()> { |
126 | | - let convolutions = grid.convolutions(); |
127 | | - |
128 | | - // TODO: the following code only works with exactly two convolutions |
129 | | - assert_eq!(convolutions.len(), 2); |
130 | | - |
131 | | - // do we have to charge-conjugate the initial states? |
132 | | - let cc1 = if let Some(pid) = convolutions[0].pid() { |
133 | | - if self.pdg1 == pid { |
134 | | - 1 |
135 | | - } else if self.pdg1 == pids::charge_conjugate_pdg_pid(pid) { |
136 | | - -1 |
137 | | - } else { |
138 | | - // TODO: return a proper error |
139 | | - return Err(()); |
140 | | - } |
141 | | - } else { |
142 | | - 0 |
143 | | - }; |
144 | | - let cc2 = if let Some(pid) = convolutions[1].pid() { |
145 | | - if self.pdg2 == pid { |
146 | | - 1 |
147 | | - } else if self.pdg2 == pids::charge_conjugate_pdg_pid(pid) { |
148 | | - -1 |
149 | | - } else { |
150 | | - // TODO: return a proper error |
151 | | - return Err(()); |
152 | | - } |
153 | | - } else { |
154 | | - 0 |
155 | | - }; |
156 | | - |
157 | | - // TODO: try to avoid calling clear |
158 | | - self.clear(); |
159 | | - |
160 | | - let mut x_grid: Vec<_> = grid |
161 | | - .subgrids() |
162 | | - .iter() |
163 | | - .filter_map(|subgrid| { |
164 | | - if subgrid.is_empty() { |
165 | | - None |
166 | | - } else { |
167 | | - let mut vec = subgrid.x1_grid().into_owned(); |
168 | | - vec.extend_from_slice(&subgrid.x2_grid()); |
169 | | - Some(vec) |
170 | | - } |
171 | | - }) |
172 | | - .flatten() |
173 | | - .collect(); |
174 | | - x_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!())); |
175 | | - x_grid.dedup(); |
176 | | - |
177 | | - let mut mur2_grid: Vec<_> = grid |
178 | | - .subgrids() |
179 | | - .iter() |
180 | | - .filter_map(|subgrid| { |
181 | | - if subgrid.is_empty() { |
182 | | - None |
183 | | - } else { |
184 | | - Some(subgrid.mu2_grid().into_owned()) |
185 | | - } |
186 | | - }) |
187 | | - .flatten() |
188 | | - .flat_map(|Mu2 { ren, .. }| { |
189 | | - xi.iter() |
190 | | - .map(|(xir, _)| xir * xir * ren) |
191 | | - .collect::<Vec<_>>() |
192 | | - }) |
193 | | - .collect(); |
194 | | - mur2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!())); |
195 | | - mur2_grid.dedup(); |
196 | | - |
197 | | - let mut muf2_grid: Vec<_> = grid |
198 | | - .subgrids() |
199 | | - .iter() |
200 | | - .filter_map(|subgrid| { |
201 | | - if subgrid.is_empty() { |
202 | | - None |
203 | | - } else { |
204 | | - Some(subgrid.mu2_grid().into_owned()) |
205 | | - } |
206 | | - }) |
207 | | - .flatten() |
208 | | - .flat_map(|Mu2 { fac, .. }| { |
209 | | - xi.iter() |
210 | | - .map(|(_, xif)| xif * xif * fac) |
211 | | - .collect::<Vec<_>>() |
212 | | - }) |
213 | | - .collect(); |
214 | | - muf2_grid.sort_by(|a, b| a.partial_cmp(b).unwrap_or_else(|| unreachable!())); |
215 | | - muf2_grid.dedup(); |
216 | | - |
217 | | - self.alphas_cache = mur2_grid.iter().map(|&mur2| (self.alphas)(mur2)).collect(); |
218 | | - |
219 | | - self.mur2_grid = mur2_grid; |
220 | | - self.muf2_grid = muf2_grid; |
221 | | - self.x_grid = x_grid; |
222 | | - self.cc1 = cc1; |
223 | | - self.cc2 = cc2; |
224 | | - |
225 | | - Ok(()) |
226 | | - } |
227 | | - |
228 | | - /// Return the PDF (multiplied with `x`) for the first initial state. |
229 | | - pub fn xfx1(&mut self, pdg_id: i32, ix1: usize, imu2: usize) -> f64 { |
230 | | - let ix1 = self.ix1[ix1]; |
231 | | - let x = self.x_grid[ix1]; |
232 | | - if self.cc1 == 0 { |
233 | | - x |
234 | | - } else { |
235 | | - let imuf2 = self.imuf2[imu2]; |
236 | | - let muf2 = self.muf2_grid[imuf2]; |
237 | | - let pid = if self.cc1 == 1 { |
238 | | - pdg_id |
239 | | - } else { |
240 | | - pids::charge_conjugate_pdg_pid(pdg_id) |
241 | | - }; |
242 | | - let (xfx, xfx_cache) = match &mut self.pdfs { |
243 | | - Pdfs::One { xfx, xfx_cache, .. } => (xfx, xfx_cache), |
244 | | - Pdfs::Two { |
245 | | - xfx1, xfx1_cache, .. |
246 | | - } => (xfx1, xfx1_cache), |
247 | | - }; |
248 | | - *xfx_cache |
249 | | - .entry((pid, ix1, imuf2)) |
250 | | - .or_insert_with(|| xfx(pid, x, muf2)) |
251 | | - } |
252 | | - } |
253 | | - |
254 | | - /// Return the PDF (multiplied with `x`) for the second initial state. |
255 | | - pub fn xfx2(&mut self, pdg_id: i32, ix2: usize, imu2: usize) -> f64 { |
256 | | - let ix2 = self.ix2[ix2]; |
257 | | - let x = self.x_grid[ix2]; |
258 | | - if self.cc2 == 0 { |
259 | | - x |
260 | | - } else { |
261 | | - let imuf2 = self.imuf2[imu2]; |
262 | | - let muf2 = self.muf2_grid[imuf2]; |
263 | | - let pid = if self.cc2 == 1 { |
264 | | - pdg_id |
265 | | - } else { |
266 | | - pids::charge_conjugate_pdg_pid(pdg_id) |
267 | | - }; |
268 | | - let (xfx, xfx_cache) = match &mut self.pdfs { |
269 | | - Pdfs::One { xfx, xfx_cache, .. } => (xfx, xfx_cache), |
270 | | - Pdfs::Two { |
271 | | - xfx2, xfx2_cache, .. |
272 | | - } => (xfx2, xfx2_cache), |
273 | | - }; |
274 | | - *xfx_cache |
275 | | - .entry((pid, ix2, imuf2)) |
276 | | - .or_insert_with(|| xfx(pid, x, muf2)) |
277 | | - } |
278 | | - } |
279 | | - |
280 | | - /// Return the strong coupling for the renormalization scale set with [`LumiCache::set_grids`], |
281 | | - /// in the grid `mu2_grid` at the index `imu2`. |
282 | | - #[must_use] |
283 | | - pub fn alphas(&self, imu2: usize) -> f64 { |
284 | | - self.alphas_cache[self.imur2[imu2]] |
285 | | - } |
286 | | - |
287 | | - /// Clears the cache. |
288 | | - pub fn clear(&mut self) { |
289 | | - self.alphas_cache.clear(); |
290 | | - self.pdfs.clear(); |
291 | | - self.mur2_grid.clear(); |
292 | | - self.muf2_grid.clear(); |
293 | | - self.x_grid.clear(); |
294 | | - } |
295 | | - |
296 | | - /// Set the grids. |
297 | | - pub fn set_grids( |
298 | | - &mut self, |
299 | | - mu2_grid: &[Mu2], |
300 | | - x1_grid: &[f64], |
301 | | - x2_grid: &[f64], |
302 | | - xir: f64, |
303 | | - xif: f64, |
304 | | - ) { |
305 | | - self.imur2 = mu2_grid |
306 | | - .iter() |
307 | | - .map(|Mu2 { ren, .. }| { |
308 | | - self.mur2_grid |
309 | | - .iter() |
310 | | - .position(|&mur2| mur2 == xir * xir * ren) |
311 | | - .unwrap_or_else(|| unreachable!()) |
312 | | - }) |
313 | | - .collect(); |
314 | | - self.imuf2 = mu2_grid |
315 | | - .iter() |
316 | | - .map(|Mu2 { fac, .. }| { |
317 | | - self.muf2_grid |
318 | | - .iter() |
319 | | - .position(|&muf2| muf2 == xif * xif * fac) |
320 | | - .unwrap_or_else(|| unreachable!()) |
321 | | - }) |
322 | | - .collect(); |
323 | | - self.ix1 = x1_grid |
324 | | - .iter() |
325 | | - .map(|x1| { |
326 | | - self.x_grid |
327 | | - .iter() |
328 | | - .position(|x| x1 == x) |
329 | | - .unwrap_or_else(|| unreachable!()) |
330 | | - }) |
331 | | - .collect(); |
332 | | - |
333 | | - self.ix2 = x2_grid |
334 | | - .iter() |
335 | | - .map(|x2| { |
336 | | - self.x_grid |
337 | | - .iter() |
338 | | - .position(|x| x2 == x) |
339 | | - .unwrap_or_else(|| unreachable!()) |
340 | | - }) |
341 | | - .collect(); |
342 | | - } |
343 | | -} |
344 | 4 |
|
345 | 5 | /// Data type that indentifies different types of convolutions. |
346 | 6 | #[derive(Debug, Eq, PartialEq)] |
|
0 commit comments