From ddbf227937af035574372d6a76a34c990dc71d93 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 27 Jun 2024 14:38:42 +0200 Subject: [PATCH] Migrate the remaining parts of `plot` to `ConvFuns` --- pineappl_cli/src/helpers.rs | 115 ++++++++++++++++++++++++------------ pineappl_cli/src/plot.rs | 63 +++++++++++--------- 2 files changed, 113 insertions(+), 65 deletions(-) diff --git a/pineappl_cli/src/helpers.rs b/pineappl_cli/src/helpers.rs index 102d03e4..d1871f19 100644 --- a/pineappl_cli/src/helpers.rs +++ b/pineappl_cli/src/helpers.rs @@ -97,23 +97,6 @@ pub fn create_conv_funs_for_set( Ok((set, conv_funs)) } -pub fn create_pdfset(pdfset: &str) -> Result<(PdfSet, Option)> { - let pdfset = pdfset.split_once('=').map_or(pdfset, |(name, _)| name); - let (pdfset, member) = pdfset - .rsplit_once('/') - .map_or((pdfset, None), |(set, member)| { - (set, Some(member.parse::().unwrap())) - }); - - Ok(( - PdfSet::new(&pdfset.parse().map_or_else( - |_| pdfset.to_owned(), - |lhaid| lhapdf::lookup_pdf(lhaid).unwrap().0, - ))?, - member, - )) -} - pub fn read_grid(input: &Path) -> Result { Grid::read(File::open(input).context(format!("unable to open '{}'", input.display()))?) .context(format!("unable to read '{}'", input.display())) @@ -381,36 +364,94 @@ pub fn convolve_limits(grid: &Grid, bins: &[usize], mode: ConvoluteMode) -> Vec< pub fn convolve_subgrid( grid: &Grid, - lhapdf: &mut Pdf, + conv_funs: &mut [Pdf], order: usize, bin: usize, lumi: usize, cfg: &GlobalConfiguration, ) -> Array3 { - // if the field 'Particle' is missing we assume it's a proton PDF - let pdf_pdg_id = lhapdf - .set() - .entry("Particle") - .map_or(Ok(2212), |string| string.parse::()) - .unwrap(); - if cfg.force_positive { - lhapdf.set_force_positive(1); + for fun in conv_funs.iter_mut() { + fun.set_force_positive(1); + } } - let x_max = lhapdf.x_max(); - let x_min = lhapdf.x_min(); - let mut pdf = |id, x, q2| { - if !cfg.allow_extrapolation && (x < x_min || x > x_max) { - 0.0 - } else { - lhapdf.xfx_q2(id, x, q2) + match conv_funs { + [fun] => { + // there's only one convolution function from which we can use the strong coupling + assert_eq!(cfg.use_alphas_from, 0); + + // if the field 'Particle' is missing we assume it's a proton PDF + let pdg_id = fun + .set() + .entry("Particle") + .map_or(Ok(2212), |string| string.parse::()) + .unwrap(); + + let x_max = fun.x_max(); + let x_min = fun.x_min(); + let mut alphas = |q2| fun.alphas_q2(q2); + let mut fun = |id, x, q2| { + if !cfg.allow_extrapolation && (x < x_min || x > x_max) { + 0.0 + } else { + fun.xfx_q2(id, x, q2) + } + }; + + let mut cache = LumiCache::with_one(pdg_id, &mut fun, &mut alphas); + + grid.convolve_subgrid(&mut cache, order, bin, lumi, 1.0, 1.0) } - }; - let mut alphas = |q2| lhapdf.alphas_q2(q2); - let mut cache = LumiCache::with_one(pdf_pdg_id, &mut pdf, &mut alphas); + [fun1, fun2] => { + let pdg_id1 = fun1 + .set() + .entry("Particle") + .map_or(Ok(2212), |string| string.parse::()) + .unwrap(); + + let pdg_id2 = fun2 + .set() + .entry("Particle") + .map_or(Ok(2212), |string| string.parse::()) + .unwrap(); + + let x_max1 = fun1.x_max(); + let x_min1 = fun1.x_min(); + let x_max2 = fun2.x_max(); + let x_min2 = fun2.x_min(); + + let mut alphas = |q2| match cfg.use_alphas_from { + 0 => fun1.alphas_q2(q2), + 1 => fun2.alphas_q2(q2), + _ => panic!( + "expected `use_alphas_from` to be `0` or `1`, is {}", + cfg.use_alphas_from + ), + }; + let mut fun1 = |id, x, q2| { + if !cfg.allow_extrapolation && (x < x_min1 || x > x_max1) { + 0.0 + } else { + fun1.xfx_q2(id, x, q2) + } + }; + + let mut fun2 = |id, x, q2| { + if !cfg.allow_extrapolation && (x < x_min2 || x > x_max2) { + 0.0 + } else { + fun2.xfx_q2(id, x, q2) + } + }; + + let mut cache = + LumiCache::with_two(pdg_id1, &mut fun1, pdg_id2, &mut fun2, &mut alphas); - grid.convolve_subgrid(&mut cache, order, bin, lumi, 1.0, 1.0) + grid.convolve_subgrid(&mut cache, order, bin, lumi, 1.0, 1.0) + } + _ => unimplemented!(), + } } pub fn parse_integer_range(range: &str) -> Result> { diff --git a/pineappl_cli/src/plot.rs b/pineappl_cli/src/plot.rs index 2bc9b5a8..87acbe1b 100644 --- a/pineappl_cli/src/plot.rs +++ b/pineappl_cli/src/plot.rs @@ -13,7 +13,6 @@ use std::fmt::Write; use std::num::NonZeroUsize; use std::path::{Path, PathBuf}; use std::process::ExitCode; -use std::slice; use std::thread; /// Creates a matplotlib script plotting the contents of the grid. @@ -529,19 +528,13 @@ impl Subcommand for Opts { metadata = format_metadata(&vector), ); } else { - let (pdfset1, pdfset2) = self - .conv_funs - .iter() - .map(|fun| { - assert_eq!(fun.lhapdf_names.len(), 1); - if let Some(member) = fun.members[0] { - format!("{}/{member}", fun.lhapdf_names[0]) - } else { - fun.lhapdf_names[0].clone() - } - }) - .collect_tuple() - .unwrap(); + // TODO: enforce two arguments with clap + assert_eq!(self.conv_funs.len(), 2); + + let (set1, mut conv_funs1) = + helpers::create_conv_funs_for_set(&self.conv_funs[0], self.conv_fun_uncert_from)?; + let (set2, mut conv_funs2) = + helpers::create_conv_funs_for_set(&self.conv_funs[1], self.conv_fun_uncert_from)?; let (order, bin, channel) = self .subgrid_pull .iter() @@ -552,17 +545,15 @@ impl Subcommand for Opts { let cl = lhapdf::CL_1_SIGMA; let grid = helpers::read_grid(&self.input)?; - let (set1, member1) = helpers::create_pdfset(&pdfset1)?; - let (set2, member2) = helpers::create_pdfset(&pdfset2)?; - let mut pdfset1 = set1.mk_pdfs()?; - let mut pdfset2 = set2.mk_pdfs()?; + let member1 = self.conv_funs[0].members[self.conv_fun_uncert_from]; + let member2 = self.conv_funs[1].members[self.conv_fun_uncert_from]; - let values1: Vec<_> = pdfset1 + let values1: Vec<_> = conv_funs1 .par_iter_mut() - .map(|pdf| { + .map(|conv_funs| { let values = helpers::convolve( &grid, - slice::from_mut(pdf), + conv_funs, &[], &[bin], &[], @@ -574,12 +565,12 @@ impl Subcommand for Opts { values[0] }) .collect(); - let values2: Vec<_> = pdfset2 + let values2: Vec<_> = conv_funs2 .par_iter_mut() - .map(|pdf| { + .map(|conv_funs| { let values = helpers::convolve( &grid, - slice::from_mut(pdf), + conv_funs, &[], &[bin], &[], @@ -613,10 +604,26 @@ impl Subcommand for Opts { unc1.hypot(unc2) }; - let res1 = helpers::convolve_subgrid(&grid, &mut pdfset1[0], order, bin, channel, cfg) - .sum_axis(Axis(0)); - let res2 = helpers::convolve_subgrid(&grid, &mut pdfset2[0], order, bin, channel, cfg) - .sum_axis(Axis(0)); + // TODO: if no member is given, the zeroth is used, but we should show the averaged + // result of all members instead + let res1 = helpers::convolve_subgrid( + &grid, + &mut conv_funs1[member1.unwrap_or(0)], + order, + bin, + channel, + cfg, + ) + .sum_axis(Axis(0)); + let res2 = helpers::convolve_subgrid( + &grid, + &mut conv_funs2[member2.unwrap_or(0)], + order, + bin, + channel, + cfg, + ) + .sum_axis(Axis(0)); let subgrid = &grid.subgrids()[[order, bin, channel]]; //let q2 = subgrid.q2_grid();