From 107841926d3e37bddea162a883fb4bb3438942fa Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Wed, 1 Apr 2026 21:26:49 +0200 Subject: [PATCH 1/8] Implement scale variations --- docs/grid-creation.md | 89 +++++++++++++++++++++++++++++++++++++++--- pineapfel/src/fill.cpp | 89 ++++++++++++++++++++++++++++++++++++------ 2 files changed, 162 insertions(+), 16 deletions(-) diff --git a/docs/grid-creation.md b/docs/grid-creation.md index 2848690..57fa76e 100644 --- a/docs/grid-creation.md +++ b/docs/grid-creation.md @@ -154,12 +154,15 @@ the same nine channel types. Note that at NNLO the \(gq\) and \(qg\) channels be | 2 | NNLO | \(G_{1,\mathrm{NS}}^{(2)}\) (\(n_f\)-dep.) | \(G_{1,gq}^{(2)}\) (\(n_f\)-dep.) | \(G_{1,qg}^{(2)}\) (\(n_f\)-dep.) | \(G_{1,gg}^{(2)}\) | \(G_{1,\mathrm{PS}}^{(2)}\) | \(G_{1,q\bar{q}}^{(2)}\) | \(G_{1,qpq_{1,2,3}}^{(2)}\) | The orders are specified in the grid card via the `Orders` field. Each order entry is a -5-element array `[alpha_s, alpha, log_xir, log_xif, log_xia]`. For pure QCD coefficient -functions, only the `alpha_s` power matters; the remaining entries should be set to 0. +5-element array `[alpha_s, alpha, log_xir, log_xif, log_xia]`. For central-scale QCD +coefficient functions set all but `alpha_s` to 0. To include renormalization-scale +logarithms (DIS/SIA only), set `log_xir` to the desired power; see +[Renormalization scale variation](#renormalization-scale-variation-dissia) below. -Each entry stores the coefficient function at that **specific** power of \(\alpha_s\), -not the cumulative sum. For a complete NNLO prediction, all three orders (LO, NLO, NNLO) -must be listed so that the grid contains separate subgrids for each perturbative contribution. +Each entry stores the coefficient function at that **specific** power of \(\alpha_s\) +(and \(\ln\xi_R^2\)), not the cumulative sum. For a complete NNLO prediction, all three +orders (LO, NLO, NNLO) must be listed so that the grid contains separate subgrids for +each perturbative contribution. !!! warning Orders beyond NNLO (`alpha_s > 2`) are silently skipped during grid filling, even @@ -755,6 +758,82 @@ values are used. --- +### Renormalization scale variation {#renormalization-scale-variation-dissia} + +PineAPFEL supports renormalization-scale variation for **DIS and SIA** grids. When one or +more orders with `log_xir > 0` are requested, `build_grid()` automatically computes the +corresponding coefficient-function contributions and stores them as separate PineAPPL +subgrids. At convolution time, PineAPPL evaluates + +$$ +F(x, Q^2;\, \xi_R) \;=\; +\sum_{n,\,m} \alpha_s(\mu_R)^n \, +\bigl[\ln\xi_R^2\bigr]^m \, +W_{n,m} +$$ + +where \(\xi_R = \mu_R / Q\) and \(W_{n,m}\) are the stored subgrid weights for order +`[n, 0, m, 0, 0]`. Setting \(\xi_R = 1\) recovers the central-scale result; varying it +around 1 gives the renormalization-scale uncertainty band. + +#### Available renorm-scale orders + +The following `log_xir` orders are derived automatically from the central-scale +coefficient functions \(C_0\), \(C_1\), \(C_2\) via the renormalization group: + +| Order entry | Label | Stored weight | Derivation | +|-------------|-------|---------------|------------| +| `[1, 0, 1, 0, 0]` | NLO × \(\ln\xi_R^2\) | \(\beta_0(n_f)\,\tfrac{1}{4\pi}\,C_0\) | \(\partial_{\ln\mu_R^2} F\big\lvert_{\mathrm{NLO}}\) | +| `[2, 0, 1, 0, 0]` | NNLO × \(\ln\xi_R^2\) | \(\beta_0(n_f)\,\tfrac{1}{4\pi}\,C_1\) | \(\partial_{\ln\mu_R^2} F\big\lvert_{\mathrm{NNLO}}\) | +| `[2, 0, 2, 0, 0]` | NNLO × \(\ln^2\xi_R^2\) | \(\tfrac{\beta_0(n_f)^2}{2}\,\tfrac{1}{(4\pi)^2}\,C_0\) | \(\tfrac{1}{2}\partial^2_{\ln\mu_R^2} F\big\lvert_{\mathrm{NNLO}}\) | + +Here \(\beta_0(n_f) = 11 - \tfrac{2}{3}n_f\) (the one-loop QCD \(\beta\)-function +coefficient) and the factors of \(\tfrac{1}{4\pi}\) absorb the difference between +APFEL++'s \((\alpha_s/4\pi)^n\) convention and PineAPPL's \(\alpha_s^n\) convention, +\(n_f\) is the number of active flavours at the Q² node. + +The central-scale orders `[0,0,0,0,0]`, `[1,0,0,0,0]`, `[2,0,0,0,0]` must also be +included in `Orders` for the corresponding \(C_0\), \(C_1\), \(C_2\) subgrids to exist +(scale-log orders are derived from them). If a base order is absent, the derived log +order will be empty. + +#### Example: DIS NLO + NNLO with renorm-scale variation + +```yaml +Process: DIS +Observable: F2 +Current: NC +PidBasis: PDG +HadronPids: [2212] +ConvolutionTypes: [UNPOL_PDF] + +Orders: + - [0, 0, 0, 0, 0] # LO (C_0) + - [1, 0, 0, 0, 0] # NLO (C_1) + - [2, 0, 0, 0, 0] # NNLO (C_2) + - [1, 0, 1, 0, 0] # NLO × ln ξ_R² + - [2, 0, 1, 0, 0] # NNLO × ln ξ_R² + - [2, 0, 2, 0, 0] # NNLO × ln² ξ_R² + +Bins: + - lower: [10.0, 0.001] + upper: [100.0, 0.01] + +Normalizations: [1.0] +``` + +The six-subgrid grid can then be convoluted at any \(\xi_R\) by passing a non-unity +scale factor to `pineappl_grid_convolve_with_one` (or the Python `Grid.convolve` +wrapper). + +!!! warning "DIS/SIA only — SIDIS and factorization logs not yet implemented" + Renormalization-scale logs are currently filled only for **DIS and SIA** processes. + SIDIS grids ignore `log_xir > 0` entries (no error is raised; the subgrid is simply + left empty). **Factorization-scale logs** (`log_xif > 0`) require convolving with + DGLAP splitting functions and are not yet implemented for any process. + +--- + ### Using the CLI The `pineapfel-build` executable creates a filled PineAPPL grid from three YAML cards: diff --git a/pineapfel/src/fill.cpp b/pineapfel/src/fill.cpp index dae277e..543d2eb 100644 --- a/pineapfel/src/fill.cpp +++ b/pineapfel/src/fill.cpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -10,6 +11,7 @@ #include #include #include +#include #include namespace pineapfel { @@ -1118,16 +1120,18 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, (int)grid_def.channels.size() - (has_gluon_channel ? 1 : 0); // 5. Precompute per-channel operators for each (Q^2, order) node. - // channel_ops[order][channel_idx] accumulates the combined operator - // from all weighted initializers. Gluon and quark channels are - // handled separately to correctly incorporate heavy-quark - // contributions in the massive scheme. + // channel_ops[{alpha_s, log_xir, log_xif}][channel_idx] accumulates + // the combined operator from all weighted initializers. Gluon and + // quark channels are handled separately to correctly incorporate + // heavy-quark contributions in the massive scheme. + // Scale-log orders ({alpha_s, >=1, 0} or {alpha_s, 0, >=1}) are + // derived from the central-scale operators after the main accumulation. + using OrdKey = std::tuple; // (alpha_s, log_xir, log_xif) struct Q2Data { - int nf; + int nf; // 6-entry EW charges (NC) or nf-entry CKM (CC): - // channel_ops[order][channel_idx] = accumulated operator - std::vector charges; - std::map> channel_ops; + std::vector charges; + std::map> channel_ops; }; std::vector q2_data(nq); @@ -1208,7 +1212,8 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, // α_s^ord the result matches APFEL++ BSF which uses // (α_s/4π)^ord. const double pert_norm = std::pow(1.0 / apfel::FourPi, ord); - auto &target = q2_data[iq].channel_ops[ord]; + const OrdKey key_ord{ord, 0, 0}; + auto &target = q2_data[iq].channel_ops[key_ord]; const auto &ops_light = C.at(1).GetObjects(); const apfel::Operator &CNS = @@ -1304,16 +1309,78 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, } } } + + // Renormalization-scale logs derived from central-scale operators. + // Populate channel_ops[{n, m, 0}] for m>0 from the accumulated + // central-scale channel_ops[{n-m, 0, 0}] entries. + // + // Coefficients follow from d/d(tR) [C(ξ_R)] with tR = ln(ξ_R²): + // {1,1,0}: b₀ · C0 · (1/4π)^1 + // {2,1,0}: b₀ · C1 · (1/4π)^2 + // {2,2,0}: (b₀²/2) · C0 · (1/4π)^2 + // where each C_n is already stored in channel_ops[{n,0,0}] with its + // (1/4π)^n factor absorbed, so the extra factor per log-slot is just + // (1/4π)^{alpha_s} / (1/4π)^{base_ord} = (1/4π)^{log_xir}. + // Only computed when the requested order set contains log_xir > 0. + { + bool need_xir = false; + for (const auto &o : grid_def.orders) + if (o.log_xir > 0) { + need_xir = true; + break; + } + + if (need_xir) { + const double b0 = apfel::beta0qcd(q2_data[iq].nf); + const double inv4pi = 1.0 / apfel::FourPi; + auto &cops = q2_data[iq].channel_ops; + + auto it0 = cops.find(OrdKey{0, 0, 0}); + auto it1 = cops.find(OrdKey{1, 0, 0}); + + // {1,1,0}: b₀ * (1/4π) * C0_channels + if (it0 != cops.end()) { + const double fac = b0 * inv4pi; + for (const auto &[ich, op] : it0->second) + add_to_channel(cops[OrdKey{1, 1, 0}], + ich, + op * fac, + 1.0); + } + // {2,1,0}: b₀ * (1/4π) * C1_channels + if (it1 != cops.end()) { + const double fac = b0 * inv4pi; + for (const auto &[ich, op] : it1->second) + add_to_channel(cops[OrdKey{2, 1, 0}], + ich, + op * fac, + 1.0); + } + // {2,2,0}: (b₀²/2) * (1/4π)² * C0_channels + if (it0 != cops.end()) { + const double fac = 0.5 * b0 * b0 * inv4pi * inv4pi; + for (const auto &[ich, op] : it0->second) + add_to_channel(cops[OrdKey{2, 2, 0}], + ich, + op * fac, + 1.0); + } + } + } } // 6. Fill subgrids for each (bin, order, channel) using pre-built operators for (std::size_t iord = 0; iord < grid_def.orders.size(); iord++) { - int alpha_s = grid_def.orders[iord].alpha_s; + const auto &odef = grid_def.orders[iord]; + int alpha_s = odef.alpha_s; if (alpha_s > 2) { std::cerr << " Warning: skipping order alpha_s^" << alpha_s << " (beyond NNLO)" << std::endl; continue; } + const OrdKey ord_key{(int)odef.alpha_s, + (int)odef.log_xir, + (int)odef.log_xif}; for (std::size_t ich = 0; ich < grid_def.channels.size(); ich++) { for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { @@ -1328,7 +1395,7 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, // [x_nodes.front(), x_nodes.back()]. if (x_center >= x_nodes.front() && x_center <= x_nodes.back()) { for (std::size_t iq = 0; iq < nq; iq++) { - auto it_ord = q2_data[iq].channel_ops.find(alpha_s); + auto it_ord = q2_data[iq].channel_ops.find(ord_key); if (it_ord == q2_data[iq].channel_ops.end()) continue; auto it_ch = it_ord->second.find((int)ich); From ba960b281d16373b6be8ef6c81a63acdf16dcb3c Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Thu, 2 Apr 2026 08:58:30 +0200 Subject: [PATCH 2/8] Major change: adopt pointwise computations --- pineapfel/src/fill.cpp | 125 +++++++------ tests/test_f2nc_lhapdf.cpp | 31 +--- tests/test_grid_vs_apfelxx.cpp | 310 +++++++++++++++------------------ 3 files changed, 211 insertions(+), 255 deletions(-) diff --git a/pineapfel/src/fill.cpp b/pineapfel/src/fill.cpp index 543d2eb..94444f5 100644 --- a/pineapfel/src/fill.cpp +++ b/pineapfel/src/fill.cpp @@ -1077,26 +1077,20 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, } const std::size_t nx = x_nodes.size(); - // DIS/SIA/etc. Q² tabulation is independent of OperatorCard SIDIS knobs. - // Those fields apply only in `build_grid_sidis`; reusing them here made - // non-SIDIS grids sensitive to SIDIS defaults and broke legacy tests. - std::vector q2_nodes = derive_q2_nodes(grid_def.bins, - theory.quark_thresholds, - /*n_intermediate=*/3, - /*include_thresholds=*/true, - /*use_bin_centers_only=*/false); - const std::size_t nq = q2_nodes.size(); - - std::cout << " Grid nodes: " << nq << " Q^2 x " << nx << " x/z points" - << std::endl; - - // Concatenated node_values: [q2_0..q2_{nq-1}, x_0..x_{nx-1}] - std::vector node_values; - node_values.reserve(nq + nx); - node_values.insert(node_values.end(), q2_nodes.begin(), q2_nodes.end()); - node_values.insert(node_values.end(), x_nodes.begin(), x_nodes.end()); + // Pointwise design: one Q² node per bin at the geometric bin centre. + // Collect unique Q²_c values to avoid redundant APFEL++ coefficient + // function calls when multiple bins share the same Q²_c. + std::vector bin_q2c(grid_def.bins.size()); + std::set unique_q2c_set; + for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { + bin_q2c[ibin] = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + unique_q2c_set.insert(bin_q2c[ibin]); + } + std::vector unique_q2c(unique_q2c_set.begin(), unique_q2c_set.end()); - std::vector shape = {nq, nx}; + std::cout << " Pointwise grid: " << unique_q2c.size() + << " unique Q² centres x " << nx << " x/z points" << std::endl; bool timelike = (grid_def.process == ProcessType::SIA); bool is_cc = (grid_def.current == Current::CC); @@ -1134,7 +1128,7 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, std::map> channel_ops; }; - std::vector q2_data(nq); + std::map q2_data_map; // Helper: accumulate op*sign into target[ich]. // NOTE: APFEL++ Operator::operator= is infinitely recursive (APFEL++ bug). @@ -1149,14 +1143,15 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, else it->second += op * sign; }; - for (std::size_t iq = 0; iq < nq; iq++) { - double Q = std::sqrt(q2_nodes[iq]); - q2_data[iq].nf = apfel::NF(Q, theory.quark_thresholds); + for (double q2_c : unique_q2c) { + Q2Data &qd = q2_data_map[q2_c]; + double Q = std::sqrt(q2_c); + qd.nf = apfel::NF(Q, theory.quark_thresholds); // Compute per-quark charges and initializer charges std::vector init_charges; if (is_cc) { - int nf = q2_data[iq].nf; + int nf = qd.nf; std::vector weights(nf, 0.0); for (int q = 1; q <= nf; q++) { bool is_down = (q % 2 == 1); @@ -1179,15 +1174,15 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, } weights[q - 1] /= 2.0; } - q2_data[iq].charges = weights; - init_charges = theory.ckm; + qd.charges = weights; + init_charges = theory.ckm; } else { // ElectroWeakCharges always returns 6 entries - q2_data[iq].charges = apfel::ElectroWeakCharges(Q, timelike); - init_charges = q2_data[iq].charges; + qd.charges = apfel::ElectroWeakCharges(Q, timelike); + init_charges = qd.charges; } - const std::vector &charges = q2_data[iq].charges; + const std::vector &charges = qd.charges; for (const auto &wi : weighted_inits) { auto FObjQ = wi.init(Q, init_charges); @@ -1195,7 +1190,7 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, // nf_local: number of "light" quarks for this init // ZM (wi.actnf<0): Q-dependent from thresholds // massive (wi.actnf>=0): fixed count of zero-mass quarks - int nf_local = (wi.actnf < 0) ? q2_data[iq].nf : wi.actnf; + int nf_local = (wi.actnf < 0) ? qd.nf : wi.actnf; // Sum of EW charges over light quarks double sum_ch_light = 0; @@ -1213,7 +1208,7 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, // (α_s/4π)^ord. const double pert_norm = std::pow(1.0 / apfel::FourPi, ord); const OrdKey key_ord{ord, 0, 0}; - auto &target = q2_data[iq].channel_ops[key_ord]; + auto &target = qd.channel_ops[key_ord]; const auto &ops_light = C.at(1).GetObjects(); const apfel::Operator &CNS = @@ -1331,9 +1326,9 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, } if (need_xir) { - const double b0 = apfel::beta0qcd(q2_data[iq].nf); + const double b0 = apfel::beta0qcd(qd.nf); const double inv4pi = 1.0 / apfel::FourPi; - auto &cops = q2_data[iq].channel_ops; + auto &cops = qd.channel_ops; auto it0 = cops.find(OrdKey{0, 0, 0}); auto it1 = cops.find(OrdKey{1, 0, 0}); @@ -1384,46 +1379,50 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, for (std::size_t ich = 0; ich < grid_def.channels.size(); ich++) { for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { - double x_lo = grid_def.bins[ibin].lower.back(); - double x_hi = grid_def.bins[ibin].upper.back(); - double x_center = std::sqrt(x_lo * x_hi); + double x_lo = grid_def.bins[ibin].lower.back(); + double x_hi = grid_def.bins[ibin].upper.back(); + double x_center = std::sqrt(x_lo * x_hi); + + // Per-bin 1-node Q² subgrid at the geometric bin centre. + double q2_c = bin_q2c[ibin]; + std::vector node_vals_bin; + node_vals_bin.reserve(1 + nx); + node_vals_bin.push_back(q2_c); + node_vals_bin.insert(node_vals_bin.end(), + x_nodes.begin(), x_nodes.end()); + std::vector shape_bin = {1, nx}; - std::vector subgrid(nq * nx, 0.0); + std::vector subgrid(nx, 0.0); // NOTE: Skip bins outside the APFEL++ interpolation range; the // operator Evaluate() has undefined behaviour for x outside // [x_nodes.front(), x_nodes.back()]. if (x_center >= x_nodes.front() && x_center <= x_nodes.back()) { - for (std::size_t iq = 0; iq < nq; iq++) { - auto it_ord = q2_data[iq].channel_ops.find(ord_key); - if (it_ord == q2_data[iq].channel_ops.end()) continue; - - auto it_ch = it_ord->second.find((int)ich); - if (it_ch == it_ord->second.end()) continue; - - apfel::Distribution dist = - it_ch->second.Evaluate(x_center); - const std::vector &vals = - dist.GetDistributionJointGrid(); - - // NOTE: APFEL++: K_j satisfy (C⊗f)(x) = Σ K_j*f(x_j) - // where f is the PDF *without* the x factor. PineAPPL - // with APPL_GRID_X divides xfx(x_j) by x_j, giving - // Σ (stored_j) * xf(x_j)/x_j. Multiplying the stored - // kernel by x_j makes PineAPPL compute - // Σ (K_j*x_j) * xf/x_j = Σ K_j * xf(x_j) = F2. - // vals is indexed over the full joint grid (nx_full). - // Map to physical (≤1) indices only. - for (std::size_t iix = 0; iix < nx; iix++) { - const std::size_t ix = phys_x_indices[iix]; - if (ix < vals.size()) - subgrid[iq * nx + iix] = - vals[ix] * x_nodes[iix]; + auto it_q2 = q2_data_map.find(q2_c); + if (it_q2 != q2_data_map.end()) { + auto it_ord = it_q2->second.channel_ops.find(ord_key); + if (it_ord != it_q2->second.channel_ops.end()) { + auto it_ch = it_ord->second.find((int)ich); + if (it_ch != it_ord->second.end()) { + apfel::Distribution dist = + it_ch->second.Evaluate(x_center); + const std::vector &vals = + dist.GetDistributionJointGrid(); + + // NOTE: APFEL++: K_j satisfy (C⊗f)(x)=Σ K_j*f(x_j) + // PineAPPL divides xfx(x_j) by x_j, so store K_j*x_j. + // vals is indexed over the full joint grid (nx_full). + for (std::size_t iix = 0; iix < nx; iix++) { + const std::size_t ix = phys_x_indices[iix]; + if (ix < vals.size()) + subgrid[iix] = vals[ix] * x_nodes[iix]; + } + } } } } - set_subgrid(grid, ibin, iord, ich, node_values, subgrid, shape); + set_subgrid(grid, ibin, iord, ich, node_vals_bin, subgrid, shape_bin); } } } diff --git a/tests/test_f2nc_lhapdf.cpp b/tests/test_f2nc_lhapdf.cpp index 01f18ef..67681ac 100644 --- a/tests/test_f2nc_lhapdf.cpp +++ b/tests/test_f2nc_lhapdf.cpp @@ -111,30 +111,6 @@ int main() { 1.0, // xi_fac predictions.data()); - const int n_intermediate = 3; - std::set q2_set; - for (const auto &bin : grid_def.bins) { - double q2_lo = bin.lower[0], q2_hi = bin.upper[0]; - q2_set.insert(q2_lo); - q2_set.insert(q2_hi); - double log_lo = std::log(q2_lo), log_hi = std::log(q2_hi); - for (int k = 1; k <= n_intermediate; k++) { - double frac = static_cast(k) / (n_intermediate + 1); - q2_set.insert(std::exp(log_lo + frac * (log_hi - log_lo))); - } - } - for (double thr : theory.quark_thresholds) { - double q2t = thr * thr; - if (!q2_set.empty() && q2t > *q2_set.begin() && q2t < *q2_set.rbegin()) - q2_set.insert(q2t); - } - std::vector q2_nodes(q2_set.begin(), q2_set.end()); - - std::printf(" Q² nodes used: %zu\n", q2_nodes.size()); - for (double q2 : q2_nodes) - std::printf(" Q²=%.4e Q=%.4f\n", q2, std::sqrt(q2)); - std::printf("\n"); - const double tolerance = 5e-2; std::printf(" %-10s %-14s %-14s %-10s %s\n", "x", @@ -149,9 +125,10 @@ int main() { double x_hi = grid_def.bins[ibin].upper.back(); double x_c = std::sqrt(x_lo * x_hi); - // Sum BSF reference over the same q2_nodes that PineAPPL sums over. - double ref = 0.0; - for (double q2 : q2_nodes) ref += F2.at(0).Evaluate(x_c, std::sqrt(q2)); + // Pointwise: evaluate BSF at the geometric Q² bin centre. + double q2_c = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + double ref = F2.at(0).Evaluate(x_c, std::sqrt(q2_c)); double pred = predictions[ibin]; double rel = (std::abs(ref) > 1e-30) diff --git a/tests/test_grid_vs_apfelxx.cpp b/tests/test_grid_vs_apfelxx.cpp index 377b894..edc82df 100644 --- a/tests/test_grid_vs_apfelxx.cpp +++ b/tests/test_grid_vs_apfelxx.cpp @@ -161,7 +161,6 @@ int main() { static_cast(tabp.n_steps), tabp.interp_degree}; - auto q2_nodes = derive_q2_nodes(grid_def.bins, theory.quark_thresholds); bool timelike = (grid_def.process == pineapfel::ProcessType::SIA); int failures = 0; @@ -198,21 +197,21 @@ int main() { grid_def.bins[ibin].upper.back()); double ref = 0; - for (double q2 : q2_nodes) { - double Q = std::sqrt(q2); - int nf = apfel::NF(Q, theory.quark_thresholds); - auto charges = apfel::ElectroWeakCharges(Q, timelike); - auto FObjQ = sf_init(Q, charges); - if (FObjQ.C0.count(1) == 0) continue; + // Pointwise: evaluate at the geometric Q² bin centre + double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + int nf = apfel::NF(Q, theory.quark_thresholds); + auto charges = apfel::ElectroWeakCharges(Q, timelike); + auto FObjQ = sf_init(Q, charges); + if (FObjQ.C0.count(1)) { auto ops = FObjQ.C0.at(1).GetObjects(); - for (std::size_t ich = 0; ich < nchan; ich++) { apfel::Operator C_ch = build_channel_operator(grid_def.channels[ich], ops, charges, nf); - const auto &ch = grid_def.channels[ich]; apfel::Distribution pdf_ch(g, [&](double const &z) -> double { @@ -227,7 +226,6 @@ int main() { } return z * sum; }); - ref += (C_ch * pdf_ch).Evaluate(x_center); } } @@ -277,66 +275,67 @@ int main() { grid_def.bins[ibin].upper.back()); double ref = 0; - for (double q2 : q2_nodes) { - double Q = std::sqrt(q2); - int nf = apfel::NF(Q, theory.quark_thresholds); - auto charges = apfel::ElectroWeakCharges(Q, timelike); - auto FObjQ = sf_init(Q, charges); - double as_val = as_tab.Evaluate(Q); - - // LO contribution - if (FObjQ.C0.count(1)) { - auto ops0 = FObjQ.C0.at(1).GetObjects(); - for (std::size_t ich = 0; ich < nchan; ich++) { - apfel::Operator C_ch = - build_channel_operator(grid_def.channels[ich], - ops0, - charges, - nf); - const auto &ch = grid_def.channels[ich]; - apfel::Distribution pdf_ch(g, - [&](double const &z) -> double { - double sum = 0; - for (std::size_t ic = 0; - ic < ch.pid_combinations.size(); - ic++) { - double f_val = 1.0; - for (int pid : ch.pid_combinations[ic]) - f_val *= toy_f(pid, z); - sum += ch.factors[ic] * f_val; - } - return z * sum; // xf convention - }); - ref += (C_ch * pdf_ch).Evaluate(x_center); - } + // Pointwise: evaluate at the geometric Q² bin centre + double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + int nf = apfel::NF(Q, theory.quark_thresholds); + auto charges = apfel::ElectroWeakCharges(Q, timelike); + auto FObjQ = sf_init(Q, charges); + double as_val = as_tab.Evaluate(Q); + + // LO contribution + if (FObjQ.C0.count(1)) { + auto ops0 = FObjQ.C0.at(1).GetObjects(); + for (std::size_t ich = 0; ich < nchan; ich++) { + apfel::Operator C_ch = + build_channel_operator(grid_def.channels[ich], + ops0, + charges, + nf); + const auto &ch = grid_def.channels[ich]; + apfel::Distribution pdf_ch(g, + [&](double const &z) -> double { + double sum = 0; + for (std::size_t ic = 0; + ic < ch.pid_combinations.size(); + ic++) { + double f_val = 1.0; + for (int pid : ch.pid_combinations[ic]) + f_val *= toy_f(pid, z); + sum += ch.factors[ic] * f_val; + } + return z * sum; // xf convention + }); + ref += (C_ch * pdf_ch).Evaluate(x_center); } + } - // NLO contribution (multiplied by alpha_s) - if (FObjQ.C1.count(1)) { - auto ops1 = FObjQ.C1.at(1).GetObjects(); - for (std::size_t ich = 0; ich < nchan; ich++) { - apfel::Operator C_ch = - build_channel_operator(grid_def.channels[ich], - ops1, - charges, - nf); - const auto &ch = grid_def.channels[ich]; - apfel::Distribution pdf_ch(g, - [&](double const &z) -> double { - double sum = 0; - for (std::size_t ic = 0; - ic < ch.pid_combinations.size(); - ic++) { - double f_val = 1.0; - for (int pid : ch.pid_combinations[ic]) - f_val *= toy_f(pid, z); - sum += ch.factors[ic] * f_val; - } - return z * sum; - }); - ref += (as_val / apfel::FourPi) * - (C_ch * pdf_ch).Evaluate(x_center); - } + // NLO contribution (multiplied by alpha_s) + if (FObjQ.C1.count(1)) { + auto ops1 = FObjQ.C1.at(1).GetObjects(); + for (std::size_t ich = 0; ich < nchan; ich++) { + apfel::Operator C_ch = + build_channel_operator(grid_def.channels[ich], + ops1, + charges, + nf); + const auto &ch = grid_def.channels[ich]; + apfel::Distribution pdf_ch(g, + [&](double const &z) -> double { + double sum = 0; + for (std::size_t ic = 0; + ic < ch.pid_combinations.size(); + ic++) { + double f_val = 1.0; + for (int pid : ch.pid_combinations[ic]) + f_val *= toy_f(pid, z); + sum += ch.factors[ic] * f_val; + } + return z * sum; + }); + ref += (as_val / apfel::FourPi) * + (C_ch * pdf_ch).Evaluate(x_center); } } @@ -380,55 +379,48 @@ int main() { grid_def.bins[ibin].upper.back()); double ref = 0; - for (double q2 : q2_nodes) { - double Q = std::sqrt(q2); - int nf = apfel::NF(Q, theory.quark_thresholds); - auto charges = apfel::ElectroWeakCharges(Q, timelike); - auto FObjQ = sf_init(Q, charges); - double as_val = as_tab.Evaluate(Q); - - // Order 0 (LO), 1 (NLO), 2 (NNLO) with alpha_s^n factors - struct { - int order; - double factor; - std::map> *Cn; - } order_info[3]; - - // We need to handle the maps carefully - for (int iord = 0; iord < 3; iord++) { - double as_power = std::pow(as_val / apfel::FourPi, iord); - std::map ops_map; - - if (iord == 0 && FObjQ.C0.count(1)) - ops_map = FObjQ.C0.at(1).GetObjects(); - else if (iord == 1 && FObjQ.C1.count(1)) - ops_map = FObjQ.C1.at(1).GetObjects(); - else if (iord == 2 && FObjQ.C2.count(1)) - ops_map = FObjQ.C2.at(1).GetObjects(); - else continue; - - for (std::size_t ich = 0; ich < nchan; ich++) { - apfel::Operator C_ch = - build_channel_operator(grid_def.channels[ich], - ops_map, - charges, - nf); - const auto &ch = grid_def.channels[ich]; - apfel::Distribution pdf_ch(g, - [&](double const &z) -> double { - double sum = 0; - for (std::size_t ic = 0; - ic < ch.pid_combinations.size(); - ic++) { - double f_val = 1.0; - for (int pid : ch.pid_combinations[ic]) - f_val *= toy_f(pid, z); - sum += ch.factors[ic] * f_val; - } - return z * sum; - }); - ref += as_power * (C_ch * pdf_ch).Evaluate(x_center); - } + // Pointwise: evaluate at the geometric Q² bin centre + double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + int nf = apfel::NF(Q, theory.quark_thresholds); + auto charges = apfel::ElectroWeakCharges(Q, timelike); + auto FObjQ = sf_init(Q, charges); + double as_val = as_tab.Evaluate(Q); + + for (int iord = 0; iord < 3; iord++) { + double as_power = std::pow(as_val / apfel::FourPi, iord); + std::map ops_map; + + if (iord == 0 && FObjQ.C0.count(1)) + ops_map = FObjQ.C0.at(1).GetObjects(); + else if (iord == 1 && FObjQ.C1.count(1)) + ops_map = FObjQ.C1.at(1).GetObjects(); + else if (iord == 2 && FObjQ.C2.count(1)) + ops_map = FObjQ.C2.at(1).GetObjects(); + else continue; + + for (std::size_t ich = 0; ich < nchan; ich++) { + apfel::Operator C_ch = + build_channel_operator(grid_def.channels[ich], + ops_map, + charges, + nf); + const auto &ch = grid_def.channels[ich]; + apfel::Distribution pdf_ch(g, + [&](double const &z) -> double { + double sum = 0; + for (std::size_t ic = 0; + ic < ch.pid_combinations.size(); + ic++) { + double f_val = 1.0; + for (int pid : ch.pid_combinations[ic]) + f_val *= toy_f(pid, z); + sum += ch.factors[ic] * f_val; + } + return z * sum; + }); + ref += as_power * (C_ch * pdf_ch).Evaluate(x_center); } } @@ -512,12 +504,12 @@ int main() { for (std::size_t ibin = 0; ibin < nbins; ibin++) { double x_center = std::sqrt(grid_def.bins[ibin].lower.back() * grid_def.bins[ibin].upper.back()); - double ref = 0; - for (double q2 : q2_nodes) { - double Q = std::sqrt(q2); - ref += F2_total.Evaluate(x_center, Q); - } + // Pointwise: evaluate BSF at the geometric Q² bin centre + double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 @@ -601,9 +593,6 @@ int main() { nfm); } - auto cc_q2_nodes = - derive_q2_nodes(cc_grid_def.bins, cc_theory.quark_thresholds); - std::vector pineappl_cc(cc_nbins, 0.0); pineappl_grid_convolve_with_one(cc_grid, 2212, @@ -619,12 +608,12 @@ int main() { for (std::size_t ibin = 0; ibin < cc_nbins; ibin++) { double x_center = std::sqrt(cc_grid_def.bins[ibin].lower.back() * cc_grid_def.bins[ibin].upper.back()); - double ref = 0; - for (double q2 : cc_q2_nodes) { - double Q = std::sqrt(q2); - ref += F2CC_total.Evaluate(x_center, Q); - } + // Pointwise: evaluate BSF at the geometric Q² bin centre + double q2 = std::sqrt(cc_grid_def.bins[ibin].lower[0] * + cc_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2CC_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 @@ -817,18 +806,15 @@ int main() { nullptr, // mu_scales pineappl_pol.data()); - auto pol_q2_nodes = - derive_q2_nodes(pol_grid_def.bins, pol_theory.quark_thresholds); - for (std::size_t ibin = 0; ibin < pol_nbins; ibin++) { double x_center = std::sqrt(pol_grid_def.bins[ibin].lower.back() * pol_grid_def.bins[ibin].upper.back()); - double ref = 0; - for (double q2 : pol_q2_nodes) { - double Q = std::sqrt(q2); - ref += g1_total.Evaluate(x_center, Q); - } + // Pointwise: evaluate BSF at the geometric Q² bin centre + double q2 = std::sqrt(pol_grid_def.bins[ibin].lower[0] * + pol_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = g1_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 @@ -993,17 +979,15 @@ int main() { nullptr, pineappl_ffn.data()); - auto ffn_q2_nodes = - derive_q2_nodes(ffn_grid_def.bins, ffn_theory.quark_thresholds); - for (std::size_t ibin = 0; ibin < ffn_nbins; ibin++) { double x_center = std::sqrt(ffn_grid_def.bins[ibin].lower.back() * ffn_grid_def.bins[ibin].upper.back()); - double ref = 0; - for (double q2 : ffn_q2_nodes) { - double Q = std::sqrt(q2); - ref += F2FFN_total.Evaluate(x_center, Q); - } + + // Pointwise: evaluate BSF at the geometric Q² bin centre + double q2 = std::sqrt(ffn_grid_def.bins[ibin].lower[0] * + ffn_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2FFN_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 ? std::abs(pineappl_ffn[ibin] - ref) / std::abs(ref) @@ -1081,17 +1065,15 @@ int main() { nullptr, pineappl_mz.data()); - auto mz_q2_nodes = - derive_q2_nodes(mz_grid_def.bins, mz_theory.quark_thresholds); - for (std::size_t ibin = 0; ibin < mz_nbins; ibin++) { double x_center = std::sqrt(mz_grid_def.bins[ibin].lower.back() * mz_grid_def.bins[ibin].upper.back()); - double ref = 0; - for (double q2 : mz_q2_nodes) { - double Q = std::sqrt(q2); - ref += F2MZ_total.Evaluate(x_center, Q); - } + + // Pointwise: evaluate BSF at the geometric Q² bin centre + double q2 = std::sqrt(mz_grid_def.bins[ibin].lower[0] * + mz_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2MZ_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 ? std::abs(pineappl_mz[ibin] - ref) / std::abs(ref) @@ -1187,20 +1169,18 @@ int main() { nullptr, pineappl_fonll.data()); - auto fonll_q2_nodes = - derive_q2_nodes(fonll_grid_def.bins, fonll_theory.quark_thresholds); - for (std::size_t ibin = 0; ibin < fonll_nbins; ibin++) { double x_center = std::sqrt(fonll_grid_def.bins[ibin].lower.back() * fonll_grid_def.bins[ibin].upper.back()); - double ref = 0; - for (double q2 : fonll_q2_nodes) { - double Q = std::sqrt(q2); - // F_FONLL = F_ZM + F_FFN - F_MZ - ref += F2ZM_tot.Evaluate(x_center, Q) + - F2FFN_tot.Evaluate(x_center, Q) - - F2MZ_tot.Evaluate(x_center, Q); - } + + // Pointwise: evaluate BSF at the geometric Q² bin centre + double q2 = std::sqrt(fonll_grid_def.bins[ibin].lower[0] * + fonll_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + // F_FONLL = F_ZM + F_FFN - F_MZ + double ref = F2ZM_tot.Evaluate(x_center, Q) + + F2FFN_tot.Evaluate(x_center, Q) - + F2MZ_tot.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 ? std::abs(pineappl_fonll[ibin] - ref) / std::abs(ref) From 20f361eab3dbc5d16cf0b669180dd6f3dce0fbe8 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 3 Apr 2026 12:04:13 +0200 Subject: [PATCH 3/8] Modify the input cards config to take pointwise kinematics --- pineapfel/include/pineapfel/cards.h | 19 +--- pineapfel/src/cards.cpp | 18 ---- pineapfel/src/fill.cpp | 147 +++++++++------------------- pineapfel/src/grid_gen.cpp | 36 +++++-- pineapfel_pyapi/src/_pineapfel.cpp | 10 +- runcards/grid_dis.yaml | 20 ++-- runcards/grid_dis_cc.yaml | 20 ++-- runcards/grid_dis_ffn.yaml | 20 ++-- runcards/grid_dis_fonll.yaml | 20 ++-- runcards/grid_dis_mz.yaml | 20 ++-- runcards/grid_dis_pol.yaml | 20 ++-- runcards/grid_dis_unp_nc.yaml | 34 +++---- runcards/grid_sia.yaml | 11 +-- runcards/grid_sidis.yaml | 21 ++-- runcards/grid_sidis_nnlo.yaml | 19 ++-- runcards/grid_sidis_pol.yaml | 21 ++-- tests/sidis_helper.cpp | 28 +++--- tests/sidis_helper.h | 8 +- tests/test_grid_vs_apfelxx.cpp | 50 +++------- 19 files changed, 189 insertions(+), 353 deletions(-) diff --git a/pineapfel/include/pineapfel/cards.h b/pineapfel/include/pineapfel/cards.h index 8b026b1..eaee61a 100644 --- a/pineapfel/include/pineapfel/cards.h +++ b/pineapfel/include/pineapfel/cards.h @@ -48,26 +48,11 @@ struct OperatorCard { std::vector xi; // Relative integration precision for SIDIS NNLO coefficient functions. // Controls the cost of InitializeSidisObjects. - double sidis_int_eps = 1e-3; + double sidis_int_eps = 1e-3; // Controls the SIDIS exact (BSF) grid filling implementation. // - "legacy": use the previous external z quadrature approach. // - "bsf_exact": use APFEL++-consistent z integration via IntInterpolant. - std::string sidis_mode = "legacy"; - // TODO: Double-check this in case mass thresholds are present. - // SIDIS fill controls for Q^2/z integration strategy. - // `q2_n_intermediate`: - // number of geometric interior Q^2 nodes per bin edge interval. - // `q2_include_thresholds`: - // insert heavy-quark threshold Q^2 points in global node set. - // `q2_use_bin_centers_only`: - // ignore edge strategy and use only bin geometric centers. - // `z_quad_subdivisions`: - // split each z-bin into this many sub-intervals, each integrated - // with the internal 8-point Gauss-Legendre rule. - int sidis_q2_n_intermediate = 3; - bool sidis_q2_include_thresholds = true; - bool sidis_q2_use_bin_centers_only = false; - int sidis_z_quad_subdivisions = 1; + std::string sidis_mode = "legacy"; }; TheoryCard load_theory_card(const std::string &path); diff --git a/pineapfel/src/cards.cpp b/pineapfel/src/cards.cpp index a57975e..b41fe00 100644 --- a/pineapfel/src/cards.cpp +++ b/pineapfel/src/cards.cpp @@ -93,24 +93,6 @@ OperatorCard load_operator_card(const std::string &path) { if (config["sidis_int_eps"]) oc.sidis_int_eps = config["sidis_int_eps"].as(); - if (config["sidis_q2_n_intermediate"]) - oc.sidis_q2_n_intermediate = - config["sidis_q2_n_intermediate"].as(); - if (config["sidis_q2_include_thresholds"]) - oc.sidis_q2_include_thresholds = - config["sidis_q2_include_thresholds"].as(); - if (config["sidis_q2_use_bin_centers_only"]) - oc.sidis_q2_use_bin_centers_only = - config["sidis_q2_use_bin_centers_only"].as(); - if (config["sidis_z_quad_subdivisions"]) - oc.sidis_z_quad_subdivisions = - config["sidis_z_quad_subdivisions"].as(); - - if (oc.sidis_q2_n_intermediate < 0) - throw std::runtime_error("sidis_q2_n_intermediate must be >= 0"); - if (oc.sidis_z_quad_subdivisions < 1) - throw std::runtime_error("sidis_z_quad_subdivisions must be >= 1"); - return oc; } diff --git a/pineapfel/src/fill.cpp b/pineapfel/src/fill.cpp index 94444f5..e6a8099 100644 --- a/pineapfel/src/fill.cpp +++ b/pineapfel/src/fill.cpp @@ -225,46 +225,6 @@ static std::vector select_initializers(ProcessType process, throw std::runtime_error("select_initializers: unhandled mass scheme"); } -// Collect unique Q^2 nodes from bin edges with geometric intermediate points. -static std::vector derive_q2_nodes(const std::vector &bins, - const std::vector &thresholds, - int n_intermediate = 3, - bool include_thresholds = true, - bool use_bin_centers_only = false) { - std::set q2_set; - - for (const auto &bin : bins) { - double q2_lo = bin.lower[0]; - double q2_hi = bin.upper[0]; - if (use_bin_centers_only) { - q2_set.insert(std::sqrt(q2_lo * q2_hi)); - continue; - } - - q2_set.insert(q2_lo); - q2_set.insert(q2_hi); - - if (n_intermediate > 0) { - double log_lo = std::log(q2_lo); - double log_hi = std::log(q2_hi); - for (int i = 1; i <= n_intermediate; i++) { - double frac = static_cast(i) / (n_intermediate + 1); - q2_set.insert(std::exp(log_lo + frac * (log_hi - log_lo))); - } - } - } - - if (include_thresholds && !q2_set.empty()) { - double q2_min = *q2_set.begin(); - double q2_max = *q2_set.rbegin(); - for (double thr : thresholds) { - double q2_thr = thr * thr; - if (q2_thr > q2_min && q2_thr < q2_max) q2_set.insert(q2_thr); - } - } - - return std::vector(q2_set.begin(), q2_set.end()); -} // Evaluate a DoubleOperator at continuous output point (x_c, z_c) and return // the 2D kernel column w[jx * nj2 + jz] for all joint-grid input nodes (jx, @@ -663,26 +623,22 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, const std::size_t nx = x_nodes.size(); // physical x nodes const std::size_t nz = z_nodes.size(); // physical z nodes - std::vector q2_nodes = derive_q2_nodes(grid_def.bins, - theory.quark_thresholds, - op_card.sidis_q2_n_intermediate, - op_card.sidis_q2_include_thresholds, - op_card.sidis_q2_use_bin_centers_only); - const std::size_t nq = q2_nodes.size(); - - std::cout << " Grid nodes: " << nq << " Q^2 x " << nx << " x x " << nz - << " z points" << std::endl; - - // node_values: [q2_0..q2_{nq-1}, x_0..x_{nx-1}, z_0..z_{nz-1}] - std::vector node_values; - node_values.reserve(nq + nx + nz); - node_values.insert(node_values.end(), q2_nodes.begin(), q2_nodes.end()); - node_values.insert(node_values.end(), x_nodes.begin(), x_nodes.end()); - node_values.insert(node_values.end(), z_nodes.begin(), z_nodes.end()); + // Pointwise Q²: each bin carries exactly one Q² node at its geometric centre. + // For the Points format lower[0]==upper[0], so sqrt gives the same value. + // Collect unique Q² values to avoid redundant APFEL++ calls. + std::vector bin_q2c(grid_def.bins.size()); + std::set unique_q2c_set; + for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { + bin_q2c[ibin] = std::sqrt(grid_def.bins[ibin].lower[0] * + grid_def.bins[ibin].upper[0]); + unique_q2c_set.insert(bin_q2c[ibin]); + } + std::vector unique_q2c(unique_q2c_set.begin(), unique_q2c_set.end()); - std::vector shape = {nq, nx, nz}; + std::cout << " Grid nodes: 1 (pointwise) Q^2 x " << nx << " x x " << nz + << " z points per bin" << std::endl; - // 5. Precompute nf per Q² node. + // 5. Precompute nf per unique Q² value. // IMPORTANT: // APFEL++ sidisbuilder uses fixed electromagnetic charge factors QCh2 // in BuildChannels (not ElectroWeakCharges(Q,false)). @@ -691,10 +647,11 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, struct Q2DataSidis { int nf; }; - std::vector q2_data(nq); - for (std::size_t iq = 0; iq < nq; iq++) { - double Q = std::sqrt(q2_nodes[iq]); - q2_data[iq].nf = apfel::NF(Q, theory.quark_thresholds); + std::map q2_data_map; + for (double q2_c : unique_q2c) { + Q2DataSidis qd; + qd.nf = apfel::NF(std::sqrt(q2_c), theory.quark_thresholds); + q2_data_map[q2_c] = qd; } // Lagrange interpolators for NNLO DoubleOperator evaluation. @@ -756,30 +713,24 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, if (type_id < 0) continue; for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { - // Bin centers: dim 0=Q², dim 1=x, dim 2=z - double x_lo = grid_def.bins[ibin].lower[1]; - double x_hi = grid_def.bins[ibin].upper[1]; - double x_center = std::sqrt(x_lo * x_hi); - - double z_lo = grid_def.bins[ibin].lower[2]; - double z_hi = grid_def.bins[ibin].upper[2]; - const bool z_point_mode = (std::abs(z_hi - z_lo) < 1e-15); - - // When `sidis_q2_use_bin_centers_only` is on, `derive_q2_nodes` - // inserts exactly one Q² per bin (sqrt(lower·upper)), merged - // into a global list. PineAPPL sums the full Q² axis at - // convolution time, so every slice must be zero except the one - // matching this bin's center — otherwise each bin picks up - // every other bin's scale (fatal for multi-bin point grids). - // With the default tabulation (lo/hi + intermediates), centers - // need not appear in `q2_nodes`; keep the legacy fill-all-iq - // path there for backward compatibility. - const double q2_bin_c = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); - const double q2_match_eps = - 1e-9 * std::max(std::abs(q2_bin_c), 1.0); - - std::vector subgrid(nq * nx * nz, 0.0); + // Bin kinematics: Q2 and x are geometric centers, + // z is an integration range [z_lo, z_hi]. + double q2_c = bin_q2c[ibin]; + double x_center = std::sqrt(grid_def.bins[ibin].lower[1] * + grid_def.bins[ibin].upper[1]); + double z_lo = grid_def.bins[ibin].lower[2]; + double z_hi = grid_def.bins[ibin].upper[2]; + const bool z_point_mode = (std::abs(z_hi - z_lo) < 1e-15); + + // Per-bin node values: [q2_c, x_0..x_{nx-1}, z_0..z_{nz-1}] + std::vector node_vals_bin; + node_vals_bin.reserve(1 + nx + nz); + node_vals_bin.push_back(q2_c); + node_vals_bin.insert(node_vals_bin.end(), x_nodes.begin(), x_nodes.end()); + node_vals_bin.insert(node_vals_bin.end(), z_nodes.begin(), z_nodes.end()); + std::vector shape_bin = {1, nx, nz}; + + std::vector subgrid(nx * nz, 0.0); int z_int_lo = 0, z_int_hi = 0; std::vector z_int_w; if (use_bsf_exact) { @@ -803,11 +754,8 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, // interpolation range. External z is handled by quadrature on // [z_lo, z_hi]; individual points must lie on the joint grid. if (x_center >= x_nodes.front() && x_center <= x_nodes.back()) { - for (std::size_t iq = 0; iq < nq; iq++) { - if (op_card.sidis_q2_use_bin_centers_only && - std::abs(q2_nodes[iq] - q2_bin_c) > q2_match_eps) - continue; - int nf = q2_data[iq].nf; + { + int nf = q2_data_map.at(q2_c).nf; static constexpr double QCh2[6] = {1. / 9., 4. / 9., 1. / 9., @@ -876,8 +824,7 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, const std::size_t iz = phys_z_indices[iiz]; const double z_node = z_nodes[iiz]; - subgrid[iq * nx * nz + iix * nz + - iiz] += + subgrid[iix * nz + iiz] += fill_weight * w[ix * nz_full + iz] * (x_node * z_node) / xz_denom; } @@ -914,8 +861,7 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, const std::size_t iz = phys_z_indices[iiz]; const double z_node = z_nodes[iiz]; - subgrid[iq * nx * nz + iix * nz + - iiz] += + subgrid[iix * nz + iiz] += fill_weight * int_w * w[ix * nz_full + iz] * (x_node * z_node) / xz_denom; @@ -943,15 +889,13 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, const std::size_t iz = phys_z_indices[iiz]; const double z_node = z_nodes[iiz]; - subgrid[iq * nx * nz + iix * nz + - iiz] += + subgrid[iix * nz + iiz] += fill_weight * w[ix * nz_full + iz] * (x_node * z_node) / xz_denom; } } } else { - const int z_nsub = - op_card.sidis_z_quad_subdivisions; + const int z_nsub = 1; for (int isub = 0; isub < z_nsub; isub++) { const double z_a = z_lo + (z_hi - z_lo) * (double)isub / @@ -990,8 +934,7 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, phys_z_indices[iiz]; const double z_node = z_nodes[iiz]; - subgrid[iq * nx * nz + - iix * nz + iiz] += + subgrid[iix * nz + iiz] += fill_weight * dz_w * w[ix * nz_full + iz] * (x_node * z_node) / @@ -1005,7 +948,7 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, } } - set_subgrid(grid, ibin, iord, ich, node_values, subgrid, shape); + set_subgrid(grid, ibin, iord, ich, node_vals_bin, subgrid, shape_bin); } } } diff --git a/pineapfel/src/grid_gen.cpp b/pineapfel/src/grid_gen.cpp index 464c709..43ea195 100644 --- a/pineapfel/src/grid_gen.cpp +++ b/pineapfel/src/grid_gen.cpp @@ -258,16 +258,36 @@ GridDef load_grid_def(const std::string &path) { } } - // Bins - for (const auto &bin : config["Bins"]) { - BinDef bdef; - bdef.lower = bin["lower"].as>(); - bdef.upper = bin["upper"].as>(); - def.bins.push_back(std::move(bdef)); + // Bins (legacy) or Points (new pointwise format) + if (config["Points"]) { + for (const auto &pt : config["Points"]) { + auto coords = pt.as>(); + BinDef bdef; + if (coords.size() == 4) { + // SIDIS: [Q2, x, z_lo, z_hi] — Q2 and x pointwise, z is range + bdef.lower = {coords[0], coords[1], coords[2]}; + bdef.upper = {coords[0], coords[1], coords[3]}; + } else { + // DIS/SIA: [Q2, x] — fully pointwise + bdef.lower = coords; + bdef.upper = coords; + } + def.bins.push_back(std::move(bdef)); + } + } else { + for (const auto &bin : config["Bins"]) { + BinDef bdef; + bdef.lower = bin["lower"].as>(); + bdef.upper = bin["upper"].as>(); + def.bins.push_back(std::move(bdef)); + } } - // Normalizations - def.normalizations = config["Normalizations"].as>(); + // Normalizations (optional — defaults to 1.0 per bin) + if (config["Normalizations"]) + def.normalizations = config["Normalizations"].as>(); + else + def.normalizations.assign(def.bins.size(), 1.0); return def; } diff --git a/pineapfel_pyapi/src/_pineapfel.cpp b/pineapfel_pyapi/src/_pineapfel.cpp index b790828..e5c9cb9 100644 --- a/pineapfel_pyapi/src/_pineapfel.cpp +++ b/pineapfel_pyapi/src/_pineapfel.cpp @@ -162,15 +162,7 @@ PYBIND11_MODULE(_pineapfel, m) { .def_readwrite("tabulation", &pineapfel::OperatorCard::tabulation) .def_readwrite("xi", &pineapfel::OperatorCard::xi) .def_readwrite("sidis_mode", &pineapfel::OperatorCard::sidis_mode) - .def_readwrite("sidis_int_eps", &pineapfel::OperatorCard::sidis_int_eps) - .def_readwrite("sidis_q2_n_intermediate", - &pineapfel::OperatorCard::sidis_q2_n_intermediate) - .def_readwrite("sidis_q2_include_thresholds", - &pineapfel::OperatorCard::sidis_q2_include_thresholds) - .def_readwrite("sidis_q2_use_bin_centers_only", - &pineapfel::OperatorCard::sidis_q2_use_bin_centers_only) - .def_readwrite("sidis_z_quad_subdivisions", - &pineapfel::OperatorCard::sidis_z_quad_subdivisions); + .def_readwrite("sidis_int_eps", &pineapfel::OperatorCard::sidis_int_eps); py::class_(m, "GridDef") .def(py::init<>()) diff --git a/runcards/grid_dis.yaml b/runcards/grid_dis.yaml index 2355487..a41a905 100644 --- a/runcards/grid_dis.yaml +++ b/runcards/grid_dis.yaml @@ -10,16 +10,10 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001] - upper: [100.0, 0.001] - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - - lower: [1000.0, 0.30] - upper: [10000.0, 0.80] - - lower: [10000.0, 0.2] - upper: [100000.0, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — geometric centers of the original bins +Points: + - [31.62278, 3.162278e-4] + - [31.62278, 3.162278e-3] + - [316.2278, 3.162278e-2] + - [3162.278, 4.898979e-1] + - [31622.78, 3.741657e-1] diff --git a/runcards/grid_dis_cc.yaml b/runcards/grid_dis_cc.yaml index 9b45431..d4bd707 100644 --- a/runcards/grid_dis_cc.yaml +++ b/runcards/grid_dis_cc.yaml @@ -11,16 +11,10 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001] - upper: [100.0, 0.001] - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - - lower: [1000.0, 0.30] - upper: [10000.0, 0.80] - - lower: [10000.0, 0.2] - upper: [100000.0, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — geometric centers of the original bins +Points: + - [31.62278, 3.162278e-4] + - [31.62278, 3.162278e-3] + - [316.2278, 3.162278e-2] + - [3162.278, 4.898979e-1] + - [31622.78, 3.741657e-1] diff --git a/runcards/grid_dis_ffn.yaml b/runcards/grid_dis_ffn.yaml index d5154bc..a0c874d 100644 --- a/runcards/grid_dis_ffn.yaml +++ b/runcards/grid_dis_ffn.yaml @@ -11,16 +11,10 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001] - upper: [100.0, 0.001] - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - - lower: [1000.0, 0.30] - upper: [10000.0, 0.80] - - lower: [10000.0, 0.2] - upper: [100000.0, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — geometric centers of the original bins +Points: + - [31.62278, 3.162278e-4] + - [31.62278, 3.162278e-3] + - [316.2278, 3.162278e-2] + - [3162.278, 4.898979e-1] + - [31622.78, 3.741657e-1] diff --git a/runcards/grid_dis_fonll.yaml b/runcards/grid_dis_fonll.yaml index dac9537..4e4b914 100644 --- a/runcards/grid_dis_fonll.yaml +++ b/runcards/grid_dis_fonll.yaml @@ -11,16 +11,10 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001] - upper: [100.0, 0.001] - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - - lower: [1000.0, 0.30] - upper: [10000.0, 0.80] - - lower: [10000.0, 0.2] - upper: [100000.0, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — geometric centers of the original bins +Points: + - [31.62278, 3.162278e-4] + - [31.62278, 3.162278e-3] + - [316.2278, 3.162278e-2] + - [3162.278, 4.898979e-1] + - [31622.78, 3.741657e-1] diff --git a/runcards/grid_dis_mz.yaml b/runcards/grid_dis_mz.yaml index a042b35..914388e 100644 --- a/runcards/grid_dis_mz.yaml +++ b/runcards/grid_dis_mz.yaml @@ -11,16 +11,10 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001] - upper: [100.0, 0.001] - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - - lower: [1000.0, 0.30] - upper: [10000.0, 0.80] - - lower: [10000.0, 0.2] - upper: [100000.0, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — geometric centers of the original bins +Points: + - [31.62278, 3.162278e-4] + - [31.62278, 3.162278e-3] + - [316.2278, 3.162278e-2] + - [3162.278, 4.898979e-1] + - [31622.78, 3.741657e-1] diff --git a/runcards/grid_dis_pol.yaml b/runcards/grid_dis_pol.yaml index ce61195..ba94010 100644 --- a/runcards/grid_dis_pol.yaml +++ b/runcards/grid_dis_pol.yaml @@ -10,16 +10,10 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001] - upper: [100.0, 0.001] - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - - lower: [1000.0, 0.30] - upper: [10000.0, 0.80] - - lower: [10000.0, 0.2] - upper: [100000.0, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — geometric centers of the original bins +Points: + - [31.62278, 3.162278e-4] + - [31.62278, 3.162278e-3] + - [316.2278, 3.162278e-2] + - [3162.278, 4.898979e-1] + - [31622.78, 3.741657e-1] diff --git a/runcards/grid_dis_unp_nc.yaml b/runcards/grid_dis_unp_nc.yaml index 9a3f7b7..4942783 100644 --- a/runcards/grid_dis_unp_nc.yaml +++ b/runcards/grid_dis_unp_nc.yaml @@ -10,26 +10,14 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Bins at Q^2=10000 and the x-values from the test -# Narrow bins but finite to ensure numerical stability -Bins: - - lower: [10000.0, 1.0e-5] - upper: [10001.0, 1.00001e-5] - - lower: [10000.0, 1.0e-4] - upper: [10001.0, 1.00001e-4] - - lower: [10000.0, 1.0e-3] - upper: [10001.0, 1.00001e-3] - - lower: [10000.0, 1.0e-2] - upper: [10001.0, 1.00001e-2] - - lower: [10000.0, 1.0e-1] - upper: [10001.0, 1.00001e-1] - - lower: [10000.0, 0.3] - upper: [10001.0, 0.30003] - - lower: [10000.0, 0.5] - upper: [10001.0, 0.50005] - - lower: [10000.0, 0.7] - upper: [10001.0, 0.70007] - - lower: [10000.0, 0.9] - upper: [10001.0, 0.90009] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x] — pointwise at fixed Q^2=10000.5 (midpoint) and selected x values +Points: + - [10000.5, 1.0e-5] + - [10000.5, 1.0e-4] + - [10000.5, 1.0e-3] + - [10000.5, 1.0e-2] + - [10000.5, 1.0e-1] + - [10000.5, 0.3] + - [10000.5, 0.5] + - [10000.5, 0.7] + - [10000.5, 0.9] diff --git a/runcards/grid_sia.yaml b/runcards/grid_sia.yaml index 4a2ffc6..026da24 100644 --- a/runcards/grid_sia.yaml +++ b/runcards/grid_sia.yaml @@ -18,10 +18,7 @@ Channels: - pids: [[21]] factors: [1.0] -Bins: - - lower: [10.0, 0.2] - upper: [100.0, 0.4] - - lower: [100.0, 0.4] - upper: [1000.0, 0.6] - -Normalizations: [1.0, 1.0] +# Points: [Q2, z] — geometric centers of the original bins +Points: + - [31.62278, 2.828427e-1] + - [316.2278, 4.898979e-1] diff --git a/runcards/grid_sidis.yaml b/runcards/grid_sidis.yaml index 239ccd3..cf0c2c1 100644 --- a/runcards/grid_sidis.yaml +++ b/runcards/grid_sidis.yaml @@ -9,16 +9,11 @@ Orders: - [0, 0, 0, 0, 0] - [1, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001, 0.2] - upper: [100.0, 0.001, 0.4] - - lower: [10.0, 0.001, 0.2] - upper: [100.0, 0.01, 0.4] - - lower: [100.0, 0.01, 0.4] - upper: [1000.0, 0.1, 0.6] - - lower: [1000.0, 0.30, 0.2] - upper: [10000.0, 0.80, 0.5] - - lower: [10000.0, 0.2, 0.3] - upper: [100000.0, 0.7, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x, z_lo, z_hi] — Q2 and x are pointwise (geometric centers), +# z is an integration range. +Points: + - [31.62278, 3.162278e-4, 0.2, 0.4] + - [31.62278, 3.162278e-3, 0.2, 0.4] + - [316.2278, 3.162278e-2, 0.4, 0.6] + - [3162.278, 4.898979e-1, 0.2, 0.5] + - [31622.78, 3.741657e-1, 0.3, 0.7] diff --git a/runcards/grid_sidis_nnlo.yaml b/runcards/grid_sidis_nnlo.yaml index 726ce76..387f9cc 100644 --- a/runcards/grid_sidis_nnlo.yaml +++ b/runcards/grid_sidis_nnlo.yaml @@ -10,16 +10,9 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# 3 log-spaced x-bins at fixed Q^2 ~10 GeV^2 and z in [0.2, 0.3]. -# Reduced from 10 to 3 bins to cut build_grid and reference computation time. -# Q range is kept very narrow so all bins share the same nf and alpha_s. -# x and z values are within the APFEL++ interpolation grid [1e-5, 1]. -Bins: - - lower: [10.0, 1.000e-4, 0.2] - upper: [10.001, 2.154e-4, 0.3] - - lower: [10.0, 2.154e-3, 0.2] - upper: [10.001, 4.642e-3, 0.3] - - lower: [10.0, 1.000e-1, 0.2] - upper: [10.001, 4.000e-1, 0.3] - -Normalizations: [1.0, 1.0, 1.0] +# 3 log-spaced x points at Q^2=10 GeV^2 and z in [0.2, 0.3]. +# Points: [Q2, x, z_lo, z_hi] — Q2 and x are pointwise, z is an integration range. +Points: + - [10.0, 1.468e-4, 0.2, 0.3] + - [10.0, 3.163e-3, 0.2, 0.3] + - [10.0, 2.000e-1, 0.2, 0.3] diff --git a/runcards/grid_sidis_pol.yaml b/runcards/grid_sidis_pol.yaml index 58cb01f..869a9ed 100644 --- a/runcards/grid_sidis_pol.yaml +++ b/runcards/grid_sidis_pol.yaml @@ -9,16 +9,11 @@ Orders: - [0, 0, 0, 0, 0] - [1, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.0001, 0.2] - upper: [100.0, 0.001, 0.4] - - lower: [10.0, 0.001, 0.2] - upper: [100.0, 0.01, 0.4] - - lower: [100.0, 0.01, 0.4] - upper: [1000.0, 0.1, 0.6] - - lower: [1000.0, 0.30, 0.2] - upper: [10000.0, 0.80, 0.5] - - lower: [10000.0, 0.2, 0.3] - upper: [100000.0, 0.7, 0.7] - -Normalizations: [1.0, 1.0, 1.0, 1.0, 1.0] +# Points: [Q2, x, z_lo, z_hi] — Q2 and x are pointwise (geometric centers), +# z is an integration range. +Points: + - [31.62278, 3.162278e-4, 0.2, 0.4] + - [31.62278, 3.162278e-3, 0.2, 0.4] + - [316.2278, 3.162278e-2, 0.4, 0.6] + - [3162.278, 4.898979e-1, 0.2, 0.5] + - [31622.78, 3.741657e-1, 0.3, 0.7] diff --git a/tests/sidis_helper.cpp b/tests/sidis_helper.cpp index b556423..dec79a7 100644 --- a/tests/sidis_helper.cpp +++ b/tests/sidis_helper.cpp @@ -151,7 +151,7 @@ static double eval_dop_at(const apfel::Grid &g, std::vector compute_sidis_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, - const std::vector &q2_nodes, + const std::vector &bin_q2, const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, const std::vector & /*charges_dummy*/, @@ -162,8 +162,8 @@ std::vector compute_sidis_reference(const apfel::Grid &g, const apfel::LagrangeInterpolator li{g}; - double q2_max_val = 0; - for (double q2 : q2_nodes) q2_max_val = std::max(q2_max_val, q2); + double q2_max_val = + *std::max_element(bin_q2.begin(), bin_q2.end()); int nf_max = apfel::NF(std::sqrt(q2_max_val), thresholds); auto channels = pineapfel::derive_channels(pineapfel::ProcessType::SIDIS, @@ -189,7 +189,8 @@ std::vector compute_sidis_reference(const apfel::Grid &g, const double z_lo = bin_z_bounds[ibin][0]; const double z_hi = bin_z_bounds[ibin][1]; - for (double q2 : q2_nodes) { + { + double q2 = bin_q2[ibin]; double Q = std::sqrt(q2); int nf = apfel::NF(Q, thresholds); double as_val = alphas_func(Q); @@ -244,7 +245,7 @@ std::vector compute_sidis_reference(const apfel::Grid &g, std::vector compute_sidis_pol_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, - const std::vector &q2_nodes, + const std::vector &bin_q2, const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, int max_alpha_s, @@ -254,8 +255,8 @@ std::vector compute_sidis_pol_reference(const apfel::Grid &g, const apfel::LagrangeInterpolator li{g}; - double q2_max_val = 0; - for (double q2 : q2_nodes) q2_max_val = std::max(q2_max_val, q2); + double q2_max_val = + *std::max_element(bin_q2.begin(), bin_q2.end()); int nf_max = apfel::NF(std::sqrt(q2_max_val), thresholds); auto channels = pineapfel::derive_channels(pineapfel::ProcessType::SIDIS, @@ -281,7 +282,8 @@ std::vector compute_sidis_pol_reference(const apfel::Grid &g, const double z_lo = bin_z_bounds[ibin][0]; const double z_hi = bin_z_bounds[ibin][1]; - for (double q2 : q2_nodes) { + { + double q2 = bin_q2[ibin]; double Q = std::sqrt(q2); int nf = apfel::NF(Q, thresholds); double as_val = alphas_func(Q); @@ -340,7 +342,7 @@ std::vector compute_sidis_nnlo_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, int nf_max, - const std::vector &q2_nodes, + const std::vector &bin_q2, const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, std::function toy_f, @@ -359,7 +361,8 @@ std::vector compute_sidis_nnlo_reference(const apfel::Grid &g, const double z_lo = bin_z_bounds[ibin][0]; const double z_hi = bin_z_bounds[ibin][1]; - for (double q2 : q2_nodes) { + { + double q2 = bin_q2[ibin]; double Q = std::sqrt(q2); int nf = apfel::NF(Q, thresholds); double as_val = alphas_func(Q); @@ -557,7 +560,7 @@ std::vector compute_sidis_g1_nnlo_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, int nf_max, - const std::vector &q2_nodes, + const std::vector &bin_q2, const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, std::function toy_f, @@ -576,7 +579,8 @@ std::vector compute_sidis_g1_nnlo_reference(const apfel::Grid &g, const double z_lo = bin_z_bounds[ibin][0]; const double z_hi = bin_z_bounds[ibin][1]; - for (double q2 : q2_nodes) { + { + double q2 = bin_q2[ibin]; double Q = std::sqrt(q2); int nf = apfel::NF(Q, thresholds); double as_val = alphas_func(Q); diff --git a/tests/sidis_helper.h b/tests/sidis_helper.h index 4d72d39..9a46351 100644 --- a/tests/sidis_helper.h +++ b/tests/sidis_helper.h @@ -13,7 +13,7 @@ std::vector compute_sidis_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, - const std::vector &q2_nodes, + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, // {lo, hi} per bin const std::vector> &bin_z_bounds, // {lo, hi} per bin const std::vector &charges_dummy, // unused placeholder @@ -29,7 +29,7 @@ std::vector compute_sidis_reference(const apfel::Grid &g, std::vector compute_sidis_pol_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, - const std::vector &q2_nodes, + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, // {lo, hi} per bin const std::vector> &bin_z_bounds, // {lo, hi} per bin int max_alpha_s, @@ -45,7 +45,7 @@ std::vector compute_sidis_nnlo_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, int nf_max, - const std::vector &q2_nodes, + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, std::function toy_f, @@ -57,7 +57,7 @@ std::vector compute_sidis_g1_nnlo_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, int nf_max, - const std::vector &q2_nodes, + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, std::function toy_f, diff --git a/tests/test_grid_vs_apfelxx.cpp b/tests/test_grid_vs_apfelxx.cpp index edc82df..97baa82 100644 --- a/tests/test_grid_vs_apfelxx.cpp +++ b/tests/test_grid_vs_apfelxx.cpp @@ -11,7 +11,6 @@ #include #include #include -#include #include // ---------------------------------------------------------------- @@ -87,33 +86,6 @@ static apfel::Operator build_channel_operator( return e_q_sq * CNS + (sum_ch / 6.0) * (CS - CNS); } -// ---------------------------------------------------------------- -// Reproduce derive_q2_nodes from fill.cpp -// ---------------------------------------------------------------- -static std::vector derive_q2_nodes( - const std::vector &bins, - const std::vector &thresholds, - int n_intermediate = 3) { - std::set q2_set; - for (const auto &bin : bins) { - double q2_lo = bin.lower[0], q2_hi = bin.upper[0]; - q2_set.insert(q2_lo); - q2_set.insert(q2_hi); - if (n_intermediate > 0) { - double log_lo = std::log(q2_lo), log_hi = std::log(q2_hi); - for (int i = 1; i <= n_intermediate; i++) { - double frac = static_cast(i) / (n_intermediate + 1); - q2_set.insert(std::exp(log_lo + frac * (log_hi - log_lo))); - } - } - } - double q2_min = *q2_set.begin(), q2_max = *q2_set.rbegin(); - for (double thr : thresholds) { - double q2_thr = thr * thr; - if (q2_thr > q2_min && q2_thr < q2_max) q2_set.insert(q2_thr); - } - return std::vector(q2_set.begin(), q2_set.end()); -} // ================================================================ int main() { @@ -675,8 +647,10 @@ int main() { pineapfel::build_grid(sidis_grid_def, sidis_theory, sidis_op_card); std::size_t sidis_nbins = pineappl_grid_bin_count(sidis_grid); - auto sidis_q2_nodes = - derive_q2_nodes(sidis_grid_def.bins, sidis_theory.quark_thresholds); + // Pointwise: one Q² per bin + std::vector sidis_bin_q2; + for (const auto &bin : sidis_grid_def.bins) + sidis_bin_q2.push_back(std::sqrt(bin.lower[0] * bin.upper[0])); // PineAPPL convolution with 2 convolution functions (PDF and FF). // Both use the same toy_f functional form. @@ -709,7 +683,7 @@ int main() { auto ref_vals = compute_sidis_reference(sidis_g, sidis_sobj, sidis_theory.quark_thresholds, - sidis_q2_nodes, + sidis_bin_q2, bin_x_bounds, bin_z_bounds, {}, @@ -858,8 +832,10 @@ int main() { sidis_op_card); std::size_t pol_sidis_nbins = pineappl_grid_bin_count(pol_sidis_grid); - auto pol_sidis_q2_nodes = derive_q2_nodes(pol_sidis_grid_def.bins, - pol_sidis_theory.quark_thresholds); + // Pointwise: one Q² per bin + std::vector pol_sidis_bin_q2; + for (const auto &bin : pol_sidis_grid_def.bins) + pol_sidis_bin_q2.push_back(bin.lower[0]); // PineAPPL convolution with 2 convolution functions (POL_PDF ⊗ // UNPOL_FF) @@ -892,7 +868,7 @@ int main() { auto pol_ref_vals = compute_sidis_pol_reference(sidis_g, sidis_sobj, pol_sidis_theory.quark_thresholds, - pol_sidis_q2_nodes, + pol_sidis_bin_q2, pol_bin_x_bounds, pol_bin_z_bounds, 1, // max_alpha_s = 1 (LO + NLO) @@ -1231,8 +1207,10 @@ int main() { std::size_t nnlo_nbins = pineappl_grid_bin_count(nnlo_grid); std::size_t nnlo_nords = pineappl_grid_order_count(nnlo_grid); - auto nnlo_q2_nodes = - derive_q2_nodes(nnlo_grid_def.bins, nnlo_theory.quark_thresholds); + // Pointwise: one Q² per bin + std::vector nnlo_q2_nodes; + for (const auto &bin : nnlo_grid_def.bins) + nnlo_q2_nodes.push_back(bin.lower[0]); // PineAPPL convolution — NNLO order only (index 2). auto nnlo_mask = std::make_unique(nnlo_nords); From 849399dcbbc2a759b7190f268142c0b15a2b1ae8 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 3 Apr 2026 19:06:00 +0200 Subject: [PATCH 4/8] Update docs according to pointwise computations --- docs/configuration-cards.md | 23 +-- docs/grid-creation.md | 312 +++++++++++------------------------- 2 files changed, 94 insertions(+), 241 deletions(-) diff --git a/docs/configuration-cards.md b/docs/configuration-cards.md index f82962a..4b19f70 100644 --- a/docs/configuration-cards.md +++ b/docs/configuration-cards.md @@ -115,37 +115,20 @@ tabulation: xi: [1.0, 1.0, 1.0] # --- Optional: SIDIS grid filling only --- -# Omit these for DIS/SIA-only workflows. DIS/SIA Q^2 tabulation in build_grid(). +# Omit these for DIS/SIA-only workflows. # # sidis_mode: legacy | bsf_exact # sidis_int_eps: # APFEL++ InitializeSidisObjects tolerance (default 1e-3) -# sidis_q2_n_intermediate: -# sidis_q2_include_thresholds: -# sidis_q2_use_bin_centers_only: -# sidis_z_quad_subdivisions: ``` #### SIDIS-only YAML fields {#operator-card-sidis-fields} These keys are parsed into `pineapfel::OperatorCard` and affect **`build_grid()` only -when the grid card has `Process: SIDIS`**. They are **ignored for \(Q^2\) node -generation** in DIS and SIA (those processes use fixed legacy \(Q^2\) rules inside -the library). +when the grid card has `Process: SIDIS`**. | YAML key | C++ default (if omitted) | Purpose | |----------|-------------------------|---------| | `sidis_mode` | `legacy` | `legacy` = Gauss–Legendre \(z\) quadrature in each \(z\) bin; `bsf_exact` = APFEL++-style \(z\) weights (`IntInterpolant`) for extended \(z\) bins. | | `sidis_int_eps` | `1e-3` | Relative tolerance for `InitializeSidisObjects` / SIDIS NNLO setup (tests often use `1e-1`). | -| `sidis_q2_n_intermediate` | `3` | Extra geometric \(Q^2\) nodes between each bin’s lower/upper edge (SIDIS global \(Q^2\) list). | -| `sidis_q2_include_thresholds` | `true` | Insert \(m_q^2\) into the SIDIS global \(Q^2\) list when in range. | -| `sidis_q2_use_bin_centers_only` | `false` | If `true`, global \(Q^2\) nodes are only \(\sqrt{Q^2_\mathrm{lo} Q^2_\mathrm{hi}}\) per bin (typical for point \(Q^2\) bins; fill activates one slice per bin). | -| `sidis_z_quad_subdivisions` | `1` | Legacy mode: number of \(z\)-subintervals per bin for internal 8-point GL quadrature. | -For **recommended values** when matching APFEL++ BSF-style SIDIS (and for avoiding -PineAPPL \(Q^2\)-axis over-counting on point grids), set `sidis_mode: bsf_exact` and - -`sidis_q2_n_intermediate: 0`, `sidis_q2_include_thresholds: false`, -`sidis_q2_use_bin_centers_only: true`, `sidis_z_quad_subdivisions: 32` - -explicitly on the operator YAML. Full rationale, workflow, and copy-paste block: -[Grid creation — Operator card: SIDIS-only settings](grid-creation.md#sidis-operator-card). +See [Grid creation — Operator card: SIDIS-only settings](grid-creation.md#sidis-operator-card) for details. diff --git a/docs/grid-creation.md b/docs/grid-creation.md index 57fa76e..5447f15 100644 --- a/docs/grid-creation.md +++ b/docs/grid-creation.md @@ -287,99 +287,23 @@ in the grid card. ### Operator card: SIDIS-only settings (`sidis_*` keys) {#sidis-operator-card} -Grid **cards** describe bins and orders; **operator** cards describe interpolation -grids and several options that only affect **SIDIS** filling. For historical -compatibility, the C++ defaults on `pineapfel::OperatorCard` keep **legacy** SIDIS -tabulation (see table below). That matches older regression tests and coarse -`operator_test.yaml`-style setups. It is **not** the same as APFEL++’s -`BuildSidis*` / `sidisbuilder.h` “BSF-exact” semantics used in modern benchmarks. - -**Therefore: for any SIDIS observable grid that you intend to compare to APFEL++ -reference builders (or to reproduce published SIDIS BSF-style tables), you should -set the following fields explicitly on the operator YAML** used in the same -`build_grid()` call as your SIDIS grid card. Do not rely on C++ defaults for those -computations. - -!!! warning "DIS and SIA ignore `sidis_q2_*` in the library" - For **DIS** and **SIA**, `build_grid()` builds the global \(Q^2\) node list with - **fixed legacy logic** (bin edges, three geometric interior points per \(Q^2\) - interval, and threshold \(Q^2\) values inside the overall range). Those nodes - **do not** read `sidis_q2_n_intermediate`, `sidis_q2_include_thresholds`, or - `sidis_q2_use_bin_centers_only` from the operator card, so shared operator cards - and SIDIS defaults cannot accidentally break DIS/SIA grids or the - `test_grid_vs_apfelxx` suite. - -Only **`build_grid_sidis()`** (invoked automatically when `Process: SIDIS` in the -grid card) consults the SIDIS-specific operator-card fields below. - -#### Recommended operator-card block (BSF-style / APFEL++-aligned SIDIS) - -Use this block (together with your `xgrid` / optional `zgrid` / `tabulation` / `xi`) -when you want SIDIS grids that align with APFEL++ BSF-style treatment and avoid -known pitfalls (duplicate \(Q^2\) slices in PineAPPL convolution for point-like -\(Q^2\) bins, etc.): +Two optional fields in the **operator card** control SIDIS-specific behaviour. +They are ignored for DIS and SIA. + +| Field | Default | Purpose | +|-------|---------|---------| +| `sidis_mode` | `legacy` | `legacy` = Gauss–Legendre \(z\) quadrature inside each \(z\) bin; `bsf_exact` = APFEL++-style interpolation weights (`IntInterpolant`) over the \(z\) range. | +| `sidis_int_eps` | `1e-3` | Relative tolerance for `InitializeSidisObjects` (SIDIS NNLO coefficient-function setup). Increase to `1e-1` for faster but coarser initialisation in tests. | + +For APFEL++-aligned (BSF-style) SIDIS grids, set `sidis_mode: bsf_exact` explicitly +in the operator YAML: ```yaml -# SIDIS fill mode: use APFEL++-consistent z-handling (IntInterpolant weights in -# binned-z mode; point-z mode unchanged). "legacy" keeps the older external -# Gauss–Legendre z-quadrature path inside each z-bin. sidis_mode: bsf_exact - -# Q^2 axis tabulation for the PineAPPL subgrid (first dimension): -# - n_intermediate: extra geometric nodes *between* each bin's Q^2_lo and Q^2_hi -# (0 = none). -# - include_thresholds: insert m_q^2 for flavours whose thresholds lie in the -# global Q^2 range spanned by all bins. -# - use_bin_centers_only: build the global Q^2 list from sqrt(Q^2_lo * Q^2_hi) -# per bin only (typical for (Q^2,x,z) *point* bins). The fill then activates -# only the matching Q^2 slice per bin so PineAPPL does not over-count when it -# sums the Q^2 dimension. -sidis_q2_n_intermediate: 0 -sidis_q2_include_thresholds: false -sidis_q2_use_bin_centers_only: true - -# In sidis_mode: legacy, each z-bin is split into this many sub-intervals; each -# piece gets the built-in 8-point Gauss–Legendre rule. Larger values increase -# smoothness of the z-integration proxy. In sidis_mode: bsf_exact, z-binned -# integration uses APFEL++ interpolation weights instead; this key still affects -# the legacy quadrature path when mode is legacy. -sidis_z_quad_subdivisions: 32 +sidis_int_eps: 1e-3 # optional; 1e-1 is enough for quick tests ``` -A complete real-world example is -[`zurich_sidis_g1_nnlo_point/operator_zurich_exact.yaml`](../zurich_sidis_g1_nnlo_point/operator_zurich_exact.yaml) -in this repository (Zurich-style point kinematics). - -#### Default C++ values when YAML omits these keys - -If the operator YAML does **not** set a field, `load_operator_card()` leaves the -`OperatorCard` defaults from `pineapfel/cards.h`: - -| Field | Default | Role | -|------|---------|------| -| `sidis_mode` | `legacy` | Legacy z-quadrature SIDIS fill | -| `sidis_q2_n_intermediate` | `3` | Legacy dense \(Q^2\) sampling between bin edges | -| `sidis_q2_include_thresholds` | `true` | Legacy: add \(m_q^2\) nodes | -| `sidis_q2_use_bin_centers_only` | `false` | Legacy: use edges + intermediates | -| `sidis_z_quad_subdivisions` | `1` | Legacy: minimal z-bin subdivision | - -Those defaults exist so existing **operator** cards without SIDIS sections keep -previous behaviour and so DIS/SIA paths remain stable. **They are not a substitute -for an explicit BSF block** when your physics goal is APFEL++-consistent SIDIS. - -#### Workflow summary - -1. **DIS / SIA grids** — use any standard operator card; SIDIS keys are ignored for - \(Q^2\) node generation in `build_grid()`. -2. **SIDIS grids, legacy regression / coarse tests** — you may omit SIDIS keys and - inherit defaults (e.g. `runcards/operator_test.yaml` used with - `test_grid_vs_apfelxx`). -3. **SIDIS grids, production or APFEL++ benchmarks** — copy the recommended block - above into the operator YAML (or a dedicated SIDIS operator card) and tune - `sidis_int_eps` separately for speed vs. accuracy of `InitializeSidisObjects`. - -See also [Configuration cards — operator card](configuration-cards.md#operator-card-sidis-fields) -for the same YAML keys in the general operator-card reference. +See also [Configuration cards — operator card](configuration-cards.md#operator-card-sidis-fields). ### The grid card @@ -465,24 +389,29 @@ Orders: # - pids: [[21]] # gluon # factors: [1.0] -# Kinematic bins. Each bin is defined by lower and upper edges -# in each dimension: -# DIS: [Q^2, x] (2 dimensions) -# SIA: [Q^2, z] (2 dimensions) -# SIDIS: [Q^2, x, z] (3 dimensions) -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - -# Bin normalisation factors (one per bin). -Normalizations: [1.0, 1.0] +# Kinematic points. Each entry is a coordinate tuple: +# DIS/SIA: [Q^2, x] — fully pointwise +# SIDIS: [Q^2, x, z_lo, z_hi] — Q^2 and x are pointwise; z is an integration range +Points: + - [10.0, 0.001] + - [100.0, 0.01] + +# Bin normalisation factors (one per bin). Optional — defaults to 1.0 per bin. +# Normalizations: [1.0, 1.0] + +# Legacy alternative to Points (still accepted for backward compatibility). +# Each bin is defined by lower and upper edges in each dimension: +# DIS/SIA: [Q^2, x] (2 dimensions) +# SIDIS: [Q^2, x, z] (3 dimensions) +# The Q^2 node is taken as the geometric centre sqrt(Q^2_lo * Q^2_hi). +# Bins: +# - lower: [10.0, 0.001] +# upper: [100.0, 0.01] ``` #### DIS example -A DIS \(F_2\) grid up to NNLO with two \((Q^2, x)\) bins. +A DIS \(F_2\) grid up to NNLO at two \((Q^2, x)\) points. ```yaml Process: DIS @@ -497,13 +426,9 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001] + - [1000.0, 0.1] ``` #### SIA example @@ -525,19 +450,16 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.2] - upper: [100.0, 0.4] - - lower: [100.0, 0.4] - upper: [1000.0, 0.6] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.2] + - [1000.0, 0.6] ``` #### SIDIS example -A SIDIS \(F_T\) grid for proton→pion semi-inclusive production up to NLO. Bins are -three-dimensional \((Q^2, x, z)\) and two convolution types are required (PDF and FF): +A SIDIS \(F_T\) grid for proton→pion semi-inclusive production up to NLO. Each point +specifies \((Q^2, x)\) exactly and a \(z\) integration range \([z_\mathrm{lo}, +z_\mathrm{hi}]\). Two convolution types are required (PDF and FF): ```yaml Process: SIDIS @@ -551,22 +473,12 @@ Orders: - [0, 0, 0, 0, 0] # LO - [1, 0, 0, 0, 0] # NLO -Bins: - - lower: [10.0, 0.001, 0.2] - upper: [100.0, 0.01, 0.4] - - lower: [100.0, 0.01, 0.4] - upper: [1000.0, 0.1, 0.6] - -Normalizations: [1.0, 1.0] +# [Q^2, x, z_lo, z_hi] — Q^2 and x are pointwise; z is an integration range. +Points: + - [10.0, 0.001, 0.2, 0.4] + - [1000.0, 0.1, 0.4, 0.6] ``` -!!! note "Operator YAML for SIDIS" - The snippet above is only the **grid** card. For APFEL++-consistent BSF-style - SIDIS (recommended for new physics grids), add the `sidis_*` block from - [Operator card: SIDIS-only settings](#sidis-operator-card) to the **operator** - YAML you pass to `build_grid()`. Legacy tests may keep a plain operator card - and inherit C++ defaults. - #### SIDIS NNLO example Adding NNLO extends the channel set from 3 per quark to 9 types (see @@ -586,19 +498,11 @@ Orders: - [1, 0, 0, 0, 0] # NLO - [2, 0, 0, 0, 0] # NNLO (exact, arXiv:2401.16281) -Bins: - - lower: [10.0, 0.001, 0.2] - upper: [100.0, 0.01, 0.4] - - lower: [100.0, 0.01, 0.4] - upper: [1000.0, 0.1, 0.6] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001, 0.2, 0.4] + - [1000.0, 0.1, 0.4, 0.6] ``` -!!! note "Operator YAML for SIDIS" - Add the recommended `sidis_*` fields on the **operator** card when you need - BSF-aligned tabulation; see [Operator card: SIDIS-only settings](#sidis-operator-card). - !!! note The NNLO SIDIS coefficient functions are computed exactly from the full 2D `DoubleExpression` classes in APFEL++ via `InitializeSidisObjects()`. This replaces @@ -624,13 +528,9 @@ Orders: - [0, 0, 0, 0, 0] - [1, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001] + - [1000.0, 0.1] ``` #### Polarized SIDIS example @@ -651,19 +551,11 @@ Orders: - [0, 0, 0, 0, 0] - [1, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.001, 0.2] - upper: [100.0, 0.01, 0.4] - - lower: [100.0, 0.01, 0.4] - upper: [1000.0, 0.1, 0.6] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001, 0.2, 0.4] + - [1000.0, 0.1, 0.4, 0.6] ``` -!!! note "Operator YAML for polarized SIDIS" - Use the [SIDIS operator-card settings](#sidis-operator-card) on the operator - YAML when you need BSF-aligned \(G_1\) grids comparable to APFEL++. - #### FFN DIS example A DIS \(F_2\) grid in the fixed-flavor-number scheme. Requires `HeavyQuarkMasses` in the @@ -685,13 +577,9 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001] + - [1000.0, 0.1] ``` #### FONLL DIS example @@ -713,13 +601,9 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001] + - [1000.0, 0.1] ``` #### CC DIS example @@ -742,13 +626,9 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - - lower: [100.0, 0.01] - upper: [1000.0, 0.1] - -Normalizations: [1.0, 1.0] +Points: + - [10.0, 0.001] + - [1000.0, 0.1] ``` The theory card should include a `CKM` field with 9 squared CKM matrix elements @@ -815,11 +695,8 @@ Orders: - [2, 0, 1, 0, 0] # NNLO × ln ξ_R² - [2, 0, 2, 0, 0] # NNLO × ln² ξ_R² -Bins: - - lower: [10.0, 0.001] - upper: [100.0, 0.01] - -Normalizations: [1.0] +Points: + - [10.0, 0.001] ``` The six-subgrid grid can then be convoluted at any \(\xi_R\) by passing a non-unity @@ -905,64 +782,56 @@ The grid nodes are defined automatically: the `xgrid` definition in the operator card. This ensures consistency between the coefficient function grid and any subsequent evolution step. -- **\(Q^2\) nodes (DIS and SIA)**: Built inside `build_grid()` with **fixed legacy rules** - that do **not** read `sidis_q2_*` from the operator card: collect each bin’s \(Q^2\) - edges, insert three geometrically spaced interior points between edges (per bin - interval), and add threshold \(Q^2 = m_q^2\) for flavours whose thresholds lie in the - global \(Q^2\) range. This keeps DIS/SIA grids stable regardless of SIDIS settings. +- **\(Q^2\) nodes (DIS and SIA)**: Each bin contributes exactly **one** \(Q^2\) node, + taken from the `Points` entry directly. All coefficient functions are evaluated at + that single \(Q^2\) value, and the subgrid has shape `[1, n_x]`. -- **\(Q^2\) nodes (SIDIS)**: Built in `build_grid_sidis()` from the **operator card** - fields `sidis_q2_n_intermediate`, `sidis_q2_include_thresholds`, and - `sidis_q2_use_bin_centers_only` (defaults are legacy-style unless you override them; - see [Operator card: SIDIS-only settings](#sidis-operator-card)). When - `sidis_q2_use_bin_centers_only: true`, only one \(Q^2\) per bin enters the global list - (geometric centre), and the fill keeps a single active slice per bin to avoid - double-counting when PineAPPL sums the \(Q^2\) axis. +- **\(Q^2\) nodes (SIDIS)**: Same pointwise strategy — each bin contributes exactly one + \(Q^2\) node equal to the \(Q^2\) coordinate from the `Points` entry (geometric centre + \(\sqrt{Q^2_\mathrm{lo} \cdot Q^2_\mathrm{hi}}\) when the old `Bins` format is used). + The subgrid has shape `[1, n_x, n_z]`. #### Subgrid layout ##### DIS and SIA Each subgrid (one per combination of bin, perturbative order, and channel) is a -two-dimensional array of shape `[n_Q2, n_x]`, stored in row-major order. The -`node_values` vector concatenates the \(Q^2\) nodes followed by the \(x\)/\(z\) nodes: +two-dimensional array of shape `[1, n_x]`, stored in row-major order. The +`node_values` vector is: ```text -node_values = [Q^2_0, ..., Q^2_{nq-1}, x_0, ..., x_{nx-1}] +node_values = [Q^2, x_0, ..., x_{nx-1}] ``` -For each \(Q^2\) node, the coefficient function operator is evaluated at the bin's -\(x\)/\(z\) centre (geometric mean of the bin edges) to produce a distribution on the -APFEL++ joint grid, which populates one row of the subgrid. +The coefficient function operator is evaluated at the bin's \(Q^2\) and \(x\) point to +produce a distribution on the APFEL++ joint grid, which fills the single row. ##### SIDIS -SIDIS subgrids are three-dimensional arrays of shape `[n_Q2, n_x, n_z]`, stored in +SIDIS subgrids are three-dimensional arrays of shape `[1, n_x, n_z]`, stored in row-major order. The same APFEL++ joint grid is used for both the \(x\) (PDF) and -\(z\) (FF) dimensions. The `node_values` vector concatenates three segments: +\(z\) (FF) dimensions. The `node_values` vector is: ```text -node_values = [Q^2_0, ..., Q^2_{nq-1}, x_0, ..., x_{nx-1}, z_0, ..., z_{nz-1}] +node_values = [Q^2, x_0, ..., x_{nx-1}, z_0, ..., z_{nz-1}] ``` Each coefficient function is stored as a `DoubleOperator` — a full 2D kernel on the -\((x, z)\) grid. For each \(Q^2\) node, the 2D operator is evaluated as a column at -the bin's \((x, z)\) centre (geometric mean of the bin edges) using -`eval_double_op_column()`, which applies 2D Lagrange interpolation weights. The -subgrid entry at `[iq, ix, iz]` accumulates the weighted kernel: +\((x, z)\) grid. The 2D operator is evaluated at the bin's \((Q^2, x)\) point and over +the bin's \(z\) range \([z_\mathrm{lo}, z_\mathrm{hi}]\) using +`eval_double_op_column()`. The subgrid entry at `[ix, iz]` accumulates the weighted +kernel: ```text -subgrid[iq, ix, iz] = fill_weight * W(x[ix], z[iz]; x_centre, z_centre) +subgrid[ix, iz] = fill_weight * W(x[ix], z[iz]; x_point, z_centre) ``` where `fill_weight` is the per-channel electroweak factor (\(e_q^2\), \(\sum_i e_i^2\), or \(e_a e_b\) depending on channel type) and \(W\) is the 2D interpolation kernel -extracted from the `DoubleOperator`. The \(z\) direction may accumulate contributions -from a single centre (\(z\) point bin), from APFEL++-style interpolation weights over -joint-grid nodes (`sidis_mode: bsf_exact` and extended \(z\) bins), or from internal -Gauss–Legendre quadrature on sub-intervals (`sidis_mode: legacy`); see -[SIDIS operator-card settings](#sidis-operator-card). The same DoubleOperator column -machinery is used for all perturbative orders (LO, NLO, NNLO). +extracted from the `DoubleOperator`. The \(z\) direction accumulates contributions +from APFEL++-style interpolation weights over joint-grid nodes (`sidis_mode: bsf_exact`) +or from internal Gauss–Legendre quadrature on the \(z\) range (`sidis_mode: legacy`). +The same DoubleOperator column machinery is used for all perturbative orders (LO, NLO, NNLO). ### Programmatic grid definition @@ -980,9 +849,10 @@ def.convolution_types = {PINEAPPL_CONV_TYPE_UNPOL_PDF}; def.orders = {{0, 0, 0, 0, 0}, {1, 0, 0, 0, 0}, {2, 0, 0, 0, 0}}; // LO + NLO + NNLO // channels are auto-derived by build_grid() — no need to set them +// Pointwise bins: lower == upper == {Q^2, x} def.bins = { - {{10.0, 0.001}, {100.0, 0.01}}, - {{100.0, 0.01}, {1000.0, 0.1}}, + {{10.0, 0.001}, {10.0, 0.001}}, + {{1000.0, 0.1}, {1000.0, 0.1}}, }; def.normalizations = {1.0, 1.0}; From df1356b34ee34038e888edb7d804a5dd40d44d9f Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 3 Apr 2026 19:08:12 +0200 Subject: [PATCH 5/8] Minor updates to the run cards --- runcards/grid_dis.yaml | 2 +- runcards/grid_dis_cc.yaml | 2 +- runcards/grid_dis_ffn.yaml | 2 +- runcards/grid_dis_fonll.yaml | 2 +- runcards/grid_dis_mz.yaml | 2 +- runcards/grid_dis_pol.yaml | 2 +- runcards/grid_dis_unp_nc.yaml | 2 +- runcards/grid_sia.yaml | 2 +- runcards/grid_sidis_nnlo.yaml | 3 ++- 9 files changed, 10 insertions(+), 9 deletions(-) diff --git a/runcards/grid_dis.yaml b/runcards/grid_dis.yaml index a41a905..6bed734 100644 --- a/runcards/grid_dis.yaml +++ b/runcards/grid_dis.yaml @@ -10,7 +10,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — geometric centers of the original bins +# Points: [Q2, x] Points: - [31.62278, 3.162278e-4] - [31.62278, 3.162278e-3] diff --git a/runcards/grid_dis_cc.yaml b/runcards/grid_dis_cc.yaml index d4bd707..b7320f9 100644 --- a/runcards/grid_dis_cc.yaml +++ b/runcards/grid_dis_cc.yaml @@ -11,7 +11,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — geometric centers of the original bins +# Points: [Q2, x] Points: - [31.62278, 3.162278e-4] - [31.62278, 3.162278e-3] diff --git a/runcards/grid_dis_ffn.yaml b/runcards/grid_dis_ffn.yaml index a0c874d..37f1f04 100644 --- a/runcards/grid_dis_ffn.yaml +++ b/runcards/grid_dis_ffn.yaml @@ -11,7 +11,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — geometric centers of the original bins +# Points: [Q2, x] Points: - [31.62278, 3.162278e-4] - [31.62278, 3.162278e-3] diff --git a/runcards/grid_dis_fonll.yaml b/runcards/grid_dis_fonll.yaml index 4e4b914..8361e47 100644 --- a/runcards/grid_dis_fonll.yaml +++ b/runcards/grid_dis_fonll.yaml @@ -11,7 +11,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — geometric centers of the original bins +# Points: [Q2, x] Points: - [31.62278, 3.162278e-4] - [31.62278, 3.162278e-3] diff --git a/runcards/grid_dis_mz.yaml b/runcards/grid_dis_mz.yaml index 914388e..f4beed4 100644 --- a/runcards/grid_dis_mz.yaml +++ b/runcards/grid_dis_mz.yaml @@ -11,7 +11,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — geometric centers of the original bins +# Points: [Q2, x] Points: - [31.62278, 3.162278e-4] - [31.62278, 3.162278e-3] diff --git a/runcards/grid_dis_pol.yaml b/runcards/grid_dis_pol.yaml index ba94010..65e11f6 100644 --- a/runcards/grid_dis_pol.yaml +++ b/runcards/grid_dis_pol.yaml @@ -10,7 +10,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — geometric centers of the original bins +# Points: [Q2, x] Points: - [31.62278, 3.162278e-4] - [31.62278, 3.162278e-3] diff --git a/runcards/grid_dis_unp_nc.yaml b/runcards/grid_dis_unp_nc.yaml index 4942783..bd7dfd5 100644 --- a/runcards/grid_dis_unp_nc.yaml +++ b/runcards/grid_dis_unp_nc.yaml @@ -10,7 +10,7 @@ Orders: - [1, 0, 0, 0, 0] - [2, 0, 0, 0, 0] -# Points: [Q2, x] — pointwise at fixed Q^2=10000.5 (midpoint) and selected x values +# Points: [Q2, x] Points: - [10000.5, 1.0e-5] - [10000.5, 1.0e-4] diff --git a/runcards/grid_sia.yaml b/runcards/grid_sia.yaml index 026da24..fefd28d 100644 --- a/runcards/grid_sia.yaml +++ b/runcards/grid_sia.yaml @@ -18,7 +18,7 @@ Channels: - pids: [[21]] factors: [1.0] -# Points: [Q2, z] — geometric centers of the original bins +# Points: [Q2, z] Points: - [31.62278, 2.828427e-1] - [316.2278, 4.898979e-1] diff --git a/runcards/grid_sidis_nnlo.yaml b/runcards/grid_sidis_nnlo.yaml index 387f9cc..6e4245c 100644 --- a/runcards/grid_sidis_nnlo.yaml +++ b/runcards/grid_sidis_nnlo.yaml @@ -11,7 +11,8 @@ Orders: - [2, 0, 0, 0, 0] # 3 log-spaced x points at Q^2=10 GeV^2 and z in [0.2, 0.3]. -# Points: [Q2, x, z_lo, z_hi] — Q2 and x are pointwise, z is an integration range. +# Points: [Q2, x, z_lo, z_hi] — Q2 and x are pointwise, z +# is an integration range. Points: - [10.0, 1.468e-4, 0.2, 0.3] - [10.0, 3.163e-3, 0.2, 0.3] From 2a1af0c1e35d284ebae7942bdadb15b789381ae3 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 3 Apr 2026 19:29:25 +0200 Subject: [PATCH 6/8] Run `pre-commit` hooks across all files --- pineapfel/src/fill.cpp | 107 +++++++++++++++++------------ pineapfel/src/grid_gen.cpp | 3 +- pineapfel_pyapi/src/_pineapfel.cpp | 3 +- tests/sidis_helper.cpp | 18 +++-- tests/sidis_helper.h | 12 ++-- tests/test_f2nc_lhapdf.cpp | 4 +- tests/test_grid_vs_apfelxx.cpp | 63 +++++++++-------- 7 files changed, 113 insertions(+), 97 deletions(-) diff --git a/pineapfel/src/fill.cpp b/pineapfel/src/fill.cpp index e6a8099..50aec42 100644 --- a/pineapfel/src/fill.cpp +++ b/pineapfel/src/fill.cpp @@ -225,7 +225,6 @@ static std::vector select_initializers(ProcessType process, throw std::runtime_error("select_initializers: unhandled mass scheme"); } - // Evaluate a DoubleOperator at continuous output point (x_c, z_c) and return // the 2D kernel column w[jx * nj2 + jz] for all joint-grid input nodes (jx, // jz). @@ -620,20 +619,21 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, phys_z_indices.push_back(i); } } - const std::size_t nx = x_nodes.size(); // physical x nodes - const std::size_t nz = z_nodes.size(); // physical z nodes + const std::size_t nx = x_nodes.size(); // physical x nodes + const std::size_t nz = z_nodes.size(); // physical z nodes - // Pointwise Q²: each bin carries exactly one Q² node at its geometric centre. - // For the Points format lower[0]==upper[0], so sqrt gives the same value. - // Collect unique Q² values to avoid redundant APFEL++ calls. + // Pointwise Q²: each bin carries exactly one Q² node at its geometric + // centre. For the Points format lower[0]==upper[0], so sqrt gives the same + // value. Collect unique Q² values to avoid redundant APFEL++ calls. std::vector bin_q2c(grid_def.bins.size()); std::set unique_q2c_set; for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { - bin_q2c[ibin] = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + bin_q2c[ibin] = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); unique_q2c_set.insert(bin_q2c[ibin]); } - std::vector unique_q2c(unique_q2c_set.begin(), unique_q2c_set.end()); + std::vector unique_q2c(unique_q2c_set.begin(), + unique_q2c_set.end()); std::cout << " Grid nodes: 1 (pointwise) Q^2 x " << nx << " x x " << nz << " z points per bin" << std::endl; @@ -650,7 +650,7 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, std::map q2_data_map; for (double q2_c : unique_q2c) { Q2DataSidis qd; - qd.nf = apfel::NF(std::sqrt(q2_c), theory.quark_thresholds); + qd.nf = apfel::NF(std::sqrt(q2_c), theory.quark_thresholds); q2_data_map[q2_c] = qd; } @@ -715,24 +715,28 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { // Bin kinematics: Q2 and x are geometric centers, // z is an integration range [z_lo, z_hi]. - double q2_c = bin_q2c[ibin]; - double x_center = std::sqrt(grid_def.bins[ibin].lower[1] * - grid_def.bins[ibin].upper[1]); - double z_lo = grid_def.bins[ibin].lower[2]; - double z_hi = grid_def.bins[ibin].upper[2]; + double q2_c = bin_q2c[ibin]; + double x_center = std::sqrt(grid_def.bins[ibin].lower[1] * + grid_def.bins[ibin].upper[1]); + double z_lo = grid_def.bins[ibin].lower[2]; + double z_hi = grid_def.bins[ibin].upper[2]; const bool z_point_mode = (std::abs(z_hi - z_lo) < 1e-15); // Per-bin node values: [q2_c, x_0..x_{nx-1}, z_0..z_{nz-1}] std::vector node_vals_bin; node_vals_bin.reserve(1 + nx + nz); node_vals_bin.push_back(q2_c); - node_vals_bin.insert(node_vals_bin.end(), x_nodes.begin(), x_nodes.end()); - node_vals_bin.insert(node_vals_bin.end(), z_nodes.begin(), z_nodes.end()); + node_vals_bin.insert(node_vals_bin.end(), + x_nodes.begin(), + x_nodes.end()); + node_vals_bin.insert(node_vals_bin.end(), + z_nodes.begin(), + z_nodes.end()); std::vector shape_bin = {1, nx, nz}; - std::vector subgrid(nx * nz, 0.0); - int z_int_lo = 0, z_int_hi = 0; - std::vector z_int_w; + std::vector subgrid(nx * nz, 0.0); + int z_int_lo = 0, z_int_hi = 0; + std::vector z_int_w; if (use_bsf_exact) { if (!z_point_mode) { const apfel::SubGrid &z_joint_sg = gz.GetJointGrid(); @@ -755,7 +759,7 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, // [z_lo, z_hi]; individual points must lie on the joint grid. if (x_center >= x_nodes.front() && x_center <= x_nodes.back()) { { - int nf = q2_data_map.at(q2_c).nf; + int nf = q2_data_map.at(q2_c).nf; static constexpr double QCh2[6] = {1. / 9., 4. / 9., 1. / 9., @@ -948,7 +952,13 @@ static pineappl_grid *build_grid_sidis(const GridDef &grid_def_in, } } - set_subgrid(grid, ibin, iord, ich, node_vals_bin, subgrid, shape_bin); + set_subgrid(grid, + ibin, + iord, + ich, + node_vals_bin, + subgrid, + shape_bin); } } } @@ -1018,7 +1028,7 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, phys_x_indices.push_back(i); } } - const std::size_t nx = x_nodes.size(); + const std::size_t nx = x_nodes.size(); // Pointwise design: one Q² node per bin at the geometric bin centre. // Collect unique Q²_c values to avoid redundant APFEL++ coefficient @@ -1026,23 +1036,24 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, std::vector bin_q2c(grid_def.bins.size()); std::set unique_q2c_set; for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { - bin_q2c[ibin] = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + bin_q2c[ibin] = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); unique_q2c_set.insert(bin_q2c[ibin]); } - std::vector unique_q2c(unique_q2c_set.begin(), unique_q2c_set.end()); + std::vector unique_q2c(unique_q2c_set.begin(), + unique_q2c_set.end()); std::cout << " Pointwise grid: " << unique_q2c.size() << " unique Q² centres x " << nx << " x/z points" << std::endl; - bool timelike = (grid_def.process == ProcessType::SIA); - bool is_cc = (grid_def.current == Current::CC); + bool timelike = (grid_def.process == ProcessType::SIA); + bool is_cc = (grid_def.current == Current::CC); // Pre-determine gluon channel index (gluon is the last channel when // present) - int gluon_ich = -1; - bool has_gluon_channel = false; - int num_quark_channels = 0; + int gluon_ich = -1; + bool has_gluon_channel = false; + int num_quark_channels = 0; for (int ich = 0; ich < (int)grid_def.channels.size(); ich++) { for (const auto &combo : grid_def.channels[ich].pid_combinations) { if (combo.size() == 1 && combo[0] == 21) { @@ -1087,9 +1098,9 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, }; for (double q2_c : unique_q2c) { - Q2Data &qd = q2_data_map[q2_c]; - double Q = std::sqrt(q2_c); - qd.nf = apfel::NF(Q, theory.quark_thresholds); + Q2Data &qd = q2_data_map[q2_c]; + double Q = std::sqrt(q2_c); + qd.nf = apfel::NF(Q, theory.quark_thresholds); // Compute per-quark charges and initializer charges std::vector init_charges; @@ -1322,20 +1333,21 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, for (std::size_t ich = 0; ich < grid_def.channels.size(); ich++) { for (std::size_t ibin = 0; ibin < grid_def.bins.size(); ibin++) { - double x_lo = grid_def.bins[ibin].lower.back(); - double x_hi = grid_def.bins[ibin].upper.back(); - double x_center = std::sqrt(x_lo * x_hi); + double x_lo = grid_def.bins[ibin].lower.back(); + double x_hi = grid_def.bins[ibin].upper.back(); + double x_center = std::sqrt(x_lo * x_hi); // Per-bin 1-node Q² subgrid at the geometric bin centre. - double q2_c = bin_q2c[ibin]; + double q2_c = bin_q2c[ibin]; std::vector node_vals_bin; node_vals_bin.reserve(1 + nx); node_vals_bin.push_back(q2_c); node_vals_bin.insert(node_vals_bin.end(), - x_nodes.begin(), x_nodes.end()); + x_nodes.begin(), + x_nodes.end()); std::vector shape_bin = {1, nx}; - std::vector subgrid(nx, 0.0); + std::vector subgrid(nx, 0.0); // NOTE: Skip bins outside the APFEL++ interpolation range; the // operator Evaluate() has undefined behaviour for x outside @@ -1352,9 +1364,10 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, const std::vector &vals = dist.GetDistributionJointGrid(); - // NOTE: APFEL++: K_j satisfy (C⊗f)(x)=Σ K_j*f(x_j) - // PineAPPL divides xfx(x_j) by x_j, so store K_j*x_j. - // vals is indexed over the full joint grid (nx_full). + // NOTE: APFEL++: K_j satisfy (C⊗f)(x)=Σ + // K_j*f(x_j) PineAPPL divides xfx(x_j) by x_j, + // so store K_j*x_j. vals is indexed over the + // full joint grid (nx_full). for (std::size_t iix = 0; iix < nx; iix++) { const std::size_t ix = phys_x_indices[iix]; if (ix < vals.size()) @@ -1365,7 +1378,13 @@ pineappl_grid *build_grid(const GridDef &grid_def_in, } } - set_subgrid(grid, ibin, iord, ich, node_vals_bin, subgrid, shape_bin); + set_subgrid(grid, + ibin, + iord, + ich, + node_vals_bin, + subgrid, + shape_bin); } } } diff --git a/pineapfel/src/grid_gen.cpp b/pineapfel/src/grid_gen.cpp index 43ea195..3b9d634 100644 --- a/pineapfel/src/grid_gen.cpp +++ b/pineapfel/src/grid_gen.cpp @@ -286,8 +286,7 @@ GridDef load_grid_def(const std::string &path) { // Normalizations (optional — defaults to 1.0 per bin) if (config["Normalizations"]) def.normalizations = config["Normalizations"].as>(); - else - def.normalizations.assign(def.bins.size(), 1.0); + else def.normalizations.assign(def.bins.size(), 1.0); return def; } diff --git a/pineapfel_pyapi/src/_pineapfel.cpp b/pineapfel_pyapi/src/_pineapfel.cpp index e5c9cb9..c471ec1 100644 --- a/pineapfel_pyapi/src/_pineapfel.cpp +++ b/pineapfel_pyapi/src/_pineapfel.cpp @@ -162,7 +162,8 @@ PYBIND11_MODULE(_pineapfel, m) { .def_readwrite("tabulation", &pineapfel::OperatorCard::tabulation) .def_readwrite("xi", &pineapfel::OperatorCard::xi) .def_readwrite("sidis_mode", &pineapfel::OperatorCard::sidis_mode) - .def_readwrite("sidis_int_eps", &pineapfel::OperatorCard::sidis_int_eps); + .def_readwrite("sidis_int_eps", + &pineapfel::OperatorCard::sidis_int_eps); py::class_(m, "GridDef") .def(py::init<>()) diff --git a/tests/sidis_helper.cpp b/tests/sidis_helper.cpp index dec79a7..aba575f 100644 --- a/tests/sidis_helper.cpp +++ b/tests/sidis_helper.cpp @@ -162,18 +162,17 @@ std::vector compute_sidis_reference(const apfel::Grid &g, const apfel::LagrangeInterpolator li{g}; - double q2_max_val = - *std::max_element(bin_q2.begin(), bin_q2.end()); - int nf_max = apfel::NF(std::sqrt(q2_max_val), thresholds); + double q2_max_val = *std::max_element(bin_q2.begin(), bin_q2.end()); + int nf_max = apfel::NF(std::sqrt(q2_max_val), thresholds); - auto channels = pineapfel::derive_channels(pineapfel::ProcessType::SIDIS, + auto channels = pineapfel::derive_channels(pineapfel::ProcessType::SIDIS, pineapfel::Observable::F2, pineapfel::Current::NC, pineapfel::CCSign::Plus, nf_max); // Key lookup helpers for FT operators (LO/NLO only here) - auto ft_key = [](int alpha_s, int channel_type) -> std::string { + auto ft_key = [](int alpha_s, int channel_type) -> std::string { if (alpha_s == 0 && channel_type == 0) return "DoubleIdentity"; if (alpha_s == 1 && channel_type == 0) return "C1TQ2Q"; if (alpha_s == 1 && channel_type == 1) return "C1TQ2G"; @@ -255,18 +254,17 @@ std::vector compute_sidis_pol_reference(const apfel::Grid &g, const apfel::LagrangeInterpolator li{g}; - double q2_max_val = - *std::max_element(bin_q2.begin(), bin_q2.end()); - int nf_max = apfel::NF(std::sqrt(q2_max_val), thresholds); + double q2_max_val = *std::max_element(bin_q2.begin(), bin_q2.end()); + int nf_max = apfel::NF(std::sqrt(q2_max_val), thresholds); - auto channels = pineapfel::derive_channels(pineapfel::ProcessType::SIDIS, + auto channels = pineapfel::derive_channels(pineapfel::ProcessType::SIDIS, pineapfel::Observable::F2, pineapfel::Current::NC, pineapfel::CCSign::Plus, nf_max); // Key lookup helpers for G1 operators (LO/NLO only here) - auto g1_key = [](int alpha_s, int channel_type) -> std::string { + auto g1_key = [](int alpha_s, int channel_type) -> std::string { if (alpha_s == 0 && channel_type == 0) return "DoubleIdentity"; if (alpha_s == 1 && channel_type == 0) return "DC1Q2Q"; if (alpha_s == 1 && channel_type == 1) return "DC1Q2G"; diff --git a/tests/sidis_helper.h b/tests/sidis_helper.h index 9a46351..c042052 100644 --- a/tests/sidis_helper.h +++ b/tests/sidis_helper.h @@ -13,7 +13,7 @@ std::vector compute_sidis_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, - const std::vector &bin_q2, // one Q² per bin (pointwise) + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, // {lo, hi} per bin const std::vector> &bin_z_bounds, // {lo, hi} per bin const std::vector &charges_dummy, // unused placeholder @@ -29,7 +29,7 @@ std::vector compute_sidis_reference(const apfel::Grid &g, std::vector compute_sidis_pol_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, - const std::vector &bin_q2, // one Q² per bin (pointwise) + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, // {lo, hi} per bin const std::vector> &bin_z_bounds, // {lo, hi} per bin int max_alpha_s, @@ -45,7 +45,7 @@ std::vector compute_sidis_nnlo_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, const std::vector &thresholds, int nf_max, - const std::vector &bin_q2, // one Q² per bin (pointwise) + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, std::function toy_f, @@ -55,9 +55,9 @@ std::vector compute_sidis_nnlo_reference(const apfel::Grid &g, // convention as pineapfel fill.cpp / compute_sidis_nnlo_reference. std::vector compute_sidis_g1_nnlo_reference(const apfel::Grid &g, const apfel::SidisNNLOObjects &sobj, - const std::vector &thresholds, - int nf_max, - const std::vector &bin_q2, // one Q² per bin (pointwise) + const std::vector &thresholds, + int nf_max, + const std::vector &bin_q2, // one Q² per bin (pointwise) const std::vector> &bin_x_bounds, const std::vector> &bin_z_bounds, std::function toy_f, diff --git a/tests/test_f2nc_lhapdf.cpp b/tests/test_f2nc_lhapdf.cpp index 67681ac..e4ff6ac 100644 --- a/tests/test_f2nc_lhapdf.cpp +++ b/tests/test_f2nc_lhapdf.cpp @@ -126,8 +126,8 @@ int main() { double x_c = std::sqrt(x_lo * x_hi); // Pointwise: evaluate BSF at the geometric Q² bin centre. - double q2_c = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + double q2_c = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); double ref = F2.at(0).Evaluate(x_c, std::sqrt(q2_c)); double pred = predictions[ibin]; diff --git a/tests/test_grid_vs_apfelxx.cpp b/tests/test_grid_vs_apfelxx.cpp index 97baa82..73bbb88 100644 --- a/tests/test_grid_vs_apfelxx.cpp +++ b/tests/test_grid_vs_apfelxx.cpp @@ -86,7 +86,6 @@ static apfel::Operator build_channel_operator( return e_q_sq * CNS + (sum_ch / 6.0) * (CS - CNS); } - // ================================================================ int main() { std::cout << "=== Grid vs APFEL++ comparison test ===" << std::endl; @@ -170,8 +169,8 @@ int main() { double ref = 0; // Pointwise: evaluate at the geometric Q² bin centre - double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + double q2 = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); double Q = std::sqrt(q2); int nf = apfel::NF(Q, theory.quark_thresholds); auto charges = apfel::ElectroWeakCharges(Q, timelike); @@ -248,8 +247,8 @@ int main() { double ref = 0; // Pointwise: evaluate at the geometric Q² bin centre - double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + double q2 = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); double Q = std::sqrt(q2); int nf = apfel::NF(Q, theory.quark_thresholds); auto charges = apfel::ElectroWeakCharges(Q, timelike); @@ -352,8 +351,8 @@ int main() { double ref = 0; // Pointwise: evaluate at the geometric Q² bin centre - double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + double q2 = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); double Q = std::sqrt(q2); int nf = apfel::NF(Q, theory.quark_thresholds); auto charges = apfel::ElectroWeakCharges(Q, timelike); @@ -478,8 +477,8 @@ int main() { grid_def.bins[ibin].upper.back()); // Pointwise: evaluate BSF at the geometric Q² bin centre - double q2 = std::sqrt(grid_def.bins[ibin].lower[0] * - grid_def.bins[ibin].upper[0]); + double q2 = std::sqrt( + grid_def.bins[ibin].lower[0] * grid_def.bins[ibin].upper[0]); double Q = std::sqrt(q2); double ref = F2_total.Evaluate(x_center, Q); @@ -582,10 +581,10 @@ int main() { cc_grid_def.bins[ibin].upper.back()); // Pointwise: evaluate BSF at the geometric Q² bin centre - double q2 = std::sqrt(cc_grid_def.bins[ibin].lower[0] * - cc_grid_def.bins[ibin].upper[0]); - double Q = std::sqrt(q2); - double ref = F2CC_total.Evaluate(x_center, Q); + double q2 = std::sqrt(cc_grid_def.bins[ibin].lower[0] * + cc_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2CC_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 @@ -645,7 +644,7 @@ int main() { pineappl_grid *sidis_grid = pineapfel::build_grid(sidis_grid_def, sidis_theory, sidis_op_card); - std::size_t sidis_nbins = pineappl_grid_bin_count(sidis_grid); + std::size_t sidis_nbins = pineappl_grid_bin_count(sidis_grid); // Pointwise: one Q² per bin std::vector sidis_bin_q2; @@ -785,10 +784,10 @@ int main() { pol_grid_def.bins[ibin].upper.back()); // Pointwise: evaluate BSF at the geometric Q² bin centre - double q2 = std::sqrt(pol_grid_def.bins[ibin].lower[0] * - pol_grid_def.bins[ibin].upper[0]); - double Q = std::sqrt(q2); - double ref = g1_total.Evaluate(x_center, Q); + double q2 = std::sqrt(pol_grid_def.bins[ibin].lower[0] * + pol_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = g1_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 @@ -960,10 +959,10 @@ int main() { ffn_grid_def.bins[ibin].upper.back()); // Pointwise: evaluate BSF at the geometric Q² bin centre - double q2 = std::sqrt(ffn_grid_def.bins[ibin].lower[0] * - ffn_grid_def.bins[ibin].upper[0]); - double Q = std::sqrt(q2); - double ref = F2FFN_total.Evaluate(x_center, Q); + double q2 = std::sqrt(ffn_grid_def.bins[ibin].lower[0] * + ffn_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2FFN_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 ? std::abs(pineappl_ffn[ibin] - ref) / std::abs(ref) @@ -1046,10 +1045,10 @@ int main() { mz_grid_def.bins[ibin].upper.back()); // Pointwise: evaluate BSF at the geometric Q² bin centre - double q2 = std::sqrt(mz_grid_def.bins[ibin].lower[0] * - mz_grid_def.bins[ibin].upper[0]); - double Q = std::sqrt(q2); - double ref = F2MZ_total.Evaluate(x_center, Q); + double q2 = std::sqrt(mz_grid_def.bins[ibin].lower[0] * + mz_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); + double ref = F2MZ_total.Evaluate(x_center, Q); double rel_diff = std::abs(ref) > 1e-30 ? std::abs(pineappl_mz[ibin] - ref) / std::abs(ref) @@ -1150,11 +1149,11 @@ int main() { fonll_grid_def.bins[ibin].upper.back()); // Pointwise: evaluate BSF at the geometric Q² bin centre - double q2 = std::sqrt(fonll_grid_def.bins[ibin].lower[0] * - fonll_grid_def.bins[ibin].upper[0]); - double Q = std::sqrt(q2); + double q2 = std::sqrt(fonll_grid_def.bins[ibin].lower[0] * + fonll_grid_def.bins[ibin].upper[0]); + double Q = std::sqrt(q2); // F_FONLL = F_ZM + F_FFN - F_MZ - double ref = F2ZM_tot.Evaluate(x_center, Q) + + double ref = F2ZM_tot.Evaluate(x_center, Q) + F2FFN_tot.Evaluate(x_center, Q) - F2MZ_tot.Evaluate(x_center, Q); double rel_diff = @@ -1204,8 +1203,8 @@ int main() { // keeps all SIDIS tests on the same 20+10-node grid. pineappl_grid *nnlo_grid = pineapfel::build_grid(nnlo_grid_def, nnlo_theory, sidis_op_card); - std::size_t nnlo_nbins = pineappl_grid_bin_count(nnlo_grid); - std::size_t nnlo_nords = pineappl_grid_order_count(nnlo_grid); + std::size_t nnlo_nbins = pineappl_grid_bin_count(nnlo_grid); + std::size_t nnlo_nords = pineappl_grid_order_count(nnlo_grid); // Pointwise: one Q² per bin std::vector nnlo_q2_nodes; From 80917bb21eb0d274d3bd8eb4b2c269e9a11c5689 Mon Sep 17 00:00:00 2001 From: Radonirinaunimi Date: Fri, 3 Apr 2026 23:45:08 +0200 Subject: [PATCH 7/8] Update Python CLI tests references --- .../tests/references/cli_dis_cc_f2.npy | Bin 168 -> 168 bytes .../tests/references/cli_dis_nc_f2.npy | Bin 168 -> 168 bytes .../tests/references/cli_dis_pol_g1.npy | Bin 168 -> 168 bytes .../tests/references/cli_fk_dis_nc_f2.npy | Bin 168 -> 168 bytes 4 files changed, 0 insertions(+), 0 deletions(-) diff --git a/pineapfel_pyapi/tests/references/cli_dis_cc_f2.npy b/pineapfel_pyapi/tests/references/cli_dis_cc_f2.npy index 825604b0df8239e972d83c5cb10c7f700817fe2d..02b949b1d5683ffa50758d6393946682c062ee53 100644 GIT binary patch delta 47 zcmV+~0MP%a0jL3xfGBBe0_O*3P(M-*Gk|k{TR)_b$?pm9mOmoO{~@xa$3L!8Wd9Z6 F%RkG^6>0zg delta 47 zcmV+~0MP%a0jL3xfGDHx(;}#Fl0V(7jXqlGo*+0nHOnp6?3P4j#zatl1 F4?r(x70Cbq diff --git a/pineapfel_pyapi/tests/references/cli_dis_nc_f2.npy b/pineapfel_pyapi/tests/references/cli_dis_nc_f2.npy index 1fbb99f24d44a284944547d7b53b7d535f4fef8b..0bdfc92dcceeae58bd3f89ab9d1793802b0685a1 100644 GIT binary patch delta 47 zcmZ3%xPoy)gNCvtXR(qni0_=USzc12}?XlnKBNSSGai2ZgR^3v5-jntK DmfsPb delta 47 zcmV+~0MP%a0jL3xfG8yp;`Nehi9aiT3iS(&l|TI&o!M^r(LX=?tnwd*0zgI5kpNWI F1wi_G6-WR8 diff --git a/pineapfel_pyapi/tests/references/cli_dis_pol_g1.npy b/pineapfel_pyapi/tests/references/cli_dis_pol_g1.npy index 516398cd53812dda1bb398ddd823b0075aa59029..1205d232ede1424c6209b4e16076e3286df01f52 100644 GIT binary patch delta 47 zcmV+~0MP%a0jL3xfGE5)cA^=>_`U|?5hI(hN59pCBCMGcgFkb@U) F1wi0)6m_2WA; Date: Wed, 8 Apr 2026 14:47:49 +0200 Subject: [PATCH 8/8] Update `mpl` file --- benchmarks/mplstyle.mplstyle | 101 +++++++++++++++++++++-------------- 1 file changed, 62 insertions(+), 39 deletions(-) diff --git a/benchmarks/mplstyle.mplstyle b/benchmarks/mplstyle.mplstyle index 68cd406..84702a8 100644 --- a/benchmarks/mplstyle.mplstyle +++ b/benchmarks/mplstyle.mplstyle @@ -1,39 +1,62 @@ -axes.grid : True -axes.axisbelow : True -axes.titlesize : medium -axes.prop_cycle : cycler(color=['1f77b4', 'd62728', '9467bd', '2ca02c', 'ff7f0e', '8c564b', 'e377c2', '7f7f7f', 'bcbd22', '17becf']) -axes.labelsize : small -axes.formatter.use_mathtext : True -axes.spines.top : True -axes.spines.right : True - -errorbar.capsize : 2 - -font.size : 12 - -figure.figsize : 5.6, 3.9 - -text.usetex : True -text.latex.preamble : \usepackage{amsmath} - -grid.color : cccccc -grid.linestyle : - -grid.linewidth : 0.05 - -legend.fontsize : small -legend.numpoints : 1 -legend.scatterpoints : 1 -legend.loc : best -legend.fancybox : True -legend.framealpha : 0.8 - -lines.markersize : 4 - -mathtext.default : regular - -xtick.labelsize : small -ytick.labelsize : small -xtick.top : False -ytick.right : False - -svg.fonttype : none +axes.grid : False +axes.axisbelow : True +axes.titlesize : 16 +axes.labelsize : 16 +axes.labelpad : 6 +axes.formatter.use_mathtext : True +axes.linewidth : 1.2 +axes.spines.top : True +axes.spines.right : True +axes.prop_cycle : cycler(color=['000000', '1f77b4', 'd62728', '2ca02c', '9467bd', 'ff7f0e', '8c564b', 'e377c2', '7f7f7f', 'bcbd22', '17becf']) + +errorbar.capsize : 2 + +font.family : sans-serif +font.sans-serif : DejaVu Sans, Helvetica, Arial, Liberation Sans +font.size : 14 + +figure.figsize : 6.0, 4.5 + +# text.usetex : True +# text.latex.preamble : \usepackage{amsmath} +text.usetex : False + +grid.color : cccccc +grid.linestyle : - +grid.linewidth : 0.2 + +legend.fontsize : 14 +legend.numpoints : 1 +legend.scatterpoints : 1 +legend.loc : best +legend.frameon : False +legend.fancybox : False +legend.borderaxespad : 0.6 +legend.handlelength : 2.0 + +lines.linewidth : 2.0 +lines.markersize : 4.0 + +mathtext.default : regular +mathtext.fontset : dejavusans + +xtick.direction : in +ytick.direction : in +xtick.top : True +ytick.right : True +xtick.minor.visible : True +ytick.minor.visible : True +xtick.major.size : 7 +ytick.major.size : 7 +xtick.minor.size : 4 +ytick.minor.size : 4 +xtick.major.width : 1.2 +ytick.major.width : 1.2 +xtick.minor.width : 1.0 +ytick.minor.width : 1.0 +xtick.labelsize : 12 +ytick.labelsize : 12 + +savefig.bbox : tight +savefig.pad_inches : 0.03 +svg.fonttype : none