C++ OOP API Examples

This example briefly demonstrates how to use the NeoPDF C++ Object Oriented (OOP) API to load and evaluate parton distributions. More examples can be found in neopdf_capi/tests.

Prerequisites

Build and install the C++ API as described in the installation guide. The C++ OOP header is needed for the following examples.

Example 1: Loading and Evaluating PDFs

This example demonstrates the use of the NeoPDF C++ OOP API to load both single and multiple PDF members, evaluate parton distributions for a range of \(x\) and \(Q^2\) values, and compare results to LHAPDF.

Technical details:

  • The NeoPDF and NeoPDFs objects manage their own memory and automatically release resources when they go out of scope (RAII).

  • The evaluation of \(x f(x, Q^2)\) and \(\alpha_s(Q^2)\) is vectorized over the input axes for efficiency.

  • The code asserts that the results from NeoPDF and LHAPDF agree within a tight tolerance, providing a robust validation.

  1#include <LHAPDF/PDF.h>
  2#include <LHAPDF/GridPDF.h>
  3#include <NeoPDF.hpp>
  4#include <cassert>
  5#include <cmath>
  6#include <iomanip>
  7#include <iostream>
  8#include <string>
  9#include <vector>
 10
 11using namespace neopdf;
 12
 13const double TOLERANCE= 1e-16;
 14
 15template<typename T>
 16std::vector<T> geomspace(T start, T stop, int num, bool endpoint = false) {
 17    std::vector<T> result(num);
 18
 19    if (num == 1) {
 20        result[0] = start;
 21        return result;
 22    }
 23
 24    T log_start = std::log(start);
 25    T log_stop = std::log(stop);
 26    T step = (log_stop - log_start) / (endpoint ? (num - 1) : num);
 27
 28    for (int i = 0; i < num; ++i) {
 29        result[i] = std::exp(log_start + i * step);
 30    }
 31
 32    return result;
 33}
 34
 35template<typename T>
 36std::vector<T> linspace(T start, T stop, int num, bool endpoint = true) {
 37    std::vector<T> result(num);
 38
 39    if (num == 1) {
 40        result[0] = start;
 41        return result;
 42    }
 43
 44    T step = (stop - start) / (endpoint ? (num - 1) : num);
 45
 46    for (int i = 0; i < num; ++i) {
 47        result[i] = start + i * step;
 48    }
 49
 50    return result;
 51}
 52
 53void test_xfxq2() {
 54    std::cout << "=== Test xfxQ2 for single PDF member ===\n";
 55
 56    // disable LHAPDF banners to guarantee deterministic output
 57    LHAPDF::setVerbosity(0);
 58
 59    std::string pdfname = "NNPDF40_nnlo_as_01180";
 60    NeoPDF neo_pdf(pdfname.c_str(), 0);
 61    const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(pdfname);
 62    const LHAPDF::GridPDF& lha_pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);
 63
 64    std::vector<int> pids = {-5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5};
 65    std::vector<double> xs = geomspace(neo_pdf.x_min(), neo_pdf.x_max(), 200);
 66    std::vector<double> q2s = geomspace(neo_pdf.q2_min(), neo_pdf.q2_max(), 200);
 67
 68    // Headers of the table to print the results
 69    std::cout << std::right
 70        << std::setw(6) << "pid"
 71        << std::setw(15) << "x"
 72        << std::setw(15) << "Q2"
 73        << std::setw(15) << "LHAPDF"
 74        << std::setw(15) << "NeoPDF"
 75        << std::setw(15) << "Rel. Diff." << "\n";
 76    std::cout << std::string(81, '-') << "\n";
 77
 78    for (const auto &pid: pids) {
 79        for (const auto &x: xs) {
 80            for (const auto &q2: q2s) {
 81                double expected = lha_pdf.xfxQ2(pid, x, q2);
 82                double result = neo_pdf.xfxQ2(pid, x, q2);
 83                double reldif = std::abs(result - expected) / expected;
 84
 85                assert(std::abs(result - expected) < TOLERANCE);
 86
 87                // Print the results as a table
 88                std::cout << std::scientific << std::setprecision(8)
 89                    << std::right
 90                    << std::setw(6)  << pid
 91                    << std::setw(15) << x
 92                    << std::setw(15) << q2
 93                    << std::right
 94                    << std::setw(15) << expected
 95                    << std::setw(15) << result
 96                    << std::setw(15) << reldif << "\n";
 97                }
 98            }
 99
100    }
101}
102
103void test_alphas_q2() {
104    std::cout << "=== Test alphasQ2 for single PDF member ===\n";
105
106    // disable LHAPDF banners to guarantee deterministic output
107    LHAPDF::setVerbosity(0);
108
109    std::string pdfname = "NNPDF40_nnlo_as_01180";
110    NeoPDF neo_pdf(pdfname.c_str(), 0);
111    const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(pdfname);
112    const LHAPDF::GridPDF& lha_pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);
113
114    std::vector<double> q2_points = linspace(4.0, 1e10, 500);
115
116    // Headers of the table to print the results
117    std::cout << std::right
118        << std::setw(15) << "Q2"
119        << std::setw(15) << "LHAPDF"
120        << std::setw(15) << "NeoPDF"
121        << std::setw(15) << "Rel. Diff." << "\n";
122    std::cout << std::string(60, '-') << "\n";
123
124    for (const auto& q2: q2_points) {
125        double expected = lha_pdf.alphasQ2(q2);
126        double result = neo_pdf.alphasQ2(q2);
127        double reldif = std::abs(result - expected) / expected;
128
129        assert(std::abs(result - expected) < TOLERANCE);
130
131        // Print the results as a table
132        std::cout << std::scientific << std::setprecision(8)
133            << std::right
134            << std::setw(15) << q2
135            << std::right
136            << std::setw(15) << expected
137            << std::setw(15) << result
138            << std::setw(15) << reldif << "\n";
139    }
140}
141
142void test_all_pdf_members() {
143    std::cout << "=== Test PDFs class (loading all members) ===\n";
144
145    // disable LHAPAPDF banners to guarantee deterministic output
146    LHAPDF::setVerbosity(0);
147
148    std::string pdfname = "NNPDF40_nnlo_as_01180";
149    NeoPDFs neo_pdfs(pdfname.c_str());
150
151    std::cout << "Loaded " << neo_pdfs.size() << " PDF members\n";
152
153    // Test case: evaluate a simple point across all members
154    int pid = 1;
155    double x = 1e-9;
156    double q2 = 1.65 * 1.65;
157
158    std::cout << std::right
159        << std::setw(8) << "Member"
160        << std::setw(15) << "LHAPDF"
161        << std::setw(15) << "NeoPDF"
162        << std::setw(15) << "Rel. Diff." << "\n";
163    std::cout << std::string(53, '-') << "\n";
164
165    // Evaluate the same point across all PDF members
166    std::vector<double> results;
167    for (size_t i = 0; i < neo_pdfs.size(); ++i) {
168        const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(pdfname, i);
169        const LHAPDF::GridPDF& lha_pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);
170
171        double expected = lha_pdf.xfxQ2(pid, x, q2);
172        double result = neo_pdfs[i].xfxQ2(pid, x, q2);
173
174        double reldif = std::abs(result - expected) / expected;
175        assert(std::abs(result - expected) < TOLERANCE);
176        results.push_back(result);
177
178        std::cout << std::right
179            << std::setw(8) << i
180            << std::scientific << std::setprecision(8)
181            << std::setw(15) << expected
182            << std::setw(15) << result
183            << std::setw(15) << reldif << "\n";
184    }
185
186    // Calculate some statistics
187    double sum = 0.0;
188    for (double result : results) {
189        sum += result;
190    }
191    double mean = sum / results.size();
192
193    double variance = 0.0;
194    for (double result : results) {
195        variance += (result - mean) * (result - mean);
196    }
197    variance /= results.size();
198    double std_dev = std::sqrt(variance);
199
200    std::cout << "\nStatistics across all members:\n";
201    std::cout << "Mean: " << std::scientific << std::setprecision(8) << mean << "\n";
202    std::cout << "Std Dev: " << std_dev << "\n";
203    std::cout << "Relative Std Dev: " << std_dev / mean << "\n";
204}
205
206void test_lazy_loading() {
207    std::cout << "=== Test NeoPDFLazy class (lazy loading) ===\n";
208
209    // Disable LHAPDF banners to guarantee deterministic output
210    LHAPDF::setVerbosity(0);
211
212    std::string pdfname = "NNPDF40_nnlo_as_01180.neopdf.lz4";
213    NeoPDFLazy lazy_pdfs(pdfname);
214
215    std::cout << "Initialized lazy loader for " << pdfname << "\n";
216
217    // Test case: evaluate a simple point across all members
218    int pid = 1;
219    double x = 1e-9;
220    double q2 = 1.65 * 1.65;
221
222    std::cout << std::right
223        << std::setw(8) << "Member"
224        << std::setw(15) << "LHAPDF"
225        << std::setw(15) << "NeoPDF"
226        << std::setw(15) << "Rel. Diff." << "\n";
227    std::cout << std::string(53, '-') << "\n";
228
229    // Evaluate the same point across all PDF members
230    std::vector<double> results;
231    int member_idx = 0;
232    while (auto neo_pdf = lazy_pdfs.next()) {
233        const LHAPAPDF::PDF* basepdf = LHAPDF::mkPDF("NNPDF40_nnlo_as_01180", member_idx);
234        const LHAPDF::GridPDF& lha_pdf = *dynamic_cast<const LHAPDF::GridPDF*>(basepdf);
235
236        double expected = lha_pdf.xfxQ2(pid, x, q2);
237        double result = neo_pdf->xfxQ2(pid, x, q2);
238
239        double reldif = std::abs(result - expected) / expected;
240        assert(std::abs(result - expected) < TOLERANCE);
241        results.push_back(result);
242
243        std::cout << std::right
244            << std::setw(8) << member_idx
245            << std::scientific << std::setprecision(8)
246            << std::setw(15) << expected
247            << std::setw(15) << result
248            << std::setw(15) << reldif << "\n";
249        member_idx++;
250    }
251
252    std::cout << "\nSuccessfully iterated through all members lazily.\n";
253}
254
255int main() {
256    // Test the computation of the PDF interpolations
257    test_xfxq2();
258
259    // Test the computation of the `alphas` interpolations
260    test_alphas_q2();
261
262    // Test the PDF interpolations by loading all the members
263    test_all_pdf_members();
264
265    // Test the lazy loading of PDF members
266    test_lazy_loading();
267
268    return EXIT_SUCCESS;
269}

Example 2: Filling and Writing a NeoPDF Grid

This example illustrates how to fill and write a NeoPDF grid using the C++ OOP API. It demonstrates the process of constructing a grid for each PDF member and serializing the collection to disk.

The filling of the PDF grid in the following example assumes no dependence in the nucleon numbers \(A\) and strong coupling \(\alpha_s\) (standard LHAPDF-like PDF). Refer to the Section below in the case the grid should explicitly depend on more parameters.

Technical details:

  • The grid axes are defined as vectors for \(x\), \(Q^2\), parton IDs, nucleons, and \(\alpha_s\) values.

  • The grid data is stored in a 6D array, with the layout [nucleons][alphas][pids][kT][xs][q2s].

  • The GridWriter class manages the collection of grids and handles compression and serialization to disk.

  • Metadata is filled in a MetaData object, which includes information about the set, axis ranges, flavors, and interpolation type. This metadata is essential for correct interpretation of the grid file.

  • All memory management is automatic; no manual deallocation is required.

  • The output file is compressed and written in the .neopdf.lz4 format, suitable for use with NeoPDF (CLI) tools and APIs.

  1#include <NeoPDF.hpp>
  2#include <cassert>
  3#include <cmath>
  4#include <cstdio>
  5#include <cstdlib>
  6#include <cstring>
  7#include <iostream>
  8#include <vector>
  9
 10using namespace neopdf;
 11
 12const double TOLERANCE= 1e-16;
 13
 14int main() {
 15    const char* pdfname = "NNPDF40_nnlo_as_01180";
 16    // Load all PDF members
 17    NeoPDFs neo_pdfs(pdfname);
 18    if (neo_pdfs.size() == 0) {
 19        std::cerr << "Failed to load any PDF members!\n";
 20        return 1;
 21    }
 22    std::cout << "Loaded " << neo_pdfs.size() << " PDF members\n";
 23
 24    // Get the first PDF as a reference for metadata
 25    NeoPDF& ref_pdf = neo_pdfs[0];
 26
 27    // Extract the PID values of the PDF set
 28    auto pids = ref_pdf.pids();
 29
 30    // Extract the number of subgrids
 31    std::size_t num_subgrids = ref_pdf.num_subgrids();
 32
 33    // Create a grid writer
 34    GridWriter writer;
 35
 36    // For each member, build a grid
 37    for (size_t m = 0; m < neo_pdfs.size(); ++m) {
 38        NeoPDF& pdf = neo_pdfs[m];
 39
 40        // Start a new grid for the current member
 41        writer.new_grid();
 42
 43        // Loop over the Subgrids
 44        for (std::size_t subgrid_idx = 0; subgrid_idx != num_subgrids; subgrid_idx++) {
 45            // Extract the knot values of the parameters for the subgrid
 46            auto xs = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_MOMENTUM, subgrid_idx);
 47            auto q2s = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_SCALE, subgrid_idx);
 48            auto alphas = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_ALPHAS, subgrid_idx);
 49            auto nucleons = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_NUCLEONS, subgrid_idx);
 50            auto kts = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_KT, subgrid_idx);
 51
 52            // Compute grid_data: [q2][x][flavor], instead of [nucleon][alphas][kt][q2][x][flavor]
 53            // NOTE: This assumes that there is no 'A' and `alphas` dependence.
 54            std::vector<double> grid_data;
 55            for (double x : xs) {
 56                for (double q2 : q2s) {
 57                    for (int pid : pids) {
 58                        double val = pdf.xfxQ2(pid, x, q2);
 59                        grid_data.push_back(val);
 60                    }
 61                }
 62            }
 63
 64            // Add subgrid
 65            writer.add_subgrid(
 66                nucleons,
 67                alphas,
 68                kts,
 69                xs,
 70                q2s,
 71                grid_data
 72            );
 73        }
 74
 75        // Finalize the Grid (inc. its subgrids) for this member.
 76        writer.push_grid(pids);
 77        std::cout << "Added grid for member " << m << "\n";
 78    }
 79
 80    // Fill the running of alphas with some random values
 81    std::vector<double> alphas_qs = {2.0};
 82    std::vector<double> alphas_vals = {0.118};
 83
 84    // Extract the ranges for the momentum x and scale Q2
 85    auto x_range = ref_pdf.param_range(NEOPDF_SUBGRID_PARAMS_MOMENTUM);
 86    auto q2_range = ref_pdf.param_range(NEOPDF_SUBGRID_PARAMS_SCALE);
 87
 88    PhysicsParameters phys_params = {
 89        .flavor_scheme = "variable",
 90        .order_qcd = 2,
 91        .alphas_order_qcd = 2,
 92        .m_w = 80.352,
 93        .m_z = 91.1876,
 94        .m_up = 0.0,
 95        .m_down = 0.0,
 96        .m_strange = 0.0,
 97        .m_charm = 1.51,
 98        .m_bottom = 4.92,
 99        .m_top = 172.5,
100        .alphas_type = "ipol",
101        .number_flavors = 4,
102    };
103
104    MetaData meta = {
105        .set_desc = "NNPDF40_nnlo_as_01180 collection",
106        .set_index = 0,
107        .num_members = (uint32_t)neo_pdfs.size(),
108        .x_min = x_range[0],
109        .x_max = x_range[1],
110        .q_min = sqrt(q2_range[0]),
111        .q_max = sqrt(q2_range[1]),
112        .flavors = pids,
113        .format = "neopdf",
114        .alphas_q_values = alphas_qs,
115        .alphas_vals = alphas_vals,
116        .polarised = false,
117        .set_type = NEOPDF_SET_TYPE_SPACE_LIKE,
118        .interpolator_type = NEOPDF_INTERPOLATOR_TYPE_LOG_BICUBIC,
119        .error_type = "replicas",
120        .hadron_pid = 2212,
121        .phys_params = phys_params,
122    };
123
124    // Check if `NEOPDF_DATA_PATH` is defined and store the Grid there.
125    const char* filename = "check-writer-oop.neopdf.lz4";
126    const char* neopdf_path = std::getenv("NEOPDF_DATA_PATH");
127    std::string output_path = neopdf_path
128        ? std::string(neopdf_path) + (std::string(neopdf_path).back() == '/' ? "" : "/") + filename
129        : filename;
130
131    // Write the PDF Grid into disk
132    try {
133        writer.compress(meta, output_path);
134        std::cout << "Compression succeeded!\n";
135    } catch (const std::runtime_error& err) {
136        std::cerr << "Compression failed: " << err.what() << "\n";
137        return EXIT_FAILURE;
138    }
139
140    // If `NEOPDF_DATA_PATH` is defined, reload the grid and check ther results.
141    if (neopdf_path) {
142        int pid_test = 21;
143        double x_test = 1e-3;
144        double q2_test1 = 1e2;
145        double q2_test2 = 1e4;
146
147        double ref1 = neo_pdfs[0].xfxQ2(pid_test, x_test, q2_test1);
148        double ref2 = neo_pdfs[0].xfxQ2(pid_test, x_test, q2_test2);
149
150        NeoPDF wpdf(filename);
151        double res1 = wpdf.xfxQ2(pid_test, x_test, q2_test1);
152        double res2 = wpdf.xfxQ2(pid_test, x_test, q2_test2);
153
154        assert(std::abs(res1 - ref1) < TOLERANCE);
155        assert(std::abs(res2 - ref2) < TOLERANCE);
156    }
157
158    // Clip the interpolated values to be zero if negatives.
159    neo_pdfs[0].set_force_positive(NEOPDF_FORCE_POSITIVE_CLIP_NEGATIVE);
160    assert(neo_pdfs[0].is_force_positive() == NEOPDF_FORCE_POSITIVE_CLIP_NEGATIVE);
161
162    // Clip all the PDF members to be positive definite
163    neo_pdfs.set_force_positive_members(NEOPDF_FORCE_POSITIVE_CLIP_SMALL);
164    assert(neo_pdfs[4].is_force_positive() == NEOPDF_FORCE_POSITIVE_CLIP_SMALL);
165
166    return EXIT_SUCCESS;
167}

Example 3: Filling TMD grids with \(k_T\) Dependence

In the following example, we are going to see how to fill TMD grids which contains a dependence on the transverse momentum \(k_T\). The following example makes use of the TMDlib library to provide the TMD distributions.

  1#include "neopdf_capi.h"
  2#include "tmdlib/TMDlib.h"
  3#include <NeoPDF.hpp>
  4#include <algorithm>
  5#include <cassert>
  6#include <cmath>
  7#include <cstddef>
  8#include <cstdlib>
  9#include <fstream>
 10#include <initializer_list>
 11#include <string>
 12#include <sstream>
 13#include <tmdlib/factories.h>
 14#include <vector>
 15
 16using namespace TMDlib;
 17using namespace neopdf;
 18
 19const double TOLERANCE = 1e-16;
 20
 21struct Kinematics {
 22    std::vector<double> xs;
 23    std::vector<double> kts;
 24    std::vector<double> qs;
 25};
 26
 27std::vector<double> compute_tmds(TMD& tmd, double x, double kt, double q2) {
 28    double xbar = 0.0;
 29    double mu = sqrt(q2);
 30    std::vector<double> pdfs = tmd.TMDpdf(x, xbar, kt, mu);
 31
 32    return pdfs;
 33}
 34
 35std::vector<double> parse_array(const std::string& line) {
 36    std::vector<double> result;
 37    size_t start = line.find('[');
 38    size_t end = line.find(']');
 39
 40    if (start == std::string::npos || end == std::string::npos) {
 41        return result;
 42    }
 43
 44    double value;
 45    std::string content = line.substr(start + 1, end - start - 1);
 46    std::replace(content.begin(), content.end(), ',', ' ');
 47    std::stringstream ss(content);
 48    while (ss >> value) { result.push_back(value); }
 49
 50    return result;
 51}
 52
 53Kinematics read_kinematics() {
 54    // Parse the Kinematics separately as it is difficult to retrieve
 55    std::ifstream input_kins("MAP22_N3LL.kinematics");
 56
 57    if (!input_kins.is_open()) {
 58        input_kins.open("raw.data");
 59        if (!input_kins.is_open()) {
 60            return {{}, {}, {}};
 61        }
 62    }
 63
 64    std::string line;
 65    std::vector<double> xs, kts, qs;
 66    while (std::getline(input_kins, line)) {
 67        if (line.find("qToQg:") != std::string::npos) { kts = parse_array(line); }
 68        else if (line.find("Qg:") != std::string::npos) { qs = parse_array(line); }
 69        else if (line.find("xg:") != std::string::npos) { xs = parse_array(line); }
 70    }
 71    input_kins.close();
 72
 73    return { xs, kts, qs };
 74}
 75
 76int main() {
 77    std::string setname = "MAP22_grids_FF_Km_N3LL";
 78
 79    TMD tmd;
 80    tmd.setVerbosity(0);
 81
 82    // Extract various informations from the TMD
 83    tmd.TMDinit(setname);
 84    std::size_t n_members = tmd.TMDgetNumMembers();
 85    double xmin = tmd.TMDgetXmin();
 86    double xmax = tmd.TMDgetXmax();
 87    double q2min = tmd.TMDgetQ2min();
 88    double q2max = tmd.TMDgetQ2max();
 89
 90    // Define the kinematics
 91    Kinematics kins = read_kinematics();
 92    std::vector<double> xs = kins.xs;
 93    std::vector<double> qs = kins.qs;
 94    std::vector<double> kts = kins.kts;
 95
 96    // Square the energy scale Q
 97    std::vector<double> q2s(qs.size());
 98    for (size_t i = 0; i < qs.size(); ++i) { q2s[i] = qs[i] * qs[i]; }
 99
100    // Define some physical parameters
101    std::vector<int> pids = {-6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6};
102    std::vector<double> nucleons = { 1.0 }; // assume to be a Proton
103    std::vector<double> alphas = { 0.118 }; // assume to be determined with as=0.118
104
105    // Instantiate NeoPDF grid writer
106    GridWriter neopdf_writer;
107
108    for (std::size_t m = 0; m != n_members; m++) {
109        tmd.TMDinit(setname, m);
110        std::cout << "Member " << m << " loaded!\n";
111
112        // Start a new grid for the current member
113        neopdf_writer.new_grid();
114
115        std::vector<double> grid_data;
116        for (double kt : kts) {
117            for (double x : xs) {
118                for (double q2 : q2s) {
119                    std::vector<double> tmd_pds = tmd.TMDpdf(x, 0.0, kt, sqrt(q2));
120                    for (std::size_t pid_idx = 0; pid_idx != pids.size(); pid_idx++) {
121                        grid_data.push_back(tmd_pds[pid_idx]);
122                    }
123                }
124            }
125        }
126
127        // Add subgrid member to the Grid
128        neopdf_writer.add_subgrid(
129            nucleons,
130            alphas,
131            kts,
132            xs,
133            q2s,
134            grid_data
135        );
136
137        // Finalize the Grid (inc. its subgrids) for this member.
138        neopdf_writer.push_grid(pids);
139    }
140
141    // Fill the running of alphas with some random values
142    std::vector<double> alphas_qs = { 91.1876 };
143    std::vector<double> alphas_vals = { 0.118 };
144
145    // Construct the Metadata
146    PhysicsParameters phys_params = {
147        .flavor_scheme = "fixed",
148        .order_qcd = 2,
149        .alphas_order_qcd = 2,
150        .m_w = 80.352,
151        .m_z = 91.1876,
152        .m_up = 0.0,
153        .m_down = 0.0,
154        .m_strange = 0.0,
155        .m_charm = 1.51,
156        .m_bottom = 4.92,
157        .m_top = 172.5,
158        .alphas_type = "ipol",
159        .number_flavors = 4,
160    };
161
162    MetaData meta = {
163        .set_desc = "NNPDF40_nnlo_as_01180 collection",
164        .set_index = 0,
165        .num_members = (uint32_t)n_members,
166        .x_min = xmin,
167        .x_max = xmax,
168        .q_min = sqrt(q2min),
169        .q_max = sqrt(q2max),
170        .flavors = pids,
171        .format = "neopdf",
172        .alphas_q_values = alphas_qs,
173        .alphas_vals = alphas_vals,
174        .polarised = false,
175        .set_type = NEOPDF_SET_TYPE_SPACE_LIKE,
176        .interpolator_type = NEOPDF_INTERPOLATOR_TYPE_LOG_TRICUBIC,
177        .error_type = "replicas",
178        .hadron_pid = 2212,
179        .phys_params = phys_params,
180    };
181
182    // Check if `NEOPDF_DATA_PATH` is defined and store the Grid there.
183    const char* filename = "check-tmds.neopdf.lz4";
184    const char* neopdf_path = std::getenv("NEOPDF_DATA_PATH");
185    std::string output_path = neopdf_path
186        ? std::string(neopdf_path) + (std::string(neopdf_path).back() == '/' ? "" : "/") + filename
187        : filename;
188
189    // Write the PDF Grid into disk
190    try {
191        neopdf_writer.compress(meta, output_path);
192        std::cout << "Compression succeeded!\n";
193    } catch (const std::runtime_error& err) {
194        std::cerr << "Compression failed: " << err.what() << "\n";
195        return EXIT_FAILURE;
196    }
197
198
199    // If `NEOPDF_DATA_PATH` is defined, reload the grid and check ther results.
200    if (neopdf_path) {
201        int irep = 12;
202        int pid_test_idx = 2;
203        double x_test = xs[20];
204        double q_test = qs[25];
205
206        // Re-load the NeoPDF and TMDLib TMD sets
207        NeoPDF wpdf(filename, irep);
208        tmd.TMDinit(setname, irep);
209
210        for (double kt : kts) {
211            std::vector<double> refs = tmd.TMDpdf(x_test, 0.0, kt, q_test);
212            double ref = refs[pid_test_idx]; // Up Quark
213
214            std::vector<double> params = { kt,x_test, q_test * q_test };
215            double res = wpdf.xfxQ2_ND(pids[pid_test_idx], params);
216
217            assert(std::abs(ref - res) < TOLERANCE);
218        }
219    }
220
221    return EXIT_SUCCESS;
222}