From da1eb897a091e8c1e87762735bad25f5fa4f14b6 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 19 Feb 2022 12:22:51 +0100 Subject: [PATCH 01/26] Make fastNLO convert part of the CLI --- pineappl_cli/Cargo.toml | 8 + pineappl_cli/src/import.rs | 74 ++++ pineappl_cli/src/import/fastnlo.cpp | 656 ++++++++++++++++++++++++++++ pineappl_cli/src/import/fastnlo.hpp | 3 + pineappl_cli/src/main.rs | 2 + 5 files changed, 743 insertions(+) create mode 100644 pineappl_cli/src/import.rs create mode 100644 pineappl_cli/src/import/fastnlo.cpp create mode 100644 pineappl_cli/src/import/fastnlo.hpp diff --git a/pineappl_cli/Cargo.toml b/pineappl_cli/Cargo.toml index e96a2e34..86faa061 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" @@ -24,6 +25,10 @@ pineappl = { path = "../pineappl", version = "0.5.0" } prettytable-rs = { default-features = false, features = ["win_crlf"], version = "0.8.0" } rayon = "1.5.1" +[build-dependencies] +cxx-build = { optional = true, version = "1.0.65" } +pkg-config = { optional = true, version = "0.3.24" } + [dev-dependencies] assert_cmd = "2.0.2" assert_fs = "1.0.6" @@ -34,3 +39,6 @@ path = "src/main.rs" [package.metadata.docs.rs] rustc-args = [ "--cfg feature=\"docs-only\"" ] + +[features] +import-fastnlo = ["cxx", "cxx-build", "pkg-config"] diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs new file mode 100644 index 00000000..e7b02d40 --- /dev/null +++ b/pineappl_cli/src/import.rs @@ -0,0 +1,74 @@ +use super::helpers::Subcommand; +use anyhow::{anyhow, Result}; +use clap::{Parser, ValueHint}; +use std::path::PathBuf; + +// the indirection of two modules instead of one is needed to hide the C++ bridge behind the +// feature gate +#[cfg(feature = "import-fastnlo")] +mod fastnlo { + #[cxx::bridge] + pub mod ffi { + unsafe extern "C++" { + include!("pineappl_cli/src/import/fastnlo.hpp"); + + unsafe fn this_would_be_main(input: *const c_char, output: *const c_char) -> i32; + } + } +} + +/// 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, +} + +impl Subcommand for Opts { + #[cfg(feature = "import-fastnlo")] + fn run(&self) -> Result<()> { + use std::ffi::CString; + use std::os::unix::ffi::OsStrExt; + + // TODO: this is a correct way of converting strings, but only works on UNIX + let input = CString::new(self.input.as_path().as_os_str().as_bytes()).unwrap(); + let output = CString::new(self.output.as_path().as_os_str().as_bytes()).unwrap(); + + let errcode = unsafe { fastnlo::ffi::this_would_be_main(input.as_ptr(), output.as_ptr()) }; + + if errcode == 0 { + Ok(()) + } else { + Err(anyhow!("something went wrong")) + } + } + + #[cfg(not(feature = "import-fastnlo"))] + fn run(&self) -> Result<()> { + Err(anyhow!( + "you need to install `pineappl` with feature `import-fastnlo`" + )) + } +} + +#[cfg(test)] +mod tests { + use assert_cmd::Command; + + const HELP_STR: &str = " +"; + + #[test] + fn help() { + Command::cargo_bin("pineappl") + .unwrap() + .args(&["import", "--help"]) + .assert() + .success() + .stdout(HELP_STR); + } +} diff --git a/pineappl_cli/src/import/fastnlo.cpp b/pineappl_cli/src/import/fastnlo.cpp new file mode 100644 index 00000000..119ae4e8 --- /dev/null +++ b/pineappl_cli/src/import/fastnlo.cpp @@ -0,0 +1,656 @@ +#include "pineappl_cli/src/import/fastnlo.hpp" + +#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 this_would_be_main(char const* in, char const* out) +{ + // 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); + pineappl_grid_delete(grids.at(0)); + + if (different) + { + return 1; + } + + return 0; +} diff --git a/pineappl_cli/src/import/fastnlo.hpp b/pineappl_cli/src/import/fastnlo.hpp new file mode 100644 index 00000000..290cda38 --- /dev/null +++ b/pineappl_cli/src/import/fastnlo.hpp @@ -0,0 +1,3 @@ +#include "rust/cxx.h" + +int this_would_be_main(char const* in, char const* out); 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), From 29aa3c568429868dedb1209d3b10616047591006 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 19 Feb 2022 14:08:34 +0100 Subject: [PATCH 02/26] Rename feature `import-fastnlo` to `fastnlo` --- pineappl_cli/Cargo.toml | 2 +- pineappl_cli/src/import.rs | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/pineappl_cli/Cargo.toml b/pineappl_cli/Cargo.toml index 86faa061..b11c12ef 100644 --- a/pineappl_cli/Cargo.toml +++ b/pineappl_cli/Cargo.toml @@ -41,4 +41,4 @@ path = "src/main.rs" rustc-args = [ "--cfg feature=\"docs-only\"" ] [features] -import-fastnlo = ["cxx", "cxx-build", "pkg-config"] +fastnlo = ["cxx", "cxx-build", "pkg-config"] diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index e7b02d40..216c181e 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -5,7 +5,7 @@ use std::path::PathBuf; // the indirection of two modules instead of one is needed to hide the C++ bridge behind the // feature gate -#[cfg(feature = "import-fastnlo")] +#[cfg(feature = "fastnlo")] mod fastnlo { #[cxx::bridge] pub mod ffi { @@ -29,7 +29,7 @@ pub struct Opts { } impl Subcommand for Opts { - #[cfg(feature = "import-fastnlo")] + #[cfg(feature = "fastnlo")] fn run(&self) -> Result<()> { use std::ffi::CString; use std::os::unix::ffi::OsStrExt; @@ -47,10 +47,10 @@ impl Subcommand for Opts { } } - #[cfg(not(feature = "import-fastnlo"))] + #[cfg(not(feature = "fastnlo"))] fn run(&self) -> Result<()> { Err(anyhow!( - "you need to install `pineappl` with feature `import-fastnlo`" + "you need to install `pineappl` with feature `fastnlo`" )) } } From 7cba20412564cb1e19b150a92ad93eefd52bf736 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 19 Feb 2022 14:09:02 +0100 Subject: [PATCH 03/26] Fix unit test --- pineappl_cli/src/import.rs | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 216c181e..f81720fe 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -59,7 +59,18 @@ impl Subcommand for Opts { mod tests { use assert_cmd::Command; - const HELP_STR: &str = " + const HELP_STR: &str = "pineappl-import +Converts fastNLO tables to PineAPPL grids + +USAGE: + pineappl import + +ARGS: + Path to the input grid + Path to the converted grid + +OPTIONS: + -h, --help Print help information "; #[test] From 589bd3812e30368c39fb215799a03d990d5b32e9 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 19 Feb 2022 14:11:48 +0100 Subject: [PATCH 04/26] Add missing `build.rs` --- pineappl_cli/build.rs | 62 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 pineappl_cli/build.rs diff --git a/pineappl_cli/build.rs b/pineappl_cli/build.rs new file mode 100644 index 00000000..07197360 --- /dev/null +++ b/pineappl_cli/build.rs @@ -0,0 +1,62 @@ +#[cfg(feature = "fastnlo")] +fn main() { + use std::process::Command; + + let mut bridge = cxx_build::bridge("src/import.rs"); + + // PineAPPL CAPI + let pineappl_capi = pkg_config::Config::new() + .atleast_version("0.5.0") + .probe("pineappl_capi") + .unwrap(); + + for include_path in pineappl_capi.include_paths { + bridge.include(include_path); + } + + for lib_path in pineappl_capi.link_paths { + println!("cargo:rustc-link-search={}", lib_path.to_str().unwrap()); + } + + for lib in pineappl_capi.libs { + println!("cargo:rustc-link-lib={}", lib); + } + + // fastNLO + let fnlo_include_path = String::from_utf8( + Command::new("fnlo-tk-config") + .arg("--incdir") + .output() + .unwrap() + .stdout, + ) + .unwrap(); + + bridge.include(fnlo_include_path); + + let fnlo_lib_path = String::from_utf8( + Command::new("fnlo-tk-config") + .arg("--libdir") + .output() + .unwrap() + .stdout, + ) + .unwrap(); + + println!("cargo:rustc-link-search={}", fnlo_lib_path); + + // TODO: why do I have to link statically? + println!("cargo:rustc-link-lib=static=fastnlotoolkit"); + println!("cargo:rustc-link-lib=z"); + + // bridging code + bridge.file("src/import/fastnlo.cpp").compile("cxx-bridge"); + + println!("cargo:rerun-if-changed=src/import.rs"); + println!("cargo:rerun-if-changed=src/import/import-ffirs"); + println!("cargo:rerun-if-changed=src/import/fastnlo.cpp"); + println!("cargo:rerun-if-changed=src/import/fastnlo.hpp"); +} + +#[cfg(not(feature = "fastnlo"))] +fn main() {} From 43298dc3338de29015d02852ea2ce92e835571f7 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 19 Feb 2022 18:35:06 +0100 Subject: [PATCH 05/26] Add preliminary fastNLO bindings --- pineappl_cli/src/import.rs | 158 ++++++++++++++++++++++++++++ pineappl_cli/src/import/fastnlo.cpp | 71 ++++++++++++- pineappl_cli/src/import/fastnlo.hpp | 29 ++++- 3 files changed, 255 insertions(+), 3 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index f81720fe..efe4ec3c 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -9,11 +9,169 @@ use std::path::PathBuf; mod fastnlo { #[cxx::bridge] pub mod ffi { + #[namespace = "fastNLO"] + #[repr(u32)] // TODO: this seems to work for me, but does it work everywhere? + enum EScaleFunctionalForm { + kScale1, + kScale2, + kQuadraticSum, + kQuadraticMean, + kQuadraticSumOver4, + kLinearMean, + kLinearSum, + kScaleMax, + kScaleMin, + kProd, + kS2plusS1half, + kPow4Sum, + kWgtAvg, + kS2plusS1fourth, + kExpProd2, + kExtern, + kConst, + } + + #[namespace = "fastNLO"] + #[repr(u32)] // TODO: this seems to work for me, but does it work everywhere? + enum ESMCalculation { + kFixedOrder = 0, + kThresholdCorrection = 1, + kElectroWeakCorrection = 2, + kNonPerturbativeCorrection = 3, + kContactInteraction = 10, + } + + #[namespace = "fastNLO"] + #[repr(u32)] // TODO: this seems to work for me, but does it work everywhere? + enum ESMOrder { + kLeading, + kNextToLeading, + kNextToNextToLeading, + } + + extern "C++" { + include!("fastnlotk/fastNLOConstants.h"); + + #[namespace = "fastNLO"] + type EScaleFunctionalForm; + + #[namespace = "fastNLO"] + type ESMCalculation; + + #[namespace = "fastNLO"] + type ESMOrder; + } + unsafe extern "C++" { include!("pineappl_cli/src/import/fastnlo.hpp"); unsafe fn this_would_be_main(input: *const c_char, output: *const c_char) -> i32; } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddBase.h"); + + type fastNLOCoeffAddBase; + + fn GetNObsBin(&self) -> i32; + fn GetNPDF(&self) -> i32; + fn GetNSubproc(&self) -> i32; + fn GetNpow(&self) -> i32; + + fn GetXNodes1(_: &fastNLOCoeffAddBase, _: i32) -> Vec; + fn GetXNodes2(_: &fastNLOCoeffAddBase, _: i32) -> Vec; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddFix.h"); + + type fastNLOCoeffAddFix; + + fn GetNPDFDim(&self) -> i32; + fn GetNevt(&self, _: i32, _: i32) -> f64; + 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; + + fn GetScaleNodes(_: &fastNLOCoeffAddFix, _: i32, _: i32) -> Vec; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOPDFLinearCombinations.h"); + + type fastNLOPDFLinearCombinations; + + fn CalcPDFLinearCombination( + _: &fastNLOPDFLinearCombinations, + _: &fastNLOCoeffAddBase, + _: &[f64], + _: &[f64], + _: bool, + ) -> Vec; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOLHAPDF.h"); + + type fastNLOLHAPDF; + + fn make_fastnlo_lhapdf_with_name_file_set( + name: &str, + lhapdf: &str, + set: i32, + ) -> UniquePtr; + } + + 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; + + fn GetCrossSection(_: Pin<&mut fastNLOReader>, _: bool) -> Vec; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOTable.h"); + + type fastNLOTable; + + fn GetCoeffTable(&self, _: i32) -> *mut fastNLOCoeffBase; + fn GetIpublunits(&self) -> i32; + fn GetNumDiffBin(&self) -> u32; + + fn GetBinSize(_: &fastNLOTable) -> Vec; + + // TODO: GetObsBinDimBounds(_: &fastNLOTable, _: u32, _: u32) -> (f64, f64); + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddFlex.h"); + + type fastNLOCoeffAddFlex; + + fn GetIXsectUnits(&self) -> i32; + fn GetNScaleDep(&self) -> i32; + fn GetPDFPDG(&self, _: i32) -> i32; + + fn GetScaleNodes1(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; + fn GetScaleNodes2(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; + + // TODO: GetSigmaTildes + } } } diff --git a/pineappl_cli/src/import/fastnlo.cpp b/pineappl_cli/src/import/fastnlo.cpp index 119ae4e8..5ce55ef1 100644 --- a/pineappl_cli/src/import/fastnlo.cpp +++ b/pineappl_cli/src/import/fastnlo.cpp @@ -1,10 +1,12 @@ -#include "pineappl_cli/src/import/fastnlo.hpp" - +#include #include #include +#include #include #include +#include +#include #include #include @@ -654,3 +656,68 @@ int this_would_be_main(char const* in, char const* out) return 0; } + +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) +{ + std::string arg0(name.begin(), name.end()); + std::string arg1(LHAPDFfile.begin(), LHAPDFfile.end()); + + return std::unique_ptr(new fastNLOLHAPDF(arg0, arg1, PDFSet)); +} + +rust::Vec GetCrossSection(fastNLOReader& reader, bool lNorm) +{ + 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)); +} diff --git a/pineappl_cli/src/import/fastnlo.hpp b/pineappl_cli/src/import/fastnlo.hpp index 290cda38..f2b52a59 100644 --- a/pineappl_cli/src/import/fastnlo.hpp +++ b/pineappl_cli/src/import/fastnlo.hpp @@ -1,3 +1,30 @@ -#include "rust/cxx.h" +#include + +#include +#include +#include +#include +#include +#include +#include +#include int this_would_be_main(char const* in, char const* out); + +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); + +std::unique_ptr make_fastnlo_lhapdf_with_name_file_set(rust::Str name, rust::Str LHAPDFfile, int PDFSet); + +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); From 15a91feaab521ab474179c4dd15089919e006118 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sun, 20 Feb 2022 10:26:24 +0100 Subject: [PATCH 06/26] Add help messages --- pineappl_cli/build.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pineappl_cli/build.rs b/pineappl_cli/build.rs index 07197360..85ac8ce1 100644 --- a/pineappl_cli/build.rs +++ b/pineappl_cli/build.rs @@ -8,7 +8,7 @@ fn main() { let pineappl_capi = pkg_config::Config::new() .atleast_version("0.5.0") .probe("pineappl_capi") - .unwrap(); + .expect("PineAPPL's C API not found, please install it"); for include_path in pineappl_capi.include_paths { bridge.include(include_path); @@ -27,7 +27,7 @@ fn main() { Command::new("fnlo-tk-config") .arg("--incdir") .output() - .unwrap() + .expect("fastNLO not found, please install it") .stdout, ) .unwrap(); From 6478f848fb15eb8a8867113a9da5b9f7f126c1ff Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 1 Mar 2022 11:01:12 +0100 Subject: [PATCH 07/26] Trim newlines from output --- pineappl_cli/build.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pineappl_cli/build.rs b/pineappl_cli/build.rs index 85ac8ce1..e1c11175 100644 --- a/pineappl_cli/build.rs +++ b/pineappl_cli/build.rs @@ -32,7 +32,7 @@ fn main() { ) .unwrap(); - bridge.include(fnlo_include_path); + bridge.include(fnlo_include_path.trim()); let fnlo_lib_path = String::from_utf8( Command::new("fnlo-tk-config") @@ -43,7 +43,7 @@ fn main() { ) .unwrap(); - println!("cargo:rustc-link-search={}", fnlo_lib_path); + println!("cargo:rustc-link-search={}", fnlo_lib_path.trim()); // TODO: why do I have to link statically? println!("cargo:rustc-link-lib=static=fastnlotoolkit"); From b712d44b65eb4a3e5a2a866bd1a029020bc1e0f9 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 1 Mar 2022 11:01:34 +0100 Subject: [PATCH 08/26] Remove superfluous line --- pineappl_cli/build.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/pineappl_cli/build.rs b/pineappl_cli/build.rs index e1c11175..55590a19 100644 --- a/pineappl_cli/build.rs +++ b/pineappl_cli/build.rs @@ -53,7 +53,6 @@ fn main() { bridge.file("src/import/fastnlo.cpp").compile("cxx-bridge"); println!("cargo:rerun-if-changed=src/import.rs"); - println!("cargo:rerun-if-changed=src/import/import-ffirs"); println!("cargo:rerun-if-changed=src/import/fastnlo.cpp"); println!("cargo:rerun-if-changed=src/import/fastnlo.hpp"); } From 0144ee5f3945c3e12fe9cd0aedf93b7b5aa04f4a Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 1 Mar 2022 11:02:11 +0100 Subject: [PATCH 09/26] Link fastNLO dynamically --- pineappl_cli/build.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pineappl_cli/build.rs b/pineappl_cli/build.rs index 55590a19..af81ff52 100644 --- a/pineappl_cli/build.rs +++ b/pineappl_cli/build.rs @@ -46,8 +46,7 @@ fn main() { println!("cargo:rustc-link-search={}", fnlo_lib_path.trim()); // TODO: why do I have to link statically? - println!("cargo:rustc-link-lib=static=fastnlotoolkit"); - println!("cargo:rustc-link-lib=z"); + println!("cargo:rustc-link-lib=fastnlotoolkit"); // bridging code bridge.file("src/import/fastnlo.cpp").compile("cxx-bridge"); From 305f40d104f34a1cc1c9313dae0f7ec4f8a4edbf Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 8 Mar 2022 09:43:18 +0100 Subject: [PATCH 10/26] Move C++ code to new crate and migrate code - the C++ code is now contained in the new crate `pineappl_fastnlo` which is a very thin wrapper - migrated most of the code to Rust, the code for the flex has yet to be writen (commented and still in C++) - unfortunately there is still a bug in the converter; the bug is most likely in `convert_coeff_add_fix`, `create_lumi` is most likely correct --- Cargo.toml | 1 + pineappl_cli/Cargo.toml | 7 +- pineappl_cli/build.rs | 60 --- pineappl_cli/src/import.rs | 710 +++++++++++++++++++++------ pineappl_cli/src/import/fastnlo.cpp | 723 ---------------------------- pineappl_cli/src/import/fastnlo.hpp | 30 -- pineappl_fastnlo/Cargo.toml | 17 + pineappl_fastnlo/build.rs | 34 ++ pineappl_fastnlo/src/fastnlo.cpp | 140 ++++++ pineappl_fastnlo/src/fastnlo.hpp | 70 +++ pineappl_fastnlo/src/lib.rs | 182 +++++++ 11 files changed, 1014 insertions(+), 960 deletions(-) delete mode 100644 pineappl_cli/build.rs delete mode 100644 pineappl_cli/src/import/fastnlo.cpp delete mode 100644 pineappl_cli/src/import/fastnlo.hpp create mode 100644 pineappl_fastnlo/Cargo.toml create mode 100644 pineappl_fastnlo/build.rs create mode 100644 pineappl_fastnlo/src/fastnlo.cpp create mode 100644 pineappl_fastnlo/src/fastnlo.hpp create mode 100644 pineappl_fastnlo/src/lib.rs 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/pineappl_cli/Cargo.toml b/pineappl_cli/Cargo.toml index b11c12ef..90e1e5e1 100644 --- a/pineappl_cli/Cargo.toml +++ b/pineappl_cli/Cargo.toml @@ -22,13 +22,10 @@ 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" -[build-dependencies] -cxx-build = { optional = true, version = "1.0.65" } -pkg-config = { optional = true, version = "0.3.24" } - [dev-dependencies] assert_cmd = "2.0.2" assert_fs = "1.0.6" @@ -41,4 +38,4 @@ path = "src/main.rs" rustc-args = [ "--cfg feature=\"docs-only\"" ] [features] -fastnlo = ["cxx", "cxx-build", "pkg-config"] +fastnlo = ["pineappl_fastnlo"] diff --git a/pineappl_cli/build.rs b/pineappl_cli/build.rs deleted file mode 100644 index af81ff52..00000000 --- a/pineappl_cli/build.rs +++ /dev/null @@ -1,60 +0,0 @@ -#[cfg(feature = "fastnlo")] -fn main() { - use std::process::Command; - - let mut bridge = cxx_build::bridge("src/import.rs"); - - // PineAPPL CAPI - let pineappl_capi = pkg_config::Config::new() - .atleast_version("0.5.0") - .probe("pineappl_capi") - .expect("PineAPPL's C API not found, please install it"); - - for include_path in pineappl_capi.include_paths { - bridge.include(include_path); - } - - for lib_path in pineappl_capi.link_paths { - println!("cargo:rustc-link-search={}", lib_path.to_str().unwrap()); - } - - for lib in pineappl_capi.libs { - println!("cargo:rustc-link-lib={}", lib); - } - - // fastNLO - let fnlo_include_path = String::from_utf8( - Command::new("fnlo-tk-config") - .arg("--incdir") - .output() - .expect("fastNLO not found, please install it") - .stdout, - ) - .unwrap(); - - bridge.include(fnlo_include_path.trim()); - - let fnlo_lib_path = String::from_utf8( - Command::new("fnlo-tk-config") - .arg("--libdir") - .output() - .unwrap() - .stdout, - ) - .unwrap(); - - println!("cargo:rustc-link-search={}", fnlo_lib_path.trim()); - - // TODO: why do I have to link statically? - println!("cargo:rustc-link-lib=fastnlotoolkit"); - - // bridging code - bridge.file("src/import/fastnlo.cpp").compile("cxx-bridge"); - - println!("cargo:rerun-if-changed=src/import.rs"); - println!("cargo:rerun-if-changed=src/import/fastnlo.cpp"); - println!("cargo:rerun-if-changed=src/import/fastnlo.hpp"); -} - -#[cfg(not(feature = "fastnlo"))] -fn main() {} diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index efe4ec3c..6c2b1547 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -1,176 +1,609 @@ -use super::helpers::Subcommand; +use super::helpers::{self, Subcommand}; use anyhow::{anyhow, Result}; use clap::{Parser, ValueHint}; use std::path::PathBuf; -// the indirection of two modules instead of one is needed to hide the C++ bridge behind the -// feature gate #[cfg(feature = "fastnlo")] mod fastnlo { - #[cxx::bridge] - pub mod ffi { - #[namespace = "fastNLO"] - #[repr(u32)] // TODO: this seems to work for me, but does it work everywhere? - enum EScaleFunctionalForm { - kScale1, - kScale2, - kQuadraticSum, - kQuadraticMean, - kQuadraticSumOver4, - kLinearMean, - kLinearSum, - kScaleMax, - kScaleMin, - kProd, - kS2plusS1half, - kPow4Sum, - kWgtAvg, - kS2plusS1fourth, - kExpProd2, - kExtern, - kConst, + use super::*; + use lhapdf::Pdf; + 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, + fastNLOPDFLinearCombinations, ESMCalculation, ESMOrder, EScaleFunctionalForm, + }; + use std::convert::{TryFrom, TryInto}; + use std::path::Path; + use std::ptr; + + fn pid_to_pdg_id(pid: i32) -> i32 { + if pid >= -6 && pid <= 6 { + if pid == 0 { + 21 + } else { + pid + } + } else { + unreachable!() } + } - #[namespace = "fastNLO"] - #[repr(u32)] // TODO: this seems to work for me, but does it work everywhere? - enum ESMCalculation { - kFixedOrder = 0, - kThresholdCorrection = 1, - kElectroWeakCorrection = 2, - kNonPerturbativeCorrection = 3, - kContactInteraction = 10, + 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)); } - #[namespace = "fastNLO"] - #[repr(u32)] // TODO: this seems to work for me, but does it work everywhere? - enum ESMOrder { - kLeading, - kNextToLeading, - kNextToNextToLeading, - } + // if the PDF coefficient vector was empty, we must reconstruct the lumi function + if lumis.is_empty() { + let nsubproc = table.GetNSubproc().try_into().unwrap(); - extern "C++" { - include!("fastnlotk/fastNLOConstants.h"); + let mut xfx1 = [0.0; 13]; + let mut xfx2 = [0.0; 13]; - #[namespace = "fastNLO"] - type EScaleFunctionalForm; + let mut entries = Vec::new(); + entries.resize(nsubproc, Vec::new()); - #[namespace = "fastNLO"] - type ESMCalculation; + for a in 0..13 { + xfx1[a] = 1.0; - #[namespace = "fastNLO"] - type ESMOrder; - } + for b in 0..13 { + xfx2[b] = 1.0; - unsafe extern "C++" { - include!("pineappl_cli/src/import/fastnlo.hpp"); + let lumi = ffi::CalcPDFLinearCombination(comb, table, &xfx1, &xfx2, false); - unsafe fn this_would_be_main(input: *const c_char, output: *const c_char) -> i32; - } + assert!(lumi.len() == nsubproc); - unsafe extern "C++" { - include!("fastnlotk/fastNLOCoeffAddBase.h"); + for i in 0..lumi.len() { + if lumi[i] == 0.0 { + continue; + } - type fastNLOCoeffAddBase; + let ap = pid_to_pdg_id(i32::try_from(a).unwrap() - 6); + let bp = pid_to_pdg_id(i32::try_from(b).unwrap() - 6); - fn GetNObsBin(&self) -> i32; - fn GetNPDF(&self) -> i32; - fn GetNSubproc(&self) -> i32; - fn GetNpow(&self) -> i32; + entries[i].push((ap, bp, lumi[i])); + } - fn GetXNodes1(_: &fastNLOCoeffAddBase, _: i32) -> Vec; - fn GetXNodes2(_: &fastNLOCoeffAddBase, _: i32) -> Vec; - } + xfx2[b] = 0.0; + } - unsafe extern "C++" { - include!("fastnlotk/fastNLOCoeffAddFix.h"); + xfx1[a] = 0.0; + } - type fastNLOCoeffAddFix; + lumis = entries.into_iter().map(LumiEntry::new).collect(); + } - fn GetNPDFDim(&self) -> i32; - fn GetNevt(&self, _: i32, _: i32) -> f64; - 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; + lumis + } - fn GetScaleNodes(_: &fastNLOCoeffAddFix, _: i32, _: i32) -> Vec; + 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 n_obs_bin = table_as_add_base.GetNObsBin(); + let n_subproc = table_as_add_base.GetNSubproc(); + let total_scalevars = table.GetTotalScalevars(); + let total_scalenodes: usize = table.GetTotalScalenodes().try_into().unwrap(); + + for obs in 0..n_obs_bin { + 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..n_subproc { + let factor = table.GetNevt(obs, subproc); + + for j in 0..total_scalevars { + // 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); + + let mut non_zero_subgrid = false; + + 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 { + non_zero_subgrid = true; + array[[mu2_slice, ix1, ix2]] = + 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 non_zero_subgrid { + grid.set_subgrid( + 0, + obs.try_into().unwrap(), + subproc.try_into().unwrap(), + ImportOnlySubgridV2::new( + array, + mu2_values, + x1_values.clone(), + x2_values.clone(), + ) + .into(), + ); + } + } + } } - unsafe extern "C++" { - include!("fastnlotk/fastNLOPDFLinearCombinations.h"); - - type fastNLOPDFLinearCombinations; + grid + } - fn CalcPDFLinearCombination( - _: &fastNLOPDFLinearCombinations, - _: &fastNLOCoeffAddBase, - _: &[f64], - _: &[f64], - _: bool, - ) -> Vec; - } + fn convert_coeff_add_flex( + _table: &fastNLOCoeffAddFlex, + _comb: &fastNLOPDFLinearCombinations, + _mur_ff: EScaleFunctionalForm, + _muf_ff: EScaleFunctionalForm, + _bins: usize, + _alpha: u32, + _ipub_units: i32, + ) -> Grid { + //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; + todo!() + } - unsafe extern "C++" { - include!("fastnlotk/fastNLOLHAPDF.h"); + pub fn convert_fastnlo_table(input: &Path, pdfset: &str) -> Result { + // TODO: read this from an argument + let alpha = 0; - type fastNLOLHAPDF; + let mut file = + ffi::make_fastnlo_lhapdf_with_name_file_set(input.to_str().unwrap(), pdfset, 0); + let file_as_reader = ffi::downcast_lhapdf_to_reader(file.as_ref().unwrap()); + let file_as_table = ffi::downcast_lhapdf_to_table(file.as_ref().unwrap()); - fn make_fastnlo_lhapdf_with_name_file_set( - name: &str, - lhapdf: &str, - set: i32, - ) -> UniquePtr; - } + let id_lo = file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kLeading); + let id_nlo = file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kNextToLeading); + let id_nnlo = + file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kNextToNextToLeading); - unsafe extern "C++" { - include!("fastnlotk/fastNLOCoeffBase.h"); + let mut ids = Vec::new(); - type fastNLOCoeffBase; + if id_lo >= 0 { + ids.push(id_lo); } - unsafe extern "C++" { - include!("fastnlotk/fastNLOReader.h"); - - type fastNLOReader; + if id_nlo >= 0 { + ids.push(id_nlo); + } - fn GetMuRFunctionalForm(&self) -> EScaleFunctionalForm; - fn GetMuFFunctionalForm(&self) -> EScaleFunctionalForm; - fn ContrId(&self, _: ESMCalculation, _: ESMOrder) -> i32; + if id_nnlo >= 0 { + ids.push(id_nnlo); + } - fn GetCrossSection(_: Pin<&mut fastNLOReader>, _: bool) -> Vec; + let normalizations = ffi::GetBinSize(file_as_table); + let bins = normalizations.len(); + + let mut grids = Vec::new(); + + for id in ids { + let coeff_table = file_as_table.GetCoeffTable(id); + assert_ne!(coeff_table, ptr::null_mut()); + 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!(); + } + } } - unsafe extern "C++" { - include!("fastnlotk/fastNLOTable.h"); + let mut result = grids.remove(0); + for grid in grids { + result.merge(grid)?; + } - type fastNLOTable; + result.scale_by_order(0.5 / f64::acos(-1.0), 1.0, 1.0, 1.0, 1.0); + result.optimize(); - fn GetCoeffTable(&self, _: i32) -> *mut fastNLOCoeffBase; - fn GetIpublunits(&self) -> i32; - fn GetNumDiffBin(&self) -> u32; + let dimensions: usize = file_as_table.GetNumDiffBin().try_into().unwrap(); + let mut limits = Vec::new(); + limits.reserve(2 * dimensions * bins); - fn GetBinSize(_: &fastNLOTable) -> Vec; + 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(), + ); - // TODO: GetObsBinDimBounds(_: &fastNLOTable, _: u32, _: u32) -> (f64, f64); + limits.push((pair.first, pair.second)); + } } - unsafe extern "C++" { - include!("fastnlotk/fastNLOCoeffAddFlex.h"); - - type fastNLOCoeffAddFlex; - - fn GetIXsectUnits(&self) -> i32; - fn GetNScaleDep(&self) -> i32; - fn GetPDFPDG(&self, _: i32) -> i32; - - fn GetScaleNodes1(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; - fn GetScaleNodes2(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; + result.set_remapper(BinRemapper::new(normalizations.clone(), limits).unwrap())?; + + let file_as_reader_mut = ffi::downcast_lhapdf_to_reader_mut(file.as_mut().unwrap()); + let results = ffi::GetCrossSection(file_as_reader_mut, false); + let pdf = pdfset + .parse() + .map_or_else(|_| Pdf::with_setname_and_member(pdfset, 0), Pdf::with_lhaid); + let other_results = helpers::convolute(&result, &pdf, &[], &[], &[], 1); + + let mut different = false; + + for (i, (one, mut two)) in results + .into_iter() + .zip(other_results.into_iter()) + .enumerate() + { + two *= normalizations[i]; + + // catches the case where both results are zero + if one == two { + continue; + } + + if (two / one - 1.0).abs() > 1e-10 { + println!( + ">>> fastNLO: {} PineAPPL: {} fN/P: {} P/fN: {}", + one, + two, + one / two, + two / one + ); + different = true; + } else { + println!(">>> Success!"); + } + } - // TODO: GetSigmaTildes + if different { + Err(anyhow!("grids are different")) + } else { + Ok(result) } } } @@ -184,25 +617,17 @@ pub struct Opts { /// 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, } impl Subcommand for Opts { #[cfg(feature = "fastnlo")] fn run(&self) -> Result<()> { - use std::ffi::CString; - use std::os::unix::ffi::OsStrExt; - - // TODO: this is a correct way of converting strings, but only works on UNIX - let input = CString::new(self.input.as_path().as_os_str().as_bytes()).unwrap(); - let output = CString::new(self.output.as_path().as_os_str().as_bytes()).unwrap(); + let grid = fastnlo::convert_fastnlo_table(&self.input, &self.pdfset)?; - let errcode = unsafe { fastnlo::ffi::this_would_be_main(input.as_ptr(), output.as_ptr()) }; - - if errcode == 0 { - Ok(()) - } else { - Err(anyhow!("something went wrong")) - } + helpers::write_grid(&self.output, &grid) } #[cfg(not(feature = "fastnlo"))] @@ -221,11 +646,12 @@ mod tests { Converts fastNLO tables to PineAPPL grids USAGE: - pineappl import + pineappl import 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: -h, --help Print help information diff --git a/pineappl_cli/src/import/fastnlo.cpp b/pineappl_cli/src/import/fastnlo.cpp deleted file mode 100644 index 5ce55ef1..00000000 --- a/pineappl_cli/src/import/fastnlo.cpp +++ /dev/null @@ -1,723 +0,0 @@ -#include -#include -#include -#include -#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 this_would_be_main(char const* in, char const* out) -{ - // 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); - pineappl_grid_delete(grids.at(0)); - - if (different) - { - return 1; - } - - return 0; -} - -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) -{ - std::string arg0(name.begin(), name.end()); - std::string arg1(LHAPDFfile.begin(), LHAPDFfile.end()); - - return std::unique_ptr(new fastNLOLHAPDF(arg0, arg1, PDFSet)); -} - -rust::Vec GetCrossSection(fastNLOReader& reader, bool lNorm) -{ - 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)); -} diff --git a/pineappl_cli/src/import/fastnlo.hpp b/pineappl_cli/src/import/fastnlo.hpp deleted file mode 100644 index f2b52a59..00000000 --- a/pineappl_cli/src/import/fastnlo.hpp +++ /dev/null @@ -1,30 +0,0 @@ -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -int this_would_be_main(char const* in, char const* out); - -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); - -std::unique_ptr make_fastnlo_lhapdf_with_name_file_set(rust::Str name, rust::Str LHAPDFfile, int PDFSet); - -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); 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..990e47cf --- /dev/null +++ b/pineappl_fastnlo/src/fastnlo.cpp @@ -0,0 +1,140 @@ +#include "pineappl_fastnlo/src/fastnlo.hpp" + +#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) +{ + std::string arg0(name.begin(), name.end()); + std::string arg1(LHAPDFfile.begin(), LHAPDFfile.end()); + + return std::unique_ptr(new fastNLOLHAPDF(arg0, arg1, PDFSet)); +} + +rust::Vec GetCrossSection(fastNLOReader& reader, bool lNorm) +{ + 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; +} + +fastNLOCoeffAddBase const& downcast_coeff_add_fix_to_base(fastNLOCoeffAddFix 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..f14d429b --- /dev/null +++ b/pineappl_fastnlo/src/fastnlo.hpp @@ -0,0 +1,70 @@ +#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 +); + +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); + +fastNLOCoeffAddBase const& downcast_coeff_add_fix_to_base(fastNLOCoeffAddFix 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..4916278b --- /dev/null +++ b/pineappl_fastnlo/src/lib.rs @@ -0,0 +1,182 @@ +#[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; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddBase.h"); + + type fastNLOCoeffAddBase; + + 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 GetNevt(&self, _: i32, _: i32) -> f64; + 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; + + fn CalcPDFLinearCombination( + _: &fastNLOPDFLinearCombinations, + _: &fastNLOCoeffAddBase, + _: &[f64], + _: &[f64], + _: bool, + ) -> Vec; + } + + 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 GetNumDiffBin(&self) -> u32; + } + + unsafe extern "C++" { + include!("fastnlotk/fastNLOCoeffAddFlex.h"); + + type fastNLOCoeffAddFlex; + + fn GetIXsectUnits(&self) -> i32; + fn GetNScaleDep(&self) -> i32; + fn GetPDFPDG(&self, _: i32) -> i32; + + // TODO: GetSigmaTildes + } + + unsafe extern "C++" { + include!("pineappl_fastnlo/src/fastnlo.hpp"); + + fn GetBinSize(_: &fastNLOTable) -> Vec; + fn GetCrossSection(_: Pin<&mut fastNLOReader>, _: bool) -> Vec; + 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 GetObsBinDimBounds(_: &fastNLOTable, _: u32, _: u32) -> pair_double_double; + + fn downcast_coeff_add_fix_to_base(_: &fastNLOCoeffAddFix) -> &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; + fn make_fastnlo_lhapdf_with_name_file_set( + name: &str, + lhapdf: &str, + set: i32, + ) -> UniquePtr; + + unsafe fn dynamic_cast_coeff_add_fix(_: *mut fastNLOCoeffBase) -> *mut fastNLOCoeffAddFix; + unsafe fn dynamic_cast_coeff_add_flex(_: *mut fastNLOCoeffBase) + -> *mut fastNLOCoeffAddFlex; + } +} From 059491d3d3613ca18318de83744afdcc465a31d5 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 12 Mar 2022 09:56:07 +0100 Subject: [PATCH 11/26] Fix bug in conversion No idea why one has to transpose the indices, something isn't properly understood --- pineappl_cli/src/import.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 6c2b1547..a446c63c 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -189,7 +189,7 @@ mod fastnlo { if value != 0.0 { non_zero_subgrid = true; - array[[mu2_slice, ix1, ix2]] = + array[[mu2_slice, ix2, ix1]] = value / factor * x1_values[ix1] * x2_values[ix2]; } From 2f1bf92da356bafa78f0e8df92c11a39c1b50c7a Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 12 Mar 2022 09:57:35 +0100 Subject: [PATCH 12/26] Streamline code a bit --- pineappl_cli/src/import.rs | 44 +++++++++++++------------------------- 1 file changed, 15 insertions(+), 29 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index a446c63c..2ad51574 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -18,6 +18,7 @@ mod fastnlo { fastNLOPDFLinearCombinations, ESMCalculation, ESMOrder, EScaleFunctionalForm, }; use std::convert::{TryFrom, TryInto}; + use std::f64::consts::TAU; use std::path::Path; use std::ptr; @@ -81,15 +82,11 @@ mod fastnlo { assert!(lumi.len() == nsubproc); - for i in 0..lumi.len() { - if lumi[i] == 0.0 { - continue; - } - + 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, lumi[i])); + entries[i].push((ap, bp, l)); } xfx2[b] = 0.0; @@ -468,33 +465,21 @@ mod fastnlo { todo!() } - pub fn convert_fastnlo_table(input: &Path, pdfset: &str) -> Result { - // TODO: read this from an argument - let alpha = 0; - + pub fn convert_fastnlo_table(input: &Path, pdfset: &str, alpha: u32) -> Result { let mut file = ffi::make_fastnlo_lhapdf_with_name_file_set(input.to_str().unwrap(), pdfset, 0); let file_as_reader = ffi::downcast_lhapdf_to_reader(file.as_ref().unwrap()); let file_as_table = ffi::downcast_lhapdf_to_table(file.as_ref().unwrap()); - let id_lo = file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kLeading); - let id_nlo = file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kNextToLeading); - let id_nnlo = - file_as_reader.ContrId(ESMCalculation::kFixedOrder, ESMOrder::kNextToNextToLeading); - - let mut ids = Vec::new(); - - if id_lo >= 0 { - ids.push(id_lo); - } - - if id_nlo >= 0 { - ids.push(id_nlo); - } - - if id_nnlo >= 0 { - ids.push(id_nnlo); - } + 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 normalizations = ffi::GetBinSize(file_as_table); let bins = normalizations.len(); @@ -544,7 +529,7 @@ mod fastnlo { result.merge(grid)?; } - result.scale_by_order(0.5 / f64::acos(-1.0), 1.0, 1.0, 1.0, 1.0); + 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(); @@ -583,6 +568,7 @@ mod fastnlo { // catches the case where both results are zero if one == two { + println!(">>> Success!"); continue; } From b8c68a403565201f4e09431714e94dd2d57f0860 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 12 Mar 2022 10:03:47 +0100 Subject: [PATCH 13/26] Allow setting the exponent of alpha --- pineappl_cli/src/import.rs | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 2ad51574..549e89f6 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -606,12 +606,15 @@ pub struct Opts { /// 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)] + alpha: u32, } impl Subcommand for Opts { #[cfg(feature = "fastnlo")] fn run(&self) -> Result<()> { - let grid = fastnlo::convert_fastnlo_table(&self.input, &self.pdfset)?; + let grid = fastnlo::convert_fastnlo_table(&self.input, &self.pdfset, self.alpha)?; helpers::write_grid(&self.output, &grid) } From c43977af17b999358458c169b27fdbf084272b2b Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 12 Mar 2022 10:27:48 +0100 Subject: [PATCH 14/26] Disable normalization if not needed --- pineappl_cli/src/import.rs | 9 ++++++--- pineappl_fastnlo/src/lib.rs | 2 ++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 549e89f6..9272f769 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -481,9 +481,7 @@ mod fastnlo { .filter(|&id| id >= 0) .collect(); - let normalizations = ffi::GetBinSize(file_as_table); - let bins = normalizations.len(); - + let bins: usize = file_as_table.GetNObsBin().try_into().unwrap(); let mut grids = Vec::new(); for id in ids { @@ -534,6 +532,11 @@ mod fastnlo { 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 { diff --git a/pineappl_fastnlo/src/lib.rs b/pineappl_fastnlo/src/lib.rs index 4916278b..08ca5446 100644 --- a/pineappl_fastnlo/src/lib.rs +++ b/pineappl_fastnlo/src/lib.rs @@ -133,7 +133,9 @@ pub mod ffi { 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++" { From 7a4b51c124f564c29577d942f143ca96927c3d1c Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Sat, 12 Mar 2022 10:29:37 +0100 Subject: [PATCH 15/26] Remove `fnlo2pine` from examples This is now part of `pineappl import ...` --- examples/fnlo2pine/fnlo2pine.cpp | 660 ------------------------------- examples/fnlo2pine/meson.build | 26 -- 2 files changed, 686 deletions(-) delete mode 100644 examples/fnlo2pine/fnlo2pine.cpp delete mode 100644 examples/fnlo2pine/meson.build 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 -) From e8eb621322c1130b82921b883696a62ac729ddc1 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 11:31:56 +0100 Subject: [PATCH 16/26] Make code more idiomatic --- pineappl_cli/src/import.rs | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 9272f769..ec06ec36 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -23,14 +23,10 @@ mod fastnlo { use std::ptr; fn pid_to_pdg_id(pid: i32) -> i32 { - if pid >= -6 && pid <= 6 { - if pid == 0 { - 21 - } else { - pid - } - } else { - unreachable!() + match pid { + -6..=-1 | 1..=6 => pid, + 0 => 21, + _ => unimplemented!(), } } @@ -121,12 +117,9 @@ mod fastnlo { SubgridParams::default(), ); - let n_obs_bin = table_as_add_base.GetNObsBin(); - let n_subproc = table_as_add_base.GetNSubproc(); - let total_scalevars = table.GetTotalScalevars(); let total_scalenodes: usize = table.GetTotalScalenodes().try_into().unwrap(); - for obs in 0..n_obs_bin { + 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? @@ -136,10 +129,10 @@ mod fastnlo { ffi::GetXNodes2(table_as_add_base, obs) }; - for subproc in 0..n_subproc { - let factor = table.GetNevt(obs, subproc); + for subproc in 0..table_as_add_base.GetNSubproc() { + let factor = table_as_add_base.GetNevt(obs, subproc); - for j in 0..total_scalevars { + 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; @@ -160,8 +153,6 @@ mod fastnlo { // TODO: figure out what the general case is supposed to be assert_eq!(j, 0); - let mut non_zero_subgrid = false; - for mu2_slice in 0..total_scalenodes { let mut ix1: usize = 0; let mut ix2: usize = 0; @@ -185,7 +176,6 @@ mod fastnlo { ); if value != 0.0 { - non_zero_subgrid = true; array[[mu2_slice, ix2, ix1]] = value / factor * x1_values[ix1] * x2_values[ix2]; } @@ -210,7 +200,7 @@ mod fastnlo { } } - if non_zero_subgrid { + if !array.is_empty() { grid.set_subgrid( 0, obs.try_into().unwrap(), From fc93e17f32d7c8413db0ceb6a5b844f16cc1daca Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 11:33:15 +0100 Subject: [PATCH 17/26] Enable conversion of flexible-scale fastNLO grids --- pineappl_cli/src/import.rs | 413 ++++++++++++++----------------- pineappl_fastnlo/src/fastnlo.cpp | 22 ++ pineappl_fastnlo/src/fastnlo.hpp | 14 ++ pineappl_fastnlo/src/lib.rs | 15 +- 4 files changed, 231 insertions(+), 233 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index ec06ec36..d1b2c73c 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -6,6 +6,7 @@ use std::path::PathBuf; #[cfg(feature = "fastnlo")] mod fastnlo { use super::*; + use itertools::Itertools; use lhapdf::Pdf; use pineappl::bin::BinRemapper; use pineappl::grid::{Grid, Order}; @@ -222,237 +223,189 @@ mod fastnlo { } fn convert_coeff_add_flex( - _table: &fastNLOCoeffAddFlex, - _comb: &fastNLOPDFLinearCombinations, - _mur_ff: EScaleFunctionalForm, - _muf_ff: EScaleFunctionalForm, - _bins: usize, - _alpha: u32, - _ipub_units: i32, + table: &fastNLOCoeffAddFlex, + comb: &fastNLOPDFLinearCombinations, + mur_ff: EScaleFunctionalForm, + muf_ff: EScaleFunctionalForm, + bins: usize, + alpha: u32, + ipub_units: i32, ) -> Grid { - //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; - todo!() + 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(input: &Path, pdfset: &str, alpha: u32) -> Result { diff --git a/pineappl_fastnlo/src/fastnlo.cpp b/pineappl_fastnlo/src/fastnlo.cpp index 990e47cf..dc34bd66 100644 --- a/pineappl_fastnlo/src/fastnlo.cpp +++ b/pineappl_fastnlo/src/fastnlo.cpp @@ -91,11 +91,33 @@ rust::Vec GetPDFCoeff(fastNLOCoeffAddBase const& coeffs, std::size 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; diff --git a/pineappl_fastnlo/src/fastnlo.hpp b/pineappl_fastnlo/src/fastnlo.hpp index f14d429b..d2f21242 100644 --- a/pineappl_fastnlo/src/fastnlo.hpp +++ b/pineappl_fastnlo/src/fastnlo.hpp @@ -45,8 +45,22 @@ 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); diff --git a/pineappl_fastnlo/src/lib.rs b/pineappl_fastnlo/src/lib.rs index 08ca5446..891b2a5e 100644 --- a/pineappl_fastnlo/src/lib.rs +++ b/pineappl_fastnlo/src/lib.rs @@ -68,6 +68,7 @@ pub mod ffi { type fastNLOCoeffAddBase; + fn GetNevt(&self, _: i32, _: i32) -> f64; fn GetNObsBin(&self) -> i32; fn GetNPDF(&self) -> i32; fn GetNSubproc(&self) -> i32; @@ -80,7 +81,6 @@ pub mod ffi { type fastNLOCoeffAddFix; fn GetNPDFDim(&self) -> i32; - fn GetNevt(&self, _: i32, _: i32) -> f64; fn GetNxmax(&self, _: i32) -> i32; fn GetNxtot2(&self, _: i32) -> i32; fn GetScaleFactor(&self, _: i32) -> f64; @@ -146,8 +146,6 @@ pub mod ffi { fn GetIXsectUnits(&self) -> i32; fn GetNScaleDep(&self) -> i32; fn GetPDFPDG(&self, _: i32) -> i32; - - // TODO: GetSigmaTildes } unsafe extern "C++" { @@ -163,8 +161,19 @@ pub mod ffi { fn GetXNodes1(_: &fastNLOCoeffAddBase, _: i32) -> Vec; fn GetXNodes2(_: &fastNLOCoeffAddBase, _: i32) -> Vec; fn GetObsBinDimBounds(_: &fastNLOTable, _: u32, _: u32) -> pair_double_double; + fn GetSigmaTilde( + _: &fastNLOCoeffAddFlex, + _: usize, + _: usize, + _: usize, + _: usize, + _: usize, + _: i32, + ) -> f64; + fn GetNx(_: &fastNLOCoeffAddFlex, _: usize) -> usize; 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; From e8a05237f2b6f128d9c5da6497c0076b3ddd78fa Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 11:39:53 +0100 Subject: [PATCH 18/26] Fix help message --- pineappl_cli/src/import.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index d1b2c73c..61721454 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -581,12 +581,13 @@ mod tests { Converts fastNLO tables to PineAPPL grids USAGE: - pineappl import + pineappl import [ALPHA] 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 + LO coupling power in alpha [default: 0] OPTIONS: -h, --help Print help information From 354c41c1194d4e7bd552ecb5fe1473fdfe446292 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 11:41:57 +0100 Subject: [PATCH 19/26] Update changelog --- CHANGELOG.md | 6 ++++++ 1 file changed, 6 insertions(+) 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 From 9b3de8ffabb06d5cc45b96a7f48d45a59da689bf Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 15:35:15 +0100 Subject: [PATCH 20/26] Make `alpha` a switch instead of a parameter --- pineappl_cli/src/import.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 61721454..c05dee42 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -553,7 +553,7 @@ pub struct Opts { #[clap(validator = helpers::validate_pdfset)] pdfset: String, /// LO coupling power in alpha. - #[clap(default_value_t = 0)] + #[clap(default_value_t = 0, long)] alpha: u32, } @@ -581,16 +581,16 @@ mod tests { Converts fastNLO tables to PineAPPL grids USAGE: - pineappl import [ALPHA] + 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 - LO coupling power in alpha [default: 0] OPTIONS: - -h, --help Print help information + --alpha LO coupling power in alpha [default: 0] + -h, --help Print help information "; #[test] From 4723c7dea0586bee82daddc1a38a4fc9e052c1e8 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 15:37:21 +0100 Subject: [PATCH 21/26] Move code and use tables to print output --- pineappl_cli/src/import.rs | 118 +++++++++++++++++++++---------------- 1 file changed, 66 insertions(+), 52 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index c05dee42..4864d8ae 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -7,7 +7,6 @@ use std::path::PathBuf; mod fastnlo { use super::*; use itertools::Itertools; - use lhapdf::Pdf; use pineappl::bin::BinRemapper; use pineappl::grid::{Grid, Order}; use pineappl::import_only_subgrid::ImportOnlySubgridV2; @@ -15,12 +14,11 @@ mod fastnlo { use pineappl::sparse_array3::SparseArray3; use pineappl::subgrid::{Mu2, SubgridParams}; use pineappl_fastnlo::ffi::{ - self, fastNLOCoeffAddBase, fastNLOCoeffAddFix, fastNLOCoeffAddFlex, + self, fastNLOCoeffAddBase, fastNLOCoeffAddFix, fastNLOCoeffAddFlex, fastNLOLHAPDF, fastNLOPDFLinearCombinations, ESMCalculation, ESMOrder, EScaleFunctionalForm, }; use std::convert::{TryFrom, TryInto}; use std::f64::consts::TAU; - use std::path::Path; use std::ptr; fn pid_to_pdg_id(pid: i32) -> i32 { @@ -408,11 +406,9 @@ mod fastnlo { grid } - pub fn convert_fastnlo_table(input: &Path, pdfset: &str, alpha: u32) -> Result { - let mut file = - ffi::make_fastnlo_lhapdf_with_name_file_set(input.to_str().unwrap(), pdfset, 0); - let file_as_reader = ffi::downcast_lhapdf_to_reader(file.as_ref().unwrap()); - let file_as_table = ffi::downcast_lhapdf_to_table(file.as_ref().unwrap()); + 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), @@ -429,7 +425,11 @@ mod fastnlo { for id in ids { let coeff_table = file_as_table.GetCoeffTable(id); - assert_ne!(coeff_table, ptr::null_mut()); + + 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); @@ -496,47 +496,7 @@ mod fastnlo { result.set_remapper(BinRemapper::new(normalizations.clone(), limits).unwrap())?; - let file_as_reader_mut = ffi::downcast_lhapdf_to_reader_mut(file.as_mut().unwrap()); - let results = ffi::GetCrossSection(file_as_reader_mut, false); - let pdf = pdfset - .parse() - .map_or_else(|_| Pdf::with_setname_and_member(pdfset, 0), Pdf::with_lhaid); - let other_results = helpers::convolute(&result, &pdf, &[], &[], &[], 1); - - let mut different = false; - - for (i, (one, mut two)) in results - .into_iter() - .zip(other_results.into_iter()) - .enumerate() - { - two *= normalizations[i]; - - // catches the case where both results are zero - if one == two { - println!(">>> Success!"); - continue; - } - - if (two / one - 1.0).abs() > 1e-10 { - println!( - ">>> fastNLO: {} PineAPPL: {} fN/P: {} P/fN: {}", - one, - two, - one / two, - two / one - ); - different = true; - } else { - println!(">>> Success!"); - } - } - - if different { - Err(anyhow!("grids are different")) - } else { - Ok(result) - } + Ok(result) } } @@ -560,9 +520,63 @@ pub struct Opts { impl Subcommand for Opts { #[cfg(feature = "fastnlo")] fn run(&self) -> Result<()> { - let grid = fastnlo::convert_fastnlo_table(&self.input, &self.pdfset, self.alpha)?; + use lhapdf::Pdf; + use pineappl_fastnlo::ffi::{self, Verbosity}; + use prettytable::{cell, row}; + + let mut file = ffi::make_fastnlo_lhapdf_with_name_file_set( + self.input.to_str().unwrap(), + &self.pdfset, + 0, + ); + let grid = fastnlo::convert_fastnlo_table(&file, self.alpha)?; - helpers::write_grid(&self.output, &grid) + 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).abs() + }; + + if rel_diff > 1e-10 { + 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"))] From 4c398e7582d9970a01125bc5ebc5dab673369af1 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 15:50:52 +0100 Subject: [PATCH 22/26] Make accuracy parameter writable from the CLI --- pineappl_cli/src/import.rs | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 4864d8ae..5ca88ba8 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -515,6 +515,9 @@ pub struct Opts { /// 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, } impl Subcommand for Opts { @@ -558,7 +561,7 @@ impl Subcommand for Opts { (two / one - 1.0).abs() }; - if rel_diff > 1e-10 { + if rel_diff > self.accuracy { different = false; } @@ -603,8 +606,10 @@ ARGS: LHAPDF id or name of the PDF set to check the converted grid with OPTIONS: - --alpha LO coupling power in alpha [default: 0] - -h, --help Print help information + --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 "; #[test] From 618c018ad09bbe4328e88ba8162d233197d9c4d8 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Mon, 14 Mar 2022 15:55:16 +0100 Subject: [PATCH 23/26] Enable silencing of fastNLO --- pineappl_cli/src/import.rs | 9 ++++++ pineappl_fastnlo/src/fastnlo.cpp | 50 ++++++++++++++++++++++++++++++-- pineappl_fastnlo/src/fastnlo.hpp | 3 +- pineappl_fastnlo/src/lib.rs | 22 ++++++++++++++ 4 files changed, 80 insertions(+), 4 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 5ca88ba8..4e7b13cc 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -518,6 +518,9 @@ pub struct Opts { /// 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 { @@ -527,10 +530,15 @@ impl Subcommand for Opts { 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)?; @@ -610,6 +618,7 @@ OPTIONS: 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] diff --git a/pineappl_fastnlo/src/fastnlo.cpp b/pineappl_fastnlo/src/fastnlo.cpp index dc34bd66..aae67437 100644 --- a/pineappl_fastnlo/src/fastnlo.cpp +++ b/pineappl_fastnlo/src/fastnlo.cpp @@ -1,5 +1,14 @@ #include "pineappl_fastnlo/src/fastnlo.hpp" +// TODO: is this portable enough? +#if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__)) +#define HAVE_UNISTD_H +#endif + +#ifdef HAVE_UNISTD_H +#include +#endif + #include #include #include @@ -41,16 +50,51 @@ 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) -{ +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()); - return std::unique_ptr(new fastNLOLHAPDF(arg0, arg1, PDFSet)); + 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)); } diff --git a/pineappl_fastnlo/src/fastnlo.hpp b/pineappl_fastnlo/src/fastnlo.hpp index d2f21242..a7004381 100644 --- a/pineappl_fastnlo/src/fastnlo.hpp +++ b/pineappl_fastnlo/src/fastnlo.hpp @@ -16,7 +16,8 @@ std::unique_ptr make_fastnlo_lhapdf_with_name_file_set( rust::Str name, rust::Str LHAPDFfile, - int PDFSet + int PDFSet, + bool silence ); rust::Vec CalcPDFLinearCombination( diff --git a/pineappl_fastnlo/src/lib.rs b/pineappl_fastnlo/src/lib.rs index 891b2a5e..8a191d4f 100644 --- a/pineappl_fastnlo/src/lib.rs +++ b/pineappl_fastnlo/src/lib.rs @@ -63,6 +63,27 @@ pub mod ffi { 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"); @@ -184,6 +205,7 @@ pub mod ffi { name: &str, lhapdf: &str, set: i32, + silence: bool, ) -> UniquePtr; unsafe fn dynamic_cast_coeff_add_fix(_: *mut fastNLOCoeffBase) -> *mut fastNLOCoeffAddFix; From 2d2c5045fe0e55a170f28e450f20ccdb20d3b8b4 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 15 Mar 2022 17:43:40 +0100 Subject: [PATCH 24/26] Test more macros --- pineappl_fastnlo/src/fastnlo.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pineappl_fastnlo/src/fastnlo.cpp b/pineappl_fastnlo/src/fastnlo.cpp index aae67437..82e7e2ad 100644 --- a/pineappl_fastnlo/src/fastnlo.cpp +++ b/pineappl_fastnlo/src/fastnlo.cpp @@ -1,7 +1,8 @@ #include "pineappl_fastnlo/src/fastnlo.hpp" // TODO: is this portable enough? -#if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__)) +#if defined (__unix__) || defined(__unix) || defined(unix) || \ + (defined (__APPLE__) && defined (__MACH__)) #define HAVE_UNISTD_H #endif From 0452104af1f324fde6a5e2e23e33f9c31f6fd2af Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Tue, 15 Mar 2022 17:44:02 +0100 Subject: [PATCH 25/26] Shuffle around helper routines --- pineappl_fastnlo/src/lib.rs | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/pineappl_fastnlo/src/lib.rs b/pineappl_fastnlo/src/lib.rs index 8a191d4f..10b8c9f2 100644 --- a/pineappl_fastnlo/src/lib.rs +++ b/pineappl_fastnlo/src/lib.rs @@ -115,14 +115,6 @@ pub mod ffi { include!("fastnlotk/fastNLOPDFLinearCombinations.h"); type fastNLOPDFLinearCombinations; - - fn CalcPDFLinearCombination( - _: &fastNLOPDFLinearCombinations, - _: &fastNLOCoeffAddBase, - _: &[f64], - _: &[f64], - _: bool, - ) -> Vec; } unsafe extern "C++" { @@ -172,8 +164,18 @@ pub mod ffi { 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; @@ -181,7 +183,6 @@ pub mod ffi { fn GetScaleNodes2(_: &fastNLOCoeffAddFlex, _: i32) -> Vec; fn GetXNodes1(_: &fastNLOCoeffAddBase, _: i32) -> Vec; fn GetXNodes2(_: &fastNLOCoeffAddBase, _: i32) -> Vec; - fn GetObsBinDimBounds(_: &fastNLOTable, _: u32, _: u32) -> pair_double_double; fn GetSigmaTilde( _: &fastNLOCoeffAddFlex, _: usize, @@ -191,7 +192,6 @@ pub mod ffi { _: usize, _: i32, ) -> f64; - fn GetNx(_: &fastNLOCoeffAddFlex, _: usize) -> usize; fn downcast_coeff_add_fix_to_base(_: &fastNLOCoeffAddFix) -> &fastNLOCoeffAddBase; fn downcast_coeff_add_flex_to_base(_: &fastNLOCoeffAddFlex) -> &fastNLOCoeffAddBase; @@ -201,15 +201,16 @@ pub mod ffi { 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; - - unsafe fn dynamic_cast_coeff_add_fix(_: *mut fastNLOCoeffBase) -> *mut fastNLOCoeffAddFix; - unsafe fn dynamic_cast_coeff_add_flex(_: *mut fastNLOCoeffBase) - -> *mut fastNLOCoeffAddFlex; } } From f236fdc7e9ab899f1aac1228b4426e875762dd53 Mon Sep 17 00:00:00 2001 From: Christopher Schwan Date: Wed, 16 Mar 2022 13:13:06 +0100 Subject: [PATCH 26/26] Show relative differences with sign --- pineappl_cli/src/import.rs | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pineappl_cli/src/import.rs b/pineappl_cli/src/import.rs index 4e7b13cc..f8e75823 100644 --- a/pineappl_cli/src/import.rs +++ b/pineappl_cli/src/import.rs @@ -563,13 +563,9 @@ impl Subcommand for Opts { .enumerate() { // catches the case where both results are zero - let rel_diff = if one == two { - 0.0 - } else { - (two / one - 1.0).abs() - }; + let rel_diff = if one == two { 0.0 } else { two / one - 1.0 }; - if rel_diff > self.accuracy { + if rel_diff.abs() > self.accuracy { different = false; }