Improve CAPI examples

.github/workflows/capi.yaml

@@ -16,14 +16,12 @@ jobs:
           cargo cinstall --verbose --prefix=/usr/local/ --manifest-path pineappl_capi/Cargo.toml
           ldconfig
 
-      - name: Test C example
+      - name: Test C++ example
         run: |
-          cd examples/basic-capi-usage
+          cd examples/cpp
           # if `make` is too old, it doesn't support the `!=` operator
           sed -i "s/\([a-zA-Z_]\+\) != \(.*\)$/echo \1 = \$(\2)/e" Makefile
-          make
-          ./dyaa
-          test -f ./DY-LO-AA.pineappl.lz4
+          make test-examples
 
       - name: Test Fortran example
         run: |
@@ -33,7 +31,7 @@ jobs:
           ./dyaa
           test -f ./DY-LO-AA.pineappl.lz4
 
-      - name: Test C++ example
+      - name: Test OO C++ example
         run: |
           cd examples/object-oriented-cpp
           sed -i "s/\([a-zA-Z_]\+\) != \(.*\)$/echo \1 = \$(\2)/e" Makefile

examples/cpp/Makefile

@@ -0,0 +1,21 @@
+CXX = c++
+CXXFLAGS = -std=c++11 -O3 -Wall -Wextra
+PINEAPPL_DEPS != pkg-config --cflags --libs pineappl_capi
+LHAPDF_DEPS != pkg-config --cflags --libs lhapdf
+
+all: convolute-grid fill-grid
+
+test-examples: convolute-grid fill-grid
+	./fill-grid
+	./convolute-grid
+
+convolute-grid: convolute-grid.cpp
+	$(CXX) $(CXXFLAGS) $< $(LHAPDF_DEPS) $(PINEAPPL_DEPS) -o $@
+
+fill-grid: fill-grid.cpp
+	$(CXX) $(CXXFLAGS) $< $(PINEAPPL_DEPS) -o $@
+
+PHONY: clean
+
+clean:
+	rm -f convolute-grid fill-grid

examples/cpp/convolute-grid.cpp

@@ -0,0 +1,93 @@
+#include <LHAPDF/PDF.h>
+#include <pineappl_capi.h>
+
+#include <cstddef>
+#include <iomanip>
+#include <ios>
+#include <iostream>
+#include <string>
+#include <vector>
+
+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::printf("Usage: %s [grid] [pdf]", argv[0]);
+    }
+
+    // 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 <LHAPDF::PDF*> (pdf)->xfxQ2(id, x, q2);
+    };
+    auto alphas = [](double q2, void* pdf) {
+        return static_cast <LHAPDF::PDF*> (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<double> 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<double> 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_convolute_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.
+    pineappl_grid_convolute_with_one(grid, 2212, xfx, alphas, pdf, order_mask,
+        channel_mask, xir, xif, dxsec.data());
+
+    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
+        std::cout << std::scientific << dxsec.at(bin) << std::defaultfloat << '\n';
+    }
+}

examples/basic-capi-usage/dyaa.cpp → examples/cpp/fill-grid.cpp

@@ -1,11 +1,9 @@
 #include <pineappl_capi.h>
-#include <LHAPDF/LHAPDF.h>
 
 #include <cmath>
 #include <cstddef>
 #include <cstdio>
 #include <random>
-#include <vector>
 
 double int_photo(double s, double t, double u) {
     double alpha0 = 1.0 / 137.03599911;
@@ -93,70 +91,97 @@ void fill_grid(pineappl_grid* grid, std::size_t calls) {
         }
 
         auto weight = jacobian * int_photo(s, t, u);
-        double  q2 = 90.0 * 90.0;
-
-        pineappl_grid_fill(grid, x1, x2, q2, 0, fabs(yll), 0, weight);
+        double q2 = 90.0 * 90.0;
+        std::size_t order = 0;
+        std::size_t channel = 0;
+
+        // fill the LO `weight` into `grid` for parton fractions `x1` and `x2`, and the (squared)
+        // renormalization/factorization scale `q2`. The parameters `order` and `channel` are
+        // indices defined from the arrays `orders` and `channel` used in creating the grid. In this
+        // case they are both `0` and denote the order #0 (leading order) and the channel #0
+        // (photon-photon channel), respectively
+        pineappl_grid_fill(grid, x1, x2, q2, order, fabs(yll), channel, weight);
     }
 }
 
 int main() {
-    // create a new luminosity function for the $\gamma\gamma$ initial state
-    auto* lumi = pineappl_lumi_new();
-    int32_t pdg_ids[] = { 22, 22 };
-    double ckm_factors[] = { 1.0 };
-    pineappl_lumi_add(lumi, 1, pdg_ids, ckm_factors);
+    // ---
+    // Create all channels
+
+    // this object will contain all channels (initial states) that we define
+    auto* channels = pineappl_lumi_new();
+
+    // photon-photon initial state, where `22` is the photon (PDG MC ids)
+    int32_t pids1[] = { 22, 22 };
+
+    // factor that each channel is multiplied with when convoluting with PDFs
+    double factors1[] = { 1.0 };
+
+    // define the channel #0
+    pineappl_lumi_add(channels, 1, pids1, factors1);
+
+    // create another channel, which we won't fill, however
+
+    // this channel is the down-type-antidown-type quark channel; here we combine down-antidown,
+    // strange-antistrange and bottom-antibottom into a single channel, which is often done if the
+    // CKM matrix is taken to be diagonal
+    int32_t pids2[] = { 1, -1, 3, -3, 5, -5 };
+
+    // for each pair of particle ids we need to give a factor; in case of a non-diagonal CKM matrix
+    // we could factor out the CKM matrix elements here
+    double factors2[] = { 1.0, 1.0, 1.0 };
+
+    // define the channel #1
+    pineappl_lumi_add(channels, 3, pids2, factors2);
+
+    // ---
+    // Specify the perturbative orders that will be filled into the grid
+
+    // in this example we only fill the LO, which has the exponents
+    // - 0 of alphas,
+    // - 2 of alpha (electroweak coupling),
+    // - 0 of log (xiR^2) (renormalization scale logarithm) and
+    // - 0 of log (xiF^2) (factorization scale logarithm)
+    uint32_t orders[] = {
+        0, 2, 0, 0, // order #0: LO
+        1, 2, 0, 0, // order #1: NLO QCD
+        1, 2, 0, 1  // order #2: NLO QCD factorization log
+    };
 
-    // only LO $\alpha_\mathrm{s}^0 \alpha^2 \log^0(\xi_\mathrm{R}) \log^0(\xi_\mathrm{F})$
-    uint32_t orders[] = { 0, 2, 0, 0 };
+    // ---
+    // Specify the bin limits
 
-    // we bin in rapidity from 0 to 2.4 in steps of 0.1
+    // Similar to many Monte Carlo integrators PineAPPL supports only one-dimensional differential
+    // distributions, and only one distribution for each grid. However, one can generate multiple
+    // grids to support multiple distributions, and since every n-dimensional distribution can be
+    // written as a one-dimensional one (by using the bin index as a new binning variable, for
+    // instance), this isn't a limitation.
+
+    // we bin the rapidity of the final-state lepton pair from 0 to 2.4 in steps of 0.1
     double bins[] = {
         0.0,
         0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2,
         1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
     };
 
+    // ---
+    // Create the grid using the previously set information about orders, bins and channels
+
     // create the PineAPPL grid with default interpolation and binning parameters
     auto* keyval = pineappl_keyval_new();
-    auto* grid = pineappl_grid_new(lumi, 1, orders, 24, bins, keyval);
+    auto* grid = pineappl_grid_new(channels, 1, orders, 24, bins, keyval);
 
     // now we no longer need `keyval` and `lumi`
     pineappl_keyval_delete(keyval);
-    pineappl_lumi_delete(lumi);
+    pineappl_lumi_delete(channels);
 
-    // fill the grid with phase-space points
+    // ---
+    // Fill the grid with phase-space points
     fill_grid(grid, 10000000);
 
-    // perform a convolution of the grid with PDFs
-    auto* pdf = LHAPDF::mkPDF("NNPDF31_nlo_as_0118_luxqed", 0);
-    auto xfx = [](int32_t id, double x, double q2, void* pdf) {
-        return static_cast <LHAPDF::PDF*> (pdf)->xfxQ2(id, x, q2);
-    };
-    auto alphas = [](double q2, void* pdf) {
-        return static_cast <LHAPDF::PDF*> (pdf)->alphasQ2(q2);
-    };
-
-    std::vector<double> dxsec(24);
-    pineappl_grid_convolute_with_one(grid, 2212, xfx, alphas, pdf, nullptr,
-        nullptr, 1.0, 1.0, dxsec.data());
-
-    // print the results
-    for (std::size_t i = 0; i != 24; ++i) {
-        std::printf("%02u %.1f %.1f %.3e\n", i, bins[i], bins[i + 1], dxsec[i]);
-    }
-
-    // store some metadata in the grid
-    pineappl_grid_set_key_value(grid, "events", "10000000");
-
-    // read out the stored value and print it on stdout
-    auto* value = pineappl_grid_key_value(grid, "events");
-    std::printf("Finished running %s events.\n", value);
-
-    // delete the allocated object
-    pineappl_string_delete(value);
-
-    // write the grid to disk - with `.lz4` suffix the grid is automatically LZ4 compressed
-    char const* filename = "DY-LO-AA.pineappl.lz4";
+    // ---
+    // Write the grid to disk - with `.lz4` suffix the grid is automatically LZ4 compressed
+    char const* filename = "drell-yan-rap-ll.pineappl.lz4";
     pineappl_grid_write(grid, filename);
 
     // destroy the object