Add 2D atmospheric flux tables and tabulated trapezoid integration by austinschneider · Pull Request #105 · Harvard-Neutrino/SIREN

austinschneider · 2026-04-28T04:29:00Z

Stack: PR 7 of 14

This PR is part of a 14-PR stack decomposing dev/HNL_DIS into reviewable chunks.

Base: feat/primary-external-distribution
Head: feat/2d-flux-tables

Merge order: feat/primary-external-distribution should land before this PR.

Squashed contents

Squashed diff for paths:

projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx
projects/distributions/private/primary/energy_direction/PrimaryEnergyDirectionDistribution.cxx
projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx
projects/distributions/public/SIREN/distributions/primary/energy/TabulatedFluxDistribution.h
projects/distributions/public/SIREN/distributions/primary/energy_direction/PrimaryEnergyDirectionDistribution.h
projects/distributions/public/SIREN/distributions/primary/energy_direction/Tabulated2DFluxDistribution.h
projects/utilities/public/SIREN/utilities/Integration.h
projects/utilities/public/SIREN/utilities/Interpolator.h
projects/distributions/private/pybindings/distributions.cxx
projects/distributions/CMakeLists.txt

Source commits on dev/HNL_DIS_clean that contributed:
21451da (Nicholas Kamp): first attempt at external list primary injector
455b073 (Nicholas Kamp): introduce primary bounded and physical vertex distributions to sample vertex in the case that the inital position is already provided from the external distribution
521695f (Nicholas Kamp): fix issues with 2D tabular integral calculation, sampling errors with very short decay lengths, DarkNews errors, and kinematic erros in Dipole portal HNL DIS
73bc119 (Nicholas Kamp): sphere build fixes
82b4fdc (Nicholas Kamp): add ability to trapezoid integrate tabulated 1D fluxes
9098b42 (Nicholas Kamp): fix weighting issue in PrimaryExternalDist
9968ff7 (Nicholas Kamp): Start implementing 2D flux table for atmospheric fluxes
9dae77a (Nicholas Kamp): spherical volume position distribution
a79d440 (Nicholas Kamp): getting HNLs to work for SINE and UNDINE
a9aebd3 (Nicholas Kamp): fix log integration
c7cff6c (Nicholas Kamp): fixing build and runtime errors on 2D flux distribution

Notes

This branch is the squashed cumulative diff of the listed source commits. Per-commit history is browsable on dev/HNL_DIS_clean.
Large .fits data files have been removed from the branch and are packaged separately.

Squashed diff for paths: - projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx - projects/distributions/private/primary/energy_direction/PrimaryEnergyDirectionDistribution.cxx - projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx - projects/distributions/public/SIREN/distributions/primary/energy/TabulatedFluxDistribution.h - projects/distributions/public/SIREN/distributions/primary/energy_direction/PrimaryEnergyDirectionDistribution.h - projects/distributions/public/SIREN/distributions/primary/energy_direction/Tabulated2DFluxDistribution.h - projects/utilities/public/SIREN/utilities/Integration.h - projects/utilities/public/SIREN/utilities/Interpolator.h - projects/distributions/private/pybindings/distributions.cxx - projects/distributions/CMakeLists.txt Source commits on dev/HNL_DIS_clean that contributed: 21451da (Nicholas Kamp): first attempt at external list primary injector 455b073 (Nicholas Kamp): introduce primary bounded and physical vertex distributions to sample vertex in the case that the inital position is already provided from the external distribution 521695f (Nicholas Kamp): fix issues with 2D tabular integral calculation, sampling errors with very short decay lengths, DarkNews errors, and kinematic erros in Dipole portal HNL DIS 73bc119 (Nicholas Kamp): sphere build fixes 82b4fdc (Nicholas Kamp): add ability to trapezoid integrate tabulated 1D fluxes 9098b42 (Nicholas Kamp): fix weighting issue in PrimaryExternalDist 9968ff7 (Nicholas Kamp): Start implementing 2D flux table for atmospheric fluxes 9dae77a (Nicholas Kamp): spherical volume position distribution a79d440 (Nicholas Kamp): getting HNLs to work for SINE and UNDINE a9aebd3 (Nicholas Kamp): fix log integration c7cff6c (Nicholas Kamp): fixing build and runtime errors on 2D flux distribution

Copilot

Pull request overview

Note

Copilot was unable to run its full agentic suite in this review.

Adds support for 2D (energy, zenith) tabulated atmospheric flux distributions and introduces tabulated/trapezoid-style integration options for flux normalization, with corresponding C++ and Python binding updates.

Changes:

Added PrimaryEnergyDirectionDistribution and Tabulated2DFluxDistribution (C++ + serialization + pybind exposure).
Added 2D Simpson integration and a simple trapezoid integration helper to utilities.
Extended TabulatedFluxDistribution to optionally normalize via trapezoid integration and updated Python constructors accordingly.

Reviewed changes

Copilot reviewed 10 out of 10 changed files in this pull request and generated 9 comments.

Show a summary per file

File	Description
projects/utilities/public/SIREN/utilities/Interpolator.h	Exposes `IsLogX/IsLogY` so callers can adapt sampling/integration to log grids.
projects/utilities/public/SIREN/utilities/Integration.h	Adds `simpsonIntegrate2D` and a simple `trapezoidIntegrate` helper for tabulated data.
projects/distributions/public/SIREN/distributions/primary/energy_direction/Tabulated2DFluxDistribution.h	Declares new 2D flux distribution API + serialization.
projects/distributions/public/SIREN/distributions/primary/energy_direction/PrimaryEnergyDirectionDistribution.h	Introduces a new distribution base for coupled energy+direction sampling.
projects/distributions/public/SIREN/distributions/primary/energy/TabulatedFluxDistribution.h	Adds `romberg` toggle for integral computation + stores PDF node values.
projects/distributions/private/pybindings/distributions.cxx	Exposes new distributions to Python and updates `TabulatedFluxDistribution` constructors.
projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx	Implements table loading, normalization, and MH sampling for (E, zenith).
projects/distributions/private/primary/energy_direction/PrimaryEnergyDirectionDistribution.cxx	Implements `Sample()` and density variables for the new base class.
projects/distributions/private/primary/energy/TabulatedFluxDistribution.cxx	Implements trapezoid-based integral option and stores PDF nodes.
projects/distributions/CMakeLists.txt	Adds new sources to the build.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Copilot · 2026-04-28T04:37:14Z

+    if (fluxTable.IsLogY()) {
+        double logZenithMin = log10(zenithMin);
+        double logZenithMax = log10(zenithMax);
+        zenith = pow(10,logZenithMin + u2 * (logZenithMax - logZenithMin));
+    }
+    else {
+        zenith = zenithMin + u2 * (zenithMax - zenithMin);
+    }


ComputeIntegral() only applies the log-Jacobian for log-X, but SampleTablePoint() also supports log-Y. If fluxTable.IsLogY() is ever true, integral (and therefore pdf()) will be mis-normalized relative to sampling. Fix by applying the corresponding log-space change of variables in ComputeIntegral() for Y as well (or remove/forbid log-Y sampling if not supported)."

Suggested change

if (fluxTable.IsLogY()) {

double logZenithMin = log10(zenithMin);

double logZenithMax = log10(zenithMax);

zenith = pow(10,logZenithMin + u2 * (logZenithMax - logZenithMin));

}

else {

zenith = zenithMin + u2 * (zenithMax - zenithMin);

}

// Keep Y sampling in linear space so it matches the current integral

// normalization. Log-Y sampling requires a corresponding Jacobian in

// ComputeIntegral(), which is not applied here.

zenith = zenithMin + u2 * (zenithMax - zenithMin);

Copilot · 2026-04-28T04:37:15Z

+    // metropolis hastings variables
+    const size_t burnin = 40; // burnin parameter for MH sampling
+    mutable size_t MH_sampled_points = 0; // to track the number of samples
+    mutable double MH_density;


MH_density is never initialized, but it is used as a divisor in Metropolis-Hastings (odds = test_density / MH_density). This is undefined behavior and can produce NaNs/inf or random acceptance. Initialize MH_density (e.g., to 0.0) in the declaration or constructors, and consider resetting MH_density/MH_sampled_points when bounds change (SetEnergyBounds/SetZenithBounds) so the chain doesn’t continue with stale state after re-normalization.

Suggested change

mutable double MH_density;

mutable double MH_density = 0.0;

Copilot · 2026-04-28T04:37:15Z

+            std::stringstream ss(buf);
+            double x, y, f;
+            ss >> x >> y >> f;


This translation unit uses std::stringstream and assert but does not include <sstream> or <cassert>. This can fail to compile depending on transitive includes. Add the missing standard headers in this file.

Copilot · 2026-04-28T04:37:15Z

+    assert(energies.size()==flux.size());
+    assert(zeniths.size()==flux.size());


This translation unit uses std::stringstream and assert but does not include <sstream> or <cassert>. This can fail to compile depending on transitive includes. Add the missing standard headers in this file.

Copilot · 2026-04-28T04:37:16Z

+double simpsonIntegrate2D(const FuncType& func, double a1, double b1, double a2, double b2,
+                          double tol=1e-3, const unsigned int N1 = 10, const unsigned int N2=10, const unsigned int maxIter=15){


The convergence test divides by estimate, which becomes a division-by-zero when the integral is 0 (or extremely small), yielding inf/NaN and preventing correct convergence. Use an absolute error check (e.g., abs(estimate - prev_estimate) < tol_abs) or a mixed abs/relative criterion that guards estimate == 0. Also, Simpson’s rule requires an even number of subintervals per dimension; currently callers can pass odd N1/N2 and the weight pattern becomes invalid—validate/enforce even N1 and N2 (or document and adjust them internally).

Copilot · 2026-04-28T04:37:16Z

+     if(std::abs((estimate-prev_estimate)/estimate) < tol) {
+      return estimate;
+     }


The convergence test divides by estimate, which becomes a division-by-zero when the integral is 0 (or extremely small), yielding inf/NaN and preventing correct convergence. Use an absolute error check (e.g., abs(estimate - prev_estimate) < tol_abs) or a mixed abs/relative criterion that guards estimate == 0. Also, Simpson’s rule requires an even number of subintervals per dimension; currently callers can pass odd N1/N2 and the weight pattern becomes invalid—validate/enforce even N1 and N2 (or document and adjust them internally).

Copilot · 2026-04-28T04:37:16Z

+inline double trapezoidIntegrate(std::vector<double>& x, std::vector<double>& y){
+   if(x.size()!=y.size())
+       throw(std::runtime_error("Integration x and y vectors must be the same size"));


trapezoidIntegrate does not modify x or y, but takes them by non-const reference, which unnecessarily restricts callers (can’t pass temporaries/const vectors) and obscures intent. Change the signature to take const std::vector<double>& for both. Additionally, consider throwing when x.size() < 2 to avoid silently returning 0 for empty/1-point inputs.

Suggested change

inline double trapezoidIntegrate(std::vector<double>& x, std::vector<double>& y){

if(x.size()!=y.size())

throw(std::runtime_error("Integration x and y vectors must be the same size"));

inline double trapezoidIntegrate(const std::vector<double>& x, const std::vector<double>& y){

if(x.size()!=y.size())

throw(std::runtime_error("Integration x and y vectors must be the same size"));

if(x.size() < 2)

throw(std::runtime_error("Integration requires at least two points"));

Copilot · 2026-04-28T04:37:16Z

+void TabulatedFluxDistribution::ComputeIntegral(bool romberg) {
    std::function<double(double)> integrand = [&] (double x) -> double {
        return unnormed_pdf(x);
    };
-    integral = siren::utilities::rombergIntegrate(integrand, energyMin, energyMax);
+    if (romberg) {
+        integral = siren::utilities::rombergIntegrate(integrand, energyMin, energyMax);
+    } else {
+        integral = siren::utilities::trapezoidIntegrate(energy_nodes, pdf_nodes);
+    }
 }


When romberg == false, the integral ignores energyMin/energyMax entirely and integrates over all energy_nodes. This produces incorrect normalization if bounds are set narrower than the full table. Fix by trapezoid-integrating only over the portion of nodes within [energyMin, energyMax] (and, if needed, interpolating endpoints). Also ensure the selected integration method remains consistent after SetEnergyBounds()—currently SetEnergyBounds() calls ComputeIntegral() without passing the romberg choice (so it will revert to the default). Persist the method as a member or plumb it through all recompute sites.

Copilot · 2026-04-28T04:37:16Z

+    .def(init<std::string, bool, bool>())
+    .def(init<double, double, std::string, bool, bool>())
+    .def(init<std::vector<double>, std::vector<double>, bool, bool>())
+    .def(init<double, double, std::vector<double>, std::vector<double>, bool, bool>())


The Python API for TabulatedFluxDistribution appears to have become strictly more required-args than before (previously (filename, has_physical_normalization) etc.). Even though C++ provides defaults, these init<...>() overloads require callers to pass both booleans explicitly, which is a breaking change for existing Python users. Consider adding overloads matching the previous signatures, or bind with py::arg defaults so existing code continues to work while still exposing the new romberg toggle.

Suggested change

.def(init<std::string, bool, bool>())

.def(init<double, double, std::string, bool, bool>())

.def(init<std::vector<double>, std::vector<double>, bool, bool>())

.def(init<double, double, std::vector<double>, std::vector<double>, bool, bool>())

.def(init<std::string, bool, bool>(),

arg("filename"),

arg("has_physical_normalization") = false,

arg("romberg") = false)

.def(init<double, double, std::string, bool, bool>(),

arg("energy_min"),

arg("energy_max"),

arg("filename"),

arg("has_physical_normalization") = false,

arg("romberg") = false)

.def(init<std::vector<double>, std::vector<double>, bool, bool>(),

arg("energy_nodes"),

arg("flux_nodes"),

arg("has_physical_normalization") = false,

arg("romberg") = false)

.def(init<double, double, std::vector<double>, std::vector<double>, bool, bool>(),

arg("energy_min"),

arg("energy_max"),

arg("energy_nodes"),

arg("flux_nodes"),

arg("has_physical_normalization") = false,

arg("romberg") = false)

austinschneider · 2026-04-28T15:32:00Z

Closing to reopen from the commit author's account so they can be added as a reviewer. Branch unchanged; will be re-linked from the new PR shortly.

Copilot AI review requested due to automatic review settings April 28, 2026 04:29

austinschneider force-pushed the feat/primary-external-distribution branch from cc05a02 to 4a078b6 Compare April 28, 2026 04:34

austinschneider force-pushed the feat/2d-flux-tables branch from b9c0966 to 496977b Compare April 28, 2026 04:34

Copilot AI reviewed Apr 28, 2026

View reviewed changes

austinschneider force-pushed the feat/primary-external-distribution branch from 4a078b6 to 4500aec Compare April 28, 2026 05:09

austinschneider force-pushed the feat/2d-flux-tables branch from 496977b to 967cc0a Compare April 28, 2026 05:09

austinschneider force-pushed the feat/primary-external-distribution branch from 4500aec to e4231d5 Compare April 28, 2026 06:04

austinschneider force-pushed the feat/2d-flux-tables branch from 967cc0a to 9204f13 Compare April 28, 2026 06:04

Copilot started reviewing on behalf of austinschneider April 28, 2026 11:41 View session

austinschneider force-pushed the feat/primary-external-distribution branch from e4231d5 to 05bc3dd Compare April 28, 2026 13:26

austinschneider force-pushed the feat/2d-flux-tables branch 2 times, most recently from b717fce to c04ebb4 Compare April 28, 2026 14:36

austinschneider force-pushed the feat/primary-external-distribution branch from 05bc3dd to c967cc7 Compare April 28, 2026 14:36

austinschneider closed this Apr 28, 2026

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Add 2D atmospheric flux tables and tabulated trapezoid integration#105

Add 2D atmospheric flux tables and tabulated trapezoid integration#105
austinschneider wants to merge 1 commit into
feat/primary-external-distributionfrom
feat/2d-flux-tables

austinschneider commented Apr 28, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

Copilot AI Apr 28, 2026

Uh oh!

austinschneider commented Apr 28, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

		assert(energies.size()==flux.size());
		assert(zeniths.size()==flux.size());

		double simpsonIntegrate2D(const FuncType& func, double a1, double b1, double a2, double b2,
		double tol=1e-3, const unsigned int N1 = 10, const unsigned int N2=10, const unsigned int maxIter=15){

-    .def(init<std::string, bool, bool>())
-    .def(init<double, double, std::string, bool, bool>())
-    .def(init<std::vector<double>, std::vector<double>, bool, bool>())
-    .def(init<double, double, std::vector<double>, std::vector<double>, bool, bool>())
+    .def(init<std::string, bool, bool>(),
+         arg("filename"),
+         arg("has_physical_normalization") = false,
+         arg("romberg") = false)
+    .def(init<double, double, std::string, bool, bool>(),
+         arg("energy_min"),
+         arg("energy_max"),
+         arg("filename"),
+         arg("has_physical_normalization") = false,
+         arg("romberg") = false)
+    .def(init<std::vector<double>, std::vector<double>, bool, bool>(),
+         arg("energy_nodes"),
+         arg("flux_nodes"),
+         arg("has_physical_normalization") = false,
+         arg("romberg") = false)
+    .def(init<double, double, std::vector<double>, std::vector<double>, bool, bool>(),
+         arg("energy_min"),
+         arg("energy_max"),
+         arg("energy_nodes"),
+         arg("flux_nodes"),
+         arg("has_physical_normalization") = false,
+         arg("romberg") = false)

Uh oh!

Conversation

austinschneider commented Apr 28, 2026

Stack: PR 7 of 14

Squashed contents

Notes

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 28, 2026

Choose a reason for hiding this comment

Uh oh!

austinschneider commented Apr 28, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants