From e7fbeb527f4f19e7cf1364c48ab7760f7db01e45 Mon Sep 17 00:00:00 2001
From: leonarski_f <filip.leonarski@psi.ch>
Date: Sat, 11 Jul 2026 10:59:15 +0200
Subject: [PATCH] docs: update analysis references to match current code

Reconcile docs/CPU_DATA_ANALYSIS.md and docs/RUGNUX.md with the current
image_analysis/ and rugnux CLI (verified section-by-section against the code):
unified profile-fit Bragg integration engine, multi-lattice indexing,
azimuthal phi binning, radial parallax/bandwidth profile with sub-pixel
centring, rot3d capture-fraction handling, automatic CC1/2 resolution cutoff,
and the new rugnux options; fix the section numbering and cross-references;
remove the never-implemented French-Wilson and still-partiality descriptions.

Delete the stale in-source design notes (ICE_RING_DETECTION,
BRAGG_INTEGRATION_ENGINE, NEXTGEN_INTEGRATOR) and fix the two code comments
that pointed at them; the BraggIntegrationEngine header no longer claims it is
'not yet wired'.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
---
 common/BraggIntegrationSettings.h             |   2 +-
 docs/CHANGELOG.md                             |   1 +
 docs/CPU_DATA_ANALYSIS.md                     | 147 +++----
 docs/RUGNUX.md                                |  30 +-
 image_analysis/azint/ICE_RING_DETECTION.md    |  71 ----
 .../BRAGG_INTEGRATION_ENGINE.md               | 107 -----
 .../BraggIntegrationEngine.h                  |   7 +-
 .../bragg_integration/NEXTGEN_INTEGRATOR.md   | 391 ------------------
 8 files changed, 97 insertions(+), 659 deletions(-)
 delete mode 100644 image_analysis/azint/ICE_RING_DETECTION.md
 delete mode 100644 image_analysis/bragg_integration/BRAGG_INTEGRATION_ENGINE.md
 delete mode 100644 image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md

diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h
index 3c3b724c..bd53b7db 100644
--- a/common/BraggIntegrationSettings.h
+++ b/common/BraggIntegrationSettings.h
@@ -9,7 +9,7 @@
 // profile-fits with a measured-width Gaussian (Kabsch-style) - more accurate intensities than the
 // classical uniform BoxSum; validated on HEWL anomalous (stronger S/Cl peaks vs box-sum). BoxSum is
 // the simpler, faster fallback. ProfileEmpirical learns the profile per resolution shell from strong
-// spots - see bragg_integration/BRAGG_INTEGRATION_ENGINE.md.
+// spots - see docs/CPU_DATA_ANALYSIS.md (Bragg integration).
 enum class IntegratorMode { BoxSum, ProfileGaussian, ProfileEmpirical };
 
 class BraggIntegrationSettings {
diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 6aaed72a..177887fe 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -6,6 +6,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
 * jfjoch_viewer: Open datasets on the WSL2/UNC filesystem (paths starting `\\`); write processing outputs next to the input file, with a Browse button and independent `_process.h5` / merged `.mtz`/`.cif` toggles; and show the determined space group in the merge-statistics window.
 * rugnux: Accept an absolute `-o` output prefix in offline processing.
 * Packaging: The self-contained Linux viewer `.tgz` now bundles cuFFT, so it runs without a system CUDA toolkit (`.deb`/`.rpm` are unchanged, distro-managed).
+* Docs: Bring the analysis references up to date with the code. `docs/CPU_DATA_ANALYSIS.md` now reflects the unified profile-fit Bragg integration engine, multi-lattice indexing, azimuthal phi binning, the radial parallax/bandwidth profile with sub-pixel centring, the rot3d capture-fraction handling and the automatic CC1/2 resolution cutoff, and drops the descriptions of features that were never implemented (French-Wilson amplitudes, the still excitation-error partiality model); `docs/RUGNUX.md` documents the new `--resolution-cutoff`/`--resolution-cc-target`/`--resolution-shells`, `--min-captured-fraction`, `--mosaicity`, `--reference-column`, the azimuthal correction toggles and the geometry-override options, and corrects the `-N` default. The outdated in-source design notes (ICE_RING_DETECTION, BRAGG_INTEGRATION_ENGINE, NEXTGEN_INTEGRATOR) are removed.
 
 ### 1.0.0-rc.157
 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use.
diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md
index 27ba5295..a3343dde 100644
--- a/docs/CPU_DATA_ANALYSIS.md
+++ b/docs/CPU_DATA_ANALYSIS.md
@@ -11,7 +11,7 @@ This document describes the crystallographic algorithms implemented in Jungfrauj
 5. Bravais lattice / centering inference,
 6. geometry and lattice refinement,
 7. reflection prediction (still and rotation),
-8. Bragg integration by either 2D summation or reference-driven profile fitting,
+8. Bragg integration by either 2D box summation or profile fitting (Kabsch, reference-free),
 9. scaling and merging,
 10. merge-level error modelling and outlier rejection,
 11. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref).
@@ -87,29 +87,30 @@ $
 $
 Jungfraujoch uses $|\Delta_\mathrm{Ewald}|$ as an operational proxy for excitation error. This appears in:
 - still prediction (accept if $|\Delta_\mathrm{Ewald}|\le \Delta_\mathrm{cut}$),
-- profile radius estimation (see §7.1),
-- still partiality option in scaling/merging (§9.3).
+- profile radius estimation (see §11.1),
+- still partiality option in scaling/merging (§10.2).
 
 ---
 
 ## 2. Azimuthal integration (radial profiles)
 
-Azimuthal integration produces a 1D radial profile $I(q)$ or $I(d)$ by histogramming pixels into radial bins. Pixels are **not split** across bins; each pixel contributes wholly to a single bin.
+Azimuthal integration produces a radial profile $I(q)$ or $I(d)$ by histogramming pixels into radial bins. Pixels are **not split** across bins; each pixel contributes wholly to a single bin. By default the profile is purely radial (a single azimuthal bin), but the azimuth can optionally be split into up to 512 $\phi$ sectors (`azim_bins`, `--azim-phi-bins`), giving a **2D $q\times\phi$ profile** that exposes azimuthal anisotropy such as detector shadowing or sample texture.
 
 ### 2.1 Histogram estimator
 
-Let bin index $b(x,y)\in\{0,\dots,B-1\}$ be precomputed from $q(x,y)$ (or equivalently from $d(x,y)$). For each bin $b$:
+Let bin index $b(x,y)$ be precomputed from $q(x,y)$ (or equivalently from $d(x,y)$) and, when $\phi$ sectors are enabled, the azimuth $\phi(x,y)$ — so $b = b_q + b_\phi B_q$. For each bin $b$:
 
-- accumulate corrected intensity:
+- accumulate corrected intensity and its square:
   $
-  S_b = \sum_{(x,y):\,b(x,y)=b} I(x,y)\,C(x,y),
+  S_b = \sum_{(x,y):\,b(x,y)=b} I(x,y)\,C(x,y),\qquad
+  S^{(2)}_b = \sum I(x,y)^2\,C(x,y)^2,
   $
 - and count:
   $
   N_b = \#\{(x,y):\,b(x,y)=b \text{ and pixel is valid}\}.
   $
 
-A simple mean profile is then $ \bar{I}_b = S_b / N_b$ (when $N_b>0$). Invalid pixels (masked, saturated, detector error codes) are excluded.
+The profile reports both the mean $\bar{I}_b = S_b / N_b$ (when $N_b>0$) and a per-bin sample standard deviation $\sigma_b = \sqrt{(S^{(2)}_b - S_b^2/N_b)/(N_b-1)}$ (a spread/error estimate for each radial point). Invalid pixels (masked, saturated, detector error codes) are excluded.
 
 ### 2.2 Corrections applied
 
@@ -117,8 +118,9 @@ Two standard corrections are available:
 
 **(i) Solid angle / geometric correction.** A commonly used approximation for flat detectors gives a $\cos^3(2\theta)$ factor:
 $
-C_\Omega(2\theta) = \cos^3(2\theta).
+C_\Omega(2\theta) = \cos^3(2\theta),
 $
+applied — like the polarization term below — as a **divisor** (intensities are scaled by $1/\cos^3(2\theta)$), so pixels at higher $2\theta$, which subtend a smaller solid angle, are boosted.
 
 **(ii) Polarization correction.** With polarization coefficient $P$ (beamline dependent) and azimuth $\phi$:
 $
@@ -129,7 +131,7 @@ applied as a divisor to intensities (i.e. scale by $1/C_\mathrm{pol}$) when enab
 
 ### 2.3 Background estimate for profiles
 
-A background estimate is derived from the integrated profile using the azimuthal integration settings (details depend on the configured estimator). This background is used for monitoring and diagnostics; it is **not** the same as local Bragg-spot background used in summation integration (§8).
+A background estimate is derived from the profile as its mean intensity over a fixed low-to-mid $Q$ window (default $2\pi/5$ to $2\pi/3$ Å$^{-1}$). This background is used for monitoring and diagnostics; it is **not** the same as the local Bragg-spot background used in summation integration (§9.2).
 
 ---
 
@@ -163,6 +165,7 @@ $
 
 A pixel is considered strong if:
 - it is above a photon/count threshold, and
+- its window contains enough valid neighbours (more than 100), so the local statistics are meaningful, and
 - $\Delta>0$, and
 - the squared deviation exceeds a scaled variance:
   $
@@ -178,7 +181,7 @@ Special cases:
 
 ### 3.2 Resolution and ice-ring handling
 
-Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, pixels in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6).
+Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, spots in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6).
 
 A single per-image **ice-ring score** is derived from the azimuthally-integrated radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins so the estimate is robust to the radial binning: at a coarse Q-spacing a local shoulder can be only ~1 bin and would double-count the ring's own edge (offline processing defaults to a fine 0.01 1/Å spacing, `--azim-q-spacing`, so the rings are well resolved). (A significance/z-score was considered but is uninformative here: with many photons any real ice ring is highly significant, so the discriminating quantity is the ice *magnitude*, i.e. this ratio.) It is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity, distinct from the merge-time ice masking, which is data-driven from the per-ring merged CC1/2.
 
@@ -211,7 +214,7 @@ Jungfraujoch supports two complementary indexing strategies:
 1. **FFT-based indexing** (Rossmann-type): does not require an a priori unit cell; suitable for unknown samples.
 2. **Fast-feedback indexing** (TORO-like): requires an approximate unit cell; optimized for speed and feedback.
 
-Both feed into a common robust refinement/selection stage which maximizes the number of inliers under an indexing tolerance.
+Both feed into a common robust refinement/selection stage which maximizes the number of inliers under an indexing tolerance, and which can return **more than one lattice** per image (multi-lattice indexing; see §5.4).
 
 ### 4.1 Indexed-spot decision (inlier test)
 
@@ -225,7 +228,7 @@ Let $(h,k,l)=(\mathrm{round}(h_f),\mathrm{round}(k_f),\mathrm{round}(l_f))$ and
 $
 \delta^2 = (h_f-h)^2 + (k_f-k)^2 + (l_f-l)^2.
 $
-A spot is indexed if $\delta^2 \le \tau^2$, where $\tau$ is the configured tolerance.
+A spot is indexed if $\delta^2 < \tau^2$, where $\tau$ is the configured tolerance.
 
 For indexed spots, the reciprocal lattice point $\mathbf{p} = h\mathbf{a}^*+k\mathbf{b}^*+l\mathbf{c}^*$ is used to compute $\Delta_\mathrm{Ewald}(\mathbf{p})$ (stored as a diagnostic and later used in profile-radius estimation).
 
@@ -250,7 +253,7 @@ where $L_\mathrm{max}$ is the maximum expected real-space unit-cell edge (Å). T
 
 ### 5.2 FFT peak picking and candidate vectors
 
-For each direction, the FFT magnitude spectrum is computed; peaks correspond to periodicities along $\mathbf{u}_d$. Each direction yields a candidate real-space length $L$ with maximum spectral magnitude (subject to $L\ge L_\mathrm{min}$).
+For each direction, the FFT magnitude spectrum is computed; peaks correspond to periodicities along $\mathbf{u}_d$. Each direction yields a candidate real-space length $L$ chosen **not** by raw magnitude but by **maximum prominence above a running-mean local background** (subtracting the broad low-frequency envelope that otherwise dominates on weak or pink-beam frames), subject to $L\ge L_\mathrm{min}$.
 
 Candidate vectors are $\mathbf{v}_d = L_d\,\mathbf{u}_d$.
 
@@ -258,22 +261,16 @@ A collinearity filter removes nearly parallel vectors (e.g. within 5°) and atte
 
 ### 5.3 Lattice reduction and cell candidates
 
-Triples of candidate vectors are combined to form candidate bases $(\mathbf{A},\mathbf{B},\mathbf{C})$. A simple reduction is applied:
-$
-\mathbf{B} \leftarrow \mathbf{B} - \mathrm{round}\!\left(\frac{\mathbf{B}\cdot\mathbf{A}}{\mathbf{A}\cdot\mathbf{A}}\right)\mathbf{A},
-$
-$
-\mathbf{C} \leftarrow \mathbf{C} - \mathrm{round}\!\left(\frac{\mathbf{C}\cdot\mathbf{A}}{\mathbf{A}\cdot\mathbf{A}}\right)\mathbf{A}
-- \mathrm{round}\!\left(\frac{\mathbf{C}\cdot\mathbf{B}}{\mathbf{B}\cdot\mathbf{B}}\right)\mathbf{B}.
-$
-
-Candidates are filtered by allowed length and angle ranges.
+Triples of candidate vectors are combined to form candidate bases $(\mathbf{A},\mathbf{B},\mathbf{C})$, each reduced to its **Niggli-reduced cell** (Gruber-vector reduction) before comparison, and filtered by allowed length and angle ranges. Two passes are run: a standard pass forms shortest-vector triples from the ~30 strongest filtered directions; if the best cell then indexes fewer than half the spots, a **widened fallback** anchors the two shortest axes and lets the third range over up to ~60 candidate vectors (deduplicated by Niggli cell), catching large, elongated or superstructure cells the first pass misses.
 
 ### 5.4 Robust refinement and best-cell selection
 
-Candidate bases are refined against observed spots using an iterative inlier‑focused least‑squares procedure (trimmed/contracting threshold). The output cell is chosen to:
-1. maximize the number of indexed spots under the tolerance $\tau$, and
-2. break ties by a refined score (smaller residual threshold/score is preferred).
+Candidate bases are refined against observed spots using an iterative inlier‑focused least‑squares procedure (trimmed/contracting threshold). Candidates are then ranked:
+1. more indexed spots wins — **unless** two candidates index within ~10 % of each other, in which case
+2. the **smaller-volume** cell is preferred (when the volumes differ by more than ~5 %), avoiding a doubled supercell, then
+3. the smaller refinement score, then the spot count again.
+
+Selection is **not limited to a single lattice**: after the best cell is accepted, further lattices are added as separate crystals provided fewer than ~40 % of their indexed spots overlap an already-accepted lattice (up to two extra by default), so split or multi-lattice crystals are indexed rather than discarded.
 
 An optional reference unit cell (if supplied) restricts acceptance to cells within a relative distance tolerance in edge lengths (permutation-invariant).
 
@@ -292,9 +289,9 @@ If not, Jungfraujoch attempts to infer the most plausible Bravais lattice type f
 The output includes:
 - a conventional cell,
 - crystal system (triclinic, monoclinic, …),
-- centering symbol $P, A, B, C, I, F, R$.
+- centering symbol (one of $P, C, I, F, R$; the $A/B$ variants are not emitted here — they are handled only later as prediction absences, §8.4).
 
-This stage provides centering information used for systematic absences in prediction (§7.3) and for reporting.
+This stage provides centering information used for systematic absences in prediction (§8.4) and for reporting.
 
 **Note.** In ambiguous or special cases, forcing space group to $P1$ (no symmetry assumptions) is recommended.
 
@@ -315,6 +312,8 @@ The refinement jointly optimizes, depending on mode and constraints:
 - crystal orientation (a global rotation),
 - unit-cell parameters, with constraints determined by inferred crystal system.
 
+By default only the beam center, unit cell and crystal orientation are refined; the detector distance, tilt angles and rotation-axis direction are held fixed unless explicitly enabled. A lighter **orientation-only** mode refines just the crystal orientation (with a weak small-rotation prior on the poorly-determined out-of-plane component), for stills whose geometry is already trusted.
+
 For higher symmetries, constraints are enforced, e.g.
 - cubic: $a=b=c,\ \alpha=\beta=\gamma=90^\circ$,
 - tetragonal: $a=b$,
@@ -341,15 +340,15 @@ For oscillation/rotation data, each image corresponds to a rotation angle $\phi$
 $
 \mathbf{s}_\mathrm{obs,ref} = R(\phi)\,\mathbf{s}_\mathrm{obs},
 $
-with $R(\phi)$ constructed from the axis-angle representation of the goniometer model.
+with $R(\phi)$ constructed from the axis-angle representation of the goniometer model. The angle $\phi$ is taken at the centre of each frame's oscillation (the frame angle plus half the oscillation width).
 
 ### 7.4 Multi-stage tightening of inlier tolerance
 
-Refinement is performed in stages with decreasing acceptance tolerance for including reflections (e.g. from coarse to fine), which stabilizes convergence when starting from imperfect indexing and approximate geometry.
+Refinement is performed in stages with decreasing acceptance tolerance for including reflections (three stages, indexing tolerance $0.3\to0.2\to0.1$), which stabilizes convergence when starting from imperfect indexing and approximate geometry.
 
 ---
 
-## 9. Bragg integration
+## 8. Reflection prediction
 
 Jungfraujoch predicts reflection positions for integration by enumerating Miller indices within a resolution cutoff and accepting those that satisfy a diffraction condition model.
 
@@ -368,6 +367,8 @@ $
 $
 Accepted reflections are projected to the detector by intersecting the diffracted direction $\mathbf{S}=\mathbf{S}_0+\mathbf{p}$ with the detector plane, using the current geometry.
 
+When the beam has a finite energy bandwidth, this window is **broadened radially per reflection**: the cutoff is combined in quadrature with a bandwidth smear, $\sqrt{\Delta_\mathrm{cut}^2 + (3\,\sigma_\mathrm{bw})^2}$, where $\sigma_\mathrm{bw}\propto|p_z|$ (the reciprocal-space depth along the beam, growing as $\sim 1/d^2$). This keeps high-resolution reflections — smeared by the bandwidth into radial streaks — from being clipped. The same $\sigma_\mathrm{bw}$ is deconvolved from the measured profile radius (§11.1), so it is not double-counted.
+
 ### 8.3 Rotation prediction (Laue equation + partiality model)
 
 For rotation/oscillation datasets, Jungfraujoch solves for rotation angles $\phi$ where the rotated reciprocal lattice point satisfies the Ewald-sphere condition. In an XDS-like notation, define:
@@ -422,7 +423,7 @@ Pixels are classified by their squared distance $r^2=(x-x_p)^2+(y-y_p)^2$:
 - **signal region:** $r^2 < r_1^2$,
 - **background annulus:** $r_2^2 \le r^2 < r_3^2$.
 
-Invalid pixels (masked/bad/saturated) are excluded from both sums.
+Invalid pixels (masked/bad/saturated) are excluded from both sums. In addition, pixels lying inside the signal disk ($r<r_2$) of any *other* predicted reflection are removed from this reflection's background annulus, so a neighbouring spot cannot leak into the background estimate.
 
 ### 9.2 Box summation (seed and fallback)
 
@@ -437,23 +438,27 @@ $
 \hat{b} = \frac{B}{n_B},\qquad
 \hat{I} = S - n_S \hat{b},
 $
-with a Poisson-like uncertainty $\sigma(\hat{I})\approx\sqrt{S}$ (floored at 1). A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum. This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below.
+with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum. This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below.
 
-The background mean is computed with a single **high-outlier reject** (drop ring pixels above $\hat{b}+3\sqrt{\hat{b}}$, recompute): a bandwidth-streaked high-resolution spot or a close neighbour can leak into the ring and bias the mean high, over-subtracting and driving weak high-resolution intensities negative. A clean Poisson background is essentially unchanged by the cut.
+For the **profile-fit path on broadband (still) data**, the background mean is additionally computed with a single high-outlier reject (drop ring pixels above $\hat{b}+3\sqrt{\hat{b}}$, then recompute): a bandwidth-streaked high-resolution spot or a close neighbour can leak into the ring and bias the mean high, over-subtracting and driving weak high-resolution intensities negative. A clean Poisson background is essentially unchanged by the cut. The reject is **not** applied to plain box summation (`--integrator boxsum`) or to monochromatic/rotation data.
 
 ### 9.3 Profile-fitted extraction (default)
 
 A fixed signal disk captures a *width-dependent* fraction of each spot, which puts a multiplicative floor on the per-observation precision of strong reflections and weights weak reflections poorly. Profile fitting removes this by extracting each intensity against a fitted spot shape, without needing reference intensities. Per frame:
 
 1. **Seed.** Box-sum every reflection (§9.2) to get a rough intensity and observed centroid, and select strong spots (significance $\ge 5$).
-2. **Build the profile.** From the strong spots, form a profile **per resolution shell**: an isotropic Gaussian of the measured second moment (the default), or an empirical averaged grid. The width is shell-dependent because spot size grows with resolution; the intrinsic spot is essentially round in the detector plane (per-detector-region and crystal-anisotropy profiles were evaluated and add nothing — the real crystal anisotropy lives in the discarded rocking direction). When a finite energy bandwidth is set, however, it smears each spot **radially** by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$ (distance from the beam centre, large at high resolution), turning high-resolution spots into radial streaks. There the profile is elongated *only along the radial direction* per reflection ($\sigma^2_\mathrm{radial}=\sigma^2_\mathrm{intrinsic}+\sigma_\mathrm{bw}^2$, tangential unchanged) on a grid grown to hold the streak, capturing it without the tangential background an isotropic widening would add.
+2. **Build the profile.** For `gaussian` (the default) the width is taken **per resolution shell** from the measured second moment of the strong spots (shell-dependent because spot size grows with resolution); the intrinsic spot is essentially round in the detector plane (per-detector-region and crystal-anisotropy profiles were evaluated and add nothing — the real crystal anisotropy lives in the discarded rocking direction). For `empirical` the profile is instead the averaged, centroid-aligned, background-subtracted pixel grid of the shell's strong spots. Either way the profile is then **rebuilt for each reflection**, centred on its **sub-pixel predicted position** (the noise-free geometric centre, not the observed centroid) and, where needed, **elongated only along the radial direction** (away from the beam centre) — because two effects stretch a spot radially but not tangentially:
+   - a finite energy **bandwidth** smears each spot by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$ ($R_\mathrm{px}$ = distance from the beam centre, large at high resolution), and
+   - sensor **parallax** — the depth over which a photon converts in a thick Si/CdTe sensor — adds a term $\propto\tan^2(2\theta)$ (material- and energy-dependent), plus, on the monochromatic path, a small fixed weak-spot capture term.
+
+   These combine as $\sigma^2_\mathrm{radial}=\sigma^2_\mathrm{intrinsic}+\sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)$ (tangential unchanged), on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add.
 3. **Fit (Kabsch).** With profile $P$, background $B$ and the shell variance model, the intensity and its uncertainty are
 $
 I = \frac{\sum P\,(c-B)/v}{\sum P^2/v},\qquad
 \sigma = \sqrt{\frac{1}{\sum P^2/v}},\qquad
 v = B + \max(I,0)\,P,
 $
-where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated. The rotation/excitation partiality is carried exactly as in the box-sum path.
+where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the variance floors the background at $1/12$ (the integer-binning pixel-variance floor). The rotation/excitation partiality is carried exactly as in the box-sum path.
 
 The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`).
 
@@ -465,7 +470,7 @@ For integrated reflections, polarization correction can be applied as a multipli
 
 ## 10. Scaling and merging
 
-After per-image integration, Jungfraujoch scales observations and merges them into unique reflections. The design is intentionally compatible with XDS/XSCALE concepts, while supporting both still and rotation partiality models.
+After per-image integration, Jungfraujoch scales observations and merges them into unique reflections. The design is intentionally compatible with XDS/XSCALE concepts, and handles both still and rotation data.
 
 ### 10.1 Observation model
 
@@ -484,40 +489,30 @@ A least-squares objective is minimized:
 $
 \sum_{ij} \left(\frac{I_{ij}^{\mathrm{pred}} - I_{ij}^{\mathrm{obs}}}{\sigma_{ij}}\right)^2
 $
-with regularization on $G_i$ and optional smoothness constraints (particularly meaningful for rotation series).
+solved by robust (Cauchy) weighted least squares, with optional post-fit smoothing of the per-frame scales for rotation series (§10.3).
 
-### 10.2 Partiality models available
+### 10.2 Partiality models
 
-Jungfraujoch supports several partiality choices:
+The partiality applied is fixed by the data type and scaling stage, not chosen from a user menu:
 
-1. **Rotation partiality** (XDS-like; see §8.3):
+1. **Rotation partiality** (XDS-like; see §8.3), used for the per-frame scaling of rotation partials:
    $
    P_{ij} = \frac{1}{2}\left[
    \mathrm{erf}\!\left(\frac{\Delta\phi_{ij}+\Delta\phi/2}{\sqrt{2}\,\sigma_{M,i}/\zeta_{ij}}\right) -
    \mathrm{erf}\!\left(\frac{\Delta\phi_{ij}-\Delta\phi/2}{\sqrt{2}\,\sigma_{M,i}/\zeta_{ij}}\right)
    \right].
    $
-   Mosaicity $\sigma_{M,i}$ can be refined per image group with bounds.
+   The mosaicity $\sigma_{M,i}$ is **measured once per image at indexing** (MLE, §11.2) and held fixed during scaling — only smoothed in frame order (§10.3), never re-refined (it is degenerate with the scale $G$; §11.2).
 
-2. **Still partiality** (excitation-error proxy):
-   $
-   P_{ij} = \exp\!\left(-\frac{\Delta_\mathrm{Ewald}^2}{R_i^2}\right),
-   $
-   where $R_i^2$ is a refined width parameter (bounded).
+2. **Unity** ($P_{ij}=1$): used for the scale-on-fulls refit (§10.6), where each observation is already a complete reflection.
 
-3. **Unity**: $P_{ij}=1$.
-
-4. **Fixed**: use the per-reflection partiality carried from prediction.
+3. **Fixed**: use the per-reflection partiality carried from prediction. Still/serial images are predicted with $P=1$, so their scaling is effectively unity/fixed — there is no excitation-error still-partiality model.
 
 Reflections below a minimum partiality can be rejected from merging to avoid unstable corrections.
 
-### 10.3 Regularization and smoothness
+### 10.3 Smoothing of per-frame scales
 
-To stabilize scale determination, a weak prior $G_i\approx 1$ is used. For rotation datasets, optional smoothness encourages slowly varying scales and mosaicity:
-$
-\log G_{i-1} - 2\log G_i + \log G_{i+1} \approx 0,
-$
-(and similarly for mosaicity), reflecting the expectation of gradual changes during a rotation scan.
+The per-frame scales $G_i$ are fit by robust (Cauchy) inverse-variance-weighted ratios; there is no explicit $G\approx1$ prior. For rotation datasets, optional smoothing enforces the expectation that scale and mosaicity vary slowly across a sweep: **after** the per-frame fit, $\log G_i$ (and the mosaicity) are replaced by a centred **moving average** over a window spanning a configurable rotation range (XDS DELPHI-like; `--smooth-g`, default 5° for rot3d, off otherwise). It is a post-fit smoothing pass, not a curvature penalty inside the least-squares objective.
 
 ### 10.4 Merging estimator
 
@@ -538,12 +533,13 @@ An internal-consistency term can inflate uncertainties when multiple observation
 ### 10.5 Merging statistics
 
 Per-shell and overall merging statistics are computed on corrected intensities, including:
-- number of observations,
-- number of unique reflections,
+- number of observations and of unique reflections, and multiplicity,
 - mean $I/\sigma(I)$,
-- an R$_\mathrm{meas}$-like quantity derived from within‑HKL deviations (shell-binned).
+- $R_\mathrm{meas}$ (the redundancy-independent Diederichs–Karplus form) from within‑HKL deviations,
+- $\mathrm{CC}_{1/2}$ (half-set correlation) and, when a reference dataset is supplied, $\mathrm{CC}_\mathrm{ref}$,
+- completeness against the enumerated reflections for the cell and symmetry.
 
-Completeness requires enumeration of possible reflections given a unit cell and symmetry; where this is not fully available, completeness may be reported as 0 or omitted.
+The error model is refined as $\sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langle I\rangle)^2$ with a systematic floor $\sigma\ge b|I|$; the asymptotic signal-to-noise $\mathrm{ISa}=1/b$ is reported and written to the output files.
 
 ### 10.6 Rotation datasets: combining partials into fulls (3D integration)
 
@@ -552,7 +548,7 @@ In a rotation scan a reflection is recorded as a series of *partials* spread acr
 The combine groups each reflection's partials into rocking events (contiguous runs of frames) and reduces each event to one full:
 
 - **De-biased weighted sum.** Partials are combined by inverse-variance weighting, where each partial's variance is its background-noise component plus the *model* signal shared across the event (Kabsch profile-fit form). Using the shared model signal rather than the individual down-fluctuating intensity stops weak partials from being over-weighted, which would otherwise inflate the merged error model. The weights depend on the full, so the estimate is iterated.
-- **Captured fraction.** The partiality summed over the event, $f=\min(1,\sum_j p_j)$, measures how completely the rocking curve was sampled; it replaces a per-partial minimum-partiality cut, because an event seen over only a few percent of its curve is unreliable however many frames it spans.
+- **Captured fraction.** The partiality summed over the event, $f=\min(1,\sum_j p_j)$, measures how completely the rocking curve was sampled. A full whose curve was captured below a threshold (`--min-captured-fraction`, default 0.7 for rotation) is dropped — an event seen over only a small fraction of its curve is unreliable however many frames it spans. (The per-partial minimum-partiality cut of §10.2 still applies upstream, in the per-frame scaling.)
 - **Capture-aware uncertainty.** A full captured incompletely ($f<1$) is extrapolated and biased high. The unobserved fraction is charged as an extra systematic uncertainty, $\sigma^2 \leftarrow \sigma^2 + \big(c\,(1-f)\,I\big)^2$, so the merge down-weights these extrapolated fulls and the error model treats their scatter as expected. It is enabled by default for the rotation path.
 
 The fulls are then re-scaled in the XDS sense — a per-image scale refit directly on the complete reflections under the unity partiality model — and merged (§10.4). Because every merged observation is now a counting-statistics-limited full rather than a partiality-divided slice, the error model reaches a far higher asymptotic $I/\sigma$.
@@ -579,7 +575,7 @@ The estimated mosaicity feeds the rotation prediction (how many frames each refl
 
 ---
 
-## 12. Wilson statistics and French–Wilson treatment
+## 12. Auxiliary statistics: ⟨I/σ(I)⟩ and Wilson plot
 
 ### 12.1 Per-shell ⟨I/σ(I)⟩
 
@@ -597,24 +593,6 @@ $
 $
 A linear regression of $\log\langle I\rangle$ vs $1/d^2$ provides an estimate of $B$, subject to basic quality checks (e.g. $R^2$ threshold).
 
-### 12.3 French–Wilson (posterior expectation of I and |F|)
-
-To mitigate negative intensities and obtain physically meaningful amplitudes, Jungfraujoch implements a French–Wilson style Bayesian treatment using per-shell mean intensity as a prior scale.
-
-For each merged observation $I_\mathrm{obs}$ with uncertainty $\sigma$, the posterior over true intensity $I\ge 0$ is:
-$
-p(I\mid I_\mathrm{obs}) \propto p(I)\,\exp\!\left(-\frac{(I_\mathrm{obs}-I)^2}{2\sigma^2}\right),
-$
-with priors differing between acentric and centric cases (standard Wilson distributions).
-
-Numerical quadrature over a scaled intensity variable is used to compute posterior moments:
-- $\langle I\rangle$,
-- $\langle |F|\rangle = \langle \sqrt{I}\rangle$,
-  and an amplitude uncertainty estimate via:
-  $
-  \sigma_F \approx \sqrt{\langle I\rangle - \langle |F|\rangle^2}.
-  $
-
 ---
 
 ## 13. Practical notes and limitations
@@ -622,4 +600,9 @@ Numerical quadrature over a scaled intensity variable is used to compute posteri
 - **Bragg integration is profile-fitted by default** (per-shell Gaussian profile, Kabsch extraction; §9.3), with plain box summation available as a fallback (`--integrator boxsum`). The profiles are built per frame from that frame's strong spots, which suits fast-feedback and serial/streaming use; a profile shared across many frames (as in full offline workflows) is not currently formed.
 - **Space-group symmetry** beyond centering absences is not necessarily enforced during prediction/integration unless the space group is supplied and used downstream.
 - **Resolution masking and ice rings** are controllable; including ice-ring spots in indexing can improve robustness for some samples but may bias refinement in others.
-- **Rotation vs still modes** differ substantially in prediction and scaling because partiality is angle-driven in rotation data and excitation-error-driven in still data.
+- **Rotation vs still modes** differ substantially in prediction and scaling: partiality is angle-driven in rotation data, while stills are predicted (within an excitation-error window) and scaled with unit partiality.
+- **Space-group determination.** When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $I(h)$ vs $I(Rh)$ plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities.
+- **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ are computed and written to the merged mmCIF as a twinning diagnostic.
+- **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default).
+- **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it.
+- **Intensities only.** The merged output carries intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`); it does not convert to amplitudes $|F|$ (no French–Wilson / truncate step) — do that downstream.
diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md
index 494702fa..3d59bed4 100644
--- a/docs/RUGNUX.md
+++ b/docs/RUGNUX.md
@@ -51,7 +51,9 @@ the first pass.
   lattice symmetry). No-reference scaling additionally emits per-iteration `<prefix>_iterN_scale.dat`.
 
 Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the
-console.
+console. By default the written resolution is trimmed automatically where CC1/2 falls off
+(`--resolution-cutoff cc-logistic`, CC1/2 target 0.30); set `--scaling-high-resolution` to fix the
+limit by hand, or `--resolution-cutoff off` to keep the full range.
 
 ## Re-scaling and re-merging (`rugnux --scale`)
 
@@ -114,7 +116,7 @@ General:
 | Option | Description |
 | --- | --- |
 | `-o, --output-prefix <txt>` | Output file prefix (default: `output`) |
-| `-N, --threads <num>` | Number of worker threads (default: 1) |
+| `-N, --threads <num>` | Number of worker threads (default: all hardware threads) |
 | `-s, --start-image <num>` | First image to process (default: 0) |
 | `-e, --end-image <num>` | Last image to process (default: all) |
 | `-t, --stride <num>` | Process every *n*-th image (default: 1) |
@@ -135,7 +137,7 @@ Spot finding:
 | `--spot-threshold <num>` | Photon-count threshold for spot finding (default: 10) |
 | `--spot-high-resolution <num>` | High-resolution limit for spot finding, Å (default: 1.5) |
 | `--max-spots <num>` | Maximum spot count (default: 250) |
-| `--detect-ice-rings[=on\|off]` | Flag ice-ring spots and exclude ice-ring reflections from scaling/merging; overrides the dataset setting (default: use the dataset value) |
+| `--detect-ice-rings[=on\|off]` | Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging; overrides the dataset/master-file setting (default: use the dataset value) |
 
 Azimuthal integration (the radial profile behind the per-image ice-ring score):
 
@@ -145,6 +147,8 @@ Azimuthal integration (the radial profile behind the per-image ice-ring score):
 | `--azim-min-q <num>` | Minimum Q, 1/Å |
 | `--azim-max-q <num>` | Maximum Q, 1/Å |
 | `--azim-phi-bins <num>` | Number of azimuthal (phi) bins (default: 1) |
+| `--polarization-correction <on\|off>` | Enable/disable the azimuthal polarization correction |
+| `--solid-angle-correction <on\|off>` | Enable/disable the azimuthal solid-angle correction |
 
 Indexing:
 
@@ -179,14 +183,20 @@ Scaling and merging:
 | `--scale-fulls` / `--no-scale-fulls` | rot3d: refit a per-frame scale on the combined fulls (XDS order, Unity model); on by default for rotation data, off for stills |
 | `--smooth-g[=deg]` | rot3d: smooth the per-frame scale *G* over a degree range before the 3D combine (XDS DELPHI-like; default 5° for rotation, 0 = off) |
 | `--capture-uncertainty <num>` | rot3d: systematic sigma on under-captured fulls, ~num·(1−captured_fraction)·I (default: 1.0 for rotation, 0 otherwise) |
-| `--scaling-high-resolution <num>` | High-resolution limit for scaling, Å (default: no limit) |
+| `--min-captured-fraction <num>` | rot3d: drop a combined full whose rocking curve was captured below this fraction — edge-of-sweep truncated fulls (default: 0.7 for rotation, 0 otherwise; 0 = off) |
+| `--scaling-high-resolution <num>` | High-resolution limit for scaling, Å — manual override (default: no limit; disables the automatic cutoff below) |
+| `--resolution-cutoff <txt>` | Automatic high-resolution cutoff for the written reflections and reported shells: `cc-logistic` \| `off` (default: `cc-logistic`; ignored when `--scaling-high-resolution` is set) |
+| `--resolution-cc-target <num>` | CC1/2 target defining the `cc-logistic` fall-off (default: 0.30) |
+| `--resolution-shells <num>` | Number of resolution shells in the reported statistics table (default: 10) |
 | `--min-partiality <num>` | Minimum partiality to accept a reflection (default: 0.02) |
 | `--reject-outliers <num>` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for `rot3d`, off otherwise) |
 | `--reject-delta-cchalf <num>` | Drop images with ΔCC1/2 below mean − N·stddev (default: off) |
 | `--min-image-cc <num>` | Per-image CC limit, percent (default: no limit) |
+| `--mosaicity <num>` | Diagnostic: fix the scaling mosaicity (°) instead of using the per-image seed |
 | `--scaling-iterations <num>` | Scaling iterations with no reference data (default: 3) |
 | `--scaling-output <txt>` | Reflection output format: `cif` (mmCIF, default) \| `mtz` \| `txt` |
 | `-z, --reference-mtz <file>` | Reference MTZ (enables reference-driven scaling) |
+| `--reference-column <label>` | Reference MTZ column to use (default: auto — F-model, else IMEAN/I/…) |
 | `--write-process-h5` | Also write the (large) `_process.h5` when merging (default: only `.mtz`/`.cif`) |
 
 Integration:
@@ -196,3 +206,15 @@ Integration:
 | `--integrator <txt>` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) |
 | `--integration-radius <r>` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` |
 | `--bandwidth <num>` | Relative X-ray bandwidth FWHM (e.g. `0.01` for a 1% DMM); default from file or 0 (monochromatic) |
+
+Geometry overrides (defaults are taken from the input file; override them to reprocess with a corrected geometry):
+
+| Option | Description |
+| --- | --- |
+| `--beam-x <num>` | Beam centre X (pixel) |
+| `--beam-y <num>` | Beam centre Y (pixel) |
+| `--detector-distance <num>` | Detector distance (mm) |
+| `--wavelength <num>` | Wavelength (Å) |
+| `--rot1 <num>` | PONI detector rotation 1 (rad) |
+| `--rot2 <num>` | PONI detector rotation 2 (rad) |
+| `--polarization <num>` | Polarization factor |
diff --git a/image_analysis/azint/ICE_RING_DETECTION.md b/image_analysis/azint/ICE_RING_DETECTION.md
deleted file mode 100644
index 1a512375..00000000
--- a/image_analysis/azint/ICE_RING_DETECTION.md
+++ /dev/null
@@ -1,71 +0,0 @@
-# Auto-gated ice-ring detection & estimation (design)
-
-_Status: designed and Python-prototyped (2026-06-30), not yet implemented in C++. The current
-`--detect-ice-rings` flag uses the fixed tabulated `ICE_RING_RES_A` with a 0.03 q half-width and must
-be enabled by hand; this document describes making it automatic and data-driven._
-
-## Motivation
-
-`--detect-ice-rings` today is a manual on/off flag over a fixed ring list. Two improvements:
-1. **Auto-gate** it — turn the handling on only when ice is actually present, so it can default to
-   `auto` without risking dropping good reflections on ice-free data.
-2. **Estimate** the real ring positions/widths from the data instead of the table (rings drift slightly
-   with detector geometry, λ and ice texture/temperature).
-
-Both fall out of the azimuthally-averaged radial **background** profile, which the `azint` engine already
-produces (`AzimuthalIntegrationProfile`) — ice rings are sharp peaks in it, distinct from Bragg spots
-(which are localised, not azimuthally complete).
-
-## Method (validated in Python on raw NXmx frames)
-
-1. Accumulate a **clipped-mean radial background profile** over ~30 frames spread across the sweep
-   (clip at ~20 counts to kill Bragg/hot pixels; mask the beamstop and bad pixels). This is exactly
-   `AzimuthalIntegrationProfile::GetResult()`.
-2. Smooth-baseline it (rolling median, ±25 bins) and take `excess = profile − baseline`;
-   `rel = excess / max(baseline, ε)`.
-3. Map radius → d → q (`θ = 0.5·atan(R·pix/dist)`, `d = λ/2sinθ`, `q = 2π/d`).
-
-### Detection statistic (self-calibrating z-score)
-
-For each tabulated ice line, `z = (max rel within ±0.03 q of the line − μ_ctrl) / σ_ctrl`, where the
-control population is all valid q at least 0.06 away from every ice line (in-image noise — so the
-threshold needs no per-detector tuning).
-
-| dataset | max line z | #lines z>5 | verdict |
-|---|---:|---:|---|
-| EP0117 (diffuse ice) | 11.8 | 6 | ICE |
-| EP0124 (sharp 3 rings) | 17.0 | 3 | ICE |
-| EP0210 (control) | 1.6 | 0 | no ice |
-| MyoB (control) | 0.8 | 0 | no ice |
-
-**Gate: enable ice handling iff ≥2 ice lines have z > 5.** Large, clean margin.
-
-### Position & width estimation
-
-- **Positions:** `ICE_RING_RES_A` is accurate for the JUNGFRAU/λ=0.953 Å — measured centroid offset
-  Δq = +0.004 ± 0.011 q (~0.25 %). No recalibration needed; the table stands.
-- **Widths:** measured ring FWHM is **0.055–0.068 q** ⇒ half-width ~0.027–0.034. The fixed 0.02 q
-  under-covered the strong low-res rings (already bumped to **0.03** as the fallback). Best is per-ring
-  `half_width = 0.5·FWHM` clamped to [0.012, 0.040], estimated from the peak shape.
-
-## Proposed wiring (small, no new data path)
-
-- New `image_analysis/azint/IceRingDetect.{h,cpp}`:
-  `IceRingModel DetectIceRings(const AzimuthalIntegrationProfile&)` → `{bool present; vector<IceRing>
-  rings;}` with per-ring `{centre_q, half_width_q}`.
-- Add an `IsOnIceRing(d, const vector<IceRing>&)` overload in `common/Definitions.h` (per-ring centre
-  + width), alongside the existing fixed-table one.
-- `tools/rugnux_cli.cpp`: change `--detect-ice-rings` to take `auto|on|off` (default `auto`); in
-  `auto`, run `DetectIceRings` on the accumulated profile and set
-  `experiment.DetectIceRings(model.present)` + feed `model.rings` to the spot-finder
-  (`MarkIceRings`) and the scaling exclusion path.
-- Online (broker): accumulate the profile over the first N frames, freeze the model, apply thereafter.
-
-Prototype scripts + the `ice_ring_detect.diff` sketch are in the session scratchpad
-(`scratchpad/ice_detect/`).
-
-## Note
-
-For the EP0117 test crystal specifically, ice is a **secondary** issue — its dominant problem is
-radiation damage across the sweep, addressed by per-observation outlier rejection (`--reject-outliers`),
-not ice handling. Auto-gated ice detection matters for genuinely ice-dominated, otherwise-good crystals.
diff --git a/image_analysis/bragg_integration/BRAGG_INTEGRATION_ENGINE.md b/image_analysis/bragg_integration/BRAGG_INTEGRATION_ENGINE.md
deleted file mode 100644
index a2f5eeed..00000000
--- a/image_analysis/bragg_integration/BRAGG_INTEGRATION_ENGINE.md
+++ /dev/null
@@ -1,107 +0,0 @@
-# BraggIntegrationEngine — GPU box + profile-fit integrator
-
-Status: **standalone, not yet wired into the pipeline.** This is a GPU (CUDA) reimplementation of
-`ProfileIntegrate2D.cpp` (the default Kabsch profile-fit integrator, ~142 ms/frame on CPU) and
-`BraggIntegrate2D.cpp` (the box-sum integrator) under one roof, following the
-`AzIntEngine` / `ROIIntegration` base + CPU + GPU pattern. It reproduces the CPU result up to
-floating-point precision and runs in **< 2 ms/frame** for realistic reflection counts.
-
-## Files (all in `image_analysis/bragg_integration/`)
-
-| file | role |
-|------|------|
-| `BraggIntegrationEngine.{h,cpp}` | base class: extracts the fixed per-experiment config (r-radii, ellipse coefficients, geometry) and owns `Finalize()` (polarization + `image_scale_corr` + `vector<Reflection>` assembly) |
-| `BraggIntegrationEngineCPU.{h,cpp}` | plain-C++ engine: the fallback **and** the numeric oracle for the GPU |
-| `BraggIntegrationEngineGPU.{h,cu}` | the CUDA engine |
-| `tests/BraggIntegrationEngineGPUTest.cpp` | GPU-vs-CPU equivalence test + hidden `[bragg_bench]` benchmark |
-
-`CMakeLists.txt` builds the base + CPU always; the `.cu` + GPU header are added only under
-`JFJOCH_CUDA_AVAILABLE` (mirrors `azint`/`roi`). No change needed to `image_analysis/CMakeLists.txt`
-(the aggregate `JFJochImageAnalysis` already links `JFJochBraggIntegration` and
-`JFJochImagePreprocessing`, so the `ImagePreprocessorBuffer` symbols resolve at the final link).
-
-## Interface
-
-```cpp
-class BraggIntegrationEngine {
-  BraggIntegrationEngine(const DiffractionExperiment&);
-  virtual std::vector<Reflection> Run(const ImagePreprocessorBuffer& image,   // preprocessed int32
-                                      const std::vector<Reflection>& predicted, size_t npredicted,
-                                      int64_t image_number) = 0;
-};
-// CPU: BraggIntegrationEngineCPU(experiment)
-// GPU: BraggIntegrationEngineGPU(experiment, std::shared_ptr<CudaStream>)
-```
-
-The integrator mode comes from `experiment.GetBraggIntegrationSettings().GetIntegrator()`:
-`BoxSum` (= BraggIntegrate2D), `ProfileGaussian` (default), `ProfileEmpirical`. The box sum is also
-Pass A / the seed of the two profile modes, so it always runs.
-
-## Key design points
-
-- **Input is the preprocessed int32 buffer**, NOT the raw `CompressedImage`. Masked/bad pixels are
-  `INT32_MIN`, saturated are `INT32_MAX` (a pixel is valid iff `v != INT32_MIN && v != INT32_MAX`).
-  This is the same buffer `AzIntEngineGPU`/`ROIIntegrationGPU` consume. **Consequence:** the `±1`
-  special/saturation band that `ProfileIntegrate2D` rejects on the lossy `CompressedImage` is handled
-  upstream by the preprocessor instead — correct for the online path; sanity-check it for offline
-  `rugnux` if that ever uses this engine on a lossy-compressed stored file.
-- Config uses raw detector dims — `xpixel = GetXPixelsNum()`, `npixel = GetPixelsNum()` — matching how
-  `MXAnalysisWithoutFPGA` sizes the `ImagePreprocessorBuffer` and the frame the predicted `predicted_x`
-  live in.
-- **GPU: one CUDA block per reflection**, shared-memory reductions. Kernels (all in the `.cu`):
-  `reset → mark_mask → boxsum → learn_profile → build_profiles → fit` (Kabsch, 4 iters). The
-  resolution shell is computed inline (`compute_shell`), so there is no separate shell pass. `BoxSum`
-  mode stops after `boxsum`.
-- **Precision:** the hot path (profile learning + Gaussian fit) is `float` (FP64 is throttled on
-  consumer GPUs; the extraction is Poisson-noise limited, so float matches the double CPU to ~1e-4).
-  `double` is kept only for the box-sum background mean / centroid, where it drives weak intensities.
-- Per frame only the predicted centres (`predicted_x/y`, `d`) are uploaded; the image is already on
-  the device via `ImagePreprocessorBufferGPU::getGPUBuffer()`. Device buffers grow-only
-  (`EnsureCapacity`). Everything runs on the passed-in `CudaStream`, one sync at the end.
-
-## Build & test
-
-```
-cd build && cmake .                       # pick up the new CMake sources
-make -j$(nproc) jfjoch_test
-cd tests
-./jfjoch_test "BraggIntegrationEngineGPU_MatchesCPU"   # GPU==CPU, 4 modes (skips if no GPU)
-./jfjoch_test "[bragg_bench]"                          # perf sweep vs reflection count
-```
-
-The equivalence test compares I / sigma / bkg within a small tolerance (float vs double + unordered
-atomic summation of the learned profile → ~1e-4, not bit-exact — unlike the integer ROI test).
-
-## Performance (measured on a throttled RTX PRO 2000 Blackwell **laptop**, WSL2 — worst case)
-
-~0.5 ms @528 refl, 1.1 ms @2600, 2.3 ms @5325 (stress). ~0.45 ms of that is fixed WSL2 kernel-launch
-latency; on a native-Linux datacenter GPU it is sub-millisecond throughout. 15–20× vs the CPU engine.
-(First cut was 10.8 ms with FP64 → 3.8 ms with float → the above after collapsing ~19 CUDA API
-calls/frame into one `reset` kernel + inline shell.)
-
-## To integrate (the binding — do this on the target machine)
-
-The engine is deliberately *not* called anywhere yet. Wiring it in:
-
-1. **Provide the preprocessed buffer + stream to the integration call.** Today
-   `MXAnalysisWithoutFPGA::Analyze` (`image_analysis/MXAnalysisWithoutFPGA.cpp`, ~line 101-104) wraps
-   `preprocessor_buffer->getBuffer()` back into a `CompressedImage` and hands it to
-   `indexer.ProcessImage(...)`. To use this engine you instead want to pass the
-   `ImagePreprocessorBuffer` (the GPU subclass carries the device image) and the `CudaStream` through
-   to where integration happens.
-2. **Swap the integrator call.** `IndexAndRefine.cpp` (~line 337-339) currently does:
-   ```cpp
-   i_outcome.reflections = mode == IntegratorMode::BoxSum
-       ? BraggIntegrate2D(...)      // CompressedImage
-       : ProfileIntegrate2D(...);
-   ```
-   Replace with a `BraggIntegrationEngine::Run(buffer, prediction.GetReflections(), nrefl, msg.number)`
-   on a member engine constructed once (GPU when `get_gpu_count() > 0`, else CPU — mirror how
-   `MXAnalysisWithoutFPGA` chooses `AzIntEngineGPU` vs `AzIntEngineCPU`). Construct with the same
-   `stream` the other GPU engines share.
-3. **Keep the old functions** (`ProfileIntegrate2D`/`BraggIntegrate2D`) until the engine is validated
-   end-to-end against them on real data (`rugnux --integrator ... --dump-observations`, A/B vs
-   `XDS_ASCII.HKL`), since the `±1`-band difference above is the one behavioural change.
-
-The output `vector<Reflection>` is identical in shape (I, sigma, bkg, partiality, d, ...), so
-`ScaleImage`/`CalcISigma`/`Combine3D`/merge downstream consume it unchanged.
diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h
index 8fe5d1dc..a9e2c6a8 100644
--- a/image_analysis/bragg_integration/BraggIntegrationEngine.h
+++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h
@@ -25,9 +25,10 @@
 //   ProfileEmpirical-> per-shell learned empirical profile fit
 // The box sum is also the seed pass (Pass A) of the two profile modes, so it always runs.
 //
-// This class is intentionally standalone: it is NOT yet wired into IndexAndRefine. It takes a
-// preprocessed image + the predicted reflections and returns the same vector<Reflection> shape
-// (I, sigma, bkg, partiality, ...) that the downstream scaling/merge consumes unchanged.
+// This is the Bragg integrator used by the pipeline (bound in MXAnalysisWithoutFPGA: the GPU
+// engine when a device is present, otherwise the CPU engine). It takes a preprocessed image +
+// the predicted reflections and returns the vector<Reflection> (I, sigma, bkg, partiality, ...)
+// that the downstream scaling/merge consumes unchanged.
 // =============================================================================
 
 #include <cmath>
diff --git a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md
deleted file mode 100644
index a81b7159..00000000
--- a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md
+++ /dev/null
@@ -1,391 +0,0 @@
-# Next-generation integrator & the XDS-quality investigation
-
-_Last updated 2026-06-25. Supersedes the 2026-06-24 design draft — `ProfileIntegrate2D` is now
-built, wired, and the default integrator; this document records the design **and** everything the
-data-quality investigation established (especially the dead ends, so they are not re-tried)._
-
----
-
-## TL;DR
-
-- **The question that started this:** jfjoch's merged data was ~3–4× worse than XDS on per-observation
-  precision (R-meas, ISa, anomalous peak height) on the HEWL test crystals. Where is the gap, and can
-  we close it **without** XDS-style 3D shoeboxes (Filip's philosophy: integrate per-frame, combine at
-  the end, keep ~300 Hz throughput)?
-- **Diagnosis (decisive):** our **mean intensities are accurate everywhere** (`CC(ours, XDS) = 0.994`,
-  ≥0.986 in every resolution shell) and **unbiased** (ANODE: our anomalous peaks sit on the same real
-  atoms as XDS). The gap is **per-observation noise**, localized to the **2D box-sum integration** of
-  **strong** reflections — a ~18% multiplicative floor that a fixed disk imposes by capturing a
-  width-dependent fraction of each spot.
-- **What moved the needle:** swapping box-sum for a **profile-fitting integrator** (`ProfileIntegrate2D`,
-  Gaussian profile). That is the single lever that ever improved intensity quality. Everything on the
-  *scaling* side was neutral or negative.
-- **Current best on the test crystal** (HEWL `fixed_master.h5`, anomalous S/Cl peak height as a fraction
-  of XDS, `1.0` = match): **0.61×**, ISa **11.0**, via Gaussian profile-fit + `-R100` + `rot3d` +
-  `--scale-fulls`. Box-sum baseline was **0.46×**, ISa 8.4. XDS is 1.0× / ISa ~28.
-- **The remaining gap is XDS's full XSCALE/CORRECT scaling machinery**, whose individual pieces are each
-  neutral on this well-behaved crystal — not a quick win, and not a bias to fix.
-
----
-
-## Philosophy & constraints (Filip)
-
-1. **Per-frame integrate, combine at the end.** Jungfraujoch integrates each frame as a 2D image and
-   combines the per-frame partials of a rotation reflection into one "full" afterwards. We deliberately
-   do **not** batch frames into XDS-style 3D shoeboxes in the online path. Push per-frame + combine as
-   far as it goes.
-2. **Throughput matters.** ~300 Hz end-to-end was demonstrated; that is an order of magnitude over XDS.
-   Judge every improvement by its **performance penalty**, not just quality (this is why global scaling
-   was removed — see below).
-3. **Modular, not bolted-on.** The integrator is swappable (`BoxSum` ↔ `ProfileGaussian` ↔
-   `ProfileEmpirical`) behind a single setting; scaling, the rot3d combine, and the merge consume a
-   plain `Reflection{I, σ, partiality, d}` and are untouched by the integrator choice. Both **stills and
-   rotation** benefit.
-4. **Simple, readable code.** No speculative/defensive machinery; neutral experiments get reverted.
-
----
-
-## The pipeline today
-
-```
-per-frame 2D integration          (BraggIntegrate2D  |  ProfileIntegrate2D)   <- the lever
-  -> per-image scaling            (ScaleOnTheFly: Ceres fit G, B_resolution, mosaicity, wedge)
-  -> rot3d combine                (Combine3D: weight-sum per-frame partials -> one full per event)
-  -> scale fulls (optional)       (ScaleOnTheFly again, Unity model, on the combined fulls)
-  -> merge + error model + stats  (MergeOnTheFly: RefineErrorModel, CC1/2, R-meas, CCref, ISa)
-```
-
-Relevant knobs (`rugnux`, mirrored in the viewer's Processing-settings tabs):
-- `--integrator boxsum|gaussian|empirical` — default **gaussian** (`IntegratorMode::ProfileGaussian`).
-- `-P fixed|rot|rot3d` — partiality model. `rot3d` = `PartialityModel::Rotation` **+** the orthogonal
-  `ScalingSettings::combine_3d` bool (see "Why rot3d is a bool" below).
-- `--scale-fulls` — second scaling pass on the combined fulls (XDS order). Implies `-M`, needs `rot3d`.
-- `-R<N>` — first-pass rotation-indexing images; default **100**.
-- `-M` — run scaling; `-A` — anomalous (keep Friedel mates separate).
-
-### Current integrator: `ProfileIntegrate2D`
-
-Drop-in for `BraggIntegrate2D`, same `vector<Reflection>` output, **no reference intensities**.
-Selected at `IndexAndRefine.cpp:347`. Per frame:
-
-- **Pass A — rough + pick.** Box-sum every predicted reflection (rough `I`, background, centroid);
-  select strong spots (significance ≥ 5).
-- **Pass B — build the profile**, per resolution shell, from the strong spots:
-  - `ProfileGaussian` (**default, the keeper**): a Gaussian of the measured second moment.
-  - `ProfileEmpirical`: the averaged, centroid-aligned, background-subtracted pixel grid.
-- **Pass C — profile-fit extract (Kabsch).** `I = Σ P(c−B)/v ÷ Σ P²/v`, de-biased variance
-  `v = B + max(I,0)·P` (iterated), `σ = √(1/Σ P²/v)`. Rotation `partiality` carried exactly as
-  `BraggIntegrate2D` (so `rot3d` works unchanged).
-
-**Why it removes the floor:** profile weighting integrates out to the spot tail (the width-varying part
-a fixed disk clips) while down-weighting the noisy edge (avoiding the wide-box background penalty). The
-floor is multiplicative and lives in strong reflections, so matched weighting + de-biased variance is
-the right lever — **not** a tighter aperture (a tight profile *lost* to a generous box for stills, the
-0.4-px centroid-undersampling floor).
-
-**Measured, do not re-try:**
-- **Anisotropy adds nothing — SUPERSEDED 2026-06-30 (round 3 below).** The original test was confounded
-  by an over-wide (contaminated) width; once the width is de-contaminated, a **radial** elongation does
-  help at high resolution. Per-detector-region profiles still add nothing.
-- **The empirical profile under-performs the Gaussian** as built (per-frame, integer-pixel-binned,
-  centred on `round(predicted)`): sub-pixel smear + too few strong spots per frame/shell. It would need
-  sub-pixel alignment + cross-frame accumulation to compete; not currently worth it.
-- _(Note: the header comment in `ProfileIntegrate2D.h` still says "DRAFT / not yet wired in" — that is
-  stale; it is wired and default.)_
-
-### The rot3d combine (`Combine3D`)
-
-Buckets per-frame partials by raw `(h,k,l)`, splits into contiguous-frame rocking events (frame gap ≤2),
-inverse-variance weight-sums each event's corrected partials into **one full**
-(`F = Σ wⱼ·I_corr,j / Σ wⱼ`, `σ = 1/√Σ wⱼ`; partiality = 1, `image_scale_corr` = 1), dropping events
-with `Σ pⱼ < min_partiality`.
-
-- **De-biased Poisson variance (kept):** weight each frame by `var = max(0, corr²(σ²−I)) + corr·max(0,F)`
-  (background noise + **model** signal, iterated 3×) instead of `σ_corr²`. Using the model signal (not
-  the down-fluctuating observed `I`) stops weak partials being over-weighted. ISa 5.3→6.4 (ref),
-  6.7→8.4 (no-ref).
-- **Validation that the grouping is structurally right:** jfjoch produces ~the **same observation count
-  as XDS** (936k vs 961k, multiplicity 25.1 vs 25.9). The residual gap is pure per-observation noise,
-  not bad event grouping.
-
-**Why `rot3d` is a `ScalingSettings::combine_3d` bool, not a 4th `PartialityModel` enum.** Many sites
-gate on `PartialityModel::Rotation` *exactly* (e.g. `GetRotationWedgeForScaling`, rotation prediction).
-A 4th enum value silently broke per-image rotation handling → NaN scales, all partials rejected. So
-`-P rot3d` sets `PartialityModel::Rotation` **and** `combine_3d = true`; the pipeline only ever sees
-`Rotation`, and the post-pass keys off `GetScalingSettings().GetCombine3D()`. **Do not reintroduce a
-rot3d enum value.**
-
-### Scale-on-fulls (`--scale-fulls`) — XDS order
-
-XDS/DIALS scale on the 3D-integrated **fulls**; jfjoch scales per-frame **partials** then combines, so
-the per-frame scale `G` is fit via `G·partiality·B·lp·Itrue − Iobs` — rocking-curve/partiality **model
-error leaks into G**, and the combined fulls inherit that partial scale (never scaled *as* fulls). Fix:
-a second `ScaleOnTheFly` pass on the combined fulls with the **Unity** model (`G·Itrue − I_full`, no
-partiality), between combine and merge. **This is the only post-integration win:** anomalous 0.53→0.57×,
-ISa 9.4→10.5 — and it improves **both** precision and accuracy (unlike outlier rejection, which trades
-accuracy for CC½). Kept.
-
----
-
-## The diagnosis — where the gap is (and is not)
-
-The investigation localized the gap with three independent instruments (see Tooling) and a per-reflection
-comparison against `XDS_ASCII.HKL` (same 4/mmm ASU). Findings:
-
-**Pipeline-wide systematic-error map (ours vs XDS, per reflection):**
-- `CC(ours, XDS) = 0.994` overall, ≥0.986 in **every** resolution shell → **mean intensities are accurate
-  everywhere.** The deficiencies are systematic *scale trends*, not random error.
-- Ratio ours/XDS by resolution: 3.0 Å **0.94 (strong, −6%)** → 1.8 Å 1.02 → 1.45 Å 1.19 → 1.2 Å
-  **1.27 (weak, +27%)**. Two distinct problems:
-  1. **Strong low-res:** mean ~right but ~6% low (mild saturation/non-linearity at high counts) **and**
-     per-observation scatter R-meas **13% vs XDS 2.9%** — averages out (CC fine) but kills **anomalous**
-     (Bijvoet differences need <1% agreement). This is the integration-quality floor.
-  2. **Weak high-res:** **+20–27% high = background under-subtraction.** Code-confirmed: the integrator
-     subtracts a **flat ring mean** (`ProfileIntegrate2D.cpp` `out.bkg = bkg_sum/n_bkg`), so an unmodeled
-     background **gradient** under the peak biases weak reflections high, worst at high res.
-
-**Anomalous accuracy (ANODE, external benchmark):** averaged anomalous peak height at the model atoms,
-ours vs XDS — SD_MET 6.05 vs 13.23, SG_CYS 5.56 vs 11.79, CL_CL 4.06 vs 9.58 (XDS ~2.2× stronger). **But
-both put their strongest peaks on the same real atoms** (<0.2 Å off; non-anomalous atoms in the noise).
-→ our data is **accurate / unbiased**; the gap is **precision (noise)**, ~2.2× on the anomalous ΔF,
-consistent with R-meas (~3×) and ISa (~3.5×). We are not systematically wrong, just ~2.2× noisier.
-
-**Where the strong-reflection floor lives:** it is in **strong** reflections (jf/XDS R-meas ratio 2.1×
-weak → 6.5× top-5%), a ~18% **multiplicative** floor, **flat across resolution**, **per-observation**
-(per-image/scaling is only 6.6% of 20.7%), **immune to box radius** (tight/default/wide all ≈ 18–24%),
-and present even in fully-captured spots (so not partiality/capture). It survives every box size →
-it is the **box-sum method** (uniform pixel weighting), confirmed by the profile-fit win.
-
----
-
-## Experiment ledger — what was tried (the dead ends are the value)
-
-Scored on the HEWL test crystal, anomalous S/Cl peak height as ×XDS (`1.0` = XDS), unless noted.
-Baseline progression that **worked**:
-
-| change | effect | verdict |
-|---|---|---|
-| box-sum → **2D Gaussian profile-fit** | 0.46 → **0.53×**, ISa 8.4 → 9.4, R-meas 18.9 → 14.6% | ✅ **the lever** (committed default) |
-| **`rot3d` combine** (+ de-biased Poisson var) | ISa 1.7 → 5.3 → 6.4; R-meas (partials 67% → fulls 18.9%) | ✅ keep |
-| **`--scale-fulls`** (XDS order) | 0.53 → **0.57×**, ISa 9.4 → 10.5 (precision **and** accuracy) | ✅ keep |
-| **`-R100`** (more first-pass indexing images) | 0.57 → **0.61×**, ISa 10.5 → 11.0 (~free; sweet spot ~100) | ✅ keep (committed default) |
-| masked-sentinel fix (`±1` band) | neutral on lyso; correctness (gap-clipping reflections elsewhere) | ✅ keep (committed) |
-
-Things that **did not work** (do not re-prototype without new evidence):
-
-| change | result | why it failed |
-|---|---|---|
-| 3D box-sum / single 3D background | R-meas identical, ISa bump is a σ-model artifact | the 3D box-sum is reconstructible from the 2D partials; recovers nothing |
-| **3D profile-fit** (`Integrate3D`, voxel Kabsch) | 0.53 → **0.56× (+6%)**, diminishing (box→2D +15%, 2D→3D +6%) | integration cannot bridge the rest; offline-only, slower, dropped |
-| **global (joint) scaling** | 0.54× (= 0.53 alternating), **+120 s** | block-coordinate loop already reaches the joint optimum; **removed (perf)** |
-| **absorption surface** (per-image / global quadratic in detector xy) | 0.54×, ISa 9.4 → 7.3 (worse) | radially degenerate with per-HKL resolution structure → "intensity sponge"; we're unbiased, a position-scale can't help |
-| **learned rocking-curve profile** (on partials) | ISa 9.4 → 6.7; with `--scale-fulls` 0.37× | partials are already divided by the model partiality; the learned `c(t)` spikes at curve edges and up-weights noise |
-| **per-batch geometry refinement** (30° batches) | 0.61 → 0.56–0.57× | batch deltas confirm real smooth drift, **but** per-image refine wants the most *stable* start; a local-but-noisy per-batch start hurts |
-| **sinusoidal geom-drift model** (beam/distance vs φ) | ~0 drift found; 0.57–0.58× | this crystal is well-centred (beam sub-pixel, distance 0.04%); physically sound, no drifty data to demonstrate |
-| background **plane fit** (a + gₓ·dx + g_y·dy over the ring) | anomalous identical; +27% weak bias persists | the weak bias is **not** a linear gradient (it's radial / profile-vs-summation); no superfluous code |
-| per-frame **B-factor** decay (`--scale-fulls -B`) | 0.57 → 0.56× | no significant radiation damage; the scalar `G` already absorbs it |
-| outlier rejection (`--reject-outliers`, `--reject-delta-cchalf`) | CC½ ↑ **but anomalous ↓** (0.53 → 0.49×) | **the accuracy trap** — improves precision, degrades accuracy; not Bijvoet-safe. Do not enable for anomalous |
-| raise `--min-partiality` (0.02 → 0.6) | 0.53 → 0.35× | low-partiality partials *feed* the scale fit; partiality model is good, post-refining it doesn't help |
-| cut high resolution | hurts (1.2 Å > 2.0 Å) | high-res helps; not over-weighted noise |
-| more scaling iterations; error-model tuning | flat | scaling converged; the error model is a **thermometer** (fits the scatter, doesn't reduce it) |
-| beam-stop masking | anomalous unchanged | only ~5 low-res reflections affected on this well-behaved crystal |
-| classical post-refinement | already done | cell (a=78.23, c=37.78 ≈ model), orientation (XtalOptimizer), and **mosaicity** (Ceres in `ScaleOnTheFly`) are already refined |
-
-**Radiation damage ruled out:** anomalous rises monotonically with dose (first 300 frames 0.24× → full
-1800 0.61×) → redundancy-dominated, not damage-limited; per-frame B fits flat.
-
-**`-R100` is refinement, not indexing:** FFT saturates by ~33 first-pass images; the 10.5→11.0 gain is
-refining the **global** geometry against more spots (30→50 populated images). **Per-image refinement is
-essential** and fits real per-frame signal: `-r none` (global only) collapses to ISa 5.1 / 0.31×;
-`-r orientation` 0.56×; default `BeamCenter` (per-image orient + beam) 11.0 / 0.61×. Finer + more DOF
-wins — the per-image beam refit is really a per-frame position correction (crystal motion /
-goniometer sphere-of-confusion).
-
----
-
-## Tooling (reusable instruments)
-
-- **`--dump-observations <file>`** (`rugnux`): `Combine3D` writes the unmerged fulls
-  (`h k l I σ d n_frames captured_fraction peak_frame`) for direct A/B vs `XDS_ASCII.HKL`.
-- **`lyso_test/anomalous_scoreboard.sh <jfjoch_anomalous.hkl> [label]`**: one command → SHELXC + ANODE
-  under `qemu-x86_64-static` (the SHELX binaries hit the legacy vsyscall page; this WSL2 kernel is
-  `vsyscall=none`, so run them under qemu) → prints `⟨d″/σ⟩` per shell + averaged anomalous peak height
-  at SD_MET/SG_CYS/CL_CL vs a cached XDS reference. ~30 s/rerun. **This is the accuracy benchmark** —
-  anomalous ΔF is tiny, so it is exquisitely sensitive to systematic error, and it caught the outlier-
-  rejection accuracy trap that CC½/R-meas hid.
-- **Per-reflection compare** (`ccp4-python` + gemmi, P43212 ASU): `CC(ours, XDS)` per shell and the
-  ours/XDS ratio-vs-resolution that produced the systematic-error map.
-- **`[timing]` log lines** in the scaling/merge post-pass (per sub-step, plus the serial `MergeAll` vs
-  parallel `ScaleAllImages` split) — to find performance bottlenecks.
-
----
-
-## Performance notes
-
-- **Scaling is already multi-threaded.** `ScaleOnTheFly` parallelizes across images (all cores); each
-  per-image Ceres solve is single-threaded (correct grain).
-- **Global (joint) scaling was removed** (2026-06-25): 122 s vs 3.5 s for the alternating loop (~35×) for
-  **no** quality gain (0.54 ≈ 0.53). Deleted everywhere (`GlobalScale.{h,cpp}`, the CLI flag, the
-  `ScalingSettings` field, the viewer checkbox).
-- **The serial post-pass bottlenecks** (the real MT targets, if pursued): `Combine3D` (~1.8 s),
-  `MergeAll` (~1.4 s × scaling iterations), the final merge (~0.5 s).
-- **Latent perf trap (open):** stats computation blows up to ~110 s when **no space group is set**
-  (huge unmerged set → `SearchSpaceGroup` auto-detection over all of it; with a SG, 0.14 s). Needs a
-  guard/cap — not yet touched.
-
----
-
-## Where the remaining gap is — honest assessment
-
-Box-sum → 2D profile (+15%) → 3D profile (+6%) shows integration improvements with **diminishing
-returns**; integration alone cannot reach 1.0×. ANODE says we are unbiased, just ~2× noisier. Every
-isolated *scaling* improvement (global fit, absorption surface, per-frame B, more iterations, error-model
-tuning) was neutral or negative on this well-behaved crystal. The conclusion: **the residual ~1.8× is
-XDS's full XSCALE/CORRECT machinery acting together** (multi-dimensional smooth scale + per-batch
-corrections + tighter error model, fit jointly), whose individual pieces are each neutral here — a large
-reimplementation with low odds of a quick win.
-
-**Candidate directions, with current expectation:**
-1. **Weak high-res background modelling** (the +20–27% bias) — radial / profile-vs-summation, **not** a
-   linear gradient (the plane fit was neutral). A genuine target, but the simple fix doesn't work.
-2. **Strong-reflection per-observation profile scatter** — the integration-quality floor; 2D profile-fit
-   already took the easy part, 3D profile gives a diminishing +6% and needs a 3D pixel shoebox (departs
-   from the per-frame philosophy; reader-blocked online).
-3. **A learnable model.** The systematic-error map is, in effect, the "what to do" answer for an ML
-   integrator: the weak background bias, the strong-reflection saturation/response, and the per-observation
-   profile scatter are all **model-learnable** corrections — none is a simple hand-coded fix.
-
-**What is settled and should not be re-litigated:** we are unbiased (ANODE); 3D summation/background is
-null; global scaling, absorption surface, per-batch geometry, learned rocking profile, per-frame B,
-outlier rejection, partiality cuts, and beam-stop masking are all dead on this data. The lever that ever
-worked is **2D profile-fit integration**; the keepers are **Gaussian profile-fit + rot3d + scale-fulls +
--R100** (≈ 0.61× anomalous, ISa 11.0).
-
----
-
-## Round 3 (2026-06-30): the profile width was wrong — de-contaminate + radial parallax ellipse
-
-Re-targeting **weak high-resolution** reflections (ligand-screening; metric = high-res CC1/2 / CCref, the
-only metric that can exceed XDS), a profile dump (per strong spot: the 13×13 background-subtracted,
-intensity-normalised grid + per-shell learned width, in the integrator's **own** units) showed the fitted
-Gaussian was **3–8× wider than the true spot**. The width came from the full-13×13 second moment, which is
-(a) **neighbour-contaminated** — the learning grid had no neighbour mask, catastrophic at low res where
-spots crowd near the beam (32× at shell 0) — and (b) inflated by **rectified background noise** in the far
-corners (lever arm `dx²+dy²` up to 72).
-
-Splitting the spot 2nd-moment into **radial vs tangential** (about the beam) separated two effects cleanly:
-the **tangential (isotropic)** width grows with resolution and saturates (mosaicity/divergence ∝ R_px),
-while the **radial excess is pure sensor parallax**, ∝ `tan²(2θ)`, giving σ_z ≈ 106 µm for 450 µm Si
-(physical: below the `t/√12 = 130 µm` uniform limit). The contaminated isotropic width was accidentally OK
-at high res (wide) but far too wide at low res — so simply narrowing it helped mid-res but *truncated* the
-radial parallax streak and hurt the weakest high-res shell.
-
-**The fix (shipped):**
-1. **De-contaminate the width** — measure the 2nd-moment over the signal disk **`r < r1`** (the same disk
-   the intensity is summed over), which inherently excludes neighbours and caps the radial tail. Gated to
-   the **monochromatic path**; the broadband (stills) path keeps the generous full-grid width (sparse
-   spots, the centroid-undersampling floor wants it). A neighbour-mask on the learning grid or an
-   iterative `r<kσ` window were considered but the disk cut is simpler and the iterative one *diverges*
-   at high res (the isotropic window chases the radial streak).
-2. **Radial parallax ellipse, always-on** — extend the (was stills-only bandwidth) ellipse with an
-   **analytic, material-aware** parallax term `c_par·tan²(2θ)`, `c_par = Var(z)/pixel²` from sensor
-   thickness + material + energy (`parallax_var_px2`; Si 450 µm → 0.66, JUNGFRAU Si 320 µm → 1.42, CdTe
-   tiny). Plus a fixed **weak-spot capture** term (`C_CAPTURE = 2.5`, monochromatic only) absorbing the
-   per-frame radial position scatter of weak high-res spots; the metric is a broad plateau (parc 2–4 ≈
-   equal) so it needs no tuning. On the broadband path the bandwidth streak supplies this, so it is omitted.
-
-**Result (HEWL rotation `fixed_master.h5` @ 1.0 Å, vs the round-2 baseline):** ISa **13.5 → 15.7**, CC1/2
-1.12 Å **91.3 → 95.9**, 1.05 Å **83.0 → 85.2**; external **CCref** band 1.00–1.19 Å **88.1 → 89.9**, CCxds
-**93.4 → 94.8** (per-shell CCref 1.00 Å 60→66, 1.05 Å 82→85.5); low/mid res flat, overall R-meas 9.4 → 8.7.
-Generalises: sharp serial stills (Jet8, JUNGFRAU 320 µm) get a small clean gain from the parallax term;
-broadband stills (Jet5, 1% DMM) are neutral (bandwidth dominates). The de-contamination is rotation-only
-because it *hurts* stills (their ellipse is on; tightening the seed starves it).
-
-This is the biggest single integration lever since box→Gaussian, and it confirms parallax is real and
-analytic.
-
-**Sub-pixel centring (added).** The profile is now built per reflection and centred on the **predicted
-sub-pixel position** (the learning/fit grid is otherwise binned about `round(predicted)`, mis-placing a
-shared profile by up to 0.5 px; the predicted position is noise-free geometry, unlike the observed
-centroid, which hurt). Small clean win on every metric — rotation ISa 15.7→16.2, CCref band 89.9→90.0,
-CCxds 94.8→95.0 (high-res 1.00 Å CCref 66.0→66.9), sharp serial stills 1.68 Å CC1/2 61.6→62.5, and the
-**anomalous S peak 0.92×→0.92× XDS** (SG_CYS/CL_CL nudge up) — i.e. it does not trade accuracy for it.
-The companion idea of **de-broadening** the learned width by the known 1/12 px²/axis integer-binning
-floor was tested and **rejected** (neutral-to-slightly-negative: the learned width was not over-broad
-enough to matter, so subtracting 1/12 just over-narrows).
-
-For the record, the combined round-3 lift on the **anomalous** accuracy benchmark (which the re-target had
-de-prioritised): S peak height **0.85× → 0.92× XDS** (pre-fix → profile-width fix → +sub-pixel), the
-closest to XDS yet — confirming the weak-high-res work raised genuine accuracy, not just precision.
-
-Open follow-up (deferred, untested): a neighbour mask on the learning accumulation as the root-cause
-alternative to the `r1` window (keeps this spot's genuine [r1,r2] tail; needs a multiplicity mask to
-separate self from neighbour in the overlap annulus).
-
-### Reciprocal-space profile width — `--reciprocal-profile` (tried and REMOVED 2026-07-02)
-
-**The idea.** The per-shell width is learned in **pixels**, so it changes ~4× with resolution (the
-geometric projection of a roughly-constant reciprocal-space relrod) and must be binned per shell, which
-**starves** at high resolution / on sparse data. Mapping the tangential 2nd-moment into **reciprocal
-space** (Jacobian `g_tan² = cos²2θ`) should remove the projection and give a transferable mosaicity law
-`σ²_q,tan = A + B|q| + C|q|²` (`C|q|²` = relrod variance `~(η|q|)²`), fit once per frame — like XDS/DIALS
-integrating in a reciprocal profile coordinate. It was shipped as an off-by-default dial after being
-tuned on the single sharp HEWL crystal, where it was metric-neutral.
-
-**Why it was removed — an 11-crystal mosaicity-stratified re-test (`/data/rotation_test`, off vs on vs a
-de-contaminated variant, plus a per-frame dump of the fitted widths) showed it is net-negative and cannot
-work per-frame.** Three linked, proven facts:
-1. **The `C|q|²` mosaicity term — the whole point — is unfittable per-frame.** The fitted curvature `a2`
-   came out ≈ 0 (often slightly *negative*) on **every** crystal, with **zero correlation** to the XDS
-   mosaicity (0.09→0.42°). Strong spots live at low `q` (q̄ ≈ 0.24–0.47) where `η²q²` is invisible; the
-   curvature only appears at high `q`, where there are ~0 strong spots to constrain it. The law degenerates
-   to a straight line.
-2. With `a2 ≈ 0`, the high-res width becomes a **blind `1/cos²2θ` extrapolation → 2–4× wider than
-   per-shell** (e.g. px²: pding4_001 1.07 flat → 2.0–3.0; cytC 1.10 → 2.7; MyoB 0.89 → 3.4). Meanwhile the
-   per-shell path *starves* at high res (n_strong<30 ⇒ flat global fallback that even **shrinks** vs
-   mid-res) — which turns out to be **accidentally correct**: weak, crowded high-res spots want a *narrow*
-   aperture, not the true (wide) spot shape.
-3. The over-wide profile pulls background into weak spots ⇒ **per-observation scatter (R-meas) rises**,
-   CC1/2 drops in reliable, high-multiplicity shells (pding4_001 mult 25, pding4_003 mult 12, EcwtCQ066,
-   MyoB); the mean intensities survive (CC-overall flat). A confirmatory cap (clamp the reciprocal width at
-   the widest well-sampled per-shell value) **recovers the regression** — direct proof the over-widening is
-   the harm. The only "improvements" seen (EcwtAL500, EP_cs_01-17 overall CC) were entirely in
-   noise-dominated shells (mult 2–3, CC<20 %, R-meas>100 %). **No crystal reliably wins.**
-
-The per-frame paradigm structurally cannot fit mosaicity, and un-starving the high-res width is net harm,
-so the dial (CLI flag, `BraggIntegrationSettings::reciprocal_profile`, and the per-frame fit block) was
-deleted. A real win would need cross-frame / rolling-window pooling to constrain `a2`, or a two-pass
-mosaicity seed from scaling — both against the per-frame / online + stills philosophy. Also dropped along
-the way: a separate **radial** relrod fit (radial is divergence + rotation-smear, not mosaicity → no gain).
-
-## Lineage: PixelRefine (removed 2026-06-25)
-
-`ProfileIntegrate2D` is the surviving half of an earlier experimental still-integrator, **PixelRefine**,
-which did three things bundled together: a reference-free **profile + extraction** (`J = Σ P(I−B)/v ÷
-Σ P²/v`, de-biased variance, plus a radial partiality) and a reference-based **per-image scaling**
-(`J ≈ G·B_DW·p·pol·I_ref`). The decomposition observation was that the first two are independent of the
-reference and identical in spirit to what a profile-fitting integrator does — so they were lifted out as
-`ProfileIntegrate2D`, while scaling stayed in the shared `ScaleOnTheFly` (which already does reference
-scaling when a reference is supplied: `IndexAndRefine::ReferenceIntensities` builds a
-`ScaleOnTheFly(experiment, reference)` and the per-image pass calls it on whatever the integrator
-produced — **integrator-agnostic**).
-
-**The decomposition won, so PixelRefine was deleted.** On the LysozymeJet5 serial stills, the default
-Gaussian path (`ProfileIntegrate2D` + the same reference scaling) **matched or beat** whole-PixelRefine
-in every per-shell CC½ (overall 95.7% vs 91.9%), ISa (1.6 vs 1.2) and R-meas (98.5% vs 175%), with
-CCref a tie — same indexing, same uniques. PixelRefine collected ~1.8× more observations (a wider Ewald
-band) but they were noisier and didn't help. Whole-PixelRefine was also incompatible with the rot3d
-combine (it emitted `image_scale_corr = 1/G` with no rotation partiality → garbage σ), which the
-decomposed integrator avoids by emitting `partiality` like `BraggIntegrate2D`.
-
-**Lessons inherited from PixelRefine (kept here so they survive the deletion):**
-- **Mean (not median) background** in both integrators — the median of a skewed background
-  under-subtracts and biases weak intensities positive (fake ⟨I/σ⟩ at CC≈0). Single biggest
-  "untrustworthy σ" cause.
-- **De-biased variance** (weight by background, not observed counts) — fixed negative high-res ⟨I/σ⟩
-  and scale collapse.
-- **A tight profile loses to a generous box for stills** — the ~0.4-px centroid-undersampling floor of
-  ~2×2 spots is a real sampling floor, not a fixable misprediction; recentring/adaptive-width add noise.
-  The lever is matched weighting + de-biased variance, never a tight aperture.
-- **Per-image R (profile-width) refinement is futile** — G↔R is intrinsically degenerate; refine G/B
-  freely but never R per image (it slides the per-image scale and wrecks the merge).