[ESSSANS] Q resolution per I(Q) bin - #766

Part of #763 (steps 2 and 3 of the proposed implementation). Refines #360.

Adds QResolution[RunType]: one σ per I(Q) bin, the NXcanSAS Qdev. The resolution is built from sums, so it stays correct when runs are merged, and later when banks and live-data chunks are merged.

For science review

What is computed

Each element i (one pixel and one wavelength bin) gets a Gaussian resolution (Mildner–Carpenter, keeping cos θ and cos 2θ):

σ_i² = (2π cosθ / λ)² · [R1²/(4 L1²) + R2²/4 · (1/L1 + cos2θ/L2)² + σ_pixel²]
     + Q_i² · σ_source(Ltotal, λ)² / λ²

σ_pixel² = [ℓ²/12 · c² + r²/4 · (1 − c²)] / L2²     # cylindrical pixel, radius r, length ℓ
c = (pixel axis) · e,   e = cos2θ ρ̂ − sin2θ k̂       # direction of increasing 2θ

σ_pixel is computed per pixel from pixel_shape and the detector transformation. The extent that matters is the one along e, which depends on the straw orientation relative to e in every bank: on the rear detector, pixels left and right of the beam have their straw along e (length ℓ ≈ 2 mm counts), pixels above and below have it perpendicular (diameter 7.5 mm counts). A scalar pixel size cannot describe this. For a flat detector perpendicular to the beam with pixel size ΔR along e it reduces to ΔR cos2θ / (√12 L2). The cos2θ in the sample-aperture term is the same foreshortening, applied to the sample aperture as seen from the pixel.

The wavelength bin width does not enter. Events are histogrammed in Q with their own wavelength, so the spread of Q in a bin is limited by the Q bin width, which var(Q_i) below already contains. A Δλ²/12 term counts it twice: on the Larmor data it made σ/Q at Q = 0.27 Å⁻¹ 4.9% with 50 wavelength bins and 3.2% with 400. Without it, σ is the same for both to four digits.

The resolution function of a Q bin is the mixture of the element Gaussians, weighted by the normalization term N_i (the existing denominator). This follows from I_bin = Σ C_i / Σ N_i. We report its standard deviation about its mean:

S0 = Σ N_i,   S1 = Σ N_i Q_i,   S2 = Σ N_i (σ_i² + Q_i²)     # over elements that reach the numerator
Q_mean = S1/S0
σ_bin² = S2/S0 − Q_mean²        # = mean(σ_i²) + var(Q_i)

Elements whose events get no wavelength, because the lookup table is masked there (relative spread above LookupTableRelativeErrorThreshold), are part of the denominator but not of the counts, so they are excluded from the sums. The mixture is exact if counts are binned like the denominator. Events are binned by their own wavelength; an event-level Monte-Carlo on the Larmor data agrees with σ_bin within 1% (σ_λ/λ = 2%) and within 2–3% (σ_λ = 0.12 Å, where linearizing in 1/λ starts to matter at high Q).

QResolution holds σ_bin as data and Q_mean as a coordinate. See #763 for the derivation and for the comparison with per-event weighting and Mantid Q1D.

Please check

• Source wavelength spread. σ_source comes from the variance stored in the essreduce wavelength lookup table, at the pixel's Ltotal. Tables built by simulation (tof, McStas ESS source with correlated emission time and wavelength, propagated through the choppers) store the right quantity. Tables built in analytical mode, which is what esslivedata uses, store half the wavelength range for a flat 0–5 ms source window. Without pulse-shaping choppers this is 2.8 times the true σ (0.345 Å instead of 0.12 Å at 28.7 m). The same flat window also puts the mean emission time at 2.5 ms instead of 1.65–1.85 ms, so analytical tables assign wavelengths that are 0.085–0.12 Å too short at 28.7 m (5.4% at 2 Å). That affects I(Q) itself and is not addressed here. See [essreduce] Analytical wavelength LUT: flat source window biases mean wavelength and spread #770 for a proposed fix (weighting the analytical table by the pulse profile), which would also give the right σ_source.
• Instrument parameters. R1, R2 and L1 are user parameters for now (SourceApertureRadius, SampleApertureRadius, CollimationLength). Should they come from NeXus? Are the LoKI apertures circular?
• Pixel size. LoKI pixels are 3.75 mm radius and about 2 mm long along the straw, per NeXus. The real position resolution along a straw (charge division) may be coarser. PixelScatteringAngleVariance[RunType] can be set directly to override the value from the pixel shape.
• Reference point. σ is measured about Q_mean. If I(Q) is reported at the bin centre, the second moment about the centre is σ_bin² + (Q_mean − Q_centre)². On the Larmor test data, Q_mean is within 5% of a bin width of the centre.
• Gravity. If CorrectForGravity is set, the gravity-corrected 2θ is used in σ_i. Gravity is not included as a separate resolution term.
• Values not yet validated on LoKI. The LoKI test NeXus file gives an all-NaN denominator, and its choppers have zero delay, so they transmit nothing. The only lookup table in the test data has no choppers and frame overlap.

Not in this PR

• Merging detector banks.
• Writing Qdev to NXcanSAS: save_background_subtracted_iofq is unchanged.
• Qx/Qy resolution.
• A better description of the non-Gaussian shape of the resolution function (follow-up, possibly via scippuncertainty bootstrap).

For developer review

• Three new IofQPart members. ResolutionZerothMoment (N), ResolutionFirstMoment (N·Q) and ResolutionSecondMoment (N·(σ² + Q²)) go through the existing generic bin_in_q, reduce_q and run merging, exactly like numerator and denominator. The [ESSSANS] Q resolution per Q bin, compatible with merging LoKI banks #763 plan used a single part with a size-2 moment dim. That does not work, because the sums have different units. The zeroth moment is needed instead of the denominator because excluded elements are in the denominator.
  • σ_i² is computed once, as DetectorQVariance[RunType] on the QDetector[RunType, Denominator] grid. All moments are derived from it. Excluded elements have NaN σ_i² and get weight 0.
  • Only the values of N are used. Denominator variances are dropped.
• Pixel term. PixelScatteringAngleVariance[RunType] is computed from DetectorPixelShape[RunType] and the detector transformation, on the pixel dims only, and broadcast over wavelength. solid_angle and this provider share normalization.pixel_cylinder. Detectors without a pixel shape can set the variance directly.
• Monitor term. mask_and_scale_resolution_moment is the counterpart of mask_and_scale_wavelength_q. It multiplies by the values of the monitor term, after binning in Q.
• Merging runs. _set_runs now also merges NormalizedQ of the three moments. q_resolution combines the sums only after all merging. Never compute σ per run and average it.
• Separate output. QResolution[RunType] is a separate output rather than variances on the Q coord of I(Q). I(Q) has bin-edge Q, and requiring the resolution in IntensityQ or in the save function would break the ISIS workflows, which have no lookup table.
• Source spread input. SourceWavelengthSpread[RunType] is where σ_source enters. By default it comes from LookupTable[RunType, NXdetector] and is NaN where the relative spread exceeds the table's error threshold, which marks the elements to exclude. Workflows without a lookup table (ISIS, LoKI at Larmor) can set it directly as a scalar or as anything broadcastable to (pixel, wavelength).
• Event data required. detector_q_variance raises if the numerator is not binned (event) data, since histogrammed data would need the wavelength-bin term.
• essreduce. New ess.reduce.unwrap.wavelength_spread(table, ltotal, wavelength):
  • For each distance row, it interpolates the stddev against that row's mean wavelength, then interpolates linearly in distance.
  • Non-finite entries are ignored, and wavelengths outside a row's range use the edge value.
  • Ltotal outside the table's distance range raises an error.
• Memory. Peak memory grows by a few arrays the size of the dense (pixel, wavelength) denominator.

Test plan

• Unit test: σ_i² against an independent numpy implementation of the formula (mixed units, two pixels).
• Unit tests: moments exclude elements with NaN variance; dense numerator raises.
• LoKI: SourceWavelengthSpread is NaN exactly where the relative spread exceeds the error threshold.
• Unit test: per-pixel 2θ variance against a numerical average over the cylinder volume, for four pixel orientations and a rotated detector transformation (agreement 1e-4).
• Larmor: QResolution equals a brute-force N-weighted second moment over the dense (pixel, wavelength) arrays, including pixel masks (rtol 1e-9).
• Larmor: σ from two runs merged by with_sample_runs equals σ from the summed per-run sums.
• Larmor: QResolution has the same dims and Q coord as I(Q), with wavelength bands and with DimsToKeep.
• LoKI: SourceWavelengthSpread from the test lookup table is finite and positive, and QResolution can be computed.
• wavelength_spread: exact for spreads linear in distance and wavelength, including unsorted (wrapped) table rows, unit conversion, NaN entries, and out-of-range errors.

🤖 Generated with Claude Code