C++ OOP API Example

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.

#include <LHAPDF/PDF.h>
#include <LHAPDF/GridPDF.h>
#include <NeoPDF.hpp>
#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>

using namespace neopdf;

const double TOLERANCE= 1e-16;

template<typename T>
std::vector<T> geomspace(T start, T stop, int num, bool endpoint = false) {
    std::vector<T> result(num);

    if (num == 1) {
        result[0] = start;
        return result;
    }

    T log_start = std::log(start);
    T log_stop = std::log(stop);
    T step = (log_stop - log_start) / (endpoint ? (num - 1) : num);

    for (int i = 0; i < num; ++i) {
        result[i] = std::exp(log_start + i * step);
    }

    return result;
}

template<typename T>
std::vector<T> linspace(T start, T stop, int num, bool endpoint = true) {
    std::vector<T> result(num);

    if (num == 1) {
        result[0] = start;
        return result;
    }

    T step = (stop - start) / (endpoint ? (num - 1) : num);

    for (int i = 0; i < num; ++i) {
        result[i] = start + i * step;
    }

    return result;
}

void test_xfxq2() {
    std::cout << "=== Test xfxQ2 for single PDF member ===\n";

    // disable LHAPDF banners to guarantee deterministic output
    LHAPDF::setVerbosity(0);

    std::string pdfname = "NNPDF40_nnlo_as_01180";
    NeoPDF neo_pdf(pdfname.c_str(), 0);
    const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(pdfname);
    const LHAPDF::GridPDF& lha_pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);

    std::vector<int> pids = {-5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5};
    std::vector<double> xs = geomspace(neo_pdf.x_min(), neo_pdf.x_max(), 200);
    std::vector<double> q2s = geomspace(neo_pdf.q2_min(), neo_pdf.q2_max(), 200);

    // Headers of the table to print the results
    std::cout << std::right
        << std::setw(6) << "pid"
        << std::setw(15) << "x"
        << std::setw(15) << "Q2"
        << std::setw(15) << "LHAPDF"
        << std::setw(15) << "NeoPDF"
        << std::setw(15) << "Rel. Diff." << "\n";
    std::cout << std::string(81, '-') << "\n";

    for (const auto &pid: pids) {
        for (const auto &x: xs) {
            for (const auto &q2: q2s) {
                double expected = lha_pdf.xfxQ2(pid, x, q2);
                double result = neo_pdf.xfxQ2(pid, x, q2);
                double reldif = std::abs(result - expected) / expected;

                assert(std::abs(result - expected) < TOLERANCE);

                // Print the results as a table
                std::cout << std::scientific << std::setprecision(8)
                    << std::right
                    << std::setw(6)  << pid
                    << std::setw(15) << x
                    << std::setw(15) << q2
                    << std::right
                    << std::setw(15) << expected
                    << std::setw(15) << result
                    << std::setw(15) << reldif << "\n";
                }
            }

    }
}

void test_alphas_q2() {
    std::cout << "=== Test alphasQ2 for single PDF member ===\n";

    // disable LHAPDF banners to guarantee deterministic output
    LHAPDF::setVerbosity(0);

    std::string pdfname = "NNPDF40_nnlo_as_01180";
    NeoPDF neo_pdf(pdfname.c_str(), 0);
    const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(pdfname);
    const LHAPDF::GridPDF& lha_pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);

    std::vector<double> q2_points = linspace(4.0, 1e10, 500);

    // Headers of the table to print the results
    std::cout << std::right
        << std::setw(15) << "Q2"
        << std::setw(15) << "LHAPDF"
        << std::setw(15) << "NeoPDF"
        << std::setw(15) << "Rel. Diff." << "\n";
    std::cout << std::string(60, '-') << "\n";

    for (const auto& q2: q2_points) {
        double expected = lha_pdf.alphasQ2(q2);
        double result = neo_pdf.alphasQ2(q2);
        double reldif = std::abs(result - expected) / expected;

        assert(std::abs(result - expected) < TOLERANCE);

        // Print the results as a table
        std::cout << std::scientific << std::setprecision(8)
            << std::right
            << std::setw(15) << q2
            << std::right
            << std::setw(15) << expected
            << std::setw(15) << result
            << std::setw(15) << reldif << "\n";
    }
}

void test_all_pdf_members() {
    std::cout << "=== Test PDFs class (loading all members) ===\n";

    // disable LHAPDF banners to guarantee deterministic output
    LHAPDF::setVerbosity(0);

    std::string pdfname = "NNPDF40_nnlo_as_01180";
    NeoPDFs neo_pdfs(pdfname.c_str());

    std::cout << "Loaded " << neo_pdfs.size() << " PDF members\n";

    // Test case: evaluate a simple point across all members
    int pid = 1;
    double x = 1e-9;
    double q2 = 1.65 * 1.65;

    std::cout << std::right
        << std::setw(8) << "Member"
        << std::setw(15) << "LHAPDF"
        << std::setw(15) << "NeoPDF"
        << std::setw(15) << "Rel. Diff." << "\n";
    std::cout << std::string(53, '-') << "\n";

    // Evaluate the same point across all PDF members
    std::vector<double> results;
    for (size_t i = 0; i < neo_pdfs.size(); ++i) {
        const LHAPDF::PDF* basepdf = LHAPDF::mkPDF(pdfname, i);
        const LHAPDF::GridPDF& lha_pdf = * dynamic_cast<const LHAPDF::GridPDF*>(basepdf);

        double expected = lha_pdf.xfxQ2(pid, x, q2);
        double result = neo_pdfs[i].xfxQ2(pid, x, q2);

        double reldif = std::abs(result - expected) / expected;
        assert(std::abs(result - expected) < TOLERANCE);
        results.push_back(result);

        std::cout << std::right
            << std::setw(8) << i
            << std::scientific << std::setprecision(8)
            << std::setw(15) << expected
            << std::setw(15) << result
            << std::setw(15) << reldif << "\n";
    }

    // Calculate some statistics
    double sum = 0.0;
    for (double result : results) {
        sum += result;
    }
    double mean = sum / results.size();

    double variance = 0.0;
    for (double result : results) {
        variance += (result - mean) * (result - mean);
    }
    variance /= results.size();
    double std_dev = std::sqrt(variance);

    std::cout << "\nStatistics across all members:\n";
    std::cout << "Mean: " << std::scientific << std::setprecision(8) << mean << "\n";
    std::cout << "Std Dev: " << std_dev << "\n";
    std::cout << "Relative Std Dev: " << std_dev / mean << "\n";
}

void test_lazy_loading() {
    std::cout << "=== Test NeoPDFLazy class (lazy loading) ===\n";

    // Disable LHAPDF banners to guarantee deterministic output
    LHAPDF::setVerbosity(0);

    std::string pdfname = "NNPDF40_nnlo_as_01180.neopdf.lz4";
    NeoPDFLazy lazy_pdfs(pdfname);

    std::cout << "Initialized lazy loader for " << pdfname << "\n";

    // Test case: evaluate a simple point across all members
    int pid = 1;
    double x = 1e-9;
    double q2 = 1.65 * 1.65;

    std::cout << std::right
        << std::setw(8) << "Member"
        << std::setw(15) << "LHAPDF"
        << std::setw(15) << "NeoPDF"
        << std::setw(15) << "Rel. Diff." << "\n";
    std::cout << std::string(53, '-') << "\n";

    // Evaluate the same point across all PDF members
    std::vector<double> results;
    int member_idx = 0;
    while (auto neo_pdf = lazy_pdfs.next()) {
        const LHAPDF::PDF* basepdf = LHAPDF::mkPDF("NNPDF40_nnlo_as_01180", member_idx);
        const LHAPDF::GridPDF& lha_pdf = *dynamic_cast<const LHAPDF::GridPDF*>(basepdf);

        double expected = lha_pdf.xfxQ2(pid, x, q2);
        double result = neo_pdf->xfxQ2(pid, x, q2);

        double reldif = std::abs(result - expected) / expected;
        assert(std::abs(result - expected) < TOLERANCE);
        results.push_back(result);

        std::cout << std::right
            << std::setw(8) << member_idx
            << std::scientific << std::setprecision(8)
            << std::setw(15) << expected
            << std::setw(15) << result
            << std::setw(15) << reldif << "\n";
        member_idx++;
    }

    std::cout << "\nSuccessfully iterated through all members lazily.\n";
}

int main() {
    // Test the computation of the PDF interpolations
    test_xfxq2();

    // Test the computation of the `alphas` interpolations
    test_alphas_q2();

    // Test the PDF interpolations by loading all the members
    test_all_pdf_members();

    // Test the lazy loading of PDF members
    test_lazy_loading();

    return EXIT_SUCCESS;
}

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.

NOTE

The following example fills the NeoPDF grid by re-computing the values of the subgrids from another set. This makes it possible to explicitly check that the filling of the grid is correct. However, this makes the codes very verbose. To easily spot the parts that actually fills the grid, some lines are highlighted.

#include <NeoPDF.hpp>
#include <cassert>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <vector>

using namespace neopdf;

const double TOLERANCE= 1e-16;

int main() {
    const char* pdfname = "NNPDF40_nnlo_as_01180";
    // Load all PDF members
    NeoPDFs neo_pdfs(pdfname);
    if (neo_pdfs.size() == 0) {
        std::cerr << "Failed to load any PDF members!\n";
        return 1;
    }
    std::cout << "Loaded " << neo_pdfs.size() << " PDF members\n";

    // Get the first PDF as a reference for metadata
    NeoPDF& ref_pdf = neo_pdfs[0];

    // Extract the PID values of the PDF set
    auto pids = ref_pdf.pids();

    // Extract the number of subgrids
    std::size_t num_subgrids = ref_pdf.num_subgrids();

    // Create a grid writer
    GridWriter writer;

    // For each member, build a grid
    for (size_t m = 0; m < neo_pdfs.size(); ++m) {
        NeoPDF& pdf = neo_pdfs[m];

        // Start a new grid for the current member
        writer.new_grid();

        // Loop over the Subgrids
        for (std::size_t subgrid_idx = 0; subgrid_idx != num_subgrids; subgrid_idx++) {
            // Extract the knot values of the parameters for the subgrid
            auto xs = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_MOMENTUM, subgrid_idx);
            auto q2s = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_SCALE, subgrid_idx);
            auto alphas = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_ALPHAS, subgrid_idx);
            auto nucleons = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_NUCLEONS, subgrid_idx);
            auto kts = pdf.subgrid_for_param(NEOPDF_SUBGRID_PARAMS_KT, subgrid_idx);

            // Compute grid_data: [q2][x][flavor], instead of [nucleon][alphas][kt][q2][x][flavor]
            // NOTE: This assumes that there is no 'A' and `alphas` dependence.
            std::vector<double> grid_data;
            for (double x : xs) {
                for (double q2 : q2s) {
                    for (int pid : pids) {
                        double val = pdf.xfxQ2(pid, x, q2);
                        grid_data.push_back(val);
                    }
                }
            }

            // Add subgrid
            writer.add_subgrid(
                nucleons,
                alphas,
                kts,
                xs,
                q2s,
                grid_data
            );
        }

        // Finalize the Grid (inc. its subgrids) for this member.
        writer.push_grid(pids);
        std::cout << "Added grid for member " << m << "\n";
    }

    // Fill the running of alphas with some random values
    std::vector<double> alphas_qs = {2.0};
    std::vector<double> alphas_vals = {0.118};

    // Extract the ranges for the momentum x and scale Q2
    auto x_range = ref_pdf.param_range(NEOPDF_SUBGRID_PARAMS_MOMENTUM);
    auto q2_range = ref_pdf.param_range(NEOPDF_SUBGRID_PARAMS_SCALE);

    PhysicsParameters phys_params = {
        .flavor_scheme = "variable",
        .order_qcd = 2,
        .alphas_order_qcd = 2,
        .m_w = 80.352,
        .m_z = 91.1876,
        .m_up = 0.0,
        .m_down = 0.0,
        .m_strange = 0.0,
        .m_charm = 1.51,
        .m_bottom = 4.92,
        .m_top = 172.5,
        .alphas_type = "ipol",
        .number_flavors = 4,
    };

    MetaData meta = {
        .set_desc = "NNPDF40_nnlo_as_01180 collection",
        .set_index = 0,
        .num_members = (uint32_t)neo_pdfs.size(),
        .x_min = x_range[0],
        .x_max = x_range[1],
        .q_min = sqrt(q2_range[0]),
        .q_max = sqrt(q2_range[1]),
        .flavors = pids,
        .format = "neopdf",
        .alphas_q_values = alphas_qs,
        .alphas_vals = alphas_vals,
        .polarised = false,
        .set_type = NEOPDF_SET_TYPE_SPACE_LIKE,
        .interpolator_type = NEOPDF_INTERPOLATOR_TYPE_LOG_BICUBIC,
        .error_type = "replicas",
        .hadron_pid = 2212,
        .phys_params = phys_params,
    };

    // Check if `NEOPDF_DATA_PATH` is defined and store the Grid there.
    const char* filename = "check-writer-oop.neopdf.lz4";
    const char* neopdf_path = std::getenv("NEOPDF_DATA_PATH");
    std::string output_path = neopdf_path
        ? std::string(neopdf_path) + (std::string(neopdf_path).back() == '/' ? "" : "/") + filename
        : filename;

    // Write the PDF Grid into disk
    try {
        writer.compress(meta, output_path);
        std::cout << "Compression succeeded!\n";
    } catch (const std::runtime_error& err) {
        std::cerr << "Compression failed: " << err.what() << "\n";
        return EXIT_FAILURE;
    }

    // If `NEOPDF_DATA_PATH` is defined, reload the grid and check ther results.
    if (neopdf_path) {
        int pid_test = 21;
        double x_test = 1e-3;
        double q2_test1 = 1e2;
        double q2_test2 = 1e4;

        double ref1 = neo_pdfs[0].xfxQ2(pid_test, x_test, q2_test1);
        double ref2 = neo_pdfs[0].xfxQ2(pid_test, x_test, q2_test2);

        NeoPDF wpdf(filename);
        double res1 = wpdf.xfxQ2(pid_test, x_test, q2_test1);
        double res2 = wpdf.xfxQ2(pid_test, x_test, q2_test2);

        assert(std::abs(res1 - ref1) < TOLERANCE);
        assert(std::abs(res2 - ref2) < TOLERANCE);
    }

    return EXIT_SUCCESS;
}

Filling Grids that depend on Multiple Parameters

In the case where the PDF grid depends on more parameters, the filling of grid_data in the above example simply now becomes:

std::vector<double> grid_data;
for (double nucleon : nucleons) {
    for (double alpha_s : alphas) {
        for (double x : xs) {
            for (double q2 : q2s) {
                for (int pid : pids) {
                    std::vector<double> params = {nucleon, alpha_s, x, q2};
                    double val = pdf.xfxQ2_ND(pid, params.data());
                    grid_data.push_back(val);
                }
            }
        }
    }
}

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.

#include "neopdf_capi.h"
#include "tmdlib/TMDlib.h"
#include <NeoPDF.hpp>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdlib>
#include <fstream>
#include <initializer_list>
#include <string>
#include <sstream>
#include <tmdlib/factories.h>
#include <vector>

using namespace TMDlib;
using namespace neopdf;

const double TOLERANCE = 1e-16;

struct Kinematics {
    std::vector<double> xs;
    std::vector<double> kts;
    std::vector<double> qs;
};

std::vector<double> compute_tmds(TMD& tmd, double x, double kt, double q2) {
    double xbar = 0.0;
    double mu = sqrt(q2);
    std::vector<double> pdfs = tmd.TMDpdf(x, xbar, kt, mu);

    return pdfs;
}

std::vector<double> parse_array(const std::string& line) {
    std::vector<double> result;
    size_t start = line.find('[');
    size_t end = line.find(']');

    if (start == std::string::npos || end == std::string::npos) {
        return result;
    }

    double value;
    std::string content = line.substr(start + 1, end - start - 1);
    std::replace(content.begin(), content.end(), ',', ' ');
    std::stringstream ss(content);
    while (ss >> value) { result.push_back(value); }

    return result;
}

Kinematics read_kinematics() {
    // Parse the Kinematics separately as it is difficult to retrieve
    std::ifstream input_kins("MAP22_N3LL.kinematics");

    if (!input_kins.is_open()) {
        input_kins.open("raw.data");
        if (!input_kins.is_open()) {
            return {{}, {}, {}};
        }
    }

    std::string line;
    std::vector<double> xs, kts, qs;
    while (std::getline(input_kins, line)) {
        if (line.find("qToQg:") != std::string::npos) { kts = parse_array(line); }
        else if (line.find("Qg:") != std::string::npos) { qs = parse_array(line); }
        else if (line.find("xg:") != std::string::npos) { xs = parse_array(line); }
    }
    input_kins.close();

    return { xs, kts, qs };
}

int main() {
    std::string setname = "MAP22_grids_FF_Km_N3LL";

    TMD tmd;
    tmd.setVerbosity(0);

    // Extract various informations from the TMD
    tmd.TMDinit(setname);
    std::size_t n_members = tmd.TMDgetNumMembers();
    double xmin = tmd.TMDgetXmin();
    double xmax = tmd.TMDgetXmax();
    double q2min = tmd.TMDgetQ2min();
    double q2max = tmd.TMDgetQ2max();

    // Define the kinematics
    Kinematics kins = read_kinematics();
    std::vector<double> xs = kins.xs;
    std::vector<double> qs = kins.qs;
    std::vector<double> kts = kins.kts;

    // Square the energy scale Q
    std::vector<double> q2s(qs.size());
    for (size_t i = 0; i < qs.size(); ++i) { q2s[i] = qs[i] * qs[i]; }

    // Define some physical parameters
    std::vector<int> pids = {-6, -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5, 6};
    std::vector<double> nucleons = { 1.0 }; // assume to be a Proton
    std::vector<double> alphas = { 0.118 }; // assume to be determined with as=0.118

    // Instantiate NeoPDF grid writer
    GridWriter neopdf_writer;

    for (std::size_t m = 0; m != n_members; m++) {
        tmd.TMDinit(setname, m);
        std::cout << "Member " << m << " loaded!" << "\n";

        // Start a new grid for the current member
        neopdf_writer.new_grid();

        std::vector<double> grid_data;
        for (double kt : kts) {
            for (double x : xs) {
                for (double q2 : q2s) {
                    std::vector<double> tmd_pds = tmd.TMDpdf(x, 0.0, kt, sqrt(q2));
                    for (std::size_t pid_idx = 0; pid_idx != pids.size(); pid_idx++) {
                        grid_data.push_back(tmd_pds[pid_idx]);
                    }
                }
            }
        }

        // Add subgrid member to the Grid
        neopdf_writer.add_subgrid(
            nucleons,
            alphas,
            kts,
            xs,
            q2s,
            grid_data
        );

        // Finalize the Grid (inc. its subgrids) for this member.
        neopdf_writer.push_grid(pids);
    }

    // Fill the running of alphas with some random values
    std::vector<double> alphas_qs = { 91.1876 };
    std::vector<double> alphas_vals = { 0.118 };

    // Construct the Metadata
    PhysicsParameters phys_params = {
        .flavor_scheme = "fixed",
        .order_qcd = 2,
        .alphas_order_qcd = 2,
        .m_w = 80.352,
        .m_z = 91.1876,
        .m_up = 0.0,
        .m_down = 0.0,
        .m_strange = 0.0,
        .m_charm = 1.51,
        .m_bottom = 4.92,
        .m_top = 172.5,
        .alphas_type = "ipol",
        .number_flavors = 4,
    };

    MetaData meta = {
        .set_desc = "NNPDF40_nnlo_as_01180 collection",
        .set_index = 0,
        .num_members = (uint32_t)n_members,
        .x_min = xmin,
        .x_max = xmax,
        .q_min = sqrt(q2min),
        .q_max = sqrt(q2max),
        .flavors = pids,
        .format = "neopdf",
        .alphas_q_values = alphas_qs,
        .alphas_vals = alphas_vals,
        .polarised = false,
        .set_type = NEOPDF_SET_TYPE_SPACE_LIKE,
        .interpolator_type = NEOPDF_INTERPOLATOR_TYPE_LOG_TRICUBIC,
        .error_type = "replicas",
        .hadron_pid = 2212,
        .phys_params = phys_params,
    };

    // Check if `NEOPDF_DATA_PATH` is defined and store the Grid there.
    const char* filename = "check-tmds.neopdf.lz4";
    const char* neopdf_path = std::getenv("NEOPDF_DATA_PATH");
    std::string output_path = neopdf_path
        ? std::string(neopdf_path) + (std::string(neopdf_path).back() == '/' ? "" : "/") + filename
        : filename;

    // Write the PDF Grid into disk
    try {
        neopdf_writer.compress(meta, output_path);
        std::cout << "Compression succeeded!\n";
    } catch (const std::runtime_error& err) {
        std::cerr << "Compression failed: " << err.what() << "\n";
        return EXIT_FAILURE;
    }


    // If `NEOPDF_DATA_PATH` is defined, reload the grid and check ther results.
    if (neopdf_path) {
        int irep = 12;
        int pid_test_idx = 2;
        double x_test = xs[20];
        double q_test = qs[25];

        // Re-load the NeoPDF and TMDLib TMD sets
        NeoPDF wpdf(filename, irep);
        tmd.TMDinit(setname, irep);

        for (double kt : kts) {
            std::vector<double> refs = tmd.TMDpdf(x_test, 0.0, kt, q_test);
            double ref = refs[pid_test_idx]; // Up Quark

            std::vector<double> params = { kt,x_test, q_test * q_test };
            double res = wpdf.xfxQ2_ND(pids[pid_test_idx], params);

            assert(std::abs(ref - res) < TOLERANCE);
        }
    }

    return EXIT_SUCCESS;
}