Fix bug in order masking

Author: Radonirinaunimi

examples/cpp/evolve-grid.cpp

@@ -50,6 +50,21 @@ int main() {
     // TODO: How to get a Grid that can be evolved??
     std::string filename = "LHCB_WP_7TEV_opt.pineappl.lz4";
 
+    // disable LHAPDF banners to guarantee deterministic output
+    LHAPDF::setVerbosity(0);
+    std::string pdfset = "NNPDF31_nlo_as_0118_luxqed";
+    auto pdf = std::unique_ptr<LHAPDF::PDF>(LHAPDF::mkPDF(pdfset, 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<LHAPDF::PDF*> pdfs = {pdf.get(), pdf.get()};
+    void** pdf_states = reinterpret_cast<void**>(pdfs.data());
+
     // read the grid from a file
     auto* grid = pineappl_grid_read(filename.c_str());
 
@@ -71,8 +86,8 @@ int main() {
 
     // Get the shape of the evolve info objects
     std::vector<std::size_t> evinfo_shape(5);
-    std::vector<uint8_t> apply_order = {3, 0};
-    pineappl_grid_evolve_info_shape(grid, apply_order.data(), evinfo_shape.data());
+    std::vector<uint8_t> max_orders = {3, 2};
+    pineappl_grid_evolve_info_shape(grid, max_orders.data(), evinfo_shape.data());
 
     // Get the values of the evolve info parameters. These contain, for example, the
     // information on the `x`-grid and `PID` used to interpolate the Grid.
@@ -82,7 +97,7 @@ int main() {
     std::vector<int> pids1(evinfo_shape[2]);
     std::vector<double> x1(evinfo_shape[3]);
     std::vector<double> ren1(evinfo_shape[4]);
-    pineappl_grid_evolve_info(grid, apply_order.data(), fac1.data(),
+    pineappl_grid_evolve_info(grid, max_orders.data(), fac1.data(),
         frg1.data(), pids1.data(), x1.data(), ren1.data());
 
     // ------------------ Construct the Operator Info ------------------
@@ -114,32 +129,22 @@ int main() {
         }
     }
 
-    // Construct the values of alphas
-    std::vector<double> alphas_table(ren1.size(), 0.118);
+    // Construct the values of alphas table
+    std::vector<double> alphas_table;
+    for (double q2 : ren1) {
+        double alpha = alphas(q2, pdf.get());
+        alphas_table.push_back(alpha);
+    }
+
     std::vector<double> xi = {1.0, 1.0, 1.0};
     std::vector<std::size_t> tensor_shape = {pids1.size(), x1.size(), pids1.size(), x1.size()};
 
     pineappl_fk_table* fktable = pineappl_grid_evolve(grid, opinfo_slices.data(),
-        apply_order.data(), op_slices.data(), x1.data(),
+        max_orders.data(), op_slices.data(), x1.data(),
         x1.data(), pids1.data(), pids1.data(),
         tensor_shape.data(), xi.data(), ren1.data(), alphas_table.data());
 
     // ------------------ Compare Grid & FK after convolution ------------------
-    // disable LHAPDF banners to guarantee deterministic output
-    LHAPDF::setVerbosity(0);
-    std::string pdfset = "NNPDF31_nlo_as_0118_luxqed";
-    auto pdf = std::unique_ptr<LHAPDF::PDF>(LHAPDF::mkPDF(pdfset, 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<LHAPDF::PDF*> pdfs = {pdf.get(), pdf.get()};
-    void** pdf_states = reinterpret_cast<void**>(pdfs.data());
-
     // how many bins does this grid have?
     std::size_t bins = pineappl_grid_bin_count(grid);
 
@@ -175,7 +180,7 @@ int main() {
     // Print the data
     std::cout << std::scientific << std::setprecision(6);
     for (size_t i = 0; i < dxsec_grid.size(); ++i) {
-        double reldiff = std::abs(dxsec_grid[i] - dxsec_fktable[i]) / (std::abs(dxsec_grid[i]));
+        double reldiff = (dxsec_fktable[i] - dxsec_grid[i]) / dxsec_grid[i];
         std::cout << std::setw(idx_width) << i
                   << std::setw(num_width) << dxsec_grid[i]
                   << std::setw(num_width) << dxsec_fktable[i]

examples/cpp/output

@@ -174,14 +174,14 @@ idx   left  right    dsig/dx     dx
         23     40626      (16, 12, 26)     (5442.31, 0.265114, 0.00287387)    -7.04493e-11
    Bin           Grid        FkTable        reldiff
   ----  -------------  -------------  -------------
-     0   7.545911e+02   7.603251e+02   7.598823e-03
-     1   6.902834e+02   6.949192e+02   6.715820e-03
-     2   6.002520e+02   6.037697e+02   5.860343e-03
-     3   4.855224e+02   4.878688e+02   4.832807e-03
-     4   3.619546e+02   3.634639e+02   4.169948e-03
-     5   2.458669e+02   2.467891e+02   3.750676e-03
-     6   1.158685e+02   1.163001e+02   3.724532e-03
-     7   2.751727e+01   2.763520e+01   4.285840e-03
+     0   7.545911e+02   7.545911e+02  -3.762829e-08
+     1   6.902834e+02   6.902834e+02  -3.537025e-08
+     2   6.002520e+02   6.002520e+02  -3.293028e-08
+     3   4.855224e+02   4.855223e+02  -3.037598e-08
+     4   3.619546e+02   3.619545e+02  -2.782860e-08
+     5   2.458669e+02   2.458669e+02  -2.526536e-08
+     6   1.158685e+02   1.158685e+02  -2.182740e-08
+     7   2.751727e+01   2.751727e+01  -1.732084e-08
 idx  p-p c#0 l#0 p-p~ c#0 l#  p-d c#0 l#0       p-d       dx
 --- ------------ ----------- ------------- ------------ ------
   0 5.263109e-01 5.263109e-01 5.263109e-01 5.263109e-01    0.1