[ESSSANS] Q resolution per Q bin, compatible with merging LoKI banks #763

Summary

Compute the Q resolution of I(Q) in esssans so that LoKI detector banks (as well as runs and live-data chunks) can be merged before normalization. Output one σ per Q bin (NXcanSAS Qdev), and optionally an approximation of the full resolution function.

This refines #360. It replaces the approach of scipp/esssans#185, which computed resolution per pixel and wavelength but did not say how to combine the values in a Q bin.

Interactive explainer with the toy model used for the numbers below: https://claude.ai/artifact/JhbFrDS2YKZMU9ZMc7i3FP

Background

• Instrument scientists suggested resolution terms per event (each event has its own pixel and wavelength), and want the 9 LoKI banks merged before normalization. Merging is linear, so it can happen at any point before the final division, e.g. in the (Q, λ) binning step.
• with_banks in ess/sans/workflow.py says banks are not merged because their resolution differs. The approach below removes that reason.

Resolution of one element

An element i is one (pixel, wavelength) combination. Mildner–Carpenter, keeping cos θ because the LoKI side banks reach 2θ ≈ 20° and more:

σ_i² = (2π cosθ / λ)² · [3(R1/L1)² + 3(R2/L3)² + ΔR_ang²] / 12
     + Q_i² · σ_λ² / λ²

1/L3 = 1/L1 + 1/L2
σ_λ² = σ_source(Ltotal, λ)² + Δλ²/12

ΔR_ang is the pixel size projected onto the scattering angle; for a flat detector perpendicular to the beam it is ΔR cos²(2θ) / L2. σ_source replaces the ISIS moderator time-spread file used in scipp/esssans#185: the essreduce wavelength lookup table already stores the wavelength variance per (distance, event_time_offset) (ess/reduce/unwrap/lut.py). It is currently only used to mask uncertain regions.

Combining elements into one σ per Q bin

I(Q) is a ratio of sums over the elements in a bin. N_i is the existing normalization term, i.e., the expected counts per unit cross section:

# E[x]: expectation value of x (average over counting statistics)
N_i      = Ω_pixel · D(λ) · M(λ) · T(λ)
E[C_i]   = N_i · ∫ R_i(Q') I(Q') dQ'          # R_i = Gaussian(Q_i, σ_i)
I_bin    = Σ C_i / Σ N_i
E[I_bin] = ∫ K_bin(Q') I(Q') dQ'   with   K_bin = Σ N_i R_i / Σ N_i

The resolution function of a bin is therefore a mixture of the element Gaussians weighted by N. This follows from the definition of I_bin; no approximation is involved. Its second moment gives σ:

S0 = Σ N_i                 # existing denominator
S1 = Σ N_i Q_i
S2 = Σ N_i (σ_i² + Q_i²)

Q_mean = S1 / S0
σ_bin² = S2 / S0 − Q_mean²   # = mean(σ_i²) + var(Q_i); the second term is the bin-width contribution

All of these are sums, so banks, runs and chunks merge by addition, exactly like numerator and denominator.

Alternatives

• Per-event weighting: the same formula with w_i = C_i (events) instead of N_i. Since C_i ≈ N_i I(Q_i), this gives more weight to the parts of a bin where I(Q) is large. It is noisy in bins with few events, gives different σ for sample and background, and needs one more float per event plus σ_λ looked up per event. Its one advantage: σ kept on each event allows re-binning in Q later.
• Mantid Q1D (Q1D2.cpp): σ_bin = Σ C_i √(σ_i² + ΔQ²/12) / Σ C_i. This averages σ instead of σ², so it comes out narrower than the mixture.

A toy model (two flat banks with L2 = 10 m and 3 m, a sphere with R = 80 Å, the ESS long pulse, 80 logarithmic Q bins) gives:

The differences are small. The argument for N-weighting rests on it being exact, deterministic, the same for sample and background, and cheaper.

Keeping the shape of the resolution function

The mixture is not Gaussian: it has a sharper peak and wider tails where short and long wavelengths contribute to the same bin (excess kurtosis ≈ 1 around Q ≈ 0.01 Å⁻¹ in the toy model). An approximation that also works with sums: group the elements of a bin into classes with fixed edges, and replace each class by one Gaussian with the same mean and variance:

class c = (Q sub-bin, σ/Q class)
# σ/Q class k: elements with f^k ≤ σ_i/Q_i < f^(k+1), f = class width factor
#   (f = 1.5 gives class edges at 1.7%, 2.6%, 3.9%, 5.9%, ...)
# Q sub-bin j: the Q bin split into n equal parts, element assigned by Q_i
W_c, S1_c, S2_c = same sums as above, restricted to class c
K_bin(Q') ≈ Σ_c (W_c / Σ W) · Gaussian(Q'; μ_c = S1_c / W_c, s_c² = S2_c / W_c − μ_c²)

Largest deviation from the exact mixture, relative to its peak (toy model):

Around 10–20 components per bin keep the error at or below about 1% at every bin checked.

NXcanSAS only standardizes Qdev (or dQw/dQl for slits); extra datasets or an NXnote are allowed but nothing reads them. sasmodels applies resolution as a weight matrix (Resolution.weight_matrix), so a custom Resolution subclass could use the mixture when fitting through the Python API.

Proposed implementation

1. Merge banks: map over NeXusDetectorName and reduce at NormalizedQ[RunType, IofQPart] by summation, as with_sample_runs does for runs. Can be a separate PR, before or after the resolution work (see below).
2. σ_source from the lookup table: an essreduce helper that returns σ_λ(Ltotal, λ) from the lookup table variance.
3. Resolution sums: compute σ_i² on the denominator grid (pixel × λ, after compute_Q), histogram S1 and S2 over Q like the denominator, then apply the same monitor-term multiplication, wavelength mask, band reduction and bank/run merging. Suggested structure: make the sums a third member of IofQPart,

   IofQPart = TypeVar('IofQPart', Numerator, Denominator, ResolutionMoments)

   carrying N·Q and N·(σ² + Q²) along a size-2 moment dimension. The generic providers (compute_Q, bin_in_q, reduce_q) then apply unchanged. New code is needed only for the provider that creates the sums and for a variant of mask_and_scale_wavelength_q that multiplies by the monitor term.
4. Output: set σ_bin² as variances of the Q coordinate. save_background_subtracted_iofq already writes these as resolutions via SASdata(..., Q_variances='resolutions'). Use the sample run's resolution for background-subtracted I(Q).
5. Optional: the mixture components as an extra output and NXcanSAS dataset.

Independence from bank merging

The resolution work does not depend on bank merging (step 1); for a single bank it gives the correct σ. Two rules keep the two compatible:

• Compute σ² = S2/S0 − (S1/S0)² only once, after all merging. Never compute σ per bank or per run and average it afterwards.
• Every merge point must also merge the resolution sums. Today this is only _set_runs (multiple sample or background runs), whose for part in (Numerator, Denominator) loop must include ResolutionMoments. The instrument scientists will probably process runs separately and merge the files instead, but the workflow supports multiple runs and must stay correct. A bank merge written as the same loop over parts includes the resolution without further changes.

Test for both rules: σ from two runs merged by the workflow equals σ computed from the summed S0, S1 and S2 of the individual runs.

Open questions

• Where do R1, R2, L1 and ΔR come from for LoKI: NeXus or user parameters?
• Report I(Q) at bin centers or at Q_mean? σ changes slightly with the choice of reference point.
• ΔR for straws at an angle to the scattered beam: use the radial extent of each pixel from the pixel shape?
• Is gravity relevant for the resolution at LoKI?
• Do the instrument scientists need σ per event (for re-binning event-mode I(Q)), or is σ per Q bin enough?