diff --git a/pineappl_cli/src/import/fastnlo.rs b/pineappl_cli/src/import/fastnlo.rs index 337fda4b..42a3bd55 100644 --- a/pineappl_cli/src/import/fastnlo.rs +++ b/pineappl_cli/src/import/fastnlo.rs @@ -1,7 +1,7 @@ use anyhow::Result; -// use itertools::Itertools; -// use ndarray::s; use float_cmp::approx_eq; +use itertools::Itertools; +use ndarray::s; use pineappl::bin::BinRemapper; use pineappl::boc::{Channel, Kinematics, Order, ScaleFuncForm, Scales}; use pineappl::convolutions::{Conv, ConvType}; @@ -16,6 +16,7 @@ use pineappl_fastnlo::ffi::{ fastNLOPDFLinearCombinations, EScaleFunctionalForm, }; use std::f64::consts::TAU; +use std::iter; use std::mem; fn pid_to_pdg_id(pid: i32) -> i32 { @@ -292,205 +293,228 @@ fn convert_coeff_add_fix( grid } +fn convert_scale_functional_form(ff: EScaleFunctionalForm) -> ScaleFuncForm { + match ff { + EScaleFunctionalForm::kScale1 => ScaleFuncForm::Scale(0), + EScaleFunctionalForm::kScale2 => ScaleFuncForm::Scale(1), + EScaleFunctionalForm::kQuadraticSum => ScaleFuncForm::QuadraticSum(0, 1), + EScaleFunctionalForm::kQuadraticMean => todo!(), + EScaleFunctionalForm::kQuadraticSumOver4 => todo!(), + EScaleFunctionalForm::kLinearMean => todo!(), + EScaleFunctionalForm::kLinearSum => todo!(), + EScaleFunctionalForm::kScaleMax => todo!(), + EScaleFunctionalForm::kScaleMin => todo!(), + EScaleFunctionalForm::kProd => todo!(), + EScaleFunctionalForm::kS2plusS1half => todo!(), + EScaleFunctionalForm::kPow4Sum => todo!(), + EScaleFunctionalForm::kWgtAvg => todo!(), + EScaleFunctionalForm::kS2plusS1fourth => todo!(), + EScaleFunctionalForm::kExpProd2 => todo!(), + _ => unimplemented!(), + } +} + fn convert_coeff_add_flex( - _table: &fastNLOCoeffAddFlex, - _comb: &fastNLOPDFLinearCombinations, - _mur_ff: EScaleFunctionalForm, - _muf_ff: EScaleFunctionalForm, - _bins: usize, - _alpha: u8, - _ipub_units: i32, + table: &fastNLOCoeffAddFlex, + comb: &fastNLOPDFLinearCombinations, + mur_ff: EScaleFunctionalForm, + muf_ff: EScaleFunctionalForm, + bins: usize, + alpha: u8, + ipub_units: i32, ) -> Grid { - todo!() - - // let table_as_add_base = ffi::downcast_coeff_add_flex_to_base(table); - - // let alphas = table_as_add_base.GetNpow().try_into().unwrap(); - // let orders: Vec<_> = [ - // Order::new(alphas, alpha, 0, 0, 0), - // Order::new(alphas, alpha, 1, 0, 0), - // Order::new(alphas, alpha, 0, 1, 0), - // Order::new(alphas, alpha, 2, 0, 0), - // Order::new(alphas, alpha, 0, 2, 0), - // Order::new(alphas, alpha, 1, 1, 0), - // ] - // .into_iter() - // .take(match table.GetNScaleDep() { - // 0..=4 => 1, - // 5 => 3, - // 6 => 4, - // 7 => 6, - // _ => unimplemented!(), - // }) - // .collect(); - // let orders_len = orders.len(); - - // let npdf = table_as_add_base.GetNPDF(); - // assert!(npdf <= 2); - - // let convolutions = (0..2) - // .map(|index| { - // if index < npdf { - // Convolution::UnpolPDF(table.GetPDFPDG(index)) - // } else { - // Convolution::None - // } - // }) - // .collect(); - - // let mut grid = Grid::new( - // PidBasis::Pdg, - // reconstruct_channels(table_as_add_base, comb, dis_pid), - // orders, - // (0..=bins) - // .map(|limit| u16::try_from(limit).unwrap().into()) - // .collect(), - // convolutions, - // // TODO: read out interpolation parameters from fastNLO - // 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::Scale(1), - // Kinematics::X1, - // Kinematics::X2, - // ], - // Scales { - // ren: todo!(), - // fac: todo!(), - // frg: ScaleFuncForm::NoScale, - // }, - // ); - - // let rescale = 0.1_f64.powi(table.GetIXsectUnits() - ipub_units); - - // for obs in 0..bins { - // let scale_nodes1 = ffi::GetScaleNodes1(table, obs.try_into().unwrap()); - // let scale_nodes2 = ffi::GetScaleNodes2(table, obs.try_into().unwrap()); - // let x1_values = ffi::GetXNodes1(table_as_add_base, obs.try_into().unwrap()); - // let x2_values = if npdf > 1 { - // ffi::GetXNodes2(table_as_add_base, obs.try_into().unwrap()) - // } else { - // vec![1.0] - // }; - - // let mu2_values: Vec<_> = scale_nodes1 - // .iter() - // .cartesian_product(scale_nodes2.iter()) - // .map(|(&s1, &s2)| Mu2 { - // ren: mur_ff.compute_scale(s1, s2), - // fac: muf_ff.compute_scale(s1, s2), - // frg: -1.0, - // }) - // .collect(); - // let nx = ffi::GetNx(table, obs); - - // for subproc in 0..table_as_add_base.GetNSubproc() { - // let factor = rescale / table_as_add_base.GetNevt(obs.try_into().unwrap(), subproc); - // let mut arrays = - // vec![ - // PackedArray::new(vec![mu2_values.len(), x1_values.len(), x2_values.len()]); - // orders_len - // ]; - - // for (mu2_slice, (is1, is2)) in (0..scale_nodes1.len()) - // .cartesian_product(0..scale_nodes2.len()) - // .enumerate() - // { - // let logmur2 = mu2_values[mu2_slice].ren.ln(); - // let logmuf2 = mu2_values[mu2_slice].fac.ln(); - // let logs00 = [ - // logmur2, - // logmuf2, - // logmur2 * logmur2, - // logmuf2 * logmuf2, - // logmur2 * logmuf2, - // ]; - // let logs10 = [2.0 * logmur2, 0.0, logmuf2]; - // let logs01 = [0.0, 2.0 * logmuf2, logmur2]; - - // for ix in 0..nx { - // // TODO: is this always correct? Isn't there a member function for it? - // let ix1 = ix % x1_values.len(); - // let ix2 = ix / x1_values.len(); - // let mut values = [0.0; 6]; - - // for (index, value) in values.iter_mut().enumerate().take(orders_len) { - // *value = ffi::GetSigmaTilde(table, index, obs, ix, is1, is2, subproc); - // } - - // values[0] += values[1..] - // .iter() - // .zip(logs00.iter()) - // .map(|(value, log)| value * log) - // .sum::(); - // values[1] += values[3..] - // .iter() - // .zip(logs10.iter()) - // .map(|(value, log)| value * log) - // .sum::(); - // values[2] += values[3..] - // .iter() - // .zip(logs01.iter()) - // .map(|(value, log)| value * log) - // .sum::(); - - // for (value, array) in values - // .iter() - // .copied() - // .zip(arrays.iter_mut()) - // .filter(|(value, _)| *value != 0.0) - // { - // array[[mu2_slice, ix1, ix2]] = - // value * factor * x1_values[ix1] * x2_values[ix2]; - // } - // } - // } - - // for (subgrid, array) in grid - // .subgrids_mut() - // .slice_mut(s![.., obs, usize::try_from(subproc).unwrap()]) - // .iter_mut() - // .zip(arrays.into_iter()) - // { - // if array.is_empty() { - // continue; - // } - - // *subgrid = ImportSubgridV1::new(array, todo!()).into(); - // } - // } - // } - - // grid + let table_as_add_base = ffi::downcast_coeff_add_flex_to_base(table); + + let alphas = table_as_add_base.GetNpow().try_into().unwrap(); + let orders: Vec<_> = [ + Order::new(alphas, alpha, 0, 0, 0), + Order::new(alphas, alpha, 1, 0, 0), + Order::new(alphas, alpha, 0, 1, 0), + Order::new(alphas, alpha, 2, 0, 0), + Order::new(alphas, alpha, 0, 2, 0), + Order::new(alphas, alpha, 1, 1, 0), + ] + .into_iter() + .take(match table.GetNScaleDep() { + 0..=4 => 1, + 5 => 3, + 6 => 4, + 7 => 6, + _ => unimplemented!(), + }) + .collect(); + let orders_len = orders.len(); + + let npdf = table_as_add_base.GetNPDF(); + assert!(npdf <= 2); + + let convolutions = (0..npdf) + .map(|index| Conv::new(ConvType::UnpolPDF, table.GetPDFPDG(index))) + .collect(); + + let npdf: usize = npdf.try_into().unwrap(); + + let mut grid = Grid::new( + PidBasis::Pdg, + reconstruct_channels(table_as_add_base, comb), + orders, + (0..=bins) + .map(|limit| u16::try_from(limit).unwrap().into()) + .collect(), + convolutions, + // TODO: read out interpolation parameters from fastNLO + iter::repeat(Interp::new( + 1e2, + 1e8, + 40, + 3, + ReweightMeth::NoReweight, + Map::ApplGridH0, + InterpMeth::Lagrange, + )) + .take(2) + .chain( + iter::repeat(Interp::new( + 2e-7, + 1.0, + 50, + 3, + ReweightMeth::ApplGridX, + Map::ApplGridF2, + InterpMeth::Lagrange, + )) + .take(npdf), + ) + .collect(), + [Kinematics::Scale(0), Kinematics::Scale(1)] + .into_iter() + .chain((0..npdf).map(Kinematics::X)) + .collect(), + Scales { + ren: convert_scale_functional_form(mur_ff), + fac: convert_scale_functional_form(muf_ff), + // TODO: does fastNLO not support fragmentation scales? + frg: ScaleFuncForm::NoScale, + }, + ); + + let rescale = 0.1_f64.powi(table.GetIXsectUnits() - ipub_units); + + for obs in 0..bins { + let scale_nodes1 = ffi::GetScaleNodes1(table, obs.try_into().unwrap()); + let scale_nodes2 = ffi::GetScaleNodes2(table, obs.try_into().unwrap()); + let x1_values = ffi::GetXNodes1(table_as_add_base, obs.try_into().unwrap()); + let x2_values = if npdf > 1 { + ffi::GetXNodes2(table_as_add_base, obs.try_into().unwrap()) + } else { + vec![1.0] + }; + + let mut dim = vec![scale_nodes1.len(), scale_nodes2.len(), x1_values.len()]; + if npdf > 1 { + dim.push(x2_values.len()); + } + + let nx = ffi::GetNx(table, obs); + + for subproc in 0..table_as_add_base.GetNSubproc() { + let factor = rescale / table_as_add_base.GetNevt(obs.try_into().unwrap(), subproc); + + let mut arrays = vec![PackedArray::new(dim.clone()); orders_len]; + + for (is1, is2) in (0..scale_nodes1.len()).cartesian_product(0..scale_nodes2.len()) { + let logmur2 = mur_ff + .compute_scale(scale_nodes1[is1], scale_nodes2[is2]) + .ln(); + let logmuf2 = muf_ff + .compute_scale(scale_nodes1[is1], scale_nodes2[is2]) + .ln(); + let logs00 = [ + logmur2, + logmuf2, + logmur2 * logmur2, + logmuf2 * logmuf2, + logmur2 * logmuf2, + ]; + let logs10 = [2.0 * logmur2, 0.0, logmuf2]; + let logs01 = [0.0, 2.0 * logmuf2, logmur2]; + + for ix in 0..nx { + // TODO: is this always correct? Isn't there a member function for it? + let ix1 = ix % x1_values.len(); + let ix2 = ix / x1_values.len(); + let mut values = [0.0; 6]; + + for (index, value) in values.iter_mut().enumerate().take(orders_len) { + *value = ffi::GetSigmaTilde(table, index, obs, ix, is1, is2, subproc); + } + + values[0] += values[1..] + .iter() + .zip(logs00.iter()) + .map(|(value, log)| value * log) + .sum::(); + values[1] += values[3..] + .iter() + .zip(logs10.iter()) + .map(|(value, log)| value * log) + .sum::(); + values[2] += values[3..] + .iter() + .zip(logs01.iter()) + .map(|(value, log)| value * log) + .sum::(); + + for (value, array) in values + .iter() + .copied() + .zip(arrays.iter_mut()) + .filter(|(value, _)| *value != 0.0) + { + if npdf == 1 { + array[[is1, is2, ix1]] = value * factor * x1_values[ix1]; + } else if npdf == 2 { + assert_eq!(is2, 0); + + array[[is1, is2, ix1, ix2]] = + value * factor * x1_values[ix1] * x2_values[ix2]; + } + } + } + } + + for (subgrid, array) in grid + .subgrids_mut() + .slice_mut(s![.., obs, usize::try_from(subproc).unwrap()]) + .iter_mut() + .zip(arrays.into_iter()) + { + if array.is_empty() { + continue; + } + + let node_values = if npdf == 1 { + vec![ + scale_nodes1.iter().map(|s| s * s).collect(), + scale_nodes2.iter().map(|s| s * s).collect(), + x1_values.clone(), + ] + } else { + vec![ + scale_nodes1.iter().map(|s| s * s).collect(), + // scale_nodes2.iter().map(|s| s * s).collect(), + vec![100000.0], + x1_values.clone(), + x2_values.clone(), + ] + }; + + *subgrid = ImportSubgridV1::new(array, node_values).into(); + } + } + } + + grid } pub fn convert_fastnlo_table(file: &fastNLOLHAPDF, alpha: u8) -> Result { diff --git a/pineappl_cli/tests/import.rs b/pineappl_cli/tests/import.rs index 2d0a169a..286fa72b 100644 --- a/pineappl_cli/tests/import.rs +++ b/pineappl_cli/tests/import.rs @@ -285,26 +285,26 @@ const IMPORT_DOUBLE_HADRONIC_FASTNLO_STR: &str = "b PineAPPL fastNLO rel. diff svmaxreldiff --+------------+------------+--------------+------------- 0 9.6382069e5 9.6382069e5 4.4408921e-16 8.3266727e-15 -1 3.7342594e5 3.7342594e5 1.7985613e-14 1.9095836e-14 -2 1.4195038e5 1.4195038e5 -1.0880186e-14 2.2648550e-14 -3 5.7043791e4 5.7043791e4 4.2188475e-15 7.9936058e-15 +1 3.7342594e5 3.7342594e5 1.7985613e-14 1.8651747e-14 +2 1.4195038e5 1.4195038e5 -1.0880186e-14 2.2870594e-14 +3 5.7043791e4 5.7043791e4 4.2188475e-15 7.7715612e-15 4 2.3327746e4 2.3327746e4 8.4376950e-15 1.2101431e-14 -5 1.0495603e4 1.0495603e4 1.3100632e-14 1.7985613e-14 -6 4.8153483e3 4.8153483e3 -1.6098234e-14 2.9753977e-14 -7 2.2957587e3 2.2957587e3 4.8849813e-15 3.0642155e-14 +5 1.0495603e4 1.0495603e4 1.3100632e-14 1.7874591e-14 +6 4.8153483e3 4.8153483e3 -1.6098234e-14 2.9531932e-14 +7 2.2957587e3 2.2957587e3 4.6629367e-15 3.0198066e-14 8 1.1142545e3 1.1142545e3 -2.4424907e-15 1.5765167e-14 -9 5.3699925e2 5.3699925e2 -6.5503158e-15 1.8429702e-14 +9 5.3699925e2 5.3699925e2 -6.7723605e-15 1.8429702e-14 10 2.5460314e2 2.5460314e2 -7.6605389e-15 1.3544721e-14 -11 1.1847638e2 1.1847638e2 1.0658141e-14 1.2656542e-14 -12 5.7567355e1 5.7567355e1 -2.9976022e-15 9.2148511e-15 -13 2.7189719e1 2.7189719e1 1.1102230e-15 1.5543122e-14 -14 1.2791922e1 1.2791922e1 -6.9944051e-15 1.2656542e-14 -15 5.8346996e0 5.8346996e0 2.8865799e-15 1.4988011e-14 -16 2.6521590e0 2.6521590e0 7.3274720e-15 1.4765966e-14 -17 1.1726035e0 1.1726035e0 1.3100632e-14 1.4432899e-14 -18 4.8823596e-1 4.8823596e-1 8.6597396e-15 1.3655743e-14 -19 1.9564964e-1 1.9564964e-1 -4.4408921e-15 1.1102230e-14 -20 2.0326950e-2 2.0326950e-2 6.6613381e-15 1.3211654e-14 +11 1.1847638e2 1.1847638e2 1.0880186e-14 1.2989609e-14 +12 5.7567355e1 5.7567355e1 -2.8865799e-15 9.2148511e-15 +13 2.7189719e1 2.7189719e1 1.3322676e-15 1.5543122e-14 +14 1.2791922e1 1.2791922e1 -6.9944051e-15 1.2878587e-14 +15 5.8346996e0 5.8346996e0 2.8865799e-15 1.4876989e-14 +16 2.6521590e0 2.6521590e0 7.1054274e-15 1.4765966e-14 +17 1.1726035e0 1.1726035e0 1.3100632e-14 1.3988810e-14 +18 4.8823596e-1 4.8823596e-1 8.6597396e-15 1.3433699e-14 +19 1.9564964e-1 1.9564964e-1 -4.6629367e-15 1.1102230e-14 +20 2.0326950e-2 2.0326950e-2 6.6613381e-15 1.2767565e-14 "; #[cfg(feature = "fastnlo")] @@ -367,7 +367,6 @@ fn import_flex_grid() { } #[test] -#[ignore] #[cfg(feature = "fastnlo")] fn import_flex_grid_scale_1() { let output = NamedTempFile::new("converted2.pineappl.lz4").unwrap(); @@ -388,7 +387,6 @@ fn import_flex_grid_scale_1() { } #[test] -#[ignore] #[cfg(feature = "fastnlo")] fn import_flex_grid_scale_2() { let output = NamedTempFile::new("converted2.pineappl.lz4").unwrap(); @@ -1289,7 +1287,6 @@ fn import_dis_applgrid() { } #[test] -#[ignore] #[cfg(feature = "fastnlo")] fn import_double_hadronic_fastnlo() { let output = NamedTempFile::new("converted9.pineappl.lz4").unwrap();