diff --git a/examples/cpp/Makefile b/examples/cpp/Makefile index 2e6db08d..39498abd 100644 --- a/examples/cpp/Makefile +++ b/examples/cpp/Makefile @@ -10,6 +10,7 @@ PROGRAMS = \ advanced-convolution \ advanced-filling \ convolve-grid \ + convolve-grid-v1 \ deprecated \ display-channels \ display-orders \ @@ -30,6 +31,9 @@ advanced-filling: advanced-filling.cpp convolve-grid: convolve-grid.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ +convolve-grid-v1: convolve-grid-v1.cpp + $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ + deprecated: deprecated.cpp $(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@ diff --git a/examples/cpp/convolve-grid-v1.cpp b/examples/cpp/convolve-grid-v1.cpp new file mode 100644 index 00000000..9e11ae43 --- /dev/null +++ b/examples/cpp/convolve-grid-v1.cpp @@ -0,0 +1,124 @@ +//////////////////////////////////////////////////////////////////////////// +// Exactly the same as `convolve-grid.cpp` but using the generalized +// convolution: `pineappl_grid_convolve`. +// +// TODO: Make it such that it does not exactly copy `convolve-grid.cpp`, +// perhaps showing as an example something with 3 Convolutions. +//////////////////////////////////////////////////////////////////////////// +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +int main(int argc, char* argv[]) { + std::string filename = "drell-yan-rap-ll.pineappl.lz4"; + std::string pdfset = "NNPDF31_nlo_as_0118_luxqed"; + + switch (argc) { + case 3: + pdfset = argv[2]; + // fall through + case 2: + filename = argv[1]; + case 1: + break; + + default: + std::cout << "Usage: " << argv[0] << " [grid] [pdf]\n"; + } + + // disable LHAPDF banners to guarantee deterministic output + LHAPDF::setVerbosity(0); + + // read the grid from a file + auto* grid = pineappl_grid_read(filename.c_str()); + + auto* pdf = LHAPDF::mkPDF(pdfset, 0); + + // define callables for the PDFs and alphas + auto xfx = [](int32_t id, double x, double q2, void* pdf) { + return static_cast (pdf)->xfxQ2(id, x, q2); + }; + auto alphas = [](double q2, void* pdf) { + return static_cast (pdf)->alphasQ2(q2); + }; + + // how many bins does this grid have? + std::size_t bins = pineappl_grid_bin_count(grid); + + // how many dimensions does each bin have? + std::size_t dims = pineappl_grid_bin_dimensions(grid); + + // allocate a vector holding the left and right bin limits for each dimension + std::vector bin_limits(2 * bins * dims); + + for (std::size_t dim = 0; dim != dims; ++dim) { + pineappl_grid_bin_limits_left(grid, dim, &bin_limits.at((2 * dim + 0) * bins)); + pineappl_grid_bin_limits_right(grid, dim, &bin_limits.at((2 * dim + 1) * bins)); + } + + // allocate a vector holding the differential cross sections + std::vector dxsec(bins); + + auto order_mask = nullptr; + auto channel_mask = nullptr; + double xir = 1.0; + double xif = 1.0; + + // perform the convolution of `grid` with the PDFs given as `xfx` and the strong coupling in + // `alphas` and the extra parameter `pdf`, which is passed to `xfx` and `alphas` as the last + // parameter. The integer `2212` is the PDG MC id for a proton and signals and `xfx` is the PDF + // of a proton. In this case we assume that both initial state hadrons' PDFs can derived from + // that of a proton. If this isn't the case, for instance for a proton-lead collision, both PDFs + // must be given separately and the function `pineappl_grid_convolve_with_two` must be used. + // The parameters `order_mask` and `channel_mask` can be used to select specific orders and + // channels, respectively. Using `xir` and `xif` the renormalization and factorization scales + // can be varied around its central values, respectively. + std::vector mu_scales = { xir, xif, 1.0 }; + using LambdaType = double(*)(int32_t, double, double, void *); + LambdaType xfxs[] = { xfx, xfx}; + pineappl_grid_convolve(grid, xfxs, alphas, pdf, order_mask, channel_mask, nullptr, 1, + mu_scales.data(), dxsec.data()); + + std::vector normalizations(bins); + + // read out the bin normalizations, which is usually the size of each bin + pineappl_grid_bin_normalizations(grid, normalizations.data()); + + // print table header + std::cout << "idx"; + for (std::size_t dim = 0; dim != dims; ++dim) { + std::cout << " left right"; + } + std::cout << " dsig/dx dx\n"; + std::cout << "---"; + for (std::size_t dim = 0; dim != dims; ++dim) { + std::cout << " ----dim #" << dim << "---"; + } + std::cout << " ------------ ------\n"; + + for (std::size_t bin = 0; bin != bins; ++bin) { + // print the bin index + std::cout << std::setw(3) << bin << ' '; + + for (std::size_t dim = 0; dim != dims; ++dim) { + double left_limit = bin_limits.at((2 * dim + 0) * bins + bin); + double right_limit = bin_limits.at((2 * dim + 1) * bins + bin); + + // print the left and right bin limit for each dimension + std::cout << std::setw(6) << left_limit << ' ' << std::setw(6) << right_limit << ' '; + } + + // print the result together with the normalization + std::cout << std::scientific << dxsec.at(bin) << std::defaultfloat << ' ' + << std::setw(6) << normalizations.at(bin) << '\n'; + } + + pineappl_grid_delete(grid); +} diff --git a/examples/cpp/fill-grid-v1.cpp b/examples/cpp/fill-grid-v1.cpp index fd0f688e..d0bc256c 100644 --- a/examples/cpp/fill-grid-v1.cpp +++ b/examples/cpp/fill-grid-v1.cpp @@ -15,7 +15,6 @@ #include #include #include -#include #include #include #include diff --git a/examples/cpp/output b/examples/cpp/output index dd0c7936..2e742d59 100644 --- a/examples/cpp/output +++ b/examples/cpp/output @@ -37,6 +37,32 @@ idx p-p c#0 l#0 p-p~ c#0 l# p-d c#0 l#0 p-d dx 22 1.967518e-02 1.967518e-02 1.967518e-02 1.967518e-02 0.1 23 5.565306e-03 5.565306e-03 5.565306e-03 5.565306e-03 0.1 idx left right dsig/dx dx +--- ----dim #0--- ------------ ------ + 0 0 0.1 5.263109e-01 0.1 + 1 0.1 0.2 5.254908e-01 0.1 + 2 0.2 0.3 5.246824e-01 0.1 + 3 0.3 0.4 5.188340e-01 0.1 + 4 0.4 0.5 5.175482e-01 0.1 + 5 0.5 0.6 5.008841e-01 0.1 + 6 0.6 0.7 4.905325e-01 0.1 + 7 0.7 0.8 4.675734e-01 0.1 + 8 0.8 0.9 4.393159e-01 0.1 + 9 0.9 1 3.992921e-01 0.1 + 10 1 1.1 3.706801e-01 0.1 + 11 1.1 1.2 3.264717e-01 0.1 + 12 1.2 1.3 2.849345e-01 0.1 + 13 1.3 1.4 2.486723e-01 0.1 + 14 1.4 1.5 2.110419e-01 0.1 + 15 1.5 1.6 1.797439e-01 0.1 + 16 1.6 1.7 1.471492e-01 0.1 + 17 1.7 1.8 1.205566e-01 0.1 + 18 1.8 1.9 9.491625e-02 0.1 + 19 1.9 2 7.255720e-02 0.1 + 20 2 2.1 5.056967e-02 0.1 + 21 2.1 2.2 3.491788e-02 0.1 + 22 2.2 2.3 1.967518e-02 0.1 + 23 2.3 2.4 5.565306e-03 0.1 +idx left right dsig/dx dx --- ----dim #0--- ------------ ------ 0 0 0.1 5.263109e-01 0.1 1 0.1 0.2 5.254908e-01 0.1