diff --git a/examples/cpp/Makefile b/examples/cpp/Makefile index 4c15c37f..2b9612ab 100644 --- a/examples/cpp/Makefile +++ b/examples/cpp/Makefile @@ -16,6 +16,7 @@ PROGRAMS = \ get-subgrids \ set-subgrids \ evolve-grid \ + evolve-grid-pure-ff \ evolve-grid-identity \ deprecated \ display-channels-deprecated \ @@ -64,6 +65,9 @@ set-subgrids: set-subgrids.cpp evolve-grid: evolve-grid.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ +evolve-grid-pure-ff: evolve-grid-pure-ff.cpp + $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ + evolve-grid-identity: evolve-grid-identity.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ diff --git a/examples/cpp/evolve-grid-pure-ff.cpp b/examples/cpp/evolve-grid-pure-ff.cpp new file mode 100644 index 00000000..2db59301 --- /dev/null +++ b/examples/cpp/evolve-grid-pure-ff.cpp @@ -0,0 +1,288 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** @brief This struct can contain arbitrary parameters that need to be passed to Evolution + * Operator Callback (`generate_identity_eko`). + */ +struct OperatorParams { + std::vector conv_types; +}; + +std::vector unravel_index(std::size_t flat_index, const std::vector& shape) { + std::size_t ndim = shape.size(); + std::vector coords(ndim); + + for (int i = static_cast(ndim) - 1; i >= 0; --i) { + coords[i] = flat_index % shape[static_cast(i)]; + flat_index /= shape[static_cast(i)]; + } + + return coords; +} + +extern "C" void generate_identity_eko( + std::size_t op_index, + double /*fac1*/, + const int* /*pids_in*/, + const double* /*x_in*/, + const int* /*pids_out*/, + const double* /*x_out*/, + const std::size_t* eko_shape, + double* eko_buffer, + void* params_state +) { + OperatorParams* op_params = static_cast(params_state); + pineappl_conv_type conv_type = op_params->conv_types[op_index]; + assert(conv_type == PINEAPPL_CONV_TYPE_UNPOL_FF); + + std::vector shape(eko_shape, eko_shape + 4); + std::size_t flat_len = std::accumulate(shape.begin(), shape.end(), 1ULL, std::multiplies()); + + for (std::size_t i = 0; i != flat_len; i++) { + std::vector coords = unravel_index(i, shape); + double delta_ik = (coords[0] == coords[2]) ? 1.0 : 0.0; + double delta_jl = (coords[1] == coords[3]) ? 1.0 : 0.0; + eko_buffer[i] = delta_ik * delta_jl; + } +} + +void print_results(const std::vector& dxsec_grid, const std::vector& dxsec_fktable) { + const int idx_width = 6; + const int num_width = 15; + const int dif_width = 15; + + std::cout << std::setw(idx_width) << "Bin" + << std::setw(num_width) << "Grid" + << std::setw(num_width) << "FkTable" + << std::setw(dif_width) << "reldiff" << std::endl; + + std::cout << std::setw(idx_width) << std::string(idx_width - 2, '-') + << std::setw(num_width) << std::string(num_width - 2, '-') + << std::setw(num_width) << std::string(num_width - 2, '-') + << std::setw(dif_width) << std::string(dif_width - 2, '-') << std::endl; + + std::cout << std::scientific << std::setprecision(6); + for (std::size_t i = 0; i < dxsec_grid.size(); ++i) { + double reldiff = (dxsec_fktable[i] - dxsec_grid[i]) / dxsec_grid[i]; + std::cout << std::setw(idx_width) << i + << std::setw(num_width) << dxsec_grid[i] + << std::setw(num_width) << dxsec_fktable[i] + << std::setw(dif_width) << reldiff + << std::endl; + } +} + +pineappl_grid* create_and_fill_pure_ff_grid() { + std::size_t nb_convolutions = 1; + auto* channels = pineappl_channels_new(nb_convolutions); + + int32_t pids1[] = { 21 }; + double factors1[] = { 1.0 }; + pineappl_channels_add(channels, 1, pids1, factors1); + + // (alphas^0, alpha^0, logxir^0, logxif^0, logxia^0) + std::vector orders = { 0, 0, 0, 0, 0 }; + + std::vector bins = { 0.0, 0.5, 1.0 }; + + pineappl_pid_basis pid_basis = PINEAPPL_PID_BASIS_PDG; + pineappl_conv convs[] = { { PINEAPPL_CONV_TYPE_UNPOL_FF, 211 } }; // hadron PID: pi+ + + pineappl_kinematics scales = { PINEAPPL_KINEMATICS_SCALE, 0 }; + pineappl_kinematics z = { PINEAPPL_KINEMATICS_X, 0 }; + pineappl_kinematics kinematics[2] = { scales, z }; + + pineappl_reweight_meth scales_reweight = PINEAPPL_REWEIGHT_METH_NO_REWEIGHT; + pineappl_reweight_meth moment_reweight = PINEAPPL_REWEIGHT_METH_APPL_GRID_X; + pineappl_map scales_mapping = PINEAPPL_MAP_APPL_GRID_H0; + pineappl_map moment_mapping = PINEAPPL_MAP_APPL_GRID_F2; + pineappl_interp_meth interpolation_meth = PINEAPPL_INTERP_METH_LAGRANGE; + pineappl_interp interpolations[2] = { + { 1e2, 1e2 * (1.0 + 1e-12), 2, 1, scales_reweight, scales_mapping, interpolation_meth }, + { 2e-7, 1.0, 50, 3, moment_reweight, moment_mapping, interpolation_meth }, + }; + + pineappl_scale_func_form scale_mu = { PINEAPPL_SCALE_FUNC_FORM_SCALE, 0 }; + pineappl_scale_func_form no_scale_mu = { PINEAPPL_SCALE_FUNC_FORM_NO_SCALE, 0 }; + pineappl_scale_func_form mu_scales[3] = { + scale_mu, // ren + no_scale_mu, // fac + scale_mu, // frg + }; + + auto* grid = pineappl_grid_new2( + bins.size() - 1, + bins.data(), + orders.size() / 5, + orders.data(), + channels, + pid_basis, + convs, + 2, + interpolations, + kinematics, + mu_scales + ); + + pineappl_channels_delete(channels); + + std::size_t order = 0; + std::size_t channel = 0; + double q2 = 100.0; // GeV^2 + for (int i = 0; i != 10; ++i) { + double zi = 0.1 + 0.08 * i; + double obs = zi; // bin observable + double weight = 1.0 + 0.1 * i; + std::vector ntuples = { q2, zi }; + pineappl_grid_fill2(grid, order, obs, channel, ntuples.data(), weight); + } + + return grid; +} + +int main() { + LHAPDF::setVerbosity(0); + + std::string ffset = "MAPFF10NLOPIsum"; + auto ff = std::unique_ptr(LHAPDF::mkPDF(ffset, 0)); + + auto xfx = [](int32_t id, double x, double q2, void* ff) { + return static_cast(ff)->xfxQ2(id, x, q2); + }; + auto alphas = [](double q2, void* ff) { + return static_cast(ff)->alphasQ2(q2); + }; + + auto* grid = create_and_fill_pure_ff_grid(); + + std::vector evinfo_shape(5); + pineappl_grid_evolve_info_shape(grid, nullptr, evinfo_shape.data()); + + std::vector fac1(evinfo_shape[0]); + std::vector frg1(evinfo_shape[1]); + std::vector pids_in(evinfo_shape[2]); + std::vector x_in(evinfo_shape[3]); + std::vector ren1(evinfo_shape[4]); + + pineappl_grid_evolve_info(grid, nullptr, fac1.data(), frg1.data(), pids_in.data(), x_in.data(), ren1.data()); + + const std::vector& proc1_all = fac1.empty() ? frg1 : fac1; + assert(!proc1_all.empty()); + + std::size_t n_convs = pineappl_grid_convolutions_len(grid); + assert(n_convs == 1); + std::vector conv_types(n_convs); + pineappl_grid_conv_types(grid, conv_types.data()); + + std::vector unique_convs; + for (std::size_t i = 0; i != n_convs; i++) { + pineappl_conv_type conv = conv_types[i]; + if (std::find(unique_convs.begin(), unique_convs.end(), conv) == unique_convs.end()) { + unique_convs.push_back(conv); + } + } + assert(unique_convs.size() == 1); + assert(unique_convs[0] == PINEAPPL_CONV_TYPE_UNPOL_FF); + + pineappl_pid_basis pid_basis = pineappl_grid_pid_basis(grid); + const double fac0 = proc1_all.front(); + std::vector opinfo_slices(unique_convs.size() * proc1_all.size()); + for (std::size_t i = 0; i != unique_convs.size(); i++) { + for (std::size_t j = 0; j != proc1_all.size(); j++) { + pineappl_operator_info opinfo = { + fac0, // fac0: common starting scale + proc1_all[j], // fac1: process scale (here: fragmentation scale) + pid_basis, + unique_convs[i], + }; + opinfo_slices[i * proc1_all.size() + j] = opinfo; + } + } + + std::vector alphas_table; + alphas_table.reserve(ren1.size()); + for (double q2 : ren1) { + alphas_table.push_back(alphas(q2, ff.get())); + } + + OperatorParams op_params; + op_params.conv_types = unique_convs; + + std::vector xi = { 1.0, 1.0, 1.0 }; + std::vector pids_out = pids_in; + std::vector tensor_shape = { pids_in.size(), x_in.size(), pids_out.size(), x_in.size() }; + + pineappl_grid* fktable = pineappl_grid_evolve( + grid, + unique_convs.size(), + generate_identity_eko, + opinfo_slices.data(), + pids_in.data(), + x_in.data(), + pids_out.data(), + x_in.data(), + tensor_shape.data(), + &op_params, + nullptr, + xi.data(), + ren1.data(), + alphas_table.data() + ); + + // Compare grid vs evolved FK by convolution + std::size_t bins = pineappl_grid_bin_count(grid); + std::vector mu_scales = { 1.0, 1.0, 1.0 }; + + std::vector ffs = { ff.get() }; + void** ff_states = reinterpret_cast(ffs.data()); + + std::vector dxsec_grid(bins); + pineappl_grid_convolve( + grid, + xfx, + alphas, + ff_states, + ff.get(), + nullptr, + nullptr, + nullptr, + 1, + mu_scales.data(), + dxsec_grid.data() + ); + + std::vector dxsec_fktable(bins); + auto as_one = [](double /*q2*/, void* /*state*/) { return 1.0; }; + pineappl_grid_convolve( + fktable, + xfx, + as_one, + ff_states, + nullptr, + nullptr, + nullptr, + nullptr, + 1, + mu_scales.data(), + dxsec_fktable.data() + ); + + print_results(dxsec_grid, dxsec_fktable); + + pineappl_grid_write(grid, "pure-ff-grid.pineappl.lz4"); + pineappl_grid_write(fktable, "pure-ff-fktable.pineappl.lz4"); + + pineappl_grid_delete(grid); + pineappl_grid_delete(fktable); +} diff --git a/examples/cpp/evolve-grid-pure-ff.output b/examples/cpp/evolve-grid-pure-ff.output new file mode 100644 index 00000000..b43054a7 --- /dev/null +++ b/examples/cpp/evolve-grid-pure-ff.output @@ -0,0 +1,4 @@ + Bin Grid FkTable reldiff + ---- ------------- ------------- ------------- + 0 5.525872e+01 5.525872e+01 0.000000e+00 + 1 2.298532e+00 2.298532e+00 1.932056e-16 diff --git a/pineappl/src/evolution.rs b/pineappl/src/evolution.rs index 389a4bdc..9452e522 100644 --- a/pineappl/src/evolution.rs +++ b/pineappl/src/evolution.rs @@ -1,6 +1,6 @@ //! Supporting classes and functions for [`Grid::evolve`]. -use super::boc::{Channel, Kinematics, Order}; +use super::boc::{Channel, Kinematics, Order, ScaleFuncForm}; use super::convolutions::ConvType; use super::error::{Error, Result}; use super::grid::Grid; @@ -317,14 +317,30 @@ fn ndarray_from_subgrid_orders_slice( .collect(); let rens = grid.scales().ren.calc(&node_values, grid.kinematics()); - let facs = grid.scales().fac.calc(&node_values, grid.kinematics()); + + // Pure FF grids (e.g. SIA) contains `fac = NoScale` and thus use `frg` as fact. scale. + let use_frg_for_process = matches!(grid.scales().fac, ScaleFuncForm::NoScale); + let proc_scales = if use_frg_for_process { + grid.scales().frg.calc(&node_values, grid.kinematics()) + } else { + grid.scales().fac.calc(&node_values, grid.kinematics()) + }; + let xi_proc_sq = if use_frg_for_process { + xia * xia + } else { + xif * xif + }; for (indices, value) in subgrid.indexed_iter() { - // TODO: implement evolution for non-zero fragmentation scales let ren = rens[grid.scales().ren.idx(&indices, &scale_dims)]; - let fac = facs[grid.scales().fac.idx(&indices, &scale_dims)]; + let proc_idx = if use_frg_for_process { + grid.scales().frg.idx(&indices, &scale_dims) + } else { + grid.scales().fac.idx(&indices, &scale_dims) + }; + let proc_mu2 = proc_scales[proc_idx]; - if !subgrid::node_value_eq(xif * xif * fac, fac1) { + if !subgrid::node_value_eq(xi_proc_sq * proc_mu2, fac1) { continue; } diff --git a/pineappl/src/fk_table.rs b/pineappl/src/fk_table.rs index 9ea820a0..58a2e18e 100644 --- a/pineappl/src/fk_table.rs +++ b/pineappl/src/fk_table.rs @@ -298,7 +298,7 @@ impl TryFrom for FkTable { type Error = Error; fn try_from(grid: Grid) -> Result { - let mut muf2 = -1.0; + let mut mu2_ref = -1.0; if grid.orders() != [Order { @@ -317,17 +317,24 @@ impl TryFrom for FkTable { continue; } - let [fac] = grid - .scales() - .fac - .calc(&subgrid.node_values(), grid.kinematics())[..] - else { - return Err(Error::General("multiple scales detected".to_owned())); + let node_values = subgrid.node_values(); + let kinematics = grid.kinematics(); + let fac_s = grid.scales().fac.calc(&node_values, kinematics); + let frg_s = grid.scales().frg.calc(&node_values, kinematics); + + let mu2 = match (fac_s.len(), frg_s.len()) { + (0, 0) => { + return Err(Error::General("No fac or frg scale in subgrid".to_owned())); + } + (1, 0) => fac_s[0], + (0, 1) => frg_s[0], + (1, 1) if subgrid::node_value_eq(fac_s[0], frg_s[0]) => fac_s[0], + _ => return Err(Error::General("multiple scales detected".to_owned())), }; - if muf2 < 0.0 { - muf2 = fac; - } else if !subgrid::node_value_eq(muf2, fac) { + if mu2_ref < 0.0 { + mu2_ref = mu2; + } else if !subgrid::node_value_eq(mu2_ref, mu2) { return Err(Error::General("multiple scales detected".to_owned())); } } diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs index 2719943c..0b96db3d 100644 --- a/pineappl_capi/src/lib.rs +++ b/pineappl_capi/src/lib.rs @@ -2336,7 +2336,7 @@ pub unsafe extern "C" fn pineappl_grid_evolve_info( /// Type alias for the operator callback pub type OperatorCallback = unsafe extern "C" fn( usize, // index which selects Evolution parameters - f64, // fac1 + f64, // squared process scale (fac or frg) *const i32, // `pids_in` *const f64, // `x_in` *const i32, // `pids_out` @@ -2412,9 +2412,11 @@ pub unsafe extern "C" fn pineappl_grid_evolve( let alphas = unsafe { slice::from_raw_parts(alphas, ren1_len) }; let xi = unsafe { slice::from_raw_parts(xi, 3) }; + let n_scale_slices = evolve_info.fac1.len().max(evolve_info.frg1.len()).max(1); + let operator_info = unsafe { - slice::from_raw_parts(operator_info, nb_slices * evolve_info.fac1.len()) - .chunks_exact(evolve_info.fac1.len()) + slice::from_raw_parts(operator_info, nb_slices * n_scale_slices) + .chunks_exact(n_scale_slices) }; let shape: (usize, usize, usize, usize) = <[usize; 4]>::try_from(eko_shape) diff --git a/pineappl_py/src/evolution.rs b/pineappl_py/src/evolution.rs index f23071c7..6ddf44cf 100644 --- a/pineappl_py/src/evolution.rs +++ b/pineappl_py/src/evolution.rs @@ -94,9 +94,16 @@ impl PyEvolveInfo { } /// Squared factorization scales of the `Grid`. + /// + /// Note: for pure time-like (FF) grids, PineAPPL stores only fragmentation scales. In that + /// case, this getter returns `frg1` as a backwards-compatible fallback. #[getter] fn fac1<'py>(&self, py: Python<'py>) -> Bound<'py, PyArray1> { - self.evolve_info.fac1.clone().into_pyarray(py) + if self.evolve_info.fac1.is_empty() && !self.evolve_info.frg1.is_empty() { + self.evolve_info.frg1.clone().into_pyarray(py) + } else { + self.evolve_info.fac1.clone().into_pyarray(py) + } } /// Squared fragmentation scales of the `Grid`.