nereids_pipeline/

calibration.rs

//! Energy calibration for TOF neutron instruments.
//!
//! Finds the flight path length (L) and TOF delay (t₀) that best align
//! a measured transmission spectrum with the ENDF resonance model.
//!
//! The energy-TOF relationship is:
//!
//!   E = C · (L / (t − t₀))²
//!
//! where C = mₙ / 2 ≈ 5.2276e-9 [eV·s²/m²].
//!
//! When L or t₀ differ from the values assumed during data reduction,
//! resonance positions shift in the energy domain, causing catastrophic
//! chi² degradation (e.g. 436 → 2.7 for a 0.3% L correction on VENUS).

use nereids_core::constants::{EV_TO_JOULES, NEUTRON_MASS_KG};
use nereids_endf::resonance::ResonanceData;
use nereids_physics::transmission::{self, InstrumentParams, SampleParams};

use crate::error::PipelineError;

/// Neutron mass constant: C = m_n / (2 · eV) ≈ 5.2276e-9 eV·s²/m².
///
/// E [eV] = C · (L [m] / t [s])²
///
/// Uses the CODATA 2018 values from `nereids_core::constants` so that
/// this calibration path, `EnergyScaleTransmissionModel`, and
/// `core::tof_to_energy` all agree to machine precision.
const NEUTRON_MASS_CONSTANT: f64 = 0.5 * NEUTRON_MASS_KG / EV_TO_JOULES;

/// Lower / upper bounds (log10) on the `n_total` (areal density,
/// atoms/barn) search interval for `calibrate_energy`.  The search
/// runs in `log10(n)` so the band is sampled with relative — rather
/// than absolute — resolution.
///
/// The internal search band is `[~5e-6, ~2e-2]` atoms/barn (a third
/// of a decade beyond each documented edge on either side).  The
/// boundary-saturation guard (`CALIBRATION_LOG10_BOUNDARY_TOL`,
/// ≈ 5 % in linear density) trims a sliver off each end, leaving
/// the *documented* user-supported interval at exactly `[1e-5,
/// 1e-2]`: the doc-stated edges are inside the tolerance window,
/// not on it.
///
/// `[1e-5, 1e-2]` covers every realistic VENUS / paper-relevant
/// density: thin diluted samples down to ~1e-5 atoms/barn (trace
/// detectability ~ Hf in matrix), the Hf calibration foil at
/// ~1e-4, and 1 mm metal foils (U, W, Ni) up to ~1e-2 atoms/barn.
/// Sample densities at the exact documented edges (`1e-5` or
/// `1e-2`) are accepted because the search band extends ~0.3
/// decades beyond them — without the buffer, a true optimum at
/// the documented edge would trip the boundary guard with a
/// "lies outside the band" diagnostic that contradicted the
/// docstring.
const CALIBRATION_LOG10_N_LO: f64 = -5.301; // log10(5e-6)
const CALIBRATION_LOG10_N_HI: f64 = -1.699; // log10(2e-2)

/// Documented lower / upper edges of the user-supported density
/// interval in `log10(n)` space (`1e-5` and `1e-2` atoms/barn).
/// Used only by the error message so the diagnostic states the
/// edges the user expects to see, not the internal buffered band.
const CALIBRATION_LOG10_N_LO_DOC: f64 = -5.0;
const CALIBRATION_LOG10_N_HI_DOC: f64 = -2.0;

/// Tolerance (in `log10(n)` space) at which the golden-section
/// iteration terminates.  `5e-5` ≈ 0.01 % relative resolution on
/// `n_total`, well below the chi² landscape's per-decade curvature
/// floor for typical SAMMY-style resonance fits.
const CALIBRATION_LOG10_N_TOL: f64 = 5e-5;

/// Tolerance (in `log10(n)` space) for the boundary-saturation
/// guard.  An optimum within `0.02` of either bound — about 5 %
/// in linear density — almost always means the true minimum lies
/// outside the supported band and the user should be told rather
/// than silently handed a railed answer.  The internal search band
/// is widened so the *documented* edges (`1e-5`, `1e-2`) remain
/// strictly inside the tolerance window even after this margin is
/// applied.
const CALIBRATION_LOG10_BOUNDARY_TOL: f64 = 0.02;

/// Golden-section search for the `n_total` that minimises
/// `chi2_of_log_n(log10(n))` on `[CALIBRATION_LOG10_N_LO,
/// CALIBRATION_LOG10_N_HI]`.
///
/// Runs in `log10(n)` so the three-decade interval gets relative
/// resolution.  Returns `(best_n, best_chi2)` — `best_n` is the
/// linear-space optimum, not the log-space value.  Uses the
/// standard two-point golden-section update: maintain `(a, b)`,
/// probe at the two golden-ratio interior points `c, d`, and shrink
/// to whichever half-interval brackets the lower value.  The non-
/// finite case (every chi² along the search returns `+inf` — e.g.
/// `SampleParams::new` rejects the entire density range at this
/// (L, t₀)) returns `(best_n, +inf)` so the outer grid search can
/// move on without latching this candidate.
fn golden_section_n_total<F>(log_lo: f64, log_hi: f64, tol: f64, mut chi2_of_log_n: F) -> (f64, f64)
where
    F: FnMut(f64) -> f64,
{
    // Golden ratio reciprocal: (√5 − 1) / 2 ≈ 0.6180.
    let phi: f64 = (5.0_f64.sqrt() - 1.0) / 2.0;

    let mut a = log_lo;
    let mut b = log_hi;
    let mut c = b - phi * (b - a);
    let mut d = a + phi * (b - a);
    let mut fc = chi2_of_log_n(c);
    let mut fd = chi2_of_log_n(d);

    // Cap iterations defensively in case `tol` is hit by NaN
    // arithmetic; for the canonical (log_lo = -5, log_hi = -2,
    // tol = 5e-5) parameters, convergence is reached in ~25 steps.
    for _ in 0..200 {
        if (b - a) <= tol {
            break;
        }
        if fc < fd {
            b = d;
            d = c;
            fd = fc;
            c = b - phi * (b - a);
            fc = chi2_of_log_n(c);
        } else {
            a = c;
            c = d;
            fc = fd;
            d = a + phi * (b - a);
            fd = chi2_of_log_n(d);
        }
    }

    // The bracket has shrunk to within `tol`; either endpoint of
    // the inner pair is within tolerance of the optimum.  Pick the
    // lower-chi² of the two final probes.
    if fc <= fd {
        (10f64.powf(c), fc)
    } else {
        (10f64.powf(d), fd)
    }
}

/// Result of energy calibration.
#[derive(Debug, Clone)]
pub struct CalibrationResult {
    /// Fitted flight path length in metres.
    pub flight_path_m: f64,
    /// Fitted TOF delay in microseconds.
    pub t0_us: f64,
    /// Fitted total areal density in atoms/barn.
    pub total_density: f64,
    /// Reduced chi-squared at the best (L, t₀, n) values.
    pub reduced_chi_squared: f64,
    /// Corrected energy grid (ascending, eV).
    pub energies_corrected: Vec<f64>,
}

/// Calibrate the energy axis of a TOF neutron measurement.
///
/// Given a measured 1D transmission spectrum and known sample composition
/// (e.g. natural Hf), finds the (L, t₀) that minimize chi² by aligning
/// the ENDF resonance positions with the measured dips.
///
/// # Search strategy
///
/// The optimisation runs as three nested coarse → fine → ultra-fine
/// grid scans on `(L, t₀)`.  At each `(L, t₀)` candidate, the third
/// parameter `n_total` (total areal density, atoms/barn) is
/// optimised by **golden-section search in `log10(n)` space** over
/// the documented user-supported interval `[1e-5, 1e-2]`
/// atoms/barn.  Searching in log space gives uniform relative
/// resolution across the three-decade band, which is necessary
/// because realistic samples span from ~1e-5 (trace) to ~1e-2
/// (1 mm metal foils).
///
/// Internally the golden section runs on a slightly wider band
/// (~5e-6 to ~2e-2) so the boundary-saturation guard's ~5 % linear
/// tolerance does not exclude a true optimum at the documented
/// edges; if the optimum lands inside the tolerance window — i.e.
/// effectively at `1e-5` or `1e-2` — the function returns
/// `Err(PipelineError::InvalidParameter)` rather than a silent
/// boundary-saturated answer, because a true minimum at the edge
/// almost always means the real optimum lies outside the supported
/// interval and the caller should supply a better initial estimate
/// or check the sample composition.
///
/// # Arguments
///
/// * `energies_nominal` — Energy grid computed with assumed L (ascending, eV)
/// * `transmission` — Measured transmission values (same length)
/// * `uncertainty` — Per-bin uncertainty (same length)
/// * `isotopes` — ENDF resonance data for each isotope
/// * `abundances` — Natural abundance fractions (same length as isotopes, sum ≤ 1)
/// * `assumed_flight_path_m` — The L used to compute `energies_nominal`
/// * `temperature_k` — Sample temperature for Doppler broadening
/// * `resolution` — Optional instrument resolution function.  When provided,
///   the forward model includes Doppler + resolution broadening, producing
///   more accurate (L, t₀) fits.  Without resolution, fitted parameters
///   absorb the missing broadening and may be biased.
///
/// # Returns
///
/// [`CalibrationResult`] with the fitted (L, t₀, n_total) and corrected energies.
#[allow(clippy::too_many_arguments)]
pub fn calibrate_energy(
    energies_nominal: &[f64],
    transmission: &[f64],
    uncertainty: &[f64],
    isotopes: &[ResonanceData],
    abundances: &[f64],
    assumed_flight_path_m: f64,
    temperature_k: f64,
    resolution: Option<&InstrumentParams>,
) -> Result<CalibrationResult, PipelineError> {
    let n = energies_nominal.len();
    if n == 0 {
        return Err(PipelineError::InvalidParameter(
            "energies_nominal must not be empty".into(),
        ));
    }
    if transmission.len() != n || uncertainty.len() != n {
        return Err(PipelineError::InvalidParameter(format!(
            "transmission ({}) and uncertainty ({}) must match energies ({})",
            transmission.len(),
            uncertainty.len(),
            n,
        )));
    }
    if isotopes.len() != abundances.len() {
        return Err(PipelineError::InvalidParameter(format!(
            "isotopes ({}) must match abundances ({})",
            isotopes.len(),
            abundances.len(),
        )));
    }

    // Validate abundance values up-front.  Without this guard, non-finite
    // or negative entries are silently multiplied into per-isotope
    // densities (`abd * n_total`), `SampleParams::new` rejects the
    // non-positive thickness, `compute_chi2` returns `INFINITY` for
    // every grid point, and the user sees "no finite chi²" or boundary
    // saturation rather than the actual cause (a bad abundance entry).
    // Equivalent guards already exist for `assumed_flight_path_m` and
    // `energies_nominal`; this closes the same gap for abundances.
    let mut total_abundance = 0.0;
    for (i, &abn) in abundances.iter().enumerate() {
        if !abn.is_finite() || abn < 0.0 {
            return Err(PipelineError::InvalidParameter(format!(
                "calibrate_energy: abundances[{i}] = {abn} is not finite and non-negative"
            )));
        }
        total_abundance += abn;
    }
    if total_abundance <= 0.0 {
        return Err(PipelineError::InvalidParameter(
            "calibrate_energy: sum of abundances must be strictly positive".into(),
        ));
    }

    // Validate scalar / array inputs up-front so the grid-search loop
    // cannot silently produce a "perfect calibration" result from
    // degenerate inputs.  Without these guards, all-NaN transmission
    // combined with the dof=1 fallback below would cause
    // chi²_reduced = 0.0 to be reported as a successful fit.
    if !assumed_flight_path_m.is_finite() || assumed_flight_path_m <= 0.0 {
        return Err(PipelineError::InvalidParameter(format!(
            "assumed_flight_path_m must be finite and positive, got {assumed_flight_path_m}",
        )));
    }
    for (i, &e) in energies_nominal.iter().enumerate() {
        if !e.is_finite() || e <= 0.0 {
            return Err(PipelineError::InvalidParameter(format!(
                "energies_nominal[{i}] must be finite and positive, got {e}",
            )));
        }
        if i > 0 && e <= energies_nominal[i - 1] {
            return Err(PipelineError::InvalidParameter(format!(
                "energies_nominal must be strictly ascending; \
                 energies_nominal[{i}]={e} <= energies_nominal[{}]={}",
                i - 1,
                energies_nominal[i - 1],
            )));
        }
    }

    // Recover TOF from nominal energies: t = L_assumed · √(C / E)
    let tof_s: Vec<f64> = energies_nominal
        .iter()
        .map(|&e| assumed_flight_path_m * (NEUTRON_MASS_CONSTANT / e).sqrt())
        .collect();

    // Pre-filter valid bins (finite T, positive sigma)
    let valid: Vec<bool> = transmission
        .iter()
        .zip(uncertainty.iter())
        .map(|(&t, &s)| t.is_finite() && s.is_finite() && s > 0.0)
        .collect();

    // Require enough valid bins to constrain the three fitted
    // parameters (L, t₀, n_total).  Previously, when every bin was
    // invalid, `compute_chi2` returned 0.0 for every grid point, the
    // first candidate latched as "best", and the dof=1 fallback turned
    // that into a reported `chi²_reduced = 0.0` — i.e. a totally
    // degenerate input was indistinguishable from a perfect calibration.
    const N_FITTED_PARAMS: usize = 3;
    let n_valid = valid.iter().filter(|&&v| v).count();
    if n_valid < N_FITTED_PARAMS {
        return Err(PipelineError::InvalidParameter(format!(
            "calibrate_energy requires at least {N_FITTED_PARAMS} bins with finite \
             transmission and positive uncertainty, got {n_valid} valid out of {n}",
        )));
    }

    // ── Phase 1: Coarse grid search over (L, t₀) ───────────────────
    // L: ±1.5 % around assumed (0.1 % steps → 31 points)
    // t₀: -5 to +10 µs (1 µs steps → 16 points)
    // n_total: at each (L, t₀), golden-section search in log10(n)
    //          on the full `[1e-5, 1e-2]` density band.  The
    //          previous implementation used a 5-point hard-coded
    //          scan `{5e-5, 1e-4, 1.5e-4, 2e-4, 3e-4}` followed by
    //          multiplicative refinements, which left the final
    //          density anchored inside a `[2.25e-5, 4.95e-4]` band
    //          — incompatible with the 1 mm metal-foil densities
    //          (U/W/Ni at ~5e-3) the paper relies on.

    let l_center = assumed_flight_path_m;
    let mut best_chi2 = f64::INFINITY;
    let mut best_l = l_center;
    let mut best_t0_us = 0.0f64;
    let mut best_n = 1e-4;

    // Coarse L: 0.2% steps, ±1.5%
    let l_steps: Vec<f64> = (-15..=15)
        .map(|i| l_center * (1.0 + i as f64 * 0.001))
        .collect();
    // Coarse t₀: 1 µs steps, -5 to +10 µs
    let t0_steps: Vec<f64> = (-5..=10).map(|i| i as f64).collect();

    for &l in &l_steps {
        for &t0 in &t0_steps {
            let t0_s = t0 * 1e-6;
            // Correct energies
            let e_corr: Vec<f64> = tof_s
                .iter()
                .map(|&t| {
                    let t_corr = t - t0_s;
                    if t_corr <= 0.0 {
                        f64::NAN
                    } else {
                        NEUTRON_MASS_CONSTANT * (l / t_corr).powi(2)
                    }
                })
                .collect();

            // Skip if any NaN
            if e_corr.iter().any(|e| !e.is_finite() || *e <= 0.0) {
                continue;
            }

            // Optimise n_total at this (L, t₀) by golden section in
            // log10(n) over the full configured search band.
            let (n_opt, chi2_opt) = golden_section_n_total(
                CALIBRATION_LOG10_N_LO,
                CALIBRATION_LOG10_N_HI,
                CALIBRATION_LOG10_N_TOL,
                |log_n| {
                    compute_chi2(
                        &e_corr,
                        transmission,
                        uncertainty,
                        isotopes,
                        abundances,
                        10f64.powf(log_n),
                        temperature_k,
                        &valid,
                        resolution,
                    )
                },
            );
            if chi2_opt < best_chi2 {
                best_chi2 = chi2_opt;
                best_l = l;
                best_t0_us = t0;
                best_n = n_opt;
            }
        }
    }

    // ── Phase 2: Fine grid search around coarse (L, t₀) best ───────
    // L: ±0.05%, 0.01% steps
    // t₀: ±2 µs, 0.25 µs steps
    // n_total: golden-section in log10(n) on the full search band
    //          at each (L, t₀) candidate, same as Phase 1.  Running
    //          the same minimiser over the full band — rather than
    //          a ±50 % window around the Phase-1 winner — guards
    //          against the chi² landscape's coupling between
    //          (L, t₀) and density: a coarse-grid winner can sit on
    //          a slightly biased density that the Phase-2 (L, t₀)
    //          refinement should be allowed to walk away from.

    let l_fine: Vec<f64> = (-5..=5)
        .map(|i| best_l * (1.0 + i as f64 * 0.0001))
        .collect();
    let t0_fine: Vec<f64> = (-8..=8).map(|i| best_t0_us + i as f64 * 0.25).collect();

    for &l in &l_fine {
        for &t0 in &t0_fine {
            let t0_s = t0 * 1e-6;
            let e_corr: Vec<f64> = tof_s
                .iter()
                .map(|&t| {
                    let t_corr = t - t0_s;
                    if t_corr <= 0.0 {
                        f64::NAN
                    } else {
                        NEUTRON_MASS_CONSTANT * (l / t_corr).powi(2)
                    }
                })
                .collect();
            if e_corr.iter().any(|e| !e.is_finite() || *e <= 0.0) {
                continue;
            }
            let (n_opt, chi2_opt) = golden_section_n_total(
                CALIBRATION_LOG10_N_LO,
                CALIBRATION_LOG10_N_HI,
                CALIBRATION_LOG10_N_TOL,
                |log_n| {
                    compute_chi2(
                        &e_corr,
                        transmission,
                        uncertainty,
                        isotopes,
                        abundances,
                        10f64.powf(log_n),
                        temperature_k,
                        &valid,
                        resolution,
                    )
                },
            );
            if chi2_opt < best_chi2 {
                best_chi2 = chi2_opt;
                best_l = l;
                best_t0_us = t0;
                best_n = n_opt;
            }
        }
    }

    // ── Phase 3: Ultra-fine refinement ──────────────────────────────
    // L: ±0.005%, 0.001% steps
    // t₀: ±0.5 µs, 0.05 µs steps
    // n_total: golden-section in log10(n) on the full search band
    //          at each (L, t₀) candidate.

    let l_ultra: Vec<f64> = (-5..=5)
        .map(|i| best_l * (1.0 + i as f64 * 0.00001))
        .collect();
    let t0_ultra: Vec<f64> = (-10..=10).map(|i| best_t0_us + i as f64 * 0.05).collect();

    for &l in &l_ultra {
        for &t0 in &t0_ultra {
            let t0_s = t0 * 1e-6;
            let e_corr: Vec<f64> = tof_s
                .iter()
                .map(|&t| {
                    let t_corr = t - t0_s;
                    if t_corr <= 0.0 {
                        f64::NAN
                    } else {
                        NEUTRON_MASS_CONSTANT * (l / t_corr).powi(2)
                    }
                })
                .collect();
            if e_corr.iter().any(|e| !e.is_finite() || *e <= 0.0) {
                continue;
            }
            let (n_opt, chi2_opt) = golden_section_n_total(
                CALIBRATION_LOG10_N_LO,
                CALIBRATION_LOG10_N_HI,
                CALIBRATION_LOG10_N_TOL,
                |log_n| {
                    compute_chi2(
                        &e_corr,
                        transmission,
                        uncertainty,
                        isotopes,
                        abundances,
                        10f64.powf(log_n),
                        temperature_k,
                        &valid,
                        resolution,
                    )
                },
            );
            if chi2_opt < best_chi2 {
                best_chi2 = chi2_opt;
                best_l = l;
                best_t0_us = t0;
                best_n = n_opt;
            }
        }
    }

    // Post-grid-search sanity check: if every `compute_chi2` call along
    // the entire 3-phase grid returned `f64::INFINITY` (e.g. because
    // `SampleParams::new` or `forward_model` failed at every candidate,
    // or because the residuals overflowed for non-finite-but-passing
    // transmission values such as 1e308), `best_chi2` remains
    // `INFINITY` and the caller would otherwise receive a
    // `CalibrationResult { reduced_chi_squared: inf, .. }` — the same
    // silent-failure class as the zero-valid-bins case the up-front
    // guard now rejects.  Reject explicitly here so calibration
    // failure is always an `Err`, never an `Ok` with a sentinel chi².
    if !best_chi2.is_finite() {
        return Err(PipelineError::InvalidParameter(format!(
            "calibrate_energy: grid search produced no finite chi² across all \
             (L, t₀, n_total) candidates — likely cause is forward-model failure \
             or non-finite residuals (e.g. wildly out-of-scale transmission); \
             best_chi2 = {best_chi2}",
        )));
    }

    // Boundary-saturation guard: if the n_total optimum lies within
    // tolerance of either density bound, the true minimum almost
    // certainly sits outside the supported range and the
    // calibration is unreliable.  Returning `Ok` with `best_n` ≈
    // boundary would silently rail the density and let the (L, t₀)
    // parameters absorb the missing density freedom by compensating
    // bias — exactly the silent-failure pattern the post-search
    // chi² guard above also defends against, but with a boundary-
    // specific diagnostic.
    //
    // The internal search band is `[~5e-6, ~2e-2]`; the documented
    // edges are `[1e-5, 1e-2]`.  Densities at the documented edges
    // lie outside the tolerance window (a true optimum at `1e-5`
    // sits ~0.3 decades above the internal lower bound, comfortably
    // past the ~0.02-log10 tolerance) so the guard fires only when
    // the optimum has actually saturated against the wider buffer.
    let log_best_n = best_n.log10();
    let n_lo = 10f64.powf(CALIBRATION_LOG10_N_LO_DOC);
    let n_hi = 10f64.powf(CALIBRATION_LOG10_N_HI_DOC);
    if (log_best_n - CALIBRATION_LOG10_N_LO).abs() < CALIBRATION_LOG10_BOUNDARY_TOL
        || (CALIBRATION_LOG10_N_HI - log_best_n).abs() < CALIBRATION_LOG10_BOUNDARY_TOL
    {
        return Err(PipelineError::InvalidParameter(format!(
            "calibrate_energy: n_total optimum {best_n:.3e} atoms/barn is at the \
             search boundary [{n_lo:.0e}, {n_hi:.0e}]; the true optimum likely lies \
             outside this band.  Provide a better initial density estimate, check \
             the sample composition / abundances, or extend the search range."
        )));
    }

    // Compute corrected energy grid at the best parameters
    let t0_best_s = best_t0_us * 1e-6;
    let energies_corrected: Vec<f64> = tof_s
        .iter()
        .map(|&t| NEUTRON_MASS_CONSTANT * (best_l / (t - t0_best_s)).powi(2))
        .collect();

    // Final chi2r (reduced).  The up-front guard ensures
    // `n_valid >= N_FITTED_PARAMS`, so we always have a non-negative
    // dof.  We still clamp to `max(1)` so that the exact-fit edge case
    // (n_valid == N_FITTED_PARAMS, dof = 0) reports a finite value
    // instead of dividing by zero.
    let dof = n_valid.saturating_sub(N_FITTED_PARAMS).max(1);
    let chi2r = best_chi2 / dof as f64;

    Ok(CalibrationResult {
        flight_path_m: best_l,
        t0_us: best_t0_us,
        total_density: best_n,
        reduced_chi_squared: chi2r,
        energies_corrected,
    })
}

/// Compute total chi² for a given (E_corrected, n_total) against measured data.
#[allow(clippy::too_many_arguments)]
fn compute_chi2(
    energies: &[f64],
    transmission: &[f64],
    uncertainty: &[f64],
    isotopes: &[ResonanceData],
    abundances: &[f64],
    n_total: f64,
    temperature_k: f64,
    valid: &[bool],
    resolution: Option<&InstrumentParams>,
) -> f64 {
    // Build (isotope, density) pairs
    let pairs: Vec<(ResonanceData, f64)> = isotopes
        .iter()
        .zip(abundances.iter())
        .map(|(iso, &abd)| (iso.clone(), abd * n_total))
        .collect();

    let sample = match SampleParams::new(temperature_k, pairs) {
        Ok(s) => s,
        Err(_) => return f64::INFINITY,
    };

    // P-5: Include resolution broadening when available.
    // Without it, fitted L and t₀ absorb the missing broadening bias.
    let model = match transmission::forward_model(energies, &sample, resolution) {
        Ok(m) => m,
        Err(_) => return f64::INFINITY,
    };

    // Chi²
    let mut chi2 = 0.0;
    for (i, (&t_data, &t_model)) in transmission.iter().zip(model.iter()).enumerate() {
        if !valid[i] {
            continue;
        }
        let residual = (t_data - t_model) / uncertainty[i];
        chi2 += residual * residual;
    }
    chi2
}

#[cfg(test)]
mod tests {
    use super::*;
    use nereids_endf::resonance::test_support::synthetic_single_resonance;

    /// Round-trip exercise of the public `calibrate_energy` API on a
    /// synthetic spectrum.  Uses `synthetic_single_resonance` from
    /// `nereids_endf::resonance::test_support` so the test does not require network access and
    /// runs in every CI invocation.
    ///
    /// Note on tolerances: the grid-search calibrator's L resolution
    /// is fundamentally limited by chi² curvature (broader resonances
    /// or sparser bins → broader minimum).  With only synthetic
    /// single-resonance isotopes on a sparse 0.2 eV grid, the chi²
    /// landscape near L=true_l is shallow on the scale of the
    /// 0.001 % ultra-fine step — Doppler-broadened ≈ 33 meV resonance
    /// width vs 200 meV grid spacing means each resonance is sampled
    /// by ≤ 1 bin, so (L, t₀) can drift across a wide band before chi²
    /// degrades enough to lock the minimum down.  We therefore (a) use
    /// a small true offset (0.05 % in L, 0.5 µs in t₀) that the
    /// calibrator can resolve, and (b) test the *physics* — a fit
    /// converged, the corrected energies are close to truth on the
    /// data-relevant range, and density recovery is in the right
    /// decade — rather than chasing exact ENDF-style L recovery.
    /// The bit-exact precision question is owned by the SAMMY parity
    /// tests in `nereids-physics`, not by this API smoke test.
    #[test]
    fn test_calibrate_round_trip_synthetic() {
        // Generate synthetic data with known L and t0, then recover them.
        // Small offsets (0.05 % in L, 0.5 µs in t₀) so the chi² minimum
        // is well inside Phase-2 fine grid (±0.05 % in L, ±2 µs in t₀).
        //
        // Setup (resonances, energy grid, forward-model transmission,
        // and the `calibrate_energy` call itself) is shared with the
        // density-band tests below via `calibrate_round_trip_at_density`;
        // the helper returns `(result, e_true, assumed_l)` so this
        // smoke test can also assert on corrected-energy accuracy.
        let true_n = 1.5e-4;
        let (result, e_true, assumed_l) =
            calibrate_round_trip_at_density(true_n).expect("Calibration failed");

        // Check recovery.  Wider tolerances than the Hf-178 fixture
        // because the synthetic chi² minimum is broader (see the
        // doc comment on the test above).  These bands still
        // distinguish a successful fit from a degenerate one (the
        // zero-valid-bins failure mode would report L = assumed_l
        // and chi² = 0.0).
        //
        // With the n_total golden-section refactor, the previously-
        // narrow density grid no longer pins (L, t₀) at a single
        // coarse-grid point; the calibrator now expresses the
        // genuine (L, t₀, n) degeneracy this sparse-grid synthetic
        // admits — 33 meV Doppler-broadened resonances vs 200 meV
        // grid spacing samples each resonance with ≤ 1 bin, so any
        // (L, t₀) that places the resonance near the same bin gives
        // an indistinguishable fit.  The L and t₀ parameters are
        // therefore not independently identifiable from this
        // synthetic, but the *corrected energy grid* — the actual
        // downstream deliverable — is, and is the right thing to
        // assert on.
        //
        // L and t₀ are still required to be inside the search grid
        // (Phase 1 L ±1.5 %, t₀ ∈ [-5, +10] µs) and density inside
        // the search band — anything outside would signal a
        // calibrator regression, not a degeneracy.
        assert!(
            (result.flight_path_m - assumed_l).abs() / assumed_l <= 0.015,
            "L drift outside Phase 1 ±1.5 % grid: got {}",
            result.flight_path_m,
        );
        assert!(
            (-5.0..=10.0).contains(&result.t0_us),
            "t0 outside Phase 1 grid: got {}",
            result.t0_us,
        );
        assert!(
            (result.total_density - true_n).abs() / true_n < 0.5,
            "n: got {}, expected {}",
            result.total_density,
            true_n,
        );

        // The corrected energy grid is the deliverable a downstream
        // fit uses; even when (L, t₀, n) drift inside the chi²
        // basin allowed by this sparse-grid synthetic, the corrected
        // energies must still track `e_true` to within ~1 % at the
        // resonance positions and a few % across the grid.  A
        // median relative error check is more robust to per-bin
        // scaling than max-over-bins on a 150-bin grid; we still
        // assert max < 5 % to catch a wholly-wrong calibration.
        let rel_errs: Vec<f64> = result
            .energies_corrected
            .iter()
            .zip(e_true.iter())
            .map(|(&ec, &et)| (ec - et).abs() / et)
            .collect();
        let mut sorted = rel_errs.clone();
        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
        let median = sorted[sorted.len() / 2];
        let max_err = sorted.last().copied().unwrap_or(0.0);
        assert!(
            median < 0.02,
            "corrected energies median rel err {median} should be < 2 %",
        );
        assert!(
            max_err < 0.05,
            "corrected energies max rel err {max_err} should be < 5 %",
        );
        // The synthetic round-trip uses noiseless data, so a grid point
        // that lands on (or arbitrarily close to) the true parameters
        // can legitimately yield chi² = 0.  Accept `>= 0.0` (the
        // physically valid range) rather than `> 0.0`, which was a
        // flake-prone strict-inequality.  The zero-valid-bins
        // regression that previously reported chi² = 0.0 is now
        // rejected up-front by the `n_valid >= N_FITTED_PARAMS` guard,
        // and the all-infinity grid-search case is rejected by the
        // post-search `best_chi2.is_finite()` guard — so finiteness
        // alone is the meaningful check here.
        assert!(
            result.reduced_chi_squared.is_finite() && result.reduced_chi_squared >= 0.0,
            "chi²_reduced must be finite and >= 0 (degenerate-input regressions \
             are caught up-front and post-search; this assertion guards against \
             chi² leaking as inf or NaN); got {}",
            result.reduced_chi_squared,
        );
    }

    // ── Degenerate-input guards ────────────────────────────────────────
    //
    // Before these guards, `compute_chi2` returned 0.0 when every bin
    // was skipped by the `valid` mask, the grid search latched the
    // first candidate as "best", and the dof=1 fallback at the end
    // turned that into a reported `chi²_reduced = 0.0` — a totally
    // degenerate input was indistinguishable from a perfect
    // calibration.

    /// `(energies_nominal, transmission, uncertainty, isotopes, abundances)`
    /// — the five array arguments to `calibrate_energy`.
    type CalibrationInputs = (Vec<f64>, Vec<f64>, Vec<f64>, Vec<ResonanceData>, Vec<f64>);

    /// Build a minimal valid input set for `calibrate_energy`, then let
    /// the caller mutate one field to drive a specific error path.
    fn minimal_calibration_inputs() -> CalibrationInputs {
        let iso = synthetic_single_resonance(72, 178, 176.4, 7.8);
        let energies: Vec<f64> = (0..50).map(|i| 5.0 + i as f64 * 0.4).collect();
        let transmission = vec![0.95; energies.len()];
        let uncertainty = vec![0.01; energies.len()];
        (energies, transmission, uncertainty, vec![iso], vec![1.0])
    }

    #[test]
    fn test_calibrate_all_nan_transmission_rejected() {
        // All-NaN transmission would previously yield zero valid bins,
        // compute_chi2() returned 0.0 for every grid point, and the
        // dof=1 fallback reported chi²_reduced = 0.0 as success.
        let (energies, mut transmission, uncertainty, isotopes, abundances) =
            minimal_calibration_inputs();
        for t in transmission.iter_mut() {
            *t = f64::NAN;
        }
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("all-NaN transmission must be rejected");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("valid"),
                    "error message should mention valid-bin count, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }

    #[test]
    fn test_calibrate_all_zero_uncertainty_rejected() {
        // All-zero uncertainty is the other path to zero valid bins
        // (sigma > 0 is required by the valid mask).
        let (energies, transmission, mut uncertainty, isotopes, abundances) =
            minimal_calibration_inputs();
        for s in uncertainty.iter_mut() {
            *s = 0.0;
        }
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("all-zero uncertainty must be rejected");
        assert!(
            matches!(err, PipelineError::InvalidParameter(_)),
            "expected InvalidParameter, got {err:?}"
        );
    }

    #[test]
    fn test_calibrate_nonfinite_flight_path_rejected() {
        // All non-finite or non-positive flight paths must produce
        // InvalidParameter, naming the offending field so the caller
        // can diagnose the source.
        let (energies, transmission, uncertainty, isotopes, abundances) =
            minimal_calibration_inputs();
        for bad_l in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY, 0.0, -1.0] {
            let result = calibrate_energy(
                &energies,
                &transmission,
                &uncertainty,
                &isotopes,
                &abundances,
                bad_l,
                293.6,
                None,
            );
            match result {
                Ok(_) => panic!("expected Err for L={bad_l}, got Ok"),
                Err(PipelineError::InvalidParameter(msg)) => {
                    assert!(
                        msg.contains("assumed_flight_path_m"),
                        "error message should name the offending field for L={bad_l}, got: {msg}"
                    );
                }
                Err(other) => panic!("expected InvalidParameter for L={bad_l}, got {other:?}"),
            }
        }
    }

    #[test]
    fn test_calibrate_nonascending_energies_rejected() {
        let (mut energies, transmission, uncertainty, isotopes, abundances) =
            minimal_calibration_inputs();
        // Introduce a non-ascending pair.
        energies[10] = energies[9];
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("non-ascending energies must be rejected");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("ascending"),
                    "error message should mention ascending, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }

    #[test]
    fn test_calibrate_all_infinite_chi2_rejected() {
        // Regression for the post-grid-search `best_chi2.is_finite()`
        // guard.  Before the guard, an input whose every `compute_chi2`
        // evaluation returned `f64::INFINITY` left `best_chi2`
        // initialised at `INFINITY`, no candidate ever beat it, and the
        // function returned `Ok(CalibrationResult { reduced_chi_squared:
        // inf, .. })` — the same silent-failure surface as the
        // zero-valid-bins case, just with infinity instead of zero.
        //
        // Driving the all-infinity path: feed finite but wildly
        // out-of-scale transmission (1e308).  Each finite uncertainty
        // (0.01) makes the residual ((1e308 − T_model) / 0.01) overflow
        // to `+inf`, residual² is `inf`, the per-bin sum is `inf`, and
        // every grid candidate returns `inf`.  Crucially, the
        // transmission values stay `finite()` so they pass the up-front
        // `t.is_finite()` mask and the new post-search guard is the
        // only line of defence.
        let (energies, _transmission, uncertainty, isotopes, abundances) =
            minimal_calibration_inputs();
        let transmission = vec![1e308; energies.len()];
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("all-infinity chi² across the grid must be rejected");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("finite chi²") || msg.contains("best_chi2"),
                    "error message should explain the all-infinity grid-search failure, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }

    #[test]
    fn test_calibrate_nonfinite_energy_rejected() {
        let (mut energies, transmission, uncertainty, isotopes, abundances) =
            minimal_calibration_inputs();
        energies[7] = f64::NAN;
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("NaN energy must be rejected");
        assert!(
            matches!(err, PipelineError::InvalidParameter(_)),
            "expected InvalidParameter, got {err:?}"
        );
    }

    #[test]
    fn test_calibrate_energy_rejects_negative_abundance() {
        // A negative abundance silently flips the sign of `abd * n_total`
        // and `SampleParams::new` rejects the non-positive thickness;
        // every grid point then returns chi² = INFINITY and the user
        // sees a "no finite chi²" boundary error.  The up-front guard
        // converts this to an actionable diagnostic that names the
        // offending index.
        let (energies, transmission, uncertainty, isotopes, mut abundances) =
            minimal_calibration_inputs();
        abundances[0] = -0.5;
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("negative abundance must be rejected");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("abundances[0]"),
                    "error message should name the offending index, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }

    #[test]
    fn test_calibrate_energy_rejects_nan_abundance() {
        // NaN bypasses naive `< 0.0` guards (NaN comparisons are always
        // false), so the up-front check must pair `is_finite()` with the
        // sign predicate.  Without this guard, `abd * n_total` is NaN
        // for every density, every chi² is NaN, and the user sees a
        // confusing boundary-saturation message.
        let (energies, transmission, uncertainty, isotopes, mut abundances) =
            minimal_calibration_inputs();
        abundances[0] = f64::NAN;
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("NaN abundance must be rejected");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("abundances[0]"),
                    "error message should name the offending index, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }

    #[test]
    fn test_calibrate_energy_rejects_all_zero_abundances() {
        // Each individual zero abundance is legal (the isotope is simply
        // not present in this sample), but the sum being zero means
        // every per-isotope density is zero and the transmission model
        // collapses to T == 1 — the calibrator has no signal to fit
        // (L, t₀, n_total) against.  Reject up-front rather than letting
        // the search bottom out at the band boundary.
        let (energies, transmission, uncertainty, isotopes, _) = minimal_calibration_inputs();
        let abundances = vec![0.0; isotopes.len()];
        let err = calibrate_energy(
            &energies,
            &transmission,
            &uncertainty,
            &isotopes,
            &abundances,
            25.0,
            293.6,
            None,
        )
        .expect_err("all-zero abundances must be rejected");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("sum of abundances"),
                    "error message should mention the zero-sum cause, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }

    // ── n_total search-band regression tests ──────────────────────────
    //
    // Before the golden-section refactor, `calibrate_energy` scanned
    // n_total at a hard-coded 5-point linear grid `{5e-5, 1e-4,
    // 1.5e-4, 2e-4, 3e-4}` and then refined multiplicatively, leaving
    // the final density anchored inside `[2.25e-5, 4.95e-4]`
    // atoms/barn — incompatible with every realistic VENUS /
    // SoftwareX paper density (1 mm metal foils at ~5e-3, trace
    // matrix densities up to ~1e-2).  The refactor replaces the
    // multi-stage density refinement with a true golden-section
    // search in log10(n) on `[1e-5, 1e-2]`, and adds a boundary-
    // saturation guard for the (now possible) case where the
    // optimum lies outside that band.
    //
    // The tests below verify recovery at three representative
    // densities that span the search range (1e-5 lower edge,
    // 1e-3 middle of band that the old code could not reach,
    // 5e-3 SoftwareX U-238 density that was a factor-10× outside
    // the old reachable max) plus the explicit boundary-failure
    // diagnostic at densities outside the band.

    /// Synthetic round-trip helper parameterised on `true_n`.  Builds
    /// data with two well-separated Hf-style resonances at the given
    /// true density, runs `calibrate_energy`, and on success returns
    /// `(result, e_true, assumed_l)` so individual tests can assert on
    /// density recovery and chi² finiteness (most callers) or on
    /// corrected-energy accuracy against `e_true` (the
    /// `test_calibrate_round_trip_synthetic` smoke test).
    fn calibrate_round_trip_at_density(
        true_n: f64,
    ) -> Result<(CalibrationResult, Vec<f64>, f64), PipelineError> {
        let true_l = 25.0125;
        let assumed_l = 25.0;
        let true_t0_us = 0.5;
        let temperature_k = 293.6;

        let iso_a = synthetic_single_resonance(72, 178, 176.4, 7.8);
        let iso_b = synthetic_single_resonance(72, 178, 176.4, 22.0);
        let isotopes = vec![iso_a, iso_b];
        let abundances = vec![0.5, 0.5];

        let e_nominal: Vec<f64> = (0..150).map(|i| 5.0 + i as f64 * 0.2).collect();
        let tof_s: Vec<f64> = e_nominal
            .iter()
            .map(|&e| assumed_l * (NEUTRON_MASS_CONSTANT / e).sqrt())
            .collect();
        let true_t0_s = true_t0_us * 1e-6;
        let e_true: Vec<f64> = tof_s
            .iter()
            .map(|&t| NEUTRON_MASS_CONSTANT * (true_l / (t - true_t0_s)).powi(2))
            .collect();

        let pairs: Vec<_> = isotopes
            .iter()
            .zip(abundances.iter())
            .map(|(iso, &abd)| (iso.clone(), abd * true_n))
            .collect();
        let sample = SampleParams::new(temperature_k, pairs).expect("SampleParams creation failed");
        let t_model =
            transmission::forward_model(&e_true, &sample, None).expect("forward_model failed");
        let sigma = vec![0.01; e_nominal.len()];

        let result = calibrate_energy(
            &e_nominal,
            &t_model,
            &sigma,
            &isotopes,
            &abundances,
            assumed_l,
            temperature_k,
            None,
        )?;
        Ok((result, e_true, assumed_l))
    }

    /// Near-lower-edge density: `true_n = 2e-5` sits just inside the
    /// `[1e-5, 1e-2]` documented user-supported interval — approximately
    /// 0.3 decades (a factor of 2) above the lower documented edge,
    /// comfortably outside the boundary guard's `5 %`-linear tolerance
    /// window.  Recovery must succeed.  Note the 30 % relative tolerance —
    /// at this low density the chi² landscape is shallow (single-resonance
    /// synthetic, weak signal), so the recovered density can drift further
    /// from the true value than at mid-band; the test still meaningfully
    /// distinguishes "we found roughly the right decade" from the old
    /// behaviour of being unable to reach the value at all.  The
    /// `test_calibrate_energy_boundary_saturation_error` test separately
    /// verifies the guard fires for genuinely-out-of-band densities.
    #[test]
    fn test_calibrate_energy_recovers_density_1e_5() {
        let true_n = 2e-5;
        let (result, _, _) = calibrate_round_trip_at_density(true_n)
            .expect("calibration at true_n=2e-5 must succeed");
        assert!(
            (result.total_density - true_n).abs() / true_n < 0.3,
            "n: got {}, expected {}",
            result.total_density,
            true_n,
        );
        assert!(result.reduced_chi_squared.is_finite());
    }

    /// Documented lower bound: `true_n = 1.0e-5` atoms/barn is exactly
    /// the lower edge promised by the `calibrate_energy` rustdoc.  Before
    /// the search-band widening the boundary-saturation guard's
    /// `~5 %`-linear tolerance trimmed a sliver off either side and
    /// rejected truly-at-the-edge optima with a "true optimum likely
    /// lies outside this band" diagnostic that contradicted the docs.
    /// With the internal band widened to `[~5e-6, ~2e-2]`, an optimum
    /// at the documented edge sits ~0.3 decades inside the buffer and
    /// is accepted — the user-facing contract here is that the call
    /// returns `Ok(_)` (no boundary-saturation error) and recovers a
    /// density close to truth in log-space, not that the recovered
    /// value is bit-exactly bounded by the documented interval: chi²
    /// minimisation can land slightly outside `[1e-5, 1e-2]` even when
    /// the true density sits at the edge, and that is correct
    /// behaviour for a smooth optimisation landscape.
    #[test]
    fn test_calibrate_energy_accepts_density_at_documented_lower_bound() {
        let true_n = 1.0e-5;
        let (result, _, _) = calibrate_round_trip_at_density(true_n)
            .expect("calibration at the documented lower edge 1e-5 must succeed");
        // Log-space tolerance because the chi² landscape is shallow at
        // this trace density (single-resonance synthetic, weak signal):
        // a recovered-vs-truth ratio of 2× corresponds to 0.3 in
        // log10(n) and is the empirically reasonable precision floor.
        let log_err = (result.total_density.log10() - true_n.log10()).abs();
        assert!(
            log_err < 0.3,
            "log10(n) error {log_err} too large; recovered {} vs truth {true_n}",
            result.total_density,
        );
        assert!(result.reduced_chi_squared.is_finite());
    }

    /// Documented upper bound: `true_n = 1.0e-2` atoms/barn is exactly
    /// the upper edge promised by `calibrate_energy`'s rustdoc — the
    /// `1 mm metal foil` use case that drives the SoftwareX paper's
    /// calibration narrative.  Sister test to
    /// `test_calibrate_energy_accepts_density_at_documented_lower_bound`;
    /// the search-band widening keeps the upper documented edge inside
    /// the boundary guard's tolerance buffer.  As with the lower-bound
    /// test, the assertion is on chi²-resolution recovery (log-space
    /// proximity to truth), not on hard-bounding the result inside
    /// `[1e-5, 1e-2]` — a smooth optimum at the edge can land just
    /// outside without indicating any defect.
    #[test]
    fn test_calibrate_energy_accepts_density_at_documented_upper_bound() {
        let true_n = 1.0e-2;
        let (result, _, _) = calibrate_round_trip_at_density(true_n)
            .expect("calibration at the documented upper edge 1e-2 must succeed");
        // Tighter log-space tolerance than the lower edge: at this
        // high density the resonance is saturated, the chi²
        // landscape is sharp and the (L, t₀, n) trade-off basin is
        // narrow.  log_err < 0.05 ≈ 12 % linear is comfortable.
        let log_err = (result.total_density.log10() - true_n.log10()).abs();
        assert!(
            log_err < 0.05,
            "log10(n) error {log_err} too large; recovered {} vs truth {true_n}",
            result.total_density,
        );
        assert!(result.reduced_chi_squared.is_finite());
    }

    /// Phase-1-grid-point density: `true_n = 1e-4` was the historical
    /// reachable-band centre.  Recovery is the easiest case for the
    /// chi² landscape and tightens the tolerance accordingly.
    #[test]
    fn test_calibrate_energy_recovers_density_1e_4() {
        let true_n = 1e-4;
        let (result, _, _) = calibrate_round_trip_at_density(true_n)
            .expect("calibration at true_n=1e-4 must succeed");
        assert!(
            (result.total_density - true_n).abs() / true_n < 0.1,
            "n: got {}, expected {}",
            result.total_density,
            true_n,
        );
        assert!(result.reduced_chi_squared.is_finite());
    }

    /// Mid-band density: `true_n = 1e-3` was **unreachable** under
    /// the previous 5-point scan + multiplicative refinement; the
    /// best the old code could return was ~4.95e-4.  After the
    /// log-space golden-section refactor this density must round-
    /// trip with full Phase-3 precision.
    #[test]
    fn test_calibrate_energy_recovers_density_1e_3() {
        let true_n = 1e-3;
        let (result, _, _) = calibrate_round_trip_at_density(true_n)
            .expect("calibration at true_n=1e-3 must succeed");
        assert!(
            (result.total_density - true_n).abs() / true_n < 0.1,
            "n: got {}, expected {} — old code saturated at ~4.95e-4",
            result.total_density,
            true_n,
        );
        assert!(result.reduced_chi_squared.is_finite());
    }

    /// SoftwareX U-238 reference density: 1 mm metal foil at ~5e-3
    /// atoms/barn.  This is the density the paper figure scripts
    /// (`gen_fig_physics.py`, `gen_fig_closed_loop.py`) use and
    /// that the paper's calibration narrative relies on; it sits a
    /// factor 10× above the old reachable maximum.
    #[test]
    fn test_calibrate_energy_recovers_density_5e_3() {
        let true_n = 5e-3;
        let (result, _, _) = calibrate_round_trip_at_density(true_n)
            .expect("calibration at true_n=5e-3 must succeed");
        assert!(
            (result.total_density - true_n).abs() / true_n < 0.1,
            "n: got {}, expected {} — old code saturated at ~4.95e-4 (10× too low)",
            result.total_density,
            true_n,
        );
        assert!(result.reduced_chi_squared.is_finite());
    }

    /// Out-of-band saturation: `true_n = 1.0` atoms/barn is two
    /// orders of magnitude above the upper search bound (`1e-2`).
    /// The golden-section minimum must land on the upper bound,
    /// and the boundary-saturation guard must turn that into an
    /// `Err(InvalidParameter)` rather than the silent railed answer
    /// the old 5-point scan would have returned (the old code
    /// would have railed to `4.95e-4`, six orders of magnitude
    /// below truth, with no diagnostic).
    #[test]
    fn test_calibrate_energy_boundary_saturation_error() {
        let true_n = 1.0;
        let err = calibrate_round_trip_at_density(true_n)
            .expect_err("density outside search band must trigger boundary guard");
        match err {
            PipelineError::InvalidParameter(msg) => {
                assert!(
                    msg.contains("search boundary") || msg.contains("boundary"),
                    "error must explain boundary saturation, got: {msg}"
                );
            }
            other => panic!("expected InvalidParameter, got {other:?}"),
        }
    }
}