diff --git a/CHANGELOG.md b/CHANGELOG.md index 2a53f6af..98e2fe69 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added + +- added new subcommand `import`, which at this stage converts fastNLO tables to + PineAPPL grids. This program was previously an example and written in C++ and + now has been removed. + ### Changed - when running `pineappl convolute ... -s 1` the scale-variation columns are no diff --git a/Cargo.toml b/Cargo.toml index 31f2fcfb..a486f1b9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,6 +3,7 @@ members = [ "pineappl", "pineappl_capi", "pineappl_cli", + "pineappl_fastnlo", ] exclude = ["pineappl_py"] diff --git a/examples/fnlo2pine/fnlo2pine.cpp b/examples/fnlo2pine/fnlo2pine.cpp deleted file mode 100644 index 0726753a..00000000 --- a/examples/fnlo2pine/fnlo2pine.cpp +++ /dev/null @@ -1,660 +0,0 @@ -#include -#include -#include -#include - -#include -#include - -std::unique_ptr pdf; - -extern "C" double xfx(int32_t pdg_id, double x, double q2, void* state) -{ - return static_cast (state)->xfxQ2(pdg_id, x, q2); -} - -extern "C" double xfx2(int32_t pdg_id, double x, double, void*) -{ - assert( pdg_id == 11 ); - return x; -} - -extern "C" double alphas(double q2, void* state) -{ - return static_cast (state)->alphasQ2(q2); -} - -int32_t convert_to_pdg_id(int id) -{ - if ((id >= -6) && (id <= 6)) - { - return (id == 0) ? 21 : id; - } - else - { - assert( false ); - } -} - -void create_lumi( - fastNLOCoeffAddBase* table, - fastNLOPDFLinearCombinations const& comb, - pineappl_lumi* lumi -) { - auto const& pdf = table->GetPDFCoeff(); - - // TODO: set this to the right value if there's only one PDF - int32_t const lepton_id = (table->GetNPDF() == 2) ? 0 : 11; - - // if there's a (non-empty) PDF coefficient vector reconstruct the luminosity function; the - // advantage is that we preserve the order of the lumi entries in the PineAPPL grid - for (auto const& pdf_entries : pdf) - { - std::vector combinations; - std::vector factors; - - for (auto const& entry : pdf_entries) - { - combinations.push_back(convert_to_pdg_id(entry.first)); - - if (lepton_id == 0) - { - combinations.push_back(convert_to_pdg_id(entry.second)); - } - else - { - combinations.push_back(11); - } - - factors.push_back(1.0); - } - - pineappl_lumi_add(lumi, factors.size(), combinations.data(), factors.data()); - } - - // if there's an empty PDF coefficient vector, we reconstruct the lumi function - if (pdf.empty()) - { - std::vector xfx1(13); - std::vector xfx2(13); - - std::vector> combinations(table->GetNSubproc()); - std::vector> factors(table->GetNSubproc()); - - for (int a = 0; a != 13; ++a) - { - xfx1.at(a) = 1.0; - - for (int b = 0; b != 13; ++b) - { - xfx2.at(b) = 1.0; - - auto const& lumi = comb.CalcPDFLinearCombination(table, xfx1, xfx2, false); - - assert( static_cast (lumi.size()) == table->GetNSubproc() ); - - for (std::size_t i = 0; i != lumi.size(); ++i) - { - if (lumi.at(i) == 0.0) - { - continue; - } - - auto const ap = convert_to_pdg_id(a - 6); - auto const bp = convert_to_pdg_id(b - 6); - - combinations.at(i).push_back(ap); - combinations.at(i).push_back(bp); - factors.at(i).push_back(lumi.at(i)); - } - - xfx2.at(b) = 0.0; - } - - xfx1.at(a) = 0.0; - } - - for (std::size_t i = 0; i != combinations.size(); ++i) - { - pineappl_lumi_add(lumi, factors.at(i).size(), combinations.at(i).data(), - factors.at(i).data()); - } - } -} - -pineappl_grid* convert_coeff_add_fix( - fastNLOCoeffAddFix* table, - fastNLOPDFLinearCombinations const& comb, - std::size_t bins, - uint32_t alpha -) { - std::vector order_params = { static_cast (table->GetNpow()), alpha, 0, 0 }; - - auto* lumi = pineappl_lumi_new(); - create_lumi(table, comb, lumi); - - std::vector bin_limits(bins + 1); - std::iota(bin_limits.begin(), bin_limits.end(), 0.0); - auto* key_vals = pineappl_keyval_new(); - auto* pgrid = pineappl_grid_new(lumi, 1, order_params.data(), bins, bin_limits.data(), - key_vals); - pineappl_keyval_delete(key_vals); - pineappl_lumi_delete(lumi); - - //auto descriptions = table->GetContributionDescription(); - - //for (auto const& desc : descriptions) - //{ - // std::cout << "-- " << desc << '\n'; - //} - - std::size_t n_obs_bin = table->GetNObsBin(); - //std::size_t n_scalevar = table->GetNScalevar(); - //std::size_t n_scale_node = table->GetNScaleNode(); - //std::size_t n_scales = table->GetNScales(); - std::size_t n_subproc = table->GetNSubproc(); - std::size_t total_scalevars = table->GetTotalScalevars(); - std::size_t total_scalenodes = table->GetTotalScalenodes(); - - //std::cout << table->GetScaleDescription(0) << '\n'; - //auto const& descr = table->GetScaleDescr(); - //std::cout << descr.at(0).at(0) << std::endl; - //std::cout << table->GetNevt() << '\n'; - - for (std::size_t obs = 0; obs != n_obs_bin; ++obs) - { - auto const& x1_values = table->GetXNodes1(obs); - - // TODO: is this the correct assumption? - auto const x2_values = (table->GetNxtot2(0) == -1) ? x1_values : table->GetXNodes2(obs); - - for (std::size_t subproc = 0; subproc != n_subproc; ++subproc) - { - double const factor = table->GetNevt(obs, subproc); - - for (std::size_t j = 0; j != total_scalevars; ++j) - { - // TODO: for the time being we only extract the central scale result - if (table->GetScaleFactor(j) != 1.0) - { - continue; - } - - auto q2_values = table->GetScaleNodes(obs, j); - std::vector mu2_values; - - // the values are the unsquared q values, correct that - for (auto& value : q2_values) - { - value *= value; - mu2_values.push_back(value); - mu2_values.push_back(value); - } - - auto* subgrid = pineappl_subgrid_new2(mu2_values.size() / 2, mu2_values.data(), - x1_values.size(), x1_values.data(), x2_values.size(), x2_values.data()); - - // TODO: figure out what the general case is supposed to be - assert( j == 0 ); - - //std::cout << table->GetScaleFactor(j) << '\n'; - - bool non_zero_subgrid = false; - - for (std::size_t mu2_slice = 0; mu2_slice != total_scalenodes; ++mu2_slice) - { - std::vector slice(x1_values.size() * x2_values.size()); - bool non_zero = false; - - std::size_t ix1 = 0; - std::size_t ix2 = 0; - - for (int ix = 0; ix != table->GetNxmax(obs); ++ix) - { - assert( table->GetXIndex(obs, ix1, ix2) == ix ); - - auto const value = table->GetSigmaTilde(obs, j, mu2_slice, ix, subproc); - - if (value != 0.0) - { - non_zero = true; - slice.at(x2_values.size() * ix2 + ix1) = value / factor - * x1_values.at(ix1) * x2_values.at(ix2); - } - - ++ix1; - - switch (table->GetNPDFDim()) - { - case 2: - if (ix1 == x1_values.size()) - { - ix1 = 0; - ++ix2; - } - break; - - case 1: - if (ix1 > ix2) - { - ix1 = 0; - ++ix2; - } - break; - - default: - // TODO: NYI - assert( false ); - } - } - - if (non_zero) - { - non_zero_subgrid = true; - pineappl_subgrid_import_mu2_slice(subgrid, mu2_slice, slice.data()); - } - } - - if (non_zero_subgrid) - { - pineappl_grid_replace_and_delete(pgrid, subgrid, 0, obs, subproc); - } - else - { - pineappl_subgrid_delete(subgrid); - } - } - } - } - - return pgrid; -} - -pineappl_grid* convert_coeff_add_flex( - fastNLOCoeffAddFlex* table, - fastNLOPDFLinearCombinations const& comb, - fastNLO::EScaleFunctionalForm mur_ff, - fastNLO::EScaleFunctionalForm muf_ff, - std::size_t bins, - uint32_t alpha, - int ipub_units -) { - std::vector order_params = { static_cast (table->GetNpow()), alpha, 0, 0 }; - - auto* lumi = pineappl_lumi_new(); - create_lumi(table, comb, lumi); - - std::vector bin_limits(bins + 1); - std::iota(bin_limits.begin(), bin_limits.end(), 0.0); - auto* key_vals = pineappl_keyval_new(); - - // flexible grids always have a hadron in initial state 1 ... - pineappl_keyval_set_string(key_vals, "initial_state_1", - std::to_string(table->GetPDFPDG(0)).c_str()); - // and something else at 2 - pineappl_keyval_set_string(key_vals, "initial_state_2", "11"); - - auto* pgrid = pineappl_grid_new(lumi, 1, order_params.data(), bins, bin_limits.data(), - key_vals); - pineappl_keyval_delete(key_vals); - pineappl_lumi_delete(lumi); - - std::size_t n_obs_bin = table->GetNObsBin(); - std::size_t n_subproc = table->GetNSubproc(); - - auto const& sigma_tildes = table->GetSigmaTildes(); - - auto const rescale = std::pow(0.1, table->GetIXsectUnits() - ipub_units); - - for (std::size_t obs = 0; obs != n_obs_bin; ++obs) - { - auto const& scale_nodes1 = table->GetScaleNodes1(obs); - auto const& scale_nodes2 = table->GetScaleNodes2(obs); - auto const& x1_values = table->GetXNodes1(obs); - - std::vector mur2_values; - - switch (mur_ff) - { - case fastNLO::kScale1: - for (auto const s1 : scale_nodes1) - { - for (std::size_t i = 0; i != scale_nodes2.size(); ++i) - { - mur2_values.push_back(s1 * s1); - } - } - break; - - case fastNLO::kScale2: - for (std::size_t i = 0; i != scale_nodes1.size(); ++i) - { - for (auto const s2 : scale_nodes2) - { - mur2_values.push_back(s2 * s2); - } - } - break; - - case fastNLO::kQuadraticSum: - for (auto const s1 : scale_nodes1) - { - for (auto const s2 : scale_nodes2) - { - mur2_values.push_back(s1 * s1 + s2 * s2); - } - } - break; - - case fastNLO::kQuadraticMean: - for (auto const s1 : scale_nodes1) - { - for (auto const s2 : scale_nodes2) - { - mur2_values.push_back(0.5 * (s1 * s1 + s2 * s2)); - } - } - break; - - default: - // TODO: NYI - assert( false ); - } - - std::vector muf2_values; - - switch (muf_ff) - { - case fastNLO::kScale1: - for (auto const s1 : scale_nodes1) - { - for (std::size_t i = 0; i != scale_nodes2.size(); ++i) - { - muf2_values.push_back(s1 * s1); - } - } - break; - - case fastNLO::kScale2: - for (std::size_t i = 0; i != scale_nodes1.size(); ++i) - { - for (auto const s2 : scale_nodes2) - { - muf2_values.push_back(s2 * s2); - } - } - break; - - case fastNLO::kQuadraticSum: - for (auto const s1 : scale_nodes1) - { - for (auto const s2 : scale_nodes2) - { - muf2_values.push_back(s1 * s1 + s2 * s2); - } - } - break; - - case fastNLO::kQuadraticMean: - for (auto const s1 : scale_nodes1) - { - for (auto const s2 : scale_nodes2) - { - muf2_values.push_back(0.5 * (s1 * s1 + s2 * s2)); - } - } - break; - - default: - // TODO: NYI - assert( false ); - } - - std::vector mu2_values; - - for (std::size_t i = 0; i != scale_nodes1.size() * scale_nodes2.size(); ++i) - { - mu2_values.push_back(mur2_values.at(i)); - mu2_values.push_back(muf2_values.at(i)); - } - - std::vector x2_values{1.0}; - - for (std::size_t subproc = 0; subproc != n_subproc; ++subproc) - { - auto* subgrid = pineappl_subgrid_new2(mu2_values.size() / 2, mu2_values.data(), - x1_values.size(), x1_values.data(), x2_values.size(), x2_values.data()); - - auto const factor = rescale / table->GetNevt(obs, subproc); - bool non_zero_subgrid = false; - - std::size_t mu2_slice = 0; - - for (std::size_t is1 = 0; is1 != scale_nodes1.size(); ++is1) - { - for (std::size_t is2 = 0; is2 != scale_nodes2.size(); ++is2) - { - std::vector slice(x1_values.size()); - bool non_zero = false; - auto const logmur2 = std::log(mu2_values.at(2 * mu2_slice + 0)); - auto const logmuf2 = std::log(mu2_values.at(2 * mu2_slice + 1)); - - // flexible scale grids only allow one initial-state hadron - for (std::size_t ix = 0; ix != sigma_tildes.at(0)->at(obs).size(); ++ix) - { - double value = sigma_tildes.at(0)->at(obs).at(ix).at(is1).at(is2) - .at(subproc); - - if (table->GetNScaleDep() >= 5) - { - // mur - value += logmur2 * - sigma_tildes.at(1)->at(obs).at(ix).at(is1).at(is2).at(subproc); - // muf - value += logmuf2 * - sigma_tildes.at(2)->at(obs).at(ix).at(is1).at(is2).at(subproc); - - if (table->GetNScaleDep() >= 6) - { - // mur mur - value += logmur2 * logmur2 * - sigma_tildes.at(3)->at(obs).at(ix).at(is1).at(is2).at(subproc); - } - - if (table->GetNScaleDep() >= 7) - { - // muf muf - value += logmuf2 * logmuf2 * - sigma_tildes.at(4)->at(obs).at(ix).at(is1).at(is2).at(subproc); - // mur muf - value += logmur2 * logmuf2 * - sigma_tildes.at(5)->at(obs).at(ix).at(is1).at(is2).at(subproc); - } - } - - if (value != 0.0) - { - non_zero = true; - slice.at(ix) = value * factor * x1_values.at(ix); - } - } - - if (non_zero) - { - non_zero_subgrid = true; - pineappl_subgrid_import_mu2_slice(subgrid, mu2_slice, slice.data()); - } - - ++mu2_slice; - } - } - - if (non_zero_subgrid) - { - pineappl_grid_replace_and_delete(pgrid, subgrid, 0, obs, subproc); - } - else - { - pineappl_subgrid_delete(subgrid); - } - } - } - - return pgrid; -} - -int main(int argc, char* argv[]) -{ - if (argc != 3) - { - return EXIT_FAILURE; - } - - std::string in(argv[1]); - std::string out(argv[2]); - - // TODO: read this from an argument - uint32_t alpha = 0; - - LHAPDF::setVerbosity(0); - pdf.reset(LHAPDF::mkPDF("NNPDF31_nlo_as_0118_luxqed", 0)); - - auto file = fastNLOLHAPDF(in, "NNPDF31_nlo_as_0118_luxqed"); - - //file.SelectProcesses({ std::make_pair(0, 0) }); - //file.ActivateContribution(fastNLO::kFixedOrder, fastNLO::kNextToLeading, false); - //file.ActivateContribution(fastNLO::kFixedOrder, fastNLO::kNextToNextToLeading, false); - - auto const id_lo = file.ContrId(fastNLO::kFixedOrder, fastNLO::kLeading); - auto const id_nlo = file.ContrId(fastNLO::kFixedOrder, fastNLO::kNextToLeading); - auto const id_nnlo = file.ContrId(fastNLO::kFixedOrder, fastNLO::kNextToNextToLeading); - - std::vector ids; - - if (id_lo >= 0) - { - ids.push_back(id_lo); - } - - if (id_nlo >= 0) - { - ids.push_back(id_nlo); - } - - if (id_nnlo >= 0) - { - ids.push_back(id_nnlo); - } - - auto const& normalizations = file.GetBinSize(); - std::size_t const bins = normalizations.size(); - - std::vector grids; - - for (auto const id : ids) - { - auto coeff_table = file.GetCoeffTable(id); - assert( coeff_table != nullptr ); - - auto converted = dynamic_cast (coeff_table); - - if (converted != nullptr) - { - grids.push_back(convert_coeff_add_fix(converted, file, bins, alpha)); - } - else - { - auto converted = dynamic_cast (coeff_table); - - if (converted != nullptr) - { - auto const mur_ff = file.GetMuRFunctionalForm(); - auto const muf_ff = file.GetMuFFunctionalForm(); - - grids.push_back(convert_coeff_add_flex(converted, file, mur_ff, muf_ff, bins, alpha, - file.GetIpublunits())); - } - else - { - // TODO: NYI - assert( false ); - } - } - } - - for (std::size_t i = 1; i < grids.size(); ++i) - { - pineappl_grid_merge_and_delete(grids.at(0), grids.at(i)); - } - - pineappl_grid_scale_by_order(grids.at(0), 0.5 / std::acos(-1.0), 1.0, 1.0, 1.0, 1.0); - pineappl_grid_optimize(grids.at(0)); - - auto const dimensions = file.GetNumDiffBin(); - std::vector limits(2 * dimensions * bins); - - for (std::size_t i = 0; i != bins; ++i) - { - for (std::size_t j = 0; j != dimensions; ++j) - { - auto const& pair = file.GetObsBinDimBounds(i, j); - - limits.at(2 * (i * dimensions + j) + 0) = pair.first; - limits.at(2 * (i * dimensions + j) + 1) = pair.second; - } - } - - pineappl_grid_set_remapper(grids.at(0), dimensions, normalizations.data(), limits.data()); - - auto const& results = file.GetCrossSection(); - std::vector other_results(results.size()); - - pineappl_grid_convolute_with_one( - grids.at(0), - 2212, - xfx, - alphas, - pdf.get(), - nullptr, - nullptr, - 1.0, - 1.0, - other_results.data() - ); - - std::cout.setf(std::ios_base::scientific, std::ios_base::floatfield); - std::cout.precision(16); - - bool different = false; - - for (std::size_t i = 0; i != results.size(); ++i) - { - auto const one = results.at(i); - auto const two = other_results.at(i) * normalizations.at(i); - - // catches the case where both results are zero - if (one == two) - { - continue; - } - - if (std::fabs(two / one - 1.0) > 1e-10) - { - std::cout << ">>> fastNLO: " << one << " PineAPPL: " << two << " fN/P: " << (one / two) - << " P/fN: " << (two / one) << '\n'; - different = true; - } - else - { - std::cout << ">>> Success!\n"; - } - } - - pineappl_grid_write(grids.at(0), out.c_str()); - pineappl_grid_delete(grids.at(0)); - - if (different) - { - return 1; - } -} diff --git a/examples/fnlo2pine/meson.build b/examples/fnlo2pine/meson.build deleted file mode 100644 index 160060a4..00000000 --- a/examples/fnlo2pine/meson.build +++ /dev/null @@ -1,26 +0,0 @@ -project( - 'fnlo2pine', - 'cpp', - version : '0.1', - license : 'GPL3', - default_options : [ 'cpp_std=c++17', 'warning_level=3' ], -) - -fastnlo_config = find_program('fnlo-tk-config') - -meson.override_dependency('fnlo-tk', declare_dependency( - compile_args : run_command(fastnlo_config, '--cxxflags').stdout().split(), - link_args : run_command(fastnlo_config, '--ldflags').stdout().split(), - version : run_command(fastnlo_config, '--version').stdout(), -)) - -fastnlo_dep = dependency('fnlo-tk', version : '>=2.3.1') -pineappl_capi_dep = dependency('pineappl_capi', version : '>=0.3') -lhapdf_dep = dependency('lhapdf') - -bin = executable( - 'fnlo2pine', - 'fnlo2pine.cpp', - dependencies : [ fastnlo_dep, lhapdf_dep, pineappl_capi_dep ], - install : true -) diff --git a/pineappl_cli/Cargo.toml b/pineappl_cli/Cargo.toml index e96a2e34..90e1e5e1 100644 --- a/pineappl_cli/Cargo.toml +++ b/pineappl_cli/Cargo.toml @@ -13,6 +13,7 @@ description = "Read, write, and query PineAPPL grids" [dependencies] anyhow = "1.0.48" clap = { features = ["derive", "unstable-replace"], version = "3.1.0" } +cxx = { optional = true, version = "1.0.65" } enum_dispatch = "0.3.7" git-version = "0.3.5" itertools = "0.10.1" @@ -21,6 +22,7 @@ lhapdf = "0.1.11" ndarray = "0.15.4" num_cpus = "1.13.0" pineappl = { path = "../pineappl", version = "0.5.0" } +pineappl_fastnlo = { optional = true, path = "../pineappl_fastnlo", version = "0.1.0" } prettytable-rs = { default-features = false, features = ["win_crlf"], version = "0.8.0" } rayon = "1.5.1" @@ -34,3 +36,6 @@ path = "src/main.rs" [package.metadata.docs.rs] rustc-args = [ "--cfg feature=\"docs-only\"" ] + +[features] +fastnlo = ["pineappl_fastnlo"] diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs new file mode 100644 index 00000000..f8e75823 --- /dev/null +++ b/pineappl_cli/src/import.rs @@ -0,0 +1,629 @@ +use super::helpers::{self, Subcommand}; +use anyhow::{anyhow, Result}; +use clap::{Parser, ValueHint}; +use std::path::PathBuf; + +#[cfg(feature = "fastnlo")] +mod fastnlo { + use super::*; + use itertools::Itertools; + use pineappl::bin::BinRemapper; + use pineappl::grid::{Grid, Order}; + use pineappl::import_only_subgrid::ImportOnlySubgridV2; + use pineappl::lumi::LumiEntry; + use pineappl::sparse_array3::SparseArray3; + use pineappl::subgrid::{Mu2, SubgridParams}; + use pineappl_fastnlo::ffi::{ + self, fastNLOCoeffAddBase, fastNLOCoeffAddFix, fastNLOCoeffAddFlex, fastNLOLHAPDF, + fastNLOPDFLinearCombinations, ESMCalculation, ESMOrder, EScaleFunctionalForm, + }; + use std::convert::{TryFrom, TryInto}; + use std::f64::consts::TAU; + use std::ptr; + + fn pid_to_pdg_id(pid: i32) -> i32 { + match pid { + -6..=-1 | 1..=6 => pid, + 0 => 21, + _ => unimplemented!(), + } + } + + fn create_lumi( + table: &fastNLOCoeffAddBase, + comb: &fastNLOPDFLinearCombinations, + ) -> Vec { + // TODO: set this to the right value if there's only one PDF + let lepton_id = if table.GetNPDF() == 2 { 0 } else { 11 }; + let mut lumis = Vec::new(); + + // if there's a (non-empty) PDF coefficient vector reconstruct the luminosity function; the + // advantage is that we preserve the order of the lumi entries in the PineAPPL grid + for pdf_entry in 0..ffi::GetPDFCoeffSize(table) { + let mut entries = Vec::new(); + + for entry in ffi::GetPDFCoeff(table, pdf_entry) { + let a = pid_to_pdg_id(entry.first); + let b = if lepton_id == 0 { + pid_to_pdg_id(entry.second) + } else { + 11 + }; + let f = 1.0; + + entries.push((a, b, f)); + } + + lumis.push(LumiEntry::new(entries)); + } + + // if the PDF coefficient vector was empty, we must reconstruct the lumi function + if lumis.is_empty() { + let nsubproc = table.GetNSubproc().try_into().unwrap(); + + let mut xfx1 = [0.0; 13]; + let mut xfx2 = [0.0; 13]; + + let mut entries = Vec::new(); + entries.resize(nsubproc, Vec::new()); + + for a in 0..13 { + xfx1[a] = 1.0; + + for b in 0..13 { + xfx2[b] = 1.0; + + let lumi = ffi::CalcPDFLinearCombination(comb, table, &xfx1, &xfx2, false); + + assert!(lumi.len() == nsubproc); + + for (i, l) in lumi.iter().copied().enumerate().filter(|(_, l)| *l != 0.0) { + let ap = pid_to_pdg_id(i32::try_from(a).unwrap() - 6); + let bp = pid_to_pdg_id(i32::try_from(b).unwrap() - 6); + + entries[i].push((ap, bp, l)); + } + + xfx2[b] = 0.0; + } + + xfx1[a] = 0.0; + } + + lumis = entries.into_iter().map(LumiEntry::new).collect(); + } + + lumis + } + + fn convert_coeff_add_fix( + table: &fastNLOCoeffAddFix, + comb: &fastNLOPDFLinearCombinations, + bins: usize, + alpha: u32, + ) -> Grid { + let table_as_add_base = ffi::downcast_coeff_add_fix_to_base(table); + + let mut grid = Grid::new( + create_lumi(table_as_add_base, comb), + vec![Order { + alphas: table_as_add_base.GetNpow().try_into().unwrap(), + alpha, + logxir: 0, + logxif: 0, + }], + (0..=bins).map(|limit| limit as f64).collect(), + SubgridParams::default(), + ); + + let total_scalenodes: usize = table.GetTotalScalenodes().try_into().unwrap(); + + for obs in 0..table_as_add_base.GetNObsBin() { + let x1_values = ffi::GetXNodes1(table_as_add_base, obs); + + // TODO: is this the correct assumption? + let x2_values = if table.GetNxtot2(0) == -1 { + x1_values.clone() + } else { + ffi::GetXNodes2(table_as_add_base, obs) + }; + + for subproc in 0..table_as_add_base.GetNSubproc() { + let factor = table_as_add_base.GetNevt(obs, subproc); + + for j in 0..table.GetTotalScalevars() { + // TODO: for the time being we only extract the central scale result + if table.GetScaleFactor(j) != 1.0 { + continue; + } + + let q_values = ffi::GetScaleNodes(table, obs, j); + let mu2_values: Vec<_> = q_values + .iter() + .map(|q| Mu2 { + ren: q * q, + fac: q * q, + }) + .collect(); + + let mut array = + SparseArray3::new(mu2_values.len(), x1_values.len(), x2_values.len()); + + // TODO: figure out what the general case is supposed to be + assert_eq!(j, 0); + + for mu2_slice in 0..total_scalenodes { + let mut ix1: usize = 0; + let mut ix2: usize = 0; + + for ix in 0..table.GetNxmax(obs) { + assert_eq!( + table.GetXIndex( + obs, + ix1.try_into().unwrap(), + ix2.try_into().unwrap() + ), + ix + ); + + let value = table.GetSigmaTilde( + obs, + j, + mu2_slice.try_into().unwrap(), + ix, + subproc, + ); + + if value != 0.0 { + array[[mu2_slice, ix2, ix1]] = + value / factor * x1_values[ix1] * x2_values[ix2]; + } + + ix1 += 1; + + match table.GetNPDFDim() { + 2 => { + if ix1 == x1_values.len() { + ix1 = 0; + ix2 += 1; + } + } + 1 => { + if ix1 > ix2 { + ix1 = 0; + ix2 += 1; + } + } + _ => unreachable!(), + } + } + } + + if !array.is_empty() { + grid.set_subgrid( + 0, + obs.try_into().unwrap(), + subproc.try_into().unwrap(), + ImportOnlySubgridV2::new( + array, + mu2_values, + x1_values.clone(), + x2_values.clone(), + ) + .into(), + ); + } + } + } + } + + grid + } + + fn convert_coeff_add_flex( + table: &fastNLOCoeffAddFlex, + comb: &fastNLOPDFLinearCombinations, + mur_ff: EScaleFunctionalForm, + muf_ff: EScaleFunctionalForm, + bins: usize, + alpha: u32, + ipub_units: i32, + ) -> Grid { + let table_as_add_base = ffi::downcast_coeff_add_flex_to_base(table); + + let mut grid = Grid::new( + create_lumi(table_as_add_base, comb), + vec![Order { + alphas: table_as_add_base.GetNpow().try_into().unwrap(), + alpha, + logxir: 0, + logxif: 0, + }], + (0..=bins).map(|limit| limit as f64).collect(), + SubgridParams::default(), + ); + + grid.set_key_value("initial_state_1", &table.GetPDFPDG(0).to_string()); + grid.set_key_value("initial_state_2", "11"); + + 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 mut mur2_values = Vec::new(); + + match mur_ff { + EScaleFunctionalForm::kScale1 => { + for s1 in &scale_nodes1 { + for _ in 0..scale_nodes2.len() { + mur2_values.push(s1 * s1); + } + } + } + EScaleFunctionalForm::kScale2 => { + for _ in 0..scale_nodes1.len() { + for s2 in &scale_nodes2 { + mur2_values.push(s2 * s2); + } + } + } + EScaleFunctionalForm::kQuadraticSum => { + for s1 in &scale_nodes1 { + for s2 in &scale_nodes2 { + mur2_values.push(s1 * s1 + s2 * s2); + } + } + } + EScaleFunctionalForm::kQuadraticMean => { + for s1 in &scale_nodes1 { + for s2 in &scale_nodes2 { + mur2_values.push(0.5 * (s1 * s1 + s2 * s2)); + } + } + } + _ => { + // TODO: NYI + unimplemented!() + } + } + + let mut muf2_values = Vec::new(); + + match muf_ff { + EScaleFunctionalForm::kScale1 => { + for s1 in &scale_nodes1 { + for _ in 0..scale_nodes2.len() { + muf2_values.push(s1 * s1); + } + } + } + EScaleFunctionalForm::kScale2 => { + for _ in 0..scale_nodes1.len() { + for s2 in &scale_nodes2 { + muf2_values.push(s2 * s2); + } + } + } + EScaleFunctionalForm::kQuadraticSum => { + for s1 in &scale_nodes1 { + for s2 in &scale_nodes2 { + muf2_values.push(s1 * s1 + s2 * s2); + } + } + } + EScaleFunctionalForm::kQuadraticMean => { + for s1 in &scale_nodes1 { + for s2 in &scale_nodes2 { + muf2_values.push(0.5 * (s1 * s1 + s2 * s2)); + } + } + } + _ => { + // TODO: NYI + unimplemented!() + } + } + + let mut mu2_values = Vec::new(); + + for i in 0..(scale_nodes1.len() * scale_nodes2.len()) { + mu2_values.push(Mu2 { + ren: mur2_values[i], + fac: muf2_values[i], + }); + } + + 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 array = SparseArray3::new(mu2_values.len(), x1_values.len(), 1); + + 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(); + + // flexible scale grids only allow one initial-state hadron + for ix in 0..nx { + let mut value = ffi::GetSigmaTilde(table, 0, obs, ix, is1, is2, subproc); + + if table.GetNScaleDep() >= 5 { + // mur + value += + logmur2 * ffi::GetSigmaTilde(table, 1, obs, ix, is1, is2, subproc); + // muf + value += + logmuf2 * ffi::GetSigmaTilde(table, 2, obs, ix, is1, is2, subproc); + + if table.GetNScaleDep() >= 6 { + // mur mur + value += logmur2 + * logmur2 + * ffi::GetSigmaTilde(table, 3, obs, ix, is1, is2, subproc); + } + + if table.GetNScaleDep() >= 7 { + // muf muf + value += logmuf2 + * logmuf2 + * ffi::GetSigmaTilde(table, 4, obs, ix, is1, is2, subproc); + // mur muf + value += logmur2 + * logmuf2 + * ffi::GetSigmaTilde(table, 5, obs, ix, is1, is2, subproc); + } + } + + if value != 0.0 { + array[[mu2_slice, ix, 0]] = value * factor * x1_values[ix]; + } + } + } + + if !array.is_empty() { + grid.set_subgrid( + 0, + obs.try_into().unwrap(), + subproc.try_into().unwrap(), + ImportOnlySubgridV2::new( + array, + mu2_values.clone(), + x1_values.clone(), + vec![1.0], + ) + .into(), + ); + } + } + } + + grid + } + + pub fn convert_fastnlo_table(file: &fastNLOLHAPDF, alpha: u32) -> Result { + let file_as_reader = ffi::downcast_lhapdf_to_reader(file); + let file_as_table = ffi::downcast_lhapdf_to_table(file); + + let ids: Vec<_> = [ + file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kLeading), + file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kNextToLeading), + file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kNextToNextToLeading), + ] + .iter() + .copied() + .filter(|&id| id >= 0) + .collect(); + + let bins: usize = file_as_table.GetNObsBin().try_into().unwrap(); + let mut grids = Vec::new(); + + for id in ids { + let coeff_table = file_as_table.GetCoeffTable(id); + + if coeff_table == ptr::null_mut() { + return Err(anyhow!("could not access coefficient table")); + } + + let linear_combinations = + ffi::downcast_reader_to_pdf_linear_combinations(file_as_reader); + + let converted = unsafe { ffi::dynamic_cast_coeff_add_fix(coeff_table) }; + + if converted != ptr::null_mut() { + grids.push(convert_coeff_add_fix( + unsafe { &*converted }, + linear_combinations, + bins, + alpha, + )); + } else { + let converted = unsafe { ffi::dynamic_cast_coeff_add_flex(coeff_table) }; + + if converted != ptr::null_mut() { + let mur_ff = file_as_reader.GetMuRFunctionalForm(); + let muf_ff = file_as_reader.GetMuFFunctionalForm(); + + grids.push(convert_coeff_add_flex( + unsafe { &*converted }, + linear_combinations, + mur_ff, + muf_ff, + bins, + alpha, + file_as_table.GetIpublunits(), + )); + } else { + // TODO: NYI + unreachable!(); + } + } + } + + let mut result = grids.remove(0); + for grid in grids { + result.merge(grid)?; + } + + result.scale_by_order(1.0 / TAU, 1.0, 1.0, 1.0, 1.0); + result.optimize(); + + let dimensions: usize = file_as_table.GetNumDiffBin().try_into().unwrap(); + let mut limits = Vec::new(); + let normalizations = if file_as_table.IsNorm() { + ffi::GetBinSize(file_as_table) + } else { + vec![1.0; bins] + }; + limits.reserve(2 * dimensions * bins); + + for i in 0..bins { + for j in 0..dimensions { + let pair = ffi::GetObsBinDimBounds( + file_as_table, + i.try_into().unwrap(), + j.try_into().unwrap(), + ); + + limits.push((pair.first, pair.second)); + } + } + + result.set_remapper(BinRemapper::new(normalizations.clone(), limits).unwrap())?; + + Ok(result) + } +} + +/// Converts fastNLO tables to PineAPPL grids. +#[derive(Parser)] +pub struct Opts { + /// Path to the input grid. + #[clap(parse(from_os_str), value_hint = ValueHint::FilePath)] + input: PathBuf, + /// Path to the converted grid. + #[clap(parse(from_os_str), value_hint = ValueHint::FilePath)] + output: PathBuf, + /// LHAPDF id or name of the PDF set to check the converted grid with. + #[clap(validator = helpers::validate_pdfset)] + pdfset: String, + /// LO coupling power in alpha. + #[clap(default_value_t = 0, long)] + alpha: u32, + /// Relative threshold between the table and the converted grid when comparison fails. + #[clap(default_value = "1e-10", long)] + accuracy: f64, + /// Prevents fastNLO from printing output. + #[clap(long = "silence-fastnlo")] + silence_fastnlo: bool, +} + +impl Subcommand for Opts { + #[cfg(feature = "fastnlo")] + fn run(&self) -> Result<()> { + use lhapdf::Pdf; + use pineappl_fastnlo::ffi::{self, Verbosity}; + use prettytable::{cell, row}; + + if self.silence_fastnlo { + ffi::SetGlobalVerbosity(Verbosity::SILENT); + } + + let mut file = ffi::make_fastnlo_lhapdf_with_name_file_set( + self.input.to_str().unwrap(), + &self.pdfset, + 0, + self.silence_fastnlo, + ); + let grid = fastnlo::convert_fastnlo_table(&file, self.alpha)?; + + let table_results = ffi::GetCrossSection( + ffi::downcast_lhapdf_to_reader_mut(file.as_mut().unwrap()), + false, + ); + let pdf = self.pdfset.parse().map_or_else( + |_| Pdf::with_setname_and_member(&self.pdfset, 0), + Pdf::with_lhaid, + ); + let results = helpers::convolute(&grid, &pdf, &[], &[], &[], 1); + + let mut different = false; + + let mut table = helpers::create_table(); + table.set_titles(row![c => "b", "PineAPPL", "fastNLO", "rel. diff"]); + + for (bin, (one, two)) in results + .into_iter() + .zip(table_results.into_iter()) + .enumerate() + { + // catches the case where both results are zero + let rel_diff = if one == two { 0.0 } else { two / one - 1.0 }; + + if rel_diff.abs() > self.accuracy { + different = false; + } + + table.add_row(row![ + bin.to_string(), + r->format!("{:.7e}", one), + r->format!("{:.7e}", two), + r->format!("{:.7e}", rel_diff) + ]); + } + + table.printstd(); + + if different { + Err(anyhow!("grids are different")) + } else { + helpers::write_grid(&self.output, &grid) + } + } + + #[cfg(not(feature = "fastnlo"))] + fn run(&self) -> Result<()> { + Err(anyhow!( + "you need to install `pineappl` with feature `fastnlo`" + )) + } +} + +#[cfg(test)] +mod tests { + use assert_cmd::Command; + + const HELP_STR: &str = "pineappl-import +Converts fastNLO tables to PineAPPL grids + +USAGE: + pineappl import [OPTIONS] + +ARGS: + Path to the input grid + Path to the converted grid + LHAPDF id or name of the PDF set to check the converted grid with + +OPTIONS: + --accuracy Relative threshold between the table and the converted grid when + comparison fails [default: 1e-10] + --alpha LO coupling power in alpha [default: 0] + -h, --help Print help information + --silence-fastnlo Prevents fastNLO from printing output +"; + + #[test] + fn help() { + Command::cargo_bin("pineappl") + .unwrap() + .args(&["import", "--help"]) + .assert() + .success() + .stdout(HELP_STR); + } +} diff --git a/pineappl_cli/src/main.rs b/pineappl_cli/src/main.rs index 34e270ae..90c14331 100644 --- a/pineappl_cli/src/main.rs +++ b/pineappl_cli/src/main.rs @@ -5,6 +5,7 @@ mod convolute; mod delete; mod diff; mod helpers; +mod import; mod info; mod merge; mod obl; @@ -54,6 +55,7 @@ enum SubcommandEnum { Convolute(convolute::Opts), Delete(delete::Opts), Diff(diff::Opts), + Import(import::Opts), Info(info::Opts), Merge(merge::Opts), Obl(obl::Opts), diff --git a/pineappl_fastnlo/Cargo.toml b/pineappl_fastnlo/Cargo.toml new file mode 100644 index 00000000..38b488e8 --- /dev/null +++ b/pineappl_fastnlo/Cargo.toml @@ -0,0 +1,17 @@ +[package] +name = "pineappl_fastnlo" +version = "0.1.0" +authors = ["Christopher Schwan "] +edition = "2018" +license = "GPL-3.0-or-later" +repository = "https://github.com/N3PDF/pineappl" +readme = "README.md" +keywords = ["high-energy-physics", "physics"] +categories = ["science"] +description = "PineAPPL's interface to fastNLO" + +[dependencies] +cxx = { version = "1.0.65" } + +[build-dependencies] +cxx-build = { version = "1.0.65" } diff --git a/pineappl_fastnlo/build.rs b/pineappl_fastnlo/build.rs new file mode 100644 index 00000000..ab82dfa7 --- /dev/null +++ b/pineappl_fastnlo/build.rs @@ -0,0 +1,34 @@ +use std::process::Command; + +fn main() { + let fnlo_lib_path = String::from_utf8( + Command::new("fnlo-tk-config") + .arg("--libdir") + .output() + .expect("did not find `fnlo-tk-config`, please install fastNLO") + .stdout, + ) + .unwrap(); + + println!("cargo:rustc-link-search={}", fnlo_lib_path.trim()); + + let fnlo_include_path = String::from_utf8( + Command::new("fnlo-tk-config") + .arg("--incdir") + .output() + .expect("did not find `fnlo-tk-config`, please install fastNLO") + .stdout, + ) + .unwrap(); + + println!("cargo:rustc-link-lib=fastnlotoolkit"); + + cxx_build::bridge("src/lib.rs") + .file("src/fastnlo.cpp") + .include(fnlo_include_path.trim()) + .compile("fnlo-bridge"); + + println!("cargo:rerun-if-changed=src/lib.rs"); + println!("cargo:rerun-if-changed=src/fastnlo.cpp"); + println!("cargo:rerun-if-changed=src/fastnlo.hpp"); +} diff --git a/pineappl_fastnlo/src/fastnlo.cpp b/pineappl_fastnlo/src/fastnlo.cpp new file mode 100644 index 00000000..82e7e2ad --- /dev/null +++ b/pineappl_fastnlo/src/fastnlo.cpp @@ -0,0 +1,207 @@ +#include "pineappl_fastnlo/src/fastnlo.hpp" + +// TODO: is this portable enough? +#if defined (__unix__) || defined(__unix) || defined(unix) || \ + (defined (__APPLE__) && defined (__MACH__)) +#define HAVE_UNISTD_H +#endif + +#ifdef HAVE_UNISTD_H +#include +#endif + +#include +#include +#include + +template +rust::Vec std_vector_to_rust_vec(std::vector vector) +{ + rust::Vec result; + result.reserve(vector.size()); + std::move(vector.begin(), vector.end(), std::back_inserter(result)); + return result; +} + +rust::Vec CalcPDFLinearCombination( + fastNLOPDFLinearCombinations const& lc, + fastNLOCoeffAddBase const& base, + rust::Slice pdfx1, + rust::Slice pdfx2, + bool pdf2IsAntiParticle +) { + std::vector fx1(pdfx1.begin(), pdfx1.end()); + std::vector fx2(pdfx2.begin(), pdfx2.end()); + + return std_vector_to_rust_vec(lc.CalcPDFLinearCombination(&base, fx1, fx2, pdf2IsAntiParticle)); +} + +rust::Vec GetScaleNodes(fastNLOCoeffAddFix const& coeffs, int iObs, int iSvar) +{ + return std_vector_to_rust_vec(coeffs.GetScaleNodes(iObs, iSvar)); +} + +rust::Vec GetXNodes1(fastNLOCoeffAddBase const& coeffs, int iObsBin) +{ + return std_vector_to_rust_vec(coeffs.GetXNodes1(iObsBin)); +} + +rust::Vec GetXNodes2(fastNLOCoeffAddBase const& coeffs, int iObsBin) +{ + return std_vector_to_rust_vec(coeffs.GetXNodes2(iObsBin)); +} + +std::unique_ptr make_fastnlo_lhapdf_with_name_file_set( + rust::Str name, + rust::Str LHAPDFfile, + int PDFSet, + bool silence +) { + std::string arg0(name.begin(), name.end()); + std::string arg1(LHAPDFfile.begin(), LHAPDFfile.end()); + + std::unique_ptr result; + +#ifdef HAVE_UNISTD_H + int backup_fd; + FILE *dev_null; + + // the constructor of `fastNLOLHAPDF` isn't completely silent even we nicely ask it to, so we + // have to resort to drastic measures + if (silence) + { + fflush(stdout); + backup_fd = dup(STDOUT_FILENO); + dev_null = fopen("/dev/null", "w"); + dup2(fileno(dev_null), STDOUT_FILENO); + } +#endif + + result.reset(new fastNLOLHAPDF(arg0, arg1, PDFSet)); + +#ifdef HAVE_UNISTD_H + if (silence) + { + fflush(stdout); + fclose(dev_null); + dup2(backup_fd, STDOUT_FILENO); + close(backup_fd); + } +#endif + + return result; +} + +rust::Vec GetCrossSection(fastNLOReader& reader, bool lNorm) +{ + // if we don't unconditionally silence LHAPDF in fastNLO its header will be shown twice + LHAPDF::setVerbosity(0); + return std_vector_to_rust_vec(reader.GetCrossSection(lNorm)); +} + +rust::Vec GetBinSize(fastNLOTable const& table) +{ + return std_vector_to_rust_vec(table.GetBinSize()); +} + +rust::Vec GetScaleNodes1(fastNLOCoeffAddFlex const& coeffs, int iObsBin) +{ + return std_vector_to_rust_vec(coeffs.GetScaleNodes1(iObsBin)); +} + +rust::Vec GetScaleNodes2(fastNLOCoeffAddFlex const& coeffs, int iObsBin) +{ + return std_vector_to_rust_vec(coeffs.GetScaleNodes2(iObsBin)); +} + +std::size_t GetPDFCoeffSize(fastNLOCoeffAddBase const& coeffs) +{ + return coeffs.GetPDFCoeff().size(); +} + +rust::Vec GetPDFCoeff(fastNLOCoeffAddBase const& coeffs, std::size_t index) +{ + auto const& cpp_result = coeffs.GetPDFCoeff().at(index); + rust::Vec result; + result.reserve(cpp_result.size()); + + for (auto const& pair : cpp_result) + { + pair_int_int res; + res.first = pair.first; + res.second = pair.second; + result.push_back(res); + } + + return result; +} + +double GetSigmaTilde( + fastNLOCoeffAddFlex const& coeffs, + std::size_t mu, + std::size_t obs, + std::size_t ix, + std::size_t is1, + std::size_t is2, + int subproc +) { + return coeffs.GetSigmaTildes().at(mu)->at(obs).at(ix).at(is1).at(is2).at(subproc); +} + +std::size_t GetNx(fastNLOCoeffAddFlex const& coeffs, std::size_t obs) +{ + return coeffs.GetSigmaTildes().at(0)->at(obs).size(); +} + +fastNLOCoeffAddBase const& downcast_coeff_add_fix_to_base(fastNLOCoeffAddFix const& coeffs) +{ + return coeffs; +} + +fastNLOCoeffAddBase const& downcast_coeff_add_flex_to_base(fastNLOCoeffAddFlex const& coeffs) +{ + return coeffs; +} + +fastNLOReader const& downcast_lhapdf_to_reader(fastNLOLHAPDF const& lhapdf) +{ + return lhapdf; +} + +fastNLOReader& downcast_lhapdf_to_reader_mut(fastNLOLHAPDF& lhapdf) +{ + return lhapdf; +} + +fastNLOTable const& downcast_lhapdf_to_table(fastNLOLHAPDF const& lhapdf) +{ + return lhapdf; +} + +fastNLOCoeffAddFix* dynamic_cast_coeff_add_fix(fastNLOCoeffBase* coeffs) +{ + return dynamic_cast (coeffs); +} + +fastNLOCoeffAddFlex* dynamic_cast_coeff_add_flex(fastNLOCoeffBase* coeffs) +{ + return dynamic_cast (coeffs); +} + +fastNLOPDFLinearCombinations const& downcast_reader_to_pdf_linear_combinations( + fastNLOReader const& reader +) { + return reader; +} + +pair_double_double GetObsBinDimBounds( + fastNLOTable const& table, + unsigned int iObs, + unsigned int iDim +) { + pair_double_double result; + auto const cpp_result = table.GetObsBinDimBounds(iObs, iDim); + result.first = cpp_result.first; + result.second = cpp_result.second; + return result; +} diff --git a/pineappl_fastnlo/src/fastnlo.hpp b/pineappl_fastnlo/src/fastnlo.hpp new file mode 100644 index 00000000..a7004381 --- /dev/null +++ b/pineappl_fastnlo/src/fastnlo.hpp @@ -0,0 +1,85 @@ +#ifndef FASTNLO_HPP +#define FASTNLO_HPP + +#include "pineappl_fastnlo/src/lib.rs.h" +#include "rust/cxx.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +std::unique_ptr make_fastnlo_lhapdf_with_name_file_set( + rust::Str name, + rust::Str LHAPDFfile, + int PDFSet, + bool silence +); + +rust::Vec CalcPDFLinearCombination( + fastNLOPDFLinearCombinations const& lc, + fastNLOCoeffAddBase const& base, + rust::Slice pdfx1, + rust::Slice pdfx2, + bool pdf2IsAntiParticle +); + +rust::Vec GetScaleNodes(fastNLOCoeffAddFix const& coeffs, int iObs, int iSvar); + +rust::Vec GetXNodes1(fastNLOCoeffAddBase const& coeffs, int iObsBin); + +rust::Vec GetXNodes2(fastNLOCoeffAddBase const& coeffs, int iObsBin); + +rust::Vec GetCrossSection(fastNLOReader& reader, bool lNorm); + +rust::Vec GetBinSize(fastNLOTable const& table); + +rust::Vec GetScaleNodes1(fastNLOCoeffAddFlex const& coeffs, int iObsBin); + +rust::Vec GetScaleNodes2(fastNLOCoeffAddFlex const& coeffs, int iObsBin); + +std::size_t GetPDFCoeffSize(fastNLOCoeffAddBase const& coeffs); + +rust::Vec GetPDFCoeff(fastNLOCoeffAddBase const& coeffs, std::size_t index); + +double GetSigmaTilde( + fastNLOCoeffAddFlex const& coeffs, + std::size_t, + std::size_t, + std::size_t, + std::size_t, + std::size_t, + int +); + +std::size_t GetNx(fastNLOCoeffAddFlex const& coeffs, std::size_t); + +fastNLOCoeffAddBase const& downcast_coeff_add_fix_to_base(fastNLOCoeffAddFix const& coeffs); + +fastNLOCoeffAddBase const& downcast_coeff_add_flex_to_base(fastNLOCoeffAddFlex const& coeffs); + +fastNLOReader const& downcast_lhapdf_to_reader(fastNLOLHAPDF const& lhapdf); + +fastNLOReader& downcast_lhapdf_to_reader_mut(fastNLOLHAPDF& lhapdf); + +fastNLOTable const& downcast_lhapdf_to_table(fastNLOLHAPDF const& lhapdf); + +fastNLOCoeffAddFix* dynamic_cast_coeff_add_fix(fastNLOCoeffBase* coeffs); + +fastNLOCoeffAddFlex* dynamic_cast_coeff_add_flex(fastNLOCoeffBase* coeffs); + +fastNLOPDFLinearCombinations const& downcast_reader_to_pdf_linear_combinations( + fastNLOReader const& reader +); + +pair_double_double GetObsBinDimBounds( + fastNLOTable const& table, + unsigned int iObs, + unsigned int iDim +); + +#endif diff --git a/pineappl_fastnlo/src/lib.rs b/pineappl_fastnlo/src/lib.rs new file mode 100644 index 00000000..10b8c9f2 --- /dev/null +++ b/pineappl_fastnlo/src/lib.rs @@ -0,0 +1,216 @@ +#[cxx::bridge] +pub mod ffi { + #[namespace = "fastNLO"] + #[repr(u32)] + enum EScaleFunctionalForm { + kScale1, + kScale2, + kQuadraticSum, + kQuadraticMean, + kQuadraticSumOver4, + kLinearMean, + kLinearSum, + kScaleMax, + kScaleMin, + kProd, + kS2plusS1half, + kPow4Sum, + kWgtAvg, + kS2plusS1fourth, + kExpProd2, + kExtern, + kConst, + } + + #[namespace = "fastNLO"] + #[repr(u32)] + enum ESMCalculation { + kFixedOrder = 0, + kThresholdCorrection = 1, + kElectroWeakCorrection = 2, + kNonPerturbativeCorrection = 3, + kContactInteraction = 10, + } + + #[namespace = "fastNLO"] + #[repr(u32)] + enum ESMOrder { + kLeading, + kNextToLeading, + kNextToNextToLeading, + } + + struct pair_double_double { + first: f64, + second: f64, + } + + struct pair_int_int { + first: i32, + second: i32, + } + + extern "C++" { + include!("fastnlotk/fastNLOConstants.h"); + + #[namespace = "fastNLO"] + type EScaleFunctionalForm; + + #[namespace = "fastNLO"] + type ESMCalculation; + + #[namespace = "fastNLO"] + type ESMOrder; + } + + #[namespace = "say"] + #[repr(i32)] + enum Verbosity { + DEBUG = -1000, + MANUAL = 2, + INFO = 0, + WARNING = 1, + ERROR = 2, + SILENT = 1000, + } + + unsafe extern "C++" { + include!("fastnlotk/speaker.h"); + + #[namespace = "say"] + type Verbosity; + + #[namespace = "say"] + fn SetGlobalVerbosity(_: Verbosity) -> i32; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddBase.h"); + + type fastNLOCoeffAddBase; + + fn GetNevt(&self, _: i32, _: i32) -> f64; + fn GetNObsBin(&self) -> i32; + fn GetNPDF(&self) -> i32; + fn GetNSubproc(&self) -> i32; + fn GetNpow(&self) -> i32; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddFix.h"); + + type fastNLOCoeffAddFix; + + fn GetNPDFDim(&self) -> i32; + fn GetNxmax(&self, _: i32) -> i32; + fn GetNxtot2(&self, _: i32) -> i32; + fn GetScaleFactor(&self, _: i32) -> f64; + fn GetSigmaTilde(&self, _: i32, _: i32, _: i32, _: i32, _: i32) -> f64; + fn GetTotalScalenodes(&self) -> i32; + fn GetTotalScalevars(&self) -> i32; + fn GetXIndex(&self, _: i32, _: i32, _: i32) -> i32; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOPDFLinearCombinations.h"); + + type fastNLOPDFLinearCombinations; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOLHAPDF.h"); + + type fastNLOLHAPDF; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffBase.h"); + + type fastNLOCoeffBase; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOReader.h"); + + type fastNLOReader; + + fn GetMuRFunctionalForm(&self) -> EScaleFunctionalForm; + fn GetMuFFunctionalForm(&self) -> EScaleFunctionalForm; + fn ContrId(&self, _: ESMCalculation, _: ESMOrder) -> i32; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOTable.h"); + + type fastNLOTable; + + fn GetCoeffTable(&self, _: i32) -> *mut fastNLOCoeffBase; + fn GetIpublunits(&self) -> i32; + fn GetNObsBin(&self) -> u32; + fn GetNumDiffBin(&self) -> u32; + fn IsNorm(&self) -> bool; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddFlex.h"); + + type fastNLOCoeffAddFlex; + + fn GetIXsectUnits(&self) -> i32; + fn GetNScaleDep(&self) -> i32; + fn GetPDFPDG(&self, _: i32) -> i32; + } + + unsafe extern "C++" { + include!("pineappl_fastnlo/src/fastnlo.hpp"); + + fn CalcPDFLinearCombination( + _: &fastNLOPDFLinearCombinations, + _: &fastNLOCoeffAddBase, + _: &[f64], + _: &[f64], + _: bool, + ) -> Vec; + + fn GetBinSize(_: &fastNLOTable) -> Vec; + fn GetCrossSection(_: Pin<&mut fastNLOReader>, _: bool) -> Vec; + fn GetNx(_: &fastNLOCoeffAddFlex, _: usize) -> usize; + fn GetObsBinDimBounds(_: &fastNLOTable, _: u32, _: u32) -> pair_double_double; + fn GetPDFCoeff(_: &fastNLOCoeffAddBase, index: usize) -> Vec; + fn GetPDFCoeffSize(_: &fastNLOCoeffAddBase) -> usize; + fn GetScaleNodes(_: &fastNLOCoeffAddFix, _: i32, _: i32) -> Vec; + fn GetScaleNodes1(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; + fn GetScaleNodes2(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; + fn GetXNodes1(_: &fastNLOCoeffAddBase, _: i32) -> Vec; + fn GetXNodes2(_: &fastNLOCoeffAddBase, _: i32) -> Vec; + fn GetSigmaTilde( + _: &fastNLOCoeffAddFlex, + _: usize, + _: usize, + _: usize, + _: usize, + _: usize, + _: i32, + ) -> f64; + + fn downcast_coeff_add_fix_to_base(_: &fastNLOCoeffAddFix) -> &fastNLOCoeffAddBase; + fn downcast_coeff_add_flex_to_base(_: &fastNLOCoeffAddFlex) -> &fastNLOCoeffAddBase; + fn downcast_lhapdf_to_reader(_: &fastNLOLHAPDF) -> &fastNLOReader; + fn downcast_lhapdf_to_reader_mut(_: Pin<&mut fastNLOLHAPDF>) -> Pin<&mut fastNLOReader>; + fn downcast_lhapdf_to_table(_: &fastNLOLHAPDF) -> &fastNLOTable; + fn downcast_reader_to_pdf_linear_combinations( + _: &fastNLOReader, + ) -> &fastNLOPDFLinearCombinations; + + unsafe fn dynamic_cast_coeff_add_fix(_: *mut fastNLOCoeffBase) -> *mut fastNLOCoeffAddFix; + unsafe fn dynamic_cast_coeff_add_flex(_: *mut fastNLOCoeffBase) + -> *mut fastNLOCoeffAddFlex; + + fn make_fastnlo_lhapdf_with_name_file_set( + name: &str, + lhapdf: &str, + set: i32, + silence: bool, + ) -> UniquePtr; + } +}