From dfd34342f9d5f12e8655c299dc2d14b099f00609 Mon Sep 17 00:00:00 2001 From: DavidA Date: Fri, 21 Aug 2026 05:08:19 +0100 Subject: [PATCH 1/2] Return discrete lines exactly in lin-lin tabular energy sampling MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ContinuousTabular, KalbachMann, and CorrelatedAngleEnergy all guard the histogram interpolation branch with k >= n_discrete, returning a sampled discrete line's energy exactly — but none of them guard the lin-lin branch. A discrete hit leaves the walk with r1 < c_k, so the lin-lin inversion runs with a negative (r1 - c_k) and shifts the sampled energy below the line (bounded by the max(0, ...) clamp at E_l_k - p_l_k/frac). In ENDF/B-VIII.0 the discrete-lines-with-lin-lin combination appears in 19,913 incident-energy tables — all secondary-photon distributions (discrete gammas from (n,2n), (n,3n), (n,n') continuum and similar), so coupled neutron-photon transport samples smeared gamma lines. No neutron-product table in the library uses the combination, so neutron-only results are unaffected. The fix mirrors the existing histogram guard in all three samplers; for Kalbach-Mann the r/a parameters use the line's tabulated values, which is what the interpolation reduces to at E_out = E_l_k. Co-Authored-By: Claude Fable 5 --- src/distribution_energy.cpp | 34 ++++++++++++++++++++-------------- src/secondary_correlated.cpp | 6 ++++++ src/secondary_kalbach.cpp | 8 ++++++++ 3 files changed, 34 insertions(+), 14 deletions(-) diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 2f8e6cf1a99..6ab6f99c3b0 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -220,20 +220,26 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const } } else if (distribution_[l].interpolation == Interpolation::lin_lin) { - // Linear-linear interpolation - double E_l_k1 = distribution_[l].e_out[k + 1]; - double p_l_k1 = distribution_[l].p[k + 1]; - - if (E_l_k != E_l_k1) { - double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); - if (frac == 0.0) { - E_out = E_l_k + (r1 - c_k) / p_l_k; - } else { - E_out = - E_l_k + - (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - - p_l_k) / - frac; + // Linear-linear interpolation. The inversion applies only to the + // continuous portion of the distribution; a sampled discrete line + // (k < n_discrete) is returned exactly, as in the histogram branch — + // otherwise (r1 - c_k) is negative and the sampled energy is shifted + // off the line. + if (k >= n_discrete) { + double E_l_k1 = distribution_[l].e_out[k + 1]; + double p_l_k1 = distribution_[l].p[k + 1]; + + if (E_l_k != E_l_k1) { + double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = + E_l_k + + (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - + p_l_k) / + frac; + } } } } else { diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 42ef67ced42..33ed8590a55 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -219,6 +219,12 @@ Distribution& CorrelatedAngleEnergy::sample_dist( E_out = E_l_k; } + } else if (k < n_discrete) { + // A sampled discrete line is returned exactly (the lin-lin inversion + // below applies only to the continuous portion; with (r1 - c_k) + // negative it would shift the sampled energy off the line) + E_out = E_l_k; + } else if (distribution_[l].interpolation == Interpolation::lin_lin) { // Linear-linear interpolation double E_l_k1 = distribution_[l].e_out[k + 1]; diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 018ce1c8a9c..5ac15483dcd 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -182,6 +182,14 @@ void KalbachMann::sample_params( km_r = distribution_[l].r[k]; km_a = distribution_[l].a[k]; + } else if (k < n_discrete) { + // A sampled discrete line is returned exactly (the lin-lin inversion + // below applies only to the continuous portion; with (r1 - c_k) + // negative it would shift the sampled energy off the line) + E_out = E_l_k; + km_r = distribution_[l].r[k]; + km_a = distribution_[l].a[k]; + } else { // Linear-linear interpolation double E_l_k1 = distribution_[l].e_out[k + 1]; From 05fcd20d0c20e78bd1b23e1afb7569a2232a59cb Mon Sep 17 00:00:00 2001 From: GuySten Date: Sat, 22 Aug 2026 20:35:11 +0300 Subject: [PATCH 2/2] simplify code a bit --- src/distribution_energy.cpp | 83 +++++++++++++++++++----------------- src/secondary_correlated.cpp | 59 ++++++++++++------------- src/secondary_kalbach.cpp | 77 +++++++++++++++++---------------- 3 files changed, 110 insertions(+), 109 deletions(-) diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 6ab6f99c3b0..7712c0d763f 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -211,21 +211,23 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const } double E_l_k = distribution_[l].e_out[k]; - double p_l_k = distribution_[l].p[k]; - double E_out = E_l_k; - if (distribution_[l].interpolation == Interpolation::histogram) { - // Histogram interpolation - if (p_l_k > 0.0 && k >= n_discrete) { - E_out = E_l_k + (r1 - c_k) / p_l_k; - } - } else if (distribution_[l].interpolation == Interpolation::lin_lin) { - // Linear-linear interpolation. The inversion applies only to the - // continuous portion of the distribution; a sampled discrete line - // (k < n_discrete) is returned exactly, as in the histogram branch — - // otherwise (r1 - c_k) is negative and the sampled energy is shifted - // off the line. - if (k >= n_discrete) { + if (k < n_discrete) { + // Discrete case + return E_l_k; + } else { + // Continuous case + double p_l_k = distribution_[l].p[k]; + double E_out; + if (distribution_[l].interpolation == Interpolation::histogram) { + // Histogram interpolation + if (p_l_k > 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = E_l_k; + } + } else if (distribution_[l].interpolation == Interpolation::lin_lin) { + // Linear-linear interpolation double E_l_k1 = distribution_[l].e_out[k + 1]; double p_l_k1 = distribution_[l].p[k + 1]; @@ -240,36 +242,39 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const p_l_k) / frac; } + } else { + E_out = E_l_k; } + } else { + throw std::runtime_error { + "Unexpected interpolation for continuous energy " + "distribution."}; } - } else { - throw std::runtime_error {"Unexpected interpolation for continuous energy " - "distribution."}; - } - // Now interpolate between incident energy bins i and i + 1 - if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) { - // Interpolation for energy E1 and EK - n_energy_out = distribution_[i].e_out.size(); - n_discrete = distribution_[i].n_discrete; - const double E_i_1 = distribution_[i].e_out[n_discrete]; - const double E_i_K = distribution_[i].e_out[n_energy_out - 1]; - - n_energy_out = distribution_[i + 1].e_out.size(); - n_discrete = distribution_[i + 1].n_discrete; - const double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; - const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; - - const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); - const double E_K = E_i_K + r * (E_i1_K - E_i_K); - - if (l == i) { - return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); + // Now interpolate between incident energy bins i and i + 1 + if (!histogram_interp && n_energy_out > 1) { + // Interpolation for energy E1 and EK + n_energy_out = distribution_[i].e_out.size(); + n_discrete = distribution_[i].n_discrete; + const double E_i_1 = distribution_[i].e_out[n_discrete]; + const double E_i_K = distribution_[i].e_out[n_energy_out - 1]; + + n_energy_out = distribution_[i + 1].e_out.size(); + n_discrete = distribution_[i + 1].n_discrete; + const double E_i1_1 = distribution_[i + 1].e_out[n_discrete]; + const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1]; + + const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1); + const double E_K = E_i_K + r * (E_i1_K - E_i_K); + + if (l == i) { + return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); + } else { + return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); + } } else { - return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1); + return E_out; } - } else { - return E_out; } } diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 33ed8590a55..642ecd82b78 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -210,40 +210,37 @@ Distribution& CorrelatedAngleEnergy::sample_dist( } double E_l_k = distribution_[l].e_out[k]; - double p_l_k = distribution_[l].p[k]; - if (distribution_[l].interpolation == Interpolation::histogram) { - // Histogram interpolation - if (p_l_k > 0.0 && k >= n_discrete) { - E_out = E_l_k + (r1 - c_k) / p_l_k; - } else { - E_out = E_l_k; - } - - } else if (k < n_discrete) { - // A sampled discrete line is returned exactly (the lin-lin inversion - // below applies only to the continuous portion; with (r1 - c_k) - // negative it would shift the sampled energy off the line) + if (k < n_discrete) { + // Discrete case E_out = E_l_k; - - } else if (distribution_[l].interpolation == Interpolation::lin_lin) { - // Linear-linear interpolation - double E_l_k1 = distribution_[l].e_out[k + 1]; - double p_l_k1 = distribution_[l].p[k + 1]; - - double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); - if (frac == 0.0) { - E_out = E_l_k + (r1 - c_k) / p_l_k; - } else { - E_out = - E_l_k + - (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - - p_l_k) / - frac; + } else { + // Continuous case + double p_l_k = distribution_[l].p[k]; + if (distribution_[l].interpolation == Interpolation::histogram) { + // Histogram interpolation + if (p_l_k > 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = E_l_k; + } + } else if (distribution_[l].interpolation == Interpolation::lin_lin) { + // Linear-linear interpolation + double E_l_k1 = distribution_[l].e_out[k + 1]; + double p_l_k1 = distribution_[l].p[k + 1]; + + double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = + E_l_k + + (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - + p_l_k) / + frac; + } } - } - // Now interpolate between incident energy bins i and i + 1 - if (k >= n_discrete) { + // Now interpolate between incident energy bins i and i + 1 if (l == i) { E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); } else { diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 5ac15483dcd..a5d9b9ad0b4 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -169,54 +169,53 @@ void KalbachMann::sample_params( } double E_l_k = distribution_[l].e_out[k]; - double p_l_k = distribution_[l].p[k]; - if (distribution_[l].interpolation == Interpolation::histogram) { - // Histogram interpolation - if (p_l_k > 0.0 && k >= n_discrete) { - E_out = E_l_k + (r1 - c_k) / p_l_k; - } else { - E_out = E_l_k; - } - - // Determine Kalbach-Mann parameters - km_r = distribution_[l].r[k]; - km_a = distribution_[l].a[k]; - - } else if (k < n_discrete) { - // A sampled discrete line is returned exactly (the lin-lin inversion - // below applies only to the continuous portion; with (r1 - c_k) - // negative it would shift the sampled energy off the line) + if (k < n_discrete) { + // Discrete case E_out = E_l_k; km_r = distribution_[l].r[k]; km_a = distribution_[l].a[k]; } else { - // Linear-linear interpolation - double E_l_k1 = distribution_[l].e_out[k + 1]; - double p_l_k1 = distribution_[l].p[k + 1]; + // Continuous case + double p_l_k = distribution_[l].p[k]; + if (distribution_[l].interpolation == Interpolation::histogram) { + // Histogram interpolation + if (p_l_k > 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = E_l_k; + } + + // Determine Kalbach-Mann parameters + km_r = distribution_[l].r[k]; + km_a = distribution_[l].a[k]; - double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); - if (frac == 0.0) { - E_out = E_l_k + (r1 - c_k) / p_l_k; } else { - E_out = - E_l_k + - (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - - p_l_k) / - frac; - } + // Linear-linear interpolation + double E_l_k1 = distribution_[l].e_out[k + 1]; + double p_l_k1 = distribution_[l].p[k + 1]; + + double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k) / p_l_k; + } else { + E_out = + E_l_k + + (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) - + p_l_k) / + frac; + } - // Determine Kalbach-Mann parameters - km_r = distribution_[l].r[k] + - (E_out - E_l_k) / (E_l_k1 - E_l_k) * - (distribution_[l].r[k + 1] - distribution_[l].r[k]); - km_a = distribution_[l].a[k] + - (E_out - E_l_k) / (E_l_k1 - E_l_k) * - (distribution_[l].a[k + 1] - distribution_[l].a[k]); - } + // Determine Kalbach-Mann parameters + km_r = distribution_[l].r[k] + + (E_out - E_l_k) / (E_l_k1 - E_l_k) * + (distribution_[l].r[k + 1] - distribution_[l].r[k]); + km_a = distribution_[l].a[k] + + (E_out - E_l_k) / (E_l_k1 - E_l_k) * + (distribution_[l].a[k + 1] - distribution_[l].a[k]); + } - // Now interpolate between incident energy bins i and i + 1 - if (k >= n_discrete) { + // Now interpolate between incident energy bins i and i + 1 if (l == i) { E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1); } else {