List-mode per-shot spectra: rotation + gap correction directly on droplet sparse arrays (no full-image reconstruction)

[lg345]

Summary

Add a registered XSpect step that produces per-shot 1D spectra directly from the droplet2photon sparse photon list, without reconstructing the full detector image. Rotation, ROI selection, dispersion/energy projection, and detector-gap correction are all done as geometry operations on the photon coordinates ("list-mode" / "event-mode" processing).

The full-image droplet_reconstruction step (XSpect analysis/droplet.py) is currently the throughput bottleneck. But dense reconstruction is only needed as an intermediate because rotation and gap-patching were written as image operations. Both are point/geometry operations that act more naturally — and more accurately — on the sparse list. Since the per-shot end product is a 1D spectrum (sum rows in an ROI, histogram along dispersion), the ~540k-pixel array is never actually required.

Motivation

• Speed: per-shot cost drops from "allocate + scatter + spline-rotate a 540k-pixel image" to "transform + histogram the handful of photons in that shot" (~100–1000× fewer ops/shot; no image allocation, no interpolation).
• Correctness / artifact removal: the polynomial gap-patch injects dead columns as linear combinations of their neighbors, which spook cannot assign independent incident-energy weight to — this produces the vertical RIXS-map stripes documented in lg345/polarization_dependent_stochastic#28. Geometric exposure normalization injects no linearly-dependent data and removes that artifact.

Method

For a shot's photons (r_i, c_i) with weights w_i (droplet photon count):

1. Rotation = exact linear transform on coordinates (O(N_photons), no interpolation):

d_i = (c_i - c0)*cosθ - (r_i - r0)*sinθ      # dispersion axis
s_i = (c_i - c0)*sinθ + (r_i - r0)*cosθ      # cross-dispersion axis

(θ sign chosen to match the existing scipy.ndimage.rotate convention.)

2. ROI + projection = filter then histogram:

m        = (s_lo <= s_i <= s_hi)
spec_raw = histogram(d_i[m], bins=dispersion_bins, weights=w_i[m])

3. Gaps = static coverage/exposure normalization (computed ONCE):

rr, cc   = nonzero(live_pixel_map within ROI)      # dead columns excluded
d_live   = (cc-c0)*cosθ - (rr-r0)*sinθ
s_live   = (cc-c0)*sinθ + (rr-r0)*cosθ
coverage = histogram(d_live[s_lo<=s_live<=s_hi], bins=dispersion_bins)
# per shot:
spec = spec_raw / coverage * coverage.mean()

Because rotation maps a dead column to a diagonal set of output bins, each output bin loses only ~1/N_rows of coverage per crossing — coverage never hits zero away from θ=0, so the correction is smooth and well-conditioned. Guard the division with a coverage floor (mask bins below ~20% of median coverage).

Optional: fold the pixel→eV map into the histogram bin edges to merge reconstruction + rotation + gap correction + energy calibration into a single pass over the photon list.

What still needs a dense image (one-time only, not the hot loop)

• Rotation-angle finding (mfx101609126_rotation_diagnostic) and gap-column / live-pixel identification (mfx101609126_pixel_patch_diagnostic) — run on the sum image over all shots, built once.

Proposed API

New registered step, e.g. sparse_project_spectrum:

• inputs: droplet_droplet2phot_sparse_* arrays, rotation angle, pivot (r0,c0), ROI cross-dispersion band [s_lo,s_hi], dispersion bins (or energy bins + calibration), live-pixel map (or explicit dead-column list)
• output: per-shot 1D spectrum (N_shots, N_bins), already gap-corrected — a drop-in replacement for droplet_reconstruction → patch_pixels → rotate_detector → reduce_detector_spatial.

Acceptance / validation

• Benchmark vs the current droplet_reconstruction pipeline on mfx101609126 run 79 (Fe foil); report speedup and peak memory.
• Compare mean per-shot spectra list-mode vs dense path (expect close agreement, NOT bit-identical — coordinate rotation ≠ image-spline rotation, exposure norm ≠ polynomial patch).
• Confirm the SEER RIXS-map gap stripes (polarization_dependent_stochastic#28) are gone / reduced with exposure normalization.
• Compare spook reconstruction residuals list-mode vs dense on Fe foil (78-82) and ferricyanide (91-95).

Caveats

• Carry the droplet photon count as w_i (a droplet may be >1 photon).
• This is a new, arguably more correct method, not a bit-match to the current path — revalidate rather than expecting corr=1.0.
• Coverage floor needed only near θ=0 or with a very narrow ROI.

Cross-references

• Downstream analysis / artifact context: lg345/polarization_dependent_stochastic#26 (Experiment mfx101609126), Revert "Merging master into documentation branch" #28 (pixel-gap patching + RIXS stripe artifact), Merging code changes from Dehe for Follmer experiment #29 (rotation angles).
• Tracked in both the XSpect project and the Polarization-Dependent XANES project.

[lg345]

Prototype results (run 79, epix100_0, ROI [270,330,400,700], −2.0°)

Prototyped list-mode projection and validated against the dense path (droplet_reconstruction → rotate_detector → reduce_detector_spatial). Scripts: mfx101609126/prototype_listmode_projection.py and prototype_listmode_gap_coverage.py in the pol_xanes repo.

Accuracy — confirmed

• Mean-spectrum correlation: 0.9996 (list-mode vs dense) with rotation sign=+1 (matches scipy.ndimage.rotate reshape=False about the array centre).
• Per-shot median correlation 0.993. Not bit-identical (coordinate rotation vs image-spline rotation), as expected.

Speed — confirmed

• End-to-end (incl. shared HDF5 read): 4.0×.
• Compute-only (photons preloaded): 7.4× (0.30 ms/shot vs 2.20 ms/shot). Pure Python + numpy histogram; a vectorized/batched implementation should widen this further, and memory per shot drops from a (60,300) image to just the photon list.

Gap correction — works, but my geometry claim above was WRONG and needs correcting

The exposure/coverage normalization is implemented and produces a smooth correction where coverage is nonzero. But at the actual geometry the coverage still hits exactly zero at the gap centers, contrary to my "coverage never hits zero away from θ=0" claim in the issue body.

Why: the diagonal smear of a dead column across the summed row band is tan(θ) × N_band_rows. Here θ=2° and the row band is 30 rows → smear ≈ 1.05 columns. The dead-column groups are 3 columns wide (e.g. 80,81,82). Since 1.05 < 3, the center column of each gap group receives no live-pixel coverage → coverage=0 at cols 80, 177, 273 (6 zero bins total), and list-mode returns 0 there.

So exposure normalization does not by itself eliminate a wide multi-column gap at a shallow angle. Options:

1. Accept small gaps (mask the ≤1–2 residual zero bins; interpolate only those, not the whole column) — still avoids the per-shot polynomial patch and drastically reduces the number of synthesized bins vs the current approach.
2. Sum over the full ROI height (60 rows) instead of the 10–40 band where signal permits: smear ≈ 2.1 cols — helps but still < 3.
3. Larger effective angle / different band only helps if physically justified (θ≥5.7° would fully cover a 3-col gap across 30 rows — not our geometry).
4. Bin finer than 1 px along dispersion so partial coverage is captured, then rebin.

Honest revised conclusion

• The speed and accuracy wins are real (7.4× compute, corr 0.9996) — rotation + projection in list mode is clearly worth doing.
• The gap story is more nuanced than I claimed: exposure normalization cleanly handles the fractional coverage deficit on gap edges, but a fully-dead multi-column gap at a shallow rotation still has zero-coverage bins that need a minimal local interpolation. The advantage over the current polynomial patch is that we synthesize far fewer bins (the true zero-coverage centers only, ~1 col/gap) instead of patching whole columns per shot — which should still reduce the spook RIXS-stripe artifact (pol_xanes#28), but likely not eliminate it entirely. To be verified by running spook on list-mode output.

Next steps

• Vectorize/batch the list-mode kernel; re-benchmark at full run scale.
• Decide gap policy (option 1 vs 2) and quantify residual zero-coverage bins per detector.
• Run spook on list-mode spectra (Fe foil + ferricyanide) and compare RIXS-map stripes vs the polynomial-patch path.

[lg345]

Performance Improvement Report

Committed copy: lg345/polarization_dependent_stochastic mfx101609126/PERFORMANCE_REPORT.md @ a84306a. Reproducible via the three prototype scripts in that dir.

Issue: #96
Date: 2026-07-03
Data: mfx101609126 run 79 (Fe foil), detector epix100_0
Geometry: ROI [270,330,400,700] → (60,300); rotation −2.0°; row band 10–40; 300 dispersion bins
Environment: ana-4.0.68-py3, single core, NumPy/SciPy, compute-only timing (photons preloaded)

1. Summary

Replacing the dense path (droplet_reconstruction → rotate_detector → reduce_detector_spatial)
with a list-mode projection — rotating photon coordinates and histogramming them directly
from the droplet sparse arrays — gives a 4.2× compute speedup at 0.9997 mean-spectrum
correlation, with no full-detector image ever allocated. Accuracy and speed wins are solid;
the detector-gap handling is improved but not fully solved by geometry alone (see §4).

2. Benchmark results

Best of 3 repeats, 2000 shots, compute-only (HDF5 read excluded):

Method	ms/shot	speedup	mean-spec corr vs dense
dense (reconstruct + rotate + ROI-sum)	2.563	1.0×	—
list-mode (per-shot loop)	0.614	4.2×	0.9997
list-mode (vectorized, flatten-all)	1.817	1.4×	0.9997

Earlier end-to-end run (400 shots, incl. shared HDF5 read) showed 4.0×; a smaller
compute-only slice showed 7.4×. The 4.2× above is the most representative figure at scale.

3. Why list-mode wins (and why the "obvious" vectorization loses)

The dominant cost in the dense path is allocating and spline-rotating a (60,300) image
per shot, regardless of how few photons it contains. List-mode cost scales with the
number of photons in the ROI, which here is tiny:

• padded slots/shot: 100,000
• live photons/shot: ~74,000
• in-ROI photons/shot: ~72 (0.1% of live)

So the per-shot loop is fast precisely because it filters to the ROI first, then does O(72)
work. The "vectorized" version that flattens all shots' photons before filtering is slower
(1.4×) because it materializes and processes all ~74k×2000 padded photons before discarding
99.9% of them — defeating the early ROI filter and thrashing memory. Lesson: filter to the
ROI as early as possible; do not flatten-then-filter on heavily zero-padded fixed-length arrays.

4. Detector-gap handling (corrected analysis)

The static exposure/coverage normalization (bin all live pixels through the same rotation+ROI
once, divide each shot's spectrum by that coverage) works where coverage is nonzero, but the
original "coverage never hits zero away from θ=0" claim was wrong at this geometry:

• diagonal smear of a dead column across the row band = tan(θ) × N_band_rows = tan(2°) × 30 ≈ 1.05 col
• dead-column groups are 3 columns wide (e.g. 80,81,82)
• since 1.05 < 3, gap centers get zero coverage → 6 zero-coverage bins (cols 80,177,273)

Consequence: exposure normalization cleanly handles the fractional deficit on gap edges,
but a fully-dead multi-column gap at a shallow angle still needs a minimal local interpolation
of the residual zero-coverage bins (~1 col/gap). This is still far better than the current
per-shot polynomial patch, which synthesizes entire columns as linear combinations of
neighbors — the root cause of the spook RIXS-map stripe artifact (pol_xanes#28). List-mode
synthesizes ~1 bin/gap instead of whole columns, so the artifact should shrink substantially,
but likely not vanish. To be confirmed by running spook on list-mode output.

5. Memory

	per-shot memory
dense (60,300) float64 image	140.6 KB
list-mode ROI photon lists (~72 photons × 3 arrays)	~1.7 KB

Additional wins not yet exploited: reading only the ROI photons from HDF5, and folding the
pixel→eV calibration into the histogram bin edges (merging reconstruction + rotation + gap
correction + energy calibration into one pass).

6. Recommendations

1. Implement list-mode as the default projection for the spook per-shot pipeline
   (register sparse_project_spectrum as proposed in List-mode per-shot spectra: rotation + gap correction directly on droplet sparse arrays (no full-image reconstruction) #96). Expect ~4× at this geometry,
   more when ROI photon fraction is low.
2. Filter to the ROI before any per-shot work; keep the per-shot loop over the ROI-filtered
   photon list. Avoid the flatten-all vectorization on padded fixed-length arrays.
3. Gap policy: use coverage normalization + minimal local interpolation of residual
   zero-coverage bins only (not whole columns). Quantify residual zero-coverage bins per
   detector/geometry.
4. Validate downstream: run spook on list-mode spectra for Fe foil (78–82) and ferricyanide
   (91–95); compare RIXS-map stripes and residuals against the polynomial-patch path.
5. Not bit-identical: coordinate rotation ≠ image-spline rotation and exposure-norm ≠
   polynomial-patch; treat this as a new, more principled method and revalidate rather than
   expecting the reconstruction step's corr=1.0 bit-match.

7. Reproduce

cd mfx101609126
python benchmark_listmode.py              # table in §2
python prototype_listmode_projection.py   # accuracy + sign calibration
python prototype_listmode_gap_coverage.py # gap coverage + compute-only speedup