diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md
index e881e111..29032db2 100644
--- a/docs/CPU_DATA_ANALYSIS.md
+++ b/docs/CPU_DATA_ANALYSIS.md
@@ -403,9 +403,9 @@ Systematic absences are applied at least at the centering level (prior to full s
 
 ---
 
-## 9. 2D summation integration (three-ring method)
+## 9. 2D Bragg integration (profile fitting over a three-ring ROI)
 
-Jungfraujoch integrates predicted reflections by **summation** (no profile fitting), using a CrystFEL-inspired “three-circle / three-ring” method in the detector plane.
+Jungfraujoch integrates each predicted reflection in the detector plane over a CrystFEL-inspired “three-ring” region of interest (§9.1). The **default** extraction is **profile fitting** (Kabsch; §9.3), which weights each pixel by a fitted spot profile and so recovers weak reflections far better than plain summation; plain box summation (§9.2) is retained as the seed for the profile and as a fallback. Both methods share the same ROI and background model, and emit the same per-reflection $(I,\sigma,\text{partiality},d)$, so scaling, the rotation combine (§10.6) and merging consume either unchanged.
 
 ### 9.1 Regions of interest
 
@@ -422,7 +422,7 @@ Pixels are classified by their squared distance $r^2=(x-x_p)^2+(y-y_p)^2$:
 
 Invalid pixels (masked/bad/saturated) are excluded from both sums.
 
-### 9.2 Background subtraction and intensity estimate
+### 9.2 Box summation (seed and fallback)
 
 Let:
 - $S = \sum I(x,y)$ over signal pixels,
@@ -430,24 +430,28 @@ Let:
 - $B = \sum I(x,y)$ over background pixels,
 - $n_B$ = number of valid background pixels.
 
-Background per pixel:
+Background per pixel and integrated intensity:
 $
-\hat{b} = \frac{B}{n_B},
-$
-integrated intensity:
-$
-\hat{I} = S - n_S \hat{b}.
+\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.
 
-A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum (to avoid unstable background estimates).
+### 9.3 Profile-fitted extraction (default)
 
-### 9.3 Uncertainty model
+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:
 
-A Poisson-like estimator is used for the raw summed counts:
+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; radial/tangential anisotropy and per-detector-region profiles were evaluated and add nothing in the 2D detector plane (the spots are essentially round there — the real anisotropy lives in the discarded rocking direction).
+3. **Fit (Kabsch).** With profile $P$, background $B$ and the shell variance model, the intensity and its uncertainty are
 $
-\sigma(\hat{I}) \approx \sqrt{S},
+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,
 $
-with a minimum $\sigma\ge 1$ to avoid singular weights. (This is a pragmatic online estimate; more elaborate models may be applied downstream.)
+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.
+
+The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`).
 
 ### 9.4 Lorentz–polarization factor handling
 
@@ -611,7 +615,7 @@ Numerical quadrature over a scaled intensity variable is used to compute posteri
 
 ## 13. Practical notes and limitations
 
-- **No profile fitting** is currently performed for Bragg integration; all integration is summation-based (§9). This is appropriate for fast feedback and many serial/streaming use cases, but differs from full profile fitting workflows.
+- **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.