Document the integrator and the XDS-quality investigation

Rewrite NEXTGEN_INTEGRATOR.md from the early design draft into a full
record: the diagnosis (our intensities are accurate and unbiased per ANODE,
~2x noisier than XDS), the experiment ledger of what worked and the many
dead ends, the reusable tooling, and an honest assessment of the remaining
gap. Fix the stale ProfileIntegrate2D.h header (it is wired and the default
integrator, not a draft).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md

@@ -1,105 +1,282 @@
-# Next-generation integrator — design draft (2026-06-24)
+# Next-generation integrator & the XDS-quality investigation
 
-## The question
+_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)._
 
-Can PixelRefine's **profile-fitted intensity extraction** be separated from its **reference-based
-scaling**, so we get a profile-fitting integrator that needs no reference intensities — and could
-PixelRefine then be rebuilt as *(profile-fit integrator) + (smarter scaling)*?
+---
 
-**Yes.** They are already algebraically separate in the code.
+## TL;DR
 
-## Why we want it (the diagnosis)
+- **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.
 
-The residual ~4× gap to XDS on R-meas (jfjoch fulls 18.9% vs XDS 4.5%) was localized to the **2D
-spot integration**, not background/scaling/partiality (see `scaling-results-roadmap` memory / the
-session experiments):
-- the gap is in **strong** reflections (jf/XDS ratio 2.1× weak → 6.5× strong), a ~18% *multiplicative*
-  per-observation floor, **flat across resolution**, **per-reflection** (only 6.6% of 20.7% is
-  per-image/scaling), and **immune to box radius** (tight/default/wide all ≈ 18–24%).
-- => it is the **box-sum method** (uniform pixel weighting over a fixed disk), `BraggIntegrate2D`.
-  A fixed disk captures a *width-dependent fraction* of each spot; spot widths vary (mosaicity,
-  divergence, resolution, detector position, rocking phase) → ~18% multiplicative scatter.
-- swapping in PixelRefine's profile-fit swung ISa 6.4→16.5 (confirming the integrator is the lever,
-  though that path is currently broken — see below).
+---
 
-## The decomposition (PixelRefine today)
+## Philosophy & constraints (Filip)
 
-PixelRefine does three things; only the third needs the reference:
+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.
 
-1. **Profile (Term 2).** Measure the tangential width `R1 = sqrt(2·<eps_t²>)` per resolution shell
-   from the intensity-weighted 2nd moment of strong spots (`PixelRefine.cpp:433`). Reference-free.
-2. **Extraction.** `J = Σ Pₚ(Iₚ−B)/vₚ / Σ Pₚ²/vₚ`, `var(J)=1/Σ Pₚ²/vₚ`, with `Pₚ` the
-   area-normalized tangential profile (Gaussian of width R1) and the **de-biased** variance
-   `vₚ = max(B,1)` (`PixelRefine.cpp:491-525`, METHODS §1). It also computes the **radial partiality**
-   `p = P_radial` and a detector-clipping **completeness**. Reference-free.
-3. **Scaling (Term 1).** Fit per-image `G,B` from `J ≈ G·B_DW·p·pol·I_ref` (`PixelRefine.cpp:365`,
-   METHODS §4). **This is the only step that uses I_ref.**
+---
 
-So **integrator = steps 1+2** (output `I=J/(p·pol)`, `σ`, `partiality`), **scaling = step 3** (or the
-existing reference-free `ScaleOnTheFly`). A clean split.
+## The pipeline today
 
-This split also **fixes the rot3d-combine incompatibility**: PixelRefine-as-a-whole emits
-`image_scale_corr=1/G` with no rotation partiality, which is why `-r pixelrefine -P rot3d` produced
-garbage σ (CCref 94→84, ⟨I/σ⟩→1.3). The *decomposed* integrator emits `partiality` like
-`BraggIntegrate2D`, so `ScaleOnTheFly` + `Combine3D` + merge consume it unchanged.
+```
+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)
+```
 
-## The new-generation integrator
+Relevant knobs (`jfjoch_process`, 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).
 
-A new module `ProfileIntegrate2D`, **swappable with `BraggIntegrate2D`**, same output
-(`vector<Reflection>` with I, σ, partiality, d, …), **no reference**, no scaling. Algorithm:
+### Current integrator: `ProfileIntegrate2D`
 
-1. **Pass A (rough + learn the profile).** Box-sum every reflection (reuse `IntegrateReflection`)
-   to get a rough I and the observed centroid; select strong spots (significance ≥ 5).
-2. **Learn an EMPIRICAL profile**, per *(detector region × resolution shell)*, by averaging the
-   background-subtracted, intensity-normalized, centroid-aligned pixel grids of the strong spots.
-   **This is the key difference from PixelRefine**, which uses a *Gaussian* of width R1×6 — a model
-   that has the *same width-mismatch* as a box (and ×6 ≈ near-box, which is why PixelRefine's stills
-   ISa was only 1.1). A *learned* profile matches the real spot — including the 0.4-px undersampling
-   smear, asymmetry, and position dependence — so the weighting is genuinely optimal. (XDS learns its
-   profile over 9 detector regions × resolution; this is that idea.)
-3. **Pass B (profile-fit extract).** For each reflection, `I = Σ Pₚ(cₚ−B)/vₚ / Σ Pₚ²/vₚ` with the
-   learned `P` for its region/shell, de-biased variance `vₚ = B + max(I,0)·Pₚ` (Kabsch; iterate a
-   couple of times), `σ = sqrt(1/Σ Pₚ²/vₚ)` corrected for the signal term. Carry the rotation
-   `partiality` exactly as `BraggIntegrate2D` does (so rot3d works).
+Drop-in for `BraggIntegrate2D`, same `vector<Reflection>` output, **no reference intensities**.
+Selected at `IndexAndRefine.cpp:347`. Per frame:
 
-**Why it removes the floor:** profile-weighting integrates *out to the tail* (capturing the
-width-varying part a fixed disk clips) while *downweighting the noisy edge* (avoiding the wide-box
-background penalty that made radius 6,9,12 worse). A *learned* profile (not a Gaussian) removes the
-model-mismatch that limited PixelRefine. The floor is multiplicative and in strong reflections → this
-is the right lever.
+- **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).
 
-**Caution (from PixelRefine):** a *tight* profile LOST to a generous box for stills (the 0.4-px
-centroid undersampling floor). An empirical profile is immune to this *because it is measured from the
-undersampled spots themselves* — the template already contains the smear. The lever is *matched
-weighting + de-biased variance*, never a tight aperture.
+**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).
 
-## Plug-in (no touch to rot3d / scaling / merge)
+**Measured, do not re-try:**
+- **Anisotropy / per-region profiles add nothing.** Radial/tangential anisotropic Gaussians and
+  per-detector-region profiles were tested → neutral-to-slightly-worse. The 2D detector-plane spots are
+  ~round; the real anisotropy is in the **rocking** direction, which per-frame 2D integration discards.
+  Keep the **isotropic per-shell** width.
+- **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.)_
 
-- New `ProfileIntegrate2D(experiment, image, predicted, …) -> vector<Reflection>`, output-compatible
-  with `BraggIntegrate2D`.
-- Select it like the geometry-refinement choice (a `BraggIntegrationSettings` enum or an
-  `IntegrationAlgorithm`), default box-sum. Everything downstream (ScaleOnTheFly, Combine3D, Merge)
-  is unchanged — they only see `Reflection{I,σ,partiality,d}`.
-- Reusable by **stills and rotation** (and the FPGA path later, if the profile is precomputable).
+### The rot3d combine (`Combine3D`)
 
-## Staging
+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`.
 
-- **v1** — standalone profile-fit with the *measured-R1 Gaussian* (= PixelRefine steps 1+2 lifted out,
-  no reference). Cheap; reproduces PixelRefine extraction and proves the decomposition + the
-  rot3d-combine compatibility. Likely only a modest ISa gain (Gaussian model).
-- **v2** — **empirical learned profile per resolution shell.** The expected real win.
-- **v3** — empirical profile per *detector region × shell* (XDS-grade); anisotropic.
+- **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.
 
-## Validation
+**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.**
 
-A/B per-observation **R-meas / ISa** vs box-sum, both with rot3d, on crystal 2 (`fixed_master.h5`)
-against `XDS_ASCII.HKL` (4.5%); repeat treating it as stills (crystal 1 jet) where box-sum's hidden
-floor is ISa 1.1. Use `--dump-observations`. Target: close the 18% strong-reflection floor toward
-XDS's ~3%.
+### 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>`** (`jfjoch_process`): `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).
+
+---
 
 ## "Could PixelRefine be rebuilt as integrator + scaling?"
 
-Yes: `ProfileIntegrate2D` (steps 1+2) + a scaling routine. The scaling routine can be the existing
-reference-free `ScaleOnTheFly`, or a reference-driven one (Term 1) for the EM/self-seed loop — both
-just consume `Reflection{I,σ,partiality}`. PixelRefine then becomes *= ProfileIntegrate2D + reference
-ScaleOnTheFly*, and the experimental coupling (and its rot3d breakage) disappears.
+Yes, and this is the clean decomposition that `ProfileIntegrate2D` realizes. PixelRefine does three
+things; only the third needs a reference:
+1. **Profile** — measure the tangential width per shell from strong spots (`PixelRefine.cpp:433`).
+   Reference-free.
+2. **Extraction** — `J = Σ P(I−B)/v ÷ Σ P²/v`, de-biased variance, plus radial partiality
+   (`PixelRefine.cpp:491-525`). Reference-free.
+3. **Scaling** — `J ≈ G·B_DW·p·pol·I_ref` (`PixelRefine.cpp:365`). **The only step using the reference.**
+
+So **integrator = steps 1+2** (= `ProfileIntegrate2D`), **scaling = step 3** (or the reference-free
+`ScaleOnTheFly`). The decomposition also **fixes the rot3d incompatibility**: whole-PixelRefine emits
+`image_scale_corr = 1/G` with no rotation partiality (which gave garbage σ with `-P rot3d`), whereas the
+decomposed integrator emits `partiality` like `BraggIntegrate2D`, so `ScaleOnTheFly` + `Combine3D` +
+merge consume it unchanged. PixelRefine then becomes `ProfileIntegrate2D + reference ScaleOnTheFly`, and
+the experimental coupling disappears.

image_analysis/bragg_integration/ProfileIntegrate2D.h

@@ -4,7 +4,7 @@
 #pragma once
 
 // =============================================================================
-// ProfileIntegrate2D — profile-fitting 2D integrator (DRAFT / not yet wired in)
+// ProfileIntegrate2D — profile-fitting 2D integrator (the DEFAULT integrator)
 // =============================================================================
 //
 // A drop-in alternative to BraggIntegrate2D that replaces uniform box summation with
@@ -34,9 +34,9 @@
 // decomposition + rot3d-combine compatibility); v2 = empirical per-shell (the expected win);
 // v3 = empirical per detector-region x shell (XDS-grade).
 //
-// EXPERIMENTAL: selected by BraggIntegrationSettings::Integrator (ProfileGaussian = v1,
-// ProfileEmpirical = v2); jfjoch_process exposes it as `--integrator gaussian|empirical` (default
-// box-sum). For A/B vs XDS_ASCII.HKL via --dump-observations. Not yet exposed in the OpenAPI / viewer.
+// Selected by BraggIntegrationSettings::Integrator: ProfileGaussian (the DEFAULT, v1) or
+// ProfileEmpirical (v2); BoxSum (BraggIntegrate2D) is the fallback. jfjoch_process exposes it as
+// `--integrator boxsum|gaussian|empirical`. For A/B vs XDS_ASCII.HKL via --dump-observations.
 // =============================================================================
 
 #include <vector>