Fix ACE law 4/44/61 sampling for discrete lines and for missing continuum

[GuySten]

Fix ACE law 4/44/61 sampling of discrete lines and of missing continuum

Two defects in the discrete-plus-continuum sampling shared by
ContinuousTabular, CorrelatedAngleEnergy and KalbachMann. Both sit in the
same block of each function, so they are fixed together.

Reference for the format and for MCNP's handling: LA-UR-12-26785,
Modification of ENDF Law 4/44/61 Sampling in MCNP6.

1. Discrete line energies were not interpolated

All three samplers picked one incident energy grid l and returned that grid's
tabulated line energy verbatim:

if (k < n_discrete) {
  return distribution_[l].e_out[k];

Correct for a line at a fixed energy, since both grids hold the same value.
Wrong for a line whose energy varies with incident energy: a single physical
line is emitted as two peaks for a monoenergetic beam, and as a comb
quantised onto the incident energy grid for a continuous one.

ENDF-102 §0.6.2.1 gives the rule — the method of corresponding points
interpolates along the line joining the kth points in the two distributions —
and names the failure mode of the alternative directly: artificial peaks appear
when the distributions shift as a function of energy. OpenMC already applies
unit-base interpolation to the continuum portion of the same table, so this
also removes an internal inconsistency.

Affected data. Scanning ENDF/B-VIII.1 (557 files) for tables where
e_out[k] differs between adjacent incident energies, above a 1e-5 relative
tolerance that excludes six-significant-digit round-off in the ENDF fixed-format
fields, gives 13 tables:

nuclide	MT	lines	dE_out/dE_in	A/(A+1)
Pt-190…198 (10 isotopes)	102	995–2276	0.9947–0.9949	0.9947–0.9950
Ta-180m	102	1592	0.9944	0.9945
H-1	102	1	0.4998	0.4998
He-3	102	1	0.7493	0.7493
Os-187	22	73	−5.5046	n/a

Every MT=102 entry matches the ENDF primary-photon relation
E_gamma = E_g + A/(A+1) * E_n to four decimals, identifying these as primary
capture gammas — lines whose energy is meant to move with incident energy.
ENDF/B-VII.1 had two such tables; VIII.1 has thirteen, several with over a
thousand lines.

The effect scales with incident energy, so it is negligible for thermal capture
and largest in the fast range where these capture cross sections are small. This
is a spectrum-shape fix for fast-capture and gamma-spectroscopy work; it will not
move k-eff.

2. Out-of-bounds read when a table has no continuum

const double E_i_1 = distribution_[i].e_out[n_discrete];

e_out[n_discrete] is the first continuum outgoing energy, so this reads past
the end when a table is entirely discrete lines. The ACE format requires the
discrete line count to match at every incident energy but does not require a
continuum at each one, and MCNP made the same assumption:

MCNP assumed that if continua are present at some incident energies, then
there must be continua present at all other incident energies. […] this need
not be true.

The read is not benign: E_i_1 enters the unit-base transform with weight
1 - r, so whenever the continuum is sampled from the other grid the outgoing
energy is contaminated by whatever lies past the array. A negative emission
energy is how LANL found this in MCNP, in Ac-226 (n,4n); the common case is a
silently wrong energy.

ContinuousTabular appeared guarded by if (!histogram_interp && n_energy_out > 1), but at that point n_energy_out still holds grid l's size and is only
reassigned to grid i's on the next line — so the guard tests the wrong table.

Fix. Each side is checked for a continuum and the interpolation factor
collapses onto whichever side has one, following MCNP6:

const double r_c = cont_i ? (cont_i1 ? r : 0.0) : 1.0;

MCNP additionally uses a 0–1 eV convention for photons; AngleEnergy::sample has
no particle type available, so the neutron treatment is applied uniformly. Happy
to thread the particle type through if reviewers prefer matching MCNP exactly.

Affected data. 4557 ENDF/B-VIII.1 tables are entirely discrete lines; those
never produce a wrong value (the discrete branch returns early) but the invalid
read happens on every call, so ASan fires on any coupled neutron-photon run with
discrete level data. The tables that produce wrong numbers — a continuum at
some incident energies but not others — reproduced Appendix D of LA-UR-12-26785
exactly for ENDF/B-VII.1. VIII.1 adds isotopes that are not minor actinides:
Cs-134, Dy-159, Os-189, Os-190, Pm-147, Ta-180m, U-233. Cs-134, Pm-147 and
Dy-159 are fission products present in any burnt fuel.

Behaviour

The discrete-line fix is bitwise identical wherever a line does not move:
e_lo[k] + r * (e_hi[k] - e_lo[k]) with e_hi[k] == e_lo[k] gives exactly
e_lo[k] for any finite r. Verified across 33 values of r. Histogram
interpolation on the incident grid passes 0.0 and likewise returns the lower
grid verbatim. Of 4591 tables with discrete lines in ENDF/B-VIII.1, only 13 can
produce a different number.

The continuum fix changes nothing where both incident energies have a continuum:
r_c == r and the bounds are identical. ContinuousTabular additionally swaps
n_energy_out > 1 for a test of the actual denominator, E_hi != E_lo, which is
what it stood in for.

Neither change alters the prn(seed) call count on any path, so the random
number stream is untouched.

Follow-up

These three functions share 89 of ~100 lines between the two secondary
distributions, and 61 across all three. Every defect here was present in some
copies and not others. A separate PR consolidating the shared CDF search,
inversion and unit-base transform is prepared and will follow.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)