From 0a66c20bc4ff1007fd329d6926aa568ec7cfaf24 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Thu, 10 Oct 2024 15:12:17 +0200 Subject: [PATCH] Fix failing doctest and make DIS FK-table subgrids 2-dimensional --- pineappl_cli/src/import/fktable.rs | 144 +++++++++++++++++++---------- pineappl_cli/tests/import.rs | 3 +- 2 files changed, 94 insertions(+), 53 deletions(-) diff --git a/pineappl_cli/src/import/fktable.rs b/pineappl_cli/src/import/fktable.rs index 46cf6b98..8b58bc6b 100644 --- a/pineappl_cli/src/import/fktable.rs +++ b/pineappl_cli/src/import/fktable.rs @@ -98,6 +98,12 @@ fn read_fktable(reader: impl BufRead) -> Result { .collect() }; + let convolutions = if hadronic { + vec![Convolution::UnpolPDF(2212); 2] + } else { + vec![Convolution::UnpolPDF(2212)] + }; + // construct `Grid` let fktable = Grid::new( PidBasis::Evol, @@ -105,43 +111,65 @@ fn read_fktable(reader: impl BufRead) -> Result { vec![Order::new(0, 0, 0, 0, 0)], (0..=ndata).map(Into::into).collect(), // legacy FK-tables only support unpolarized proton PDFs + convolutions.clone(), + // TODO: what are sensible parameters for FK-tables? if hadronic { - vec![Convolution::UnpolPDF(2212); 2] + vec![ + Interp::new( + 1e2, + 1e8, + 40, + 3, + ReweightMeth::NoReweight, + Map::ApplGridH0, + InterpMeth::Lagrange, + ), + Interp::new( + 2e-7, + 1.0, + 50, + 3, + ReweightMeth::ApplGridX, + Map::ApplGridF2, + InterpMeth::Lagrange, + ), + Interp::new( + 2e-7, + 1.0, + 50, + 3, + ReweightMeth::ApplGridX, + Map::ApplGridF2, + InterpMeth::Lagrange, + ), + ] } else { - vec![Convolution::UnpolPDF(2212)] + vec![ + Interp::new( + 1e2, + 1e8, + 40, + 3, + ReweightMeth::NoReweight, + Map::ApplGridH0, + InterpMeth::Lagrange, + ), + Interp::new( + 2e-7, + 1.0, + 50, + 3, + ReweightMeth::ApplGridX, + Map::ApplGridF2, + InterpMeth::Lagrange, + ), + ] + }, + if hadronic { + vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2] + } else { + vec![Kinematics::Scale(0), Kinematics::X1] }, - // TODO: what are sensible parameters for FK-tables? - vec![ - Interp::new( - 1e2, - 1e8, - 40, - 3, - ReweightMeth::NoReweight, - Map::ApplGridH0, - InterpMeth::Lagrange, - ), - Interp::new( - 2e-7, - 1.0, - 50, - 3, - ReweightMeth::ApplGridX, - Map::ApplGridF2, - InterpMeth::Lagrange, - ), - Interp::new( - 2e-7, - 1.0, - 50, - 3, - ReweightMeth::ApplGridX, - Map::ApplGridF2, - InterpMeth::Lagrange, - ), - ], - // TODO: change kinematics for DIS - vec![Kinematics::Scale(0), Kinematics::X1, Kinematics::X2], // TODO: is this correct? Scales { ren: ScaleFuncForm::NoScale, @@ -152,9 +180,13 @@ fn read_fktable(reader: impl BufRead) -> Result { grid = Some(fktable); - arrays = iter::repeat(PackedArray::new(vec![1, nx1, nx2])) - .take(flavor_mask.iter().filter(|&&value| value).count()) - .collect(); + arrays = iter::repeat(PackedArray::new(if hadronic { + vec![1, nx1, nx2] + } else { + vec![1, nx1] + })) + .take(flavor_mask.iter().filter(|&&value| value).count()) + .collect(); } _ => match section { FkTableSection::GridInfo => { @@ -216,22 +248,29 @@ fn read_fktable(reader: impl BufRead) -> Result { { *subgrid = PackedQ1X2SubgridV1::new( array, - vec![ - NodeValues::UseThese(vec![q0 * q0]), - NodeValues::UseThese(x_grid.clone()), - NodeValues::UseThese(if hadronic { - x_grid.clone() - } else { - vec![1.0] - }), - ], + if hadronic { + vec![ + NodeValues::UseThese(vec![q0 * q0]), + NodeValues::UseThese(x_grid.clone()), + NodeValues::UseThese(x_grid.clone()), + ] + } else { + vec![ + NodeValues::UseThese(vec![q0 * q0]), + NodeValues::UseThese(x_grid.clone()), + ] + }, ) .into(); } - arrays = iter::repeat(PackedArray::new(vec![1, nx1, nx2])) - .take(flavor_mask.iter().filter(|&&value| value).count()) - .collect(); + arrays = iter::repeat(PackedArray::new(if hadronic { + vec![1, nx1, nx2] + } else { + vec![1, nx1] + })) + .take(flavor_mask.iter().filter(|&&value| value).count()) + .collect(); last_bin = bin; } @@ -257,8 +296,11 @@ fn read_fktable(reader: impl BufRead) -> Result { .zip(grid_values.iter()) .filter(|(_, value)| **value != 0.0) { - array[[0, x1, x2]] = - x_grid[x1] * if hadronic { x_grid[x2] } else { 1.0 } * value; + if hadronic { + array[[0, x1, x2]] = x_grid[x1] * x_grid[x2] * value; + } else { + array[[0, x1]] = x_grid[x1] * value; + } } } _ => {} diff --git a/pineappl_cli/tests/import.rs b/pineappl_cli/tests/import.rs index 83a8d23f..25b0d1a0 100644 --- a/pineappl_cli/tests/import.rs +++ b/pineappl_cli/tests/import.rs @@ -822,13 +822,12 @@ fn import_dis_fktable() { 0.9126417795942889, 0.9416284084927907, 0.9707498946430192, - 1.0 ] ); let table: Array3 = fk_table.table().into_dimensionality().unwrap(); - assert_eq!(table.dim(), (20, 9, 101)); + assert_eq!(table.dim(), (20, 9, 100)); assert_eq!( table .indexed_iter()