Jungfraujoch/docs/CPU_DATA_ANALYSIS_INTEGRATION.md

Data analysis: integration, scaling and merging (§8–§12)

Part of the CPU/GPU data-analysis reference; the section numbers are continuous across its four parts.

:local:
:depth: 2

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.

8.1 Enumerating reciprocal lattice points

For a maximum resolution d_\mathrm{min}, accept (h,k,l) such that: $ \lVert \mathbf{p}(h,k,l)\rVert^2 = \lVert h\mathbf{a}^* + k\mathbf{b}^* + l\mathbf{c}^*\rVert^2 \le \left(\frac{1}{d_\mathrm{min}}\right)^2. $

8.2 Still prediction (excitation-error cutoff)

For still images, the diffracting condition is approximated by an excitation-error cutoff: $ \left|\Delta_\mathrm{Ewald}(\mathbf{p})\right| \le \Delta_\mathrm{cut}. $ 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:

• rotation axis unit vector \mathbf{m}_2,
• \mathbf{S}_0 incident vector,
• \mathbf{S}(\phi)=\mathbf{S}_0+\mathbf{p}(\phi).

A key quantity is: $ \zeta = \left|\mathbf{m}_2\cdot \mathbf{e}_1\right|,\quad \mathbf{e}_1 = \frac{\mathbf{S}\times \mathbf{S}_0}{\lVert \mathbf{S}\times \mathbf{S}_0\rVert}, $ which also appears in XDS as the Lorentz component linked to the rotation axis.

A Gaussian mosaicity model yields a partiality fraction over an oscillation width \Delta\phi:

P(\phi;\sigma_M,\zeta,\Delta\phi) = \frac{1}{2}\left[\mathrm{erf}\!\left(\frac{\phi+\Delta\phi/2}{\sqrt{2}\,\sigma_M/\zeta}\right) - \mathrm{erf}\!\left(\frac{\phi-\Delta\phi/2}{\sqrt{2}\,\sigma_M/\zeta}\right)\right],

with mosaicity \sigma_M in radians.

Reflections are predicted if they meet minimum \zeta and mosaicity-window criteria, and their predicted detector coordinates fall on the active detector area.

8.4 Systematic absences (centering)

Systematic absences are applied at the centering level (prior to full space-group symmetry) when the space group is supplied by the user. With no user-fixed space group, prediction runs in P regardless of the centering the lattice search inferred: the centering-absent reflections are integrated so that the space-group search (§13) can confirm or disprove the centering from the measured intensities, and so that a missed superstructure shows up. For centering symbol C:

• I: absent if h+k+l odd,
• A: absent if k+l odd,
• B: absent if h+l odd,
• C: absent if h+k odd,
• F: absent if any of h+k, h+l, k+l is odd,
• R: absent if (-h+k+l)\bmod 3 \ne 0,
• P: no centering absences.

---

9. 2D Bragg integration (profile fitting over a three-ring ROI)

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

For each predicted reflection at (x_p,y_p), define three radii:

• r_1: inner signal radius,
• r_2: inner background radius,
• r_3: outer background radius.

The defaults are 4,6,13 px for rotation data and 6,8,14 px for stills, which have a sparser pattern and can afford the wider ring. --integration-radius sets them by hand; on rotation data r_1 is otherwise measured from the crystal's own spots (§9.5).

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. 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. (Both the annulus and that exclusion become ellipses when the option below is used; with it off, which is the default, they are the circles just described.)

Radially elongated background ring (opt-in, --integration-stencil <k>, default 0). The three radii above are one triple for the whole run, identical for every reflection at every resolution. A reflection is not round, though: a finite bandwidth streaks it radially by \sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}, with R_\mathrm{px} the distance from the beam centre — the same physical smear as §8.2's and §11.1's \sigma_\mathrm{bw}, expressed here in detector pixels where those sections use reciprocal units; the two forms are never mixed in one formula. Throughout, \text{bandwidth} is the rms relative energy spread: the user-facing --bandwidth takes a FWHM (a DMM's usual specification) and it is divided by 2.355 on input. On a radially smeared spot the fixed 6\ldots13 px ring therefore sits only $\approx1.3$–2.2 radial \sigma from the centre — on the reflection's own tails, which it then measures as background.

With k>0 the background ring becomes an ellipse, elongated along the beam→reflection direction by k\sigma_\mathrm{bw}. The radial semi-axes become r_2+k\sigma_\mathrm{bw} and r_3+k\sigma_\mathrm{bw}; the tangential half-widths stay r_2 and r_3; and the growth is capped at 2r_3, which bounds what a mis-declared bandwidth can do to the bounding box. Pixels are then classified as

• signal region: r^2 < r_1^2 — a circle, unchanged,
• background ring: r^2-q_\mathrm{in}\rho^2 \ge r_2^2 and r^2-q_\mathrm{out}\rho^2 < r_3^2,

where \rho is the pixel's radial offset (its projection on the beam→reflection direction), g=\min(k\sigma_\mathrm{bw},\,2r_3) is the capped growth, and q=1-\big(r/(r+g)\big)^2 for the boundary concerned. Written this way k=0 gives q=0 and both tests collapse onto r^2 exactly in floating point, so the default classifies every pixel exactly as the circular stencil did. The neighbour exclusion above follows: each neighbour's inner ellipse, taken in that neighbour's own radial frame, is what is masked out of this reflection's ring.

The width is the bandwidth streak alone, and deliberately not the profile's full radial variance of §9.3, which also carries the sensor parallax and weak-spot capture terms. Those two are the only terms there are on a monochromatic beam, and widening the ring by them was measured on the rotation battery: it neither helped the crystals with clean high-resolution shells nor left the weak ones alone. The bandwidth streak, by contrast, is a measured elongation of the recorded spot — principal axis along the radius to within a couple of degrees, and azimuth-independent. Keeping only it also makes the option exactly inert on a monochromatic beam, where \sigma_\mathrm{bw} is zero.

Growing the ring also grows the neighbour exclusion, so on a crowded pattern fewer background pixels survive; a reflection left with too few is rejected outright. On the data this was measured on the loss is under 0.1% of reflections, but it is not structurally zero.

Only the ring moves. The signal disk r_1 stays circular, deliberately: it sets n_S, it sets \mathrm{var}(\hat b), it is the domain the profile width is learned over (§9.3), and with --integrator boxsum it drives the all-or-nothing “every signal pixel valid” acceptance gate (§9.2), so growing it would reject any box sum carrying a single bad pixel anywhere along a long streak. In the default gaussian mode r_1 does not set the intensity at all — the fit grid, \lceil r_2\rceil, does.

What a circular r_1 loses is flux, and that loss is not a function of resolution alone: measured per reflection, it carries a directional component worth several Ų with a definite principal axis, on top of the isotropic part. Nor is there anything in the merge to absorb it. There is no per-shell scale, and there cannot usefully be one: every scale in §10 is fitted against a reference built from a reflection's own symmetry equivalents, and equivalents share s^2 exactly, so any function of s^2 lies in the exact null space of the whole scaling model — a per-shell parameter would have zero residual to fit against. (XDS and DIALS have the same null space, for the same reason.) The isotropic part of the loss is instead degenerate with the overall Wilson B and is silently reported as part of it, so the reported WILSON_B / _reflns.B_iso_Wilson_estimate carries an $r_1$-dependent contribution: measured across a constant-ring-area radius sweep it falls monotonically as the disk grows, by 0.5 Ų on sharp strong data and by up to ~10 Ų on weak wide-spot data. What this costs the data is much less than what it costs the flux, because most of the loss is matched by a proportional \sigma: it moves no CC$_{1/2}$ and no R_\text{meas}, and — to within a few hundredths of an ångström — no resolution cut.

9.2 Box summation (seed and fallback)

Let:

• S = \sum I(x,y) over signal pixels,
• n_S = number of valid signal pixels,
• B = \sum I(x,y) over background pixels,
• n_B = number of valid background pixels.

Background per pixel and integrated intensity: $ \hat{b} = \frac{B}{n_B},\qquad \hat{I} = S - n_S \hat{b}, $ with a Poisson-like uncertainty \sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big), i.e. \sqrt{S} floored both at 1 count (pixel values are photon counts) and at a small fraction r_\sigma of the intensity. The second term under the root is the uncertainty of the background estimate itself: \hat b is measured from a finite number of ring pixels, \mathrm{var}(\hat b)=\hat b/n_B, and it is subtracted n_S times over, so it enters squared. Omitting it understates the variance by 1+n_S/n_B — 1.11 with the shipped circular stencil (n_S = 45, n_B = 408) — and so understates \sigma by up to \sqrt{1+n_S/n_B} \approx 1.05, a bound attained on background-limited (weak) reflections and falling towards 1 on strong ones, where S dominates; with an elongated ring n_B grows with resolution, so the factor is no longer one number for a run. The same term is carried into the profile fit (§9.3), where it adds \big(\sum P/v \,\big/ \sum P^2/v\big)^2\,\mathrm{var}(\hat b) — the square of \partial I/\partial\hat b for that fit; n_B is the count of pixels behind the final background value, so a clip or trim that discards ring pixels raises it. A box sum is accepted as “observed” only if all signal pixels were valid and n_B exceeds a minimum — it measures what is in the disk with no model of what should be there, so it cannot renormalise a disk it has lost pixels out of. The profile modes can, and do (§9.3). This box sum is the classical estimator; it is used directly with --integrator boxsum, and otherwise seeds the profile fit below, where S and n_S then count only the pixels that were actually read.

High-side clipped background (default on). Because \hat{I}=S-n_S\hat{b} is a small difference of large numbers for weak reflections, a per-pixel background bias \delta\hat{b} becomes a fractional intensity bias \approx n_S\,\delta\hat{b}/\hat{I} that grows as \hat{I} shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above \hat{b}+n\sqrt{\hat{b}} are rejected and the mean recomputed, with n=4 (--background-clip; n=0 disables), lowered by rugnux to n=3 on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. That is only a default — the flag sets n whatever the bandwidth is. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at 4\sigma), while a 40-pixel neighbour core at +100 counts shifts the estimate by +0.009 ct/px.

The clip cuts only the high tail, which matters: the symmetric trimmed mean it replaced (drop the lowest and highest fraction f of ring pixels, f=0.10; still reachable with --background-trim, which switches the clip off) is not a consistent estimator of the mean of a right-skewed Poisson sample. It sits \approx0.1 ct/px below the true mean at every level, and with n_S\approx50 signal pixels in the r_1 disk that under-estimate adds \approx5 counts to every partial — negligible at low resolution, but a large fraction of a partial in the outermost shell. The trim also collapses once contamination exceeds \approx10\,\% of the ring, where the clip does not. Note that removing a positive background bias lowers \langle I/\sigma\rangle and raises edge R_\text{meas}, because both are inflated by information-free counts — so neither may be read as evidence against the change.

Both estimators are computed in the shared background pass, but only the trim reaches plain box summation: the high-side clip is skipped for --integrator boxsum, which therefore uses the plain ring mean unless --background-trim is given.

Radial background correction (opt-in). A ring mean estimates the background under the signal disk correctly only if the background is flat there. The signal disk and the ring are concentric, so for a background that is linear in position \langle B\rangle_\mathrm{ring}=\langle B\rangle_\mathrm{disk} identically — a plane or gradient fit buys exactly nothing. The leading error is the curvature of the radial background, which is negligible on a smooth background but reaches tens of counts on a single reflection sitting on a sharp powder ring. That error is a kernel over radial offset, $ \delta \hat b ;=; \textstyle\sum_k \kappa_k, \bar B(r_0+k), $ with \kappa the annulus-minus-disk histogram of the stencil over radial offset, averaged over azimuth, and \bar B(r) the image's own radial background curve. With the fixed circular stencil (k=0, §9.1) that single kernel serves every reflection. An elongated ring does not: its radial-offset histogram depends on how far that particular reflection's ring was grown, so \kappa becomes a small table of kernels, indexed by the growth rounded to whole pixels. The azimuthal average survives the change unaltered, because the stencil is rebuilt in the reflection's own frame at each azimuth and so stays radially aligned: what is averaged over is the sub-pixel phase of the detector grid against the radius, which is what genuinely differs between reflections. Applying it costs one short dot product per reflection and no extra pixel reads; correcting the background scalar means the box sum, the profile fit and the variance all pick it up. The curve is accumulated from the same annulus pixels the background pass already reads (a pixel's radius is the reflection's radius plus the pixel's projection on the beam→reflection direction, so no per-pixel square root is needed) and specifically from the clipped pixels, or it would carry neighbour tails and zingers — which is why the correction is inert under --integrator boxsum, that path having no clip pass.

The model is a function of radius alone, so it is applied only where that is true of the background. --background-radial takes on, off or auto. It is off by default; under auto each image's peak-excluded ice score (§3.3) is taken after spot detection and before integration, and the correction is applied to that image when the score reaches the same --ice-min-score gate. Smooth powder ice is a radial feature and is corrected; ice made of discrete crystallite spots — which the profile channel is blind to and the spot channel catches — leaves no smooth ring to model, and correcting it makes matters worse. Measured against a fixed atomic model, comparing ice bands with resolution-matched decoy bands carrying no ice: on a crystal with pure smooth ice the correction removes 43 % of the bands' excess amplitude, and the improvement is 7× larger inside the bands than outside, which is its stated mechanism; on a crystal whose ice is textured the same correction increased the excess amplitude by half; on a clean crystal it is inert to four decimal places. Auto engages only where a peak-excluded score exists (adaptive spot finding, §3.2) — a plain azimuthal profile carries the Bragg peaks and cannot support an absolute threshold, so without one auto leaves the correction off.

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. For gaussian (the default) the width is taken per resolution shell from the measured second moments of the strong spots (shell-dependent because spot size grows with resolution). The moments are anisotropic: each strong spot's pixels are rotated into its own radial/tangential frame before being accumulated, giving \sigma^2_r and \sigma^2_t separately. Stacking the spots in the detector frame instead — they sit at every azimuth — averages the two directions away, leaving only \sigma_r^2+\sigma_t^2, so radial smearing is read back as a wider tangential spot. For empirical the profile is instead the averaged, background-subtracted pixel grid of the shell's strong spots, accumulated in the detector frame on their rounded predicted positions. For gaussian only, 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 a small fixed weak-spot capture term.

   The two enter as a floor on the measured radial excess: \sigma^2_\mathrm{radial}=\sigma^2_t+\max\!\left(\sigma^2_r-\sigma^2_t,\ \sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)\right), tangential unchanged at \sigma^2_t. The measured excess is what the signal disk can resolve; the analytic term takes over for a streak too long to be measured there. The Gaussian is built on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add. The empirical profile keeps the fixed per-shell grid and gets none of this.

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 = \max!\left(B + I,P,\ \tfrac{1}{2}B\right), $ 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). The plug-in I enters as it is: half-wave rectifying it, v=B+\max(I,0)P, lets v — and with it the reported 1/\sum P^2/v — respond only to upward fluctuations of a noisy estimate, which adds \approx0.4\,\sigma\sum P^3/(\sum P^2)^2 to every \sigma whatever the count rate. That offset is invisible on strong reflections and a large fractional inflation on weak ones; the \tfrac12 B clamp keeps v positive without reintroducing it. 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 background term is floored at 0.01 ct/px — enough to keep P^2/v finite when the ring mean reads exactly zero, which a ring of n_B pixels cannot distinguish from any background below \approx1/n_B. The rotation/excitation partiality is carried exactly as in the box-sum path.

Pixels the fit cannot use (MINPK). A profile fit is the amplitude of a normalised profile, so a pixel left out of the sum renormalises the estimator by construction: it costs information — \sum P^2/v shrinks and \sigma grows — but biases nothing. That is what keeps a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload: those pixels are simply not read, and the fit is taken over the rest, exactly as the shared pixels of a crowded reflection are (--overlap exclude). The reflection is kept only while enough of the expected profile survives — at least --overlap-minpk of the profile mass that falls on the detector at all, default 0.75, which is XDS's MINPK and dials' valid_foreground_threshold. The complete reflections alone teach the profile, its resolution shells and their widths. --integrator boxsum has no profile to renormalise with and keeps the all-or-nothing rule of §9.2.

"Biases nothing" holds only while the profile model is exact. Lose the peak and the amplitude is set by the wings alone, so the result stops being a measurement of the reflection and becomes a measurement of how well the fitted shape describes it. The worst case is a pixel invalidated by the flux it saw — a detector's per-frame overload marker: that pixel goes missing because the reflection was bright, so the loss concentrates on the strong low-resolution reflections that are the largest terms of R_\mathrm{meas}, where the fit reads -50\% against the symmetry mates. MINPK cannot catch it, because it cuts on profile mass and the peak of a broad spot is a few percent of the mass. So a second condition applies alongside it, on any unreadable pixel whatever made it unreadable: no unreadable pixel may carry more than 0.9 of the profile's own peak value. As a fraction of the peak rather than a radius in pixels, that scales with the spot — for a Gaussian it is a cut at \sqrt{-2\ln f}\,\sigma = 0.46\sigma, the peak pixel alone where \sigma is 0.8 px and the crest of the ridge where the profile is a bandwidth streak — and it costs well under 0.1 % of the recovered observations.

The integrator is selected by --integrator boxsum|gaussian|empirical (default gaussian).

9.4 Lorentz–polarization factor handling

For integrated reflections, polarization correction can be applied as a multiplicative correction to the reflection scale via the geometry-based polarization term (§2.2). A Lorentz-like factor is carried as rlp in predictions, and used during scaling/merging (§10).

---

9.5 Choosing the signal radius from the crystal's own spots (rotation)

The three radii are one triple for the whole run, but on rotation data they are no longer a fixed constant: r_1 is measured from how wide this crystal's spots actually are (--adaptive-integration-radius, on by default for rotation, off for stills, ignored when --integration-radius is given).

Why r_1 matters even though it does not set the intensity. In the default gaussian mode the intensity is a profile-fit amplitude over the grid \lceil r_2\rceil (§9.3), so r_1 is not the integration domain. It is the aperture the profile width is learned over, and a second moment taken over a disk of radius a saturates at a^2/4. At r_1 = 4 the learned \sigma can therefore never exceed 2 px, and a crystal whose spots are broader than that is fitted with a profile the model cannot represent.

The measurement is independent of the integrator. It is made in the pre-scan, on the frames the beam-stop projection already reads, so it costs no extra frame reads and there is no feedback loop. On the spots the spot finder has already found, a spot is used only if it is clear of the detector edge and of the direct beam, has no neighbouring spot within 28 px, is one of the 40 strongest in its resolution band, sits on a fully readable disk, and reaches a signal-to-noise of 15 with its centroid within 2 px of the found position. For each surviving spot the background-subtracted encircled-flux curve is accumulated in 1-px annuli out to a fixed 14 px aperture and normalised at 8 px — an aperture that owes nothing to r_1, r_2 or r_3.

Pooling. Spots are stratified into five resolution bands (2–3, 3–4.5, 4.5–7, 7–12, 12–30 Å), because a weak crystal's strongest spots sit at high angle and a strong one's at low angle. Each band with enough members contributes the radius at which its median curve reaches 0.80 of its normalised flux — r_{80} — at the band's median d. Those points are fitted by weighted least squares against 1/d (the mosaic contribution to the detector footprint grows as 1/d) and evaluated at a common 5 Å, then clamped to the range the bands actually measured so the fit never extrapolates.

The radius.

r_1 = \mathrm{clamp}\!\left(\mathrm{round}(2\,r_{80}),\ 4,\ 6\right),\qquad r_2 = r_1 + 2,\qquad r_3 = \sqrt{r_2^2 + 133}

The factor 2 is not fitted: for a Gaussian r_{80} = 1.794\,\sigma, so r_1 = 2r_{80} = 3.59\,\sigma, where the truncated second moment recovers 0.990 of \sigma^2. The expression for r_3 holds the background-ring area constant at its value for the shipped 4,6,13 (13^2 - 6^2 = 133) — a ring that shrank with the disk is what makes a bare --integration-radius worse than the default it replaces. The floor of 4 is that shipped default; the ceiling of 6 is pattern density, since r_2 also drives the neighbour-ownership radius and the ring's inner edge. At r_1 = 4 the triple is bit-for-bit the shipped default, so a crystal with ordinary spots is left exactly where it was.

The sample grows until the answer settles. The frames are measured in tiers of stride 8, 4, 2, 1, each tier's sample strictly containing the previous one, and the pooling is redone after each. Measuring stops when the new r_{80} is within 0.40 px of what the smaller sample said and is at least 0.25 px clear of both radii at which the rounding in r_1 changes answer. Both conditions are load-bearing: clearance alone lets a small sample settle across a switch, and the step test alone lets it settle on one. Every frame of the sample is still read — the beam-stop mask and the beam centre are unchanged; what the tiers save is the decompression, preprocessing and spot finding the width measurement adds on top of the read.

It applies to the final pass only. A rotation run integrates twice (§7.5), and the widened radius is handed to the canonical second pass, not to the geometry pre-pass. The reason is that post-refinement takes its observed positions from the integrator, and an observed position is a first moment over the signal disk with the background still in it: a flat background adds nothing to the numerator but adds n\,b to the denominator, so every measured offset is pulled toward its prediction by I/(I + n b), and n nearly doubles between r_1 = 4 and r_1 = 6. A wider disk therefore under-corrects the geometry — enough, on a crystal whose metric is half a degree off orthorhombic, to flip the second pass's de-novo Bravais choice.

The density guard. Widening r_1 pushes r_2, the ring's inner edge, into the neighbours; a reflection whose ring is left with five or fewer clean pixels has no background and is dropped whole. The integrator counts these, and separates the ones lost to neighbouring reflections from the ones lost to the detector itself (module gaps, the beam stop, the resolution mask) — a floor that reaches a couple of percent on some geometries and does not move with r_1. Where the neighbour-driven loss exceeds 1.13 % of the predicted reflections, the pattern is too dense for the widened radius and the final pass is integrated again at the fixed 4,6,13, reported as pass 3 of 3 with the reason in PASS_DECISION.

Every integration pass, adaptive or not, now logs the radii it used together with the fraction of predicted reflections that lost their background ring, split into the neighbour and detector parts, and the profile-fit fallback rate.

---

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, and handles both still and rotation data.

10.1 Observation model

For an observation j of a unique reflection h on image (or image group) i, the predicted measured intensity is modeled as: $ I_{ij} \approx G_i , L_{ij}, P_{ij}, I_h, $ where:

• G_i is the image scale factor,
• L_{ij} is a Lorentz-like / geometry factor; predictions carry its reciprocal as rlp, so L = 1/\texttt{rlp} and the correction below is applied as a multiplication by rlp,
• P_{ij} is a partiality term (model-dependent),
• I_h is the merged (true) intensity parameter for that unique reflection.

A least-squares objective is minimized: $ \sum_{ij} \left(\frac{I_{ij}^{\mathrm{pred}} - I_{ij}^{\mathrm{obs}}}{\sigma_{ij}}\right)^2 $ 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

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), 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]. $ Here \Delta\phi_{ij} is observation $j$'s rocking offset from its exact Bragg angle on image i, and the unsubscripted \Delta\phi is the oscillation width per frame — two different quantities that share a letter. 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. Unity (P_{ij}=1): used for the scale-on-fulls refit (§10.6), where each observation is already a complete reflection.

3. Fixed: use the per-reflection partiality carried from prediction. Still/serial images are predicted with P=1, so a single-pass stills scale is effectively unity/fixed — which is exactly what --simple-stills keeps. By default the stills path instead post-refines a physical partiality: a small per-crystal orientation tilt (\delta\psi_x,\delta\psi_y) about the two axes perpendicular to the beam is refined against the running merge, and every reflection's partiality is then recomputed analytically from the refined lattice through its excitation error \Delta_\mathrm{Ewald}=\big|\,|\mathbf{q}+\mathbf{S}_0|-1/\lambda\,\big| and a Gaussian width \sigma^2=\gamma_0^2+(\gamma_e d^*)^2+(\mathrm{bw}\,|q_z|)^2 — the reciprocal-lattice point's own radius (resolution-independent), the mosaic/divergence spread, and the bandwidth smear along the beam, in quadrature. The fit typically drives \gamma_e\to0, leaving the resolution-independent \gamma_0 as the effective width. A tilt moves reflections on opposite sides of the Ewald sphere in opposite directions, so it reshapes the spatial pattern of partialities — a degree of freedom the per-image scale G does not have, and the reason the tilt is refined rather than a scalar partiality width, which would be degenerate with G. Nothing is re-integrated (the integrated intensities are fixed); the tilt is hard-bounded at about 1° and held by a soft prior, so it stays inert on sparse or weak crystals. The cycle is merge → per-crystal tilt refinement (with G profiled out by the same robust Cauchy IRLS used for the per-frame scales, §10.3) → recompute P → re-merge, repeated a few times.

Reflections below a minimum partiality can be rejected from merging to avoid unstable corrections.

10.3 Smoothing of per-frame scales

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.

The crystal orientation is smoothed the same way, and for the same reason. Geometry is re-refined independently on every frame against that frame's spots alone — as few as a dozen on a sparse crystal — so the per-frame orientation carries a real slow drift (crystal slippage, up to ~1.3° across a sweep) on top of fit noise that scales with spots per frame. Before scaling, the per-frame lattices are de-rotated to a common reference, averaged in frame order, rotated back, and every partial's \Delta\phi — hence its partiality — is recomputed from the smoothed lattice. The window is chosen per dataset by leave-one-out cross-validation (does a frame's neighbours predict its orientation?) rather than fixed, because drift and noise both vary by two orders of magnitude between crystals; it is capped, because the per-frame fit also absorbs a real per-frame systematic that smoothing too wide destroys. Only frames that actually indexed take part: a frame that did not carries an all-zero lattice, which is finite and so passes a validity check written as a finite test, and would otherwise be both averaged into its neighbours' orientation and scored in the cross-validation that picks the window. Refining less is not an alternative: with per-image refinement off the space group is lost on several crystals.

A per-frame scale enters every intensity as 1/G, so a frame whose fit is not determined by its data can amplify it without bound — and \sigma is amplified by the same factor, which makes it invisible to any $\sigma$-based outlier test. A fitted G far below the run's median is therefore treated as undetermined rather than as a successful fit, both here and in the separate refit on the combined fulls (§10.6). The bound is a ratio to the run's own median because G is not gauge-fixed: it and the merged means have an exact global multiplicative degeneracy, so no absolute value is meaningful.

10.4 Merging estimator

After refinement, corrected observations are formed: $ I^{\mathrm{corr}}{ij} = \frac{I^{\mathrm{obs}}{ij}}{G_i L_{ij} P_{ij}},\qquad \sigma^{\mathrm{corr}}{ij} = \frac{\sigma^{\mathrm{obs}}{ij}}{G_i L_{ij} P_{ij}}. $

Unique intensities are merged by inverse-variance weighted mean: $ I_h = \frac{\sum_j w_j I^{\mathrm{corr}}{ij}}{\sum_j w_j},\qquad w_j = \frac{1}{(\sigma^{\mathrm{corr}}{ij})^2}. $

The weights use an expected variance: the Poisson signal part of each \sigma^{\mathrm{corr}}_{ij} is rebuilt at the reflection's merged \langle I\rangle rather than at that observation's own intensity. Weighting by an observation's own \sigma^2 biases the inverse-variance mean low below about one photon, because an up-fluctuated observation gets a larger sigma and is then down-weighted too hard. The rotation combine already does this; for stills it is on by default, and --no-expected-variance-merge restores the observed-sigma weighting.

An internal-consistency term can inflate uncertainties when multiple observations are present, in the spirit of XSCALE.

10.5 Merging statistics

The shells are nine bins of equal width in $1/d^2$, laid between the lowest- and the highest-resolution reflection the merge actually kept — XDS's rule and XDS's count, so at the same resolution limits the two programs' tables have the same shell boundaries and can be read row for row. --resolution-shells changes the count; the binning rule does not change with it.

Per-shell and overall merging statistics are computed on corrected intensities, including:

• number of observations and of unique reflections, and multiplicity,
• mean I/\sigma(I),
• 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,
• the anomalous signal-to-noise \mathrm{SigAno} (below).

The error model is refined as \sigma_\mathrm{corr}^2 = a\,\sigma^2 + (b\,\langle I\rangle)^2, with a set by the scatter of weak (counting-limited) reflections and b the intensity-proportional systematic scatter of the strong ones. On the rotation path, ISa is the asymptotic (I\to\infty) signal-to-noise — by definition the reproducibility limit of the strongest reflections (Diederichs, Acta Cryst. D66 (2010) 733) — and is read directly from the strong symmetry equivalents as the counting-subtracted fractional scatter of well-measured reflection groups (a robust median over strong groups; the I/\sigma threshold is relaxed on weak or radiation-damaged data that has few strong reflections), rather than as 1/b of the whole-range fit, whose b is raised slightly by an intermediate-intensity excess and so understates the limit. The asymptotic value is report-only — nothing downstream reads it, and the merged \sigma is not floored at b|I| (that floor was removed). The per-observation \sigma_\mathrm{corr} (the merge weights) uses the whole-range a,b. The stills path has no asymptotic estimate and reports \mathrm{ISa}=1/b directly.

a and b are reported in XDS's convention, which is \sigma^2 = a(\sigma_0^2 + b I^2) with \mathrm{ISa}=1/\sqrt{ab}, so the printed pair can be read straight against a CORRECT.LP. The internal fit keeps the form above; only the report converts, as b_\mathrm{XDS} = b^2/a. Note that a is the same in both conventions and that the two ISa expressions are the same number, 1/\sqrt{a\cdot b^2/a} = 1/b — so the rotation log prints two ISa, the whole-range 1/b (XDS's meaning) and the strong-reflection asymptote beside it, which can only ever be the more optimistic of the two. The mmCIF follows the same split: _reflns.jfjoch_diffrn_ISa is the whole-range value, directly comparable with a CORRECT.LP, and the asymptote is written separately as _reflns.jfjoch_diffrn_ISa_asymptotic, with _reflns.jfjoch_error_model_a and _b alongside so the number can be re-derived. Note that a file written before this change carries the asymptote under the plain ISa name. A third, unrelated b appears in the space-group search (§13.1); it is fitted with the \sigma^2 coefficient held at 1 and its gate constants are calibrated in that convention.

Anomalous signal-to-noise (SigAno). The strength of the anomalous signal is reported per shell and overall as \mathrm{SigAno}=\langle|\Delta I|\rangle / \langle\sigma(\Delta I)\rangle, where \Delta I = I(+)-I(-) over acentric reflections measured in both Bijvoet hands and \sigma(\Delta I)=\sqrt{\sigma_+^2+\sigma_-^2}. It is computed from the full-multiplicity inverse-variance I(+)/I(-) split (the same one written to the output), i.e. from all observations rather than a half-set. For pure noise \mathrm{SigAno} approaches the half-normal value \sqrt{2/\pi}\approx0.8, and it rises above 1 once a real anomalous difference is present. A half-set anomalous correlation ("$\mathrm{CC}_\mathrm{anom}$") is not reported. Its two half estimates \Delta I_0,\Delta I_1 are complementary partitions of one observation pool (\Delta I_0+\Delta I_1=2\,\Delta I_\mathrm{full}), and subtracting the two Bijvoet hands cancels the large common intensity that keeps \mathrm{CC}_{1/2} non-negative, leaving the small anomalous signal against the per-half split noise; below an anomalous signal-to-noise of 1 per half that correlation tends towards -1 rather than 0. \mathrm{SigAno} has no such floor. It is emitted only when an anomalous split was made, using the standard PDBx items _reflns.pdbx_absDiff_over_sigma_anomalous (overall) and _reflns_shell.pdbx_absDiff_over_sigma_anomalous (per shell), and appears as the SigAno column of the printed merge-statistics table.

10.6 Rotation datasets: combining partials into fulls (3D integration)

In a rotation scan a reflection is recorded as a series of partials spread across the frames its rocking curve crosses. Merging those partials directly would force the merge error model to absorb the rocking-curve slicing as if it were measurement noise, capping the achievable I/\sigma. For rotation data Jungfraujoch instead combines each reflection's partials into a single full intensity first, then scales and merges the fulls — a 3D integration over the rocking curve.

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. 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.)
• Per-image rejection (opt-in). A frame whose observations correlate poorly with the merged reference is not measuring the crystal being merged — it may be off-crystal, or on a different crystal where two lattices occupy separate regions of the sample. --min-image-cc drops such frames. It has no default: the per-frame correlation measures data quality as much as frame validity, and its typical level varies widely between datasets, so no single absolute bound is generally valid.
• 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.

After scale-fulls, four correction surfaces are fitted on the combined fulls (rotation path, on by default; disable all with --no-scaling-corrections), each an alternating multiplicative refinement of the per-full scale against the merged reference:

• Decay. Radiation damage weakens later frames more at higher resolution — a resolution×time (Debye–Waller) systematic the resolution-flat per-image scale cannot capture. A single global relative-B rate is fitted, \ln(I_\mathrm{ref}/I_\mathrm{obs}) = 2\,(\mathrm{d}B/\mathrm{d}n)\,(n-\bar n)\,s^2 (frame n, s^2 = 1/4d^2), and folded into the scale. It engages only when the total relative-B over the run exceeds a physical floor (2 Ų); below that the decay is negligible and "correcting" it only spreads symmetry equivalents (same s^2, different frames). An optional per-batch relative-$B$ (--relative-b[=deg], off unless requested; 10°-of-rotation batches by default) extends the single global rate to a smooth B(n) curve — the same $s^2$-weighted decay fit solved independently over short frame batches, curvature-penalized so it cannot over-fit and cross-validated like the surfaces below — for crystals whose decay is non-linear in dose. Its cross-validation splits on ASU-group parity, not the frame parity the surfaces below use: a per-batch parameter owns whole frames and so cannot be scored on a held-out frame, whereas splitting the symmetry equivalents tests whether a batch's B generalises to reflections it was not fitted on.
• Absorption. A smooth multiplicative factor over the diffracted-beam direction expressed in the goniometer (crystal) frame: each full's predicted detector position gives the lab diffracted direction, de-rotated by the spindle so a fixed crystal-frame direction is sampled at many rotation angles and its grid cell is well-determined. Negligible at hard X-rays / thin crystals; it matters at low photon energy.
• Modulation (detector-plane flat-field). A smooth multiplicative factor over where each reflection lands on the detector (predicted x,y): symmetry-equivalents land at different positions as the crystal rotates, over-determining the surface. It absorbs detector-response and geometric systematics that inflate R_\mathrm{meas}.
• Time-dependent absorption. The same surface as Absorption, but over (rotation angle × detector position) instead of the crystal-frame direction alone. The two agree while the illuminated volume stays put — the incident path then depends only on the spindle angle, which the per-image scale already takes, and the exit path is fixed in the crystal frame. Once the diffracting volume drifts through the beam the exit path becomes a function of the spindle angle as well, and nothing time-independent describes it. Fitted last, on 12 rotation bins × a 10×10 equal-occupancy detector grid, so the two time-independent surfaces get first claim on what they can explain.

Each surface is cross-validated: fitted on even-numbered frames and kept only if it improves the held-out odd-frame agreement by a clear margin (and vice versa), scored by a σ-independent, $R_\mathrm{meas}$-like fractional agreement \sum|I_s-I_\mathrm{ref}|/\sum|I_\mathrm{ref}| rather than a studentized \chi^2 — so a surface cannot "pass" by reshaping the sigmas instead of tightening the intensities. A surface fitted to noise where its systematic is absent does not generalize and is discarded — a correction never adds scatter.

Radiation-damage report (rotation, report-only). Independently of whether any decay correction is applied, rugnux measures and reports the relative Debye–Waller B across the sweep: the per-image scale's correlation to the merge and the per-image mosaicity versus frame (dose), together with a per-batch relative-B curve whose first→last change is a single headline number (measured before any decay correction, against the least-damaged early wedge). It is written to the log and to the merged mmCIF as a data-quality-vs-dose diagnostic and never alters the merged intensities — distinct from the decay correction above, which does fold into the scale.

Each batch's B is fitted on resolution-shell means, not on single observations: \ln(I_\mathrm{ref}/I_\mathrm{obs}) of one weak observation is unbounded and biased downwards — the observation appears in the response and in its own weight, and the logarithm needs I_\mathrm{obs} > 0, which keeps only the upward half of the noise — and on decayed data that bias grows with dose until it reverses the sign of the answer. The shells are laid inside the range the run actually diffracted to, and the fit carries an intercept as well as a slope, so a batch that is merely dimmer than the run (an attenuated beam, a mis-fitted frame scale) is not reported as damage. A batch whose shells are too weak to fit, or whose solved value reaches the bound the smoothing solve clamps to, is reported as absent rather than as a number. The first→last headline is reported only where a straight line describes the curve: radiation damage is progressive, so a curve that dips and recovers is a disturbance rather than dose, and is left to the sweep-quality report (docs/RUGNUX_REPORT.md) to name.

10.7 R-free test-set flags

A fraction of the unique reflections (rfree_fraction, default 0.05) is flagged as a free (test) set, written to the output (MTZ FreeR_flag, mmCIF _refln.status = f, a text-HKL column) for model validation (§14) and for downstream refinement. The flag is a pure function of the reflection's Friedel-merged (Laue) ASU index, which gives three properties:

• all symmetry- and Friedel-equivalent reflections share one flag — in particular a Bijvoet pair I(+)/I(-), kept as two separate merged rows in anomalous mode, is never split across the work and free sets (which would bias R-free);
• the free/work decision is a deterministic hash of that key, so the same reflection always lands in the same set — reproducible run-to-run and independent of the order in which observations were merged;
• the hash depends only on the reflection index, not on this dataset's resolution range or which reflections it happens to contain, so a uniform draw takes ~rfree_fraction of the distinct reflections free and — crucially — every dataset of one crystal form gets the same free set. That cross-dataset consistency is what a multi-dataset campaign (ensemble refinement, PanDDA) requires; a per-shell stratification tied to each dataset's own d_\mathrm{min} would break it.

On small data, where rfree_fraction (default 0.05) would give too few test reflections for a statistically stable R-free (Brünger's ~500–2000 rule), the fraction is floored so at least ~500 distinct reflections are free — capped at 10 % so a large test set never steals working data. For ordinary data this floor is inactive and the fraction stays flat at rfree_fraction, preserving the cross-dataset-identical property above; it only lifts the fraction on genuinely small datasets, where per-dataset R-free stability outweighs cross-dataset identity (and a shared reference set is the way to keep exact identity there).

When a reference MTZ (--reference-mtz) carries a FreeR_flag column, its test set is imported instead: every merged reflection whose Laue-ASU index matches the reference takes the reference's flag (reflections absent from the reference keep the hash flag). This lets a whole fragment-screening campaign inherit one shared free set from the apo/reference dataset. The CCP4/refmac convention (test set = flag 0, including the historical 0–19 form) is assumed, with the complement taken automatically if flag 0 would be the majority (a phenix-style file where 1 marks free).

10.8 French–Wilson amplitudes

The last step of the merge estimates a Bayesian structure-factor amplitude |F| for each unique reflection from its intensity I and error \sigma, so the output carries amplitudes alongside intensities (a naïve \sqrt{\max(I,0)} turns every weak or negative measurement into a biased — or zero — amplitude). With the Wilson prior for the true intensity J\ge 0 at that resolution,

$ P_\mathrm{acentric}(J) \propto e^{-J/\Sigma},\qquad P_\mathrm{centric}(J) \propto J^{-1/2},e^{-J/2\Sigma}, $

and a Gaussian likelihood \mathcal{N}(I;J,\sigma^2), the posterior mean amplitude and its uncertainty are

$ \langle |F|\rangle = \frac{\int_0^\infty \sqrt{J},\mathcal{N}(I;J,\sigma^2),P(J),\mathrm{d}J}{\int_0^\infty \mathcal{N}(I;J,\sigma^2),P(J),\mathrm{d}J},\qquad \sigma_F = \sqrt{\langle J\rangle - \langle|F|\rangle^2}. $

The prior mean is \Sigma = \varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell}, where \varepsilon is the reflection's epsilon (symmetry-enhancement) multiplicity and \langle I/\varepsilon\rangle is the Wilson mean in its resolution shell (so reflections on symmetry elements, and each shell, are treated correctly). Strong reflections (I>4\sigma) short-circuit to |F|=\sqrt{I}, where the French–Wilson bias is negligible; a reflection with an unusable I/\sigma falls back to \sqrt{\max(I,0)}. The integral is evaluated numerically with a log-shift for stability.

Amplitudes are written as MTZ F/SIGF, mmCIF _refln.F_meas_au/F_meas_sigma_au, and appended to the text HKL, alongside the intensity columns. The same |F| feed the model-validation step (§14), so the reflection file and the maps use one consistent set of amplitudes.

10.9 Reference data: fixing the space group and resolving the indexing ambiguity

A reference dataset (--reference-mtz) supplies known intensities for the same crystal form, and is used in two ways.

Fix the space group and cell. Unless overridden on the command line (-S for the space group, -C for the cell), the reference's space group is adopted and its cell is used as the soft reference cell — indexing may still drift the cell within tolerance, so a small mismatch between reference and data is absorbed rather than rejected. This applies to both stills and rotation data.

Resolve the indexing (merohedral) ambiguity. When the lattice symmetry is higher than the crystal's Laue symmetry (e.g. P3, P4, P6, C2), more than one indexing of the same lattice is geometrically valid, and the two solutions produce different merged intensities that a self-consistent scale cannot tell apart — only an external reference can. The candidate reindexings are the identity together with the twin-law cosets of the metric symmetry (from the unit-cell metric and the Laue group); each is scored by the intensity correlation \mathrm{CC}_\mathrm{ref} of the reindexed merge against the reference, and the data are re-merged in the best-correlating indexing. The reindex is metric-preserving — only the hkl labels change, the cell is unchanged — and it is a no-op for a holohedral crystal, which has no twin laws (the lattice and Laue symmetry coincide). For rotation data this is done once, after the space group is determined, and the whole merge is then repeated in the chosen indexing. For stills it has to be done per image, at integration time: each crystal is indexed independently, so a run resolves the ambiguity image by image (the image's partiality/Lorentz-corrected intensities are correlated with the reference under each candidate operator, which is scale-invariant, and the best-correlating one is adopted once and for good) — otherwise the merge would average reflections that are not symmetry mates. In neither workflow is the reference a scale target: both scale against their own data (§10.2), so \mathrm{ISa} and the merging statistics come from the data alone and no cross-dataset systematic is imported. Because the stills choice is made at integration time, a later re-merge of stored reflections cannot repair a dataset integrated without a reference. Where there is no reference dataset but there is a model (--model), the reference intensities are computed from it instead - |F_\mathrm{model}|^2 from the atomic structure factors with a flat bulk-solvent contribution at the standard constants (k_\mathrm{sol}=0.35, B_\mathrm{sol}=46 Å^2), which are not fitted because there are no observations yet. Nothing is scaled against them; they serve only to rank the candidate indexings, and the correlation that does the ranking is scale-invariant. This needs the cell and space group up front (-C / -S, as serial indexing wants anyway); on rotation data the same job is done after the merge, in §14.5, where a merge exists to fit the model to properly.

10.10 Ice rings at the scale and merge stages

Where the gate of §3.3 has found ice, reflections falling within \pm w in q = 2\pi/d (§1.2) of a hexagonal-ice band (w=0.03 Å^{-1} offline, about the measured ring half-width) are marked. Marked reflections are excluded where a model is fitted — the per-frame scale G, the per-image correlation, and the P1 merge the space-group search runs on — because ice contamination is a positive bias, not extra scatter, and a least-squares scale absorbs it into G and into the error-model b, where it damages every other reflection on the same frame. They are kept in the final merge, which is also what the established scaling programs do by default, so the affected shells keep their completeness.

Nothing on an ice band is deleted from the merged output. Deleting the bands was implemented, measured against an external arbiter rather than against the merge's own statistics, and removed: on the one rotation-battery crystal where a band was both dead by its own merged \mathrm{CC}_{1/2} and scorable by anomalous peak height, dropping it changed the mean anomalous density at the known sites by -0.001\pm0.018\,\sigma — about 2 % of the site height — while removing 1149 unique reflections whose mean I/\sigma was 3.62 against the dataset's own 3.05, i.e. better-than-average data, and costing 6 to 8 points of completeness in the affected shell.

---

11. Mosaicity and “profile radius” monitoring

11.1 Profile radius (intrinsic excitation-error width)

The “profile radius” is the intrinsic angular width of a reflection — crystal mosaicity plus beam divergence — estimated from the spread of \Delta_\mathrm{Ewald} over indexed spots, $ R \approx \sqrt{\tfrac{1}{N}\sum_i \Delta_{\mathrm{Ewald},i}^2}. $ When the beam has a finite energy bandwidth, that bandwidth smears each reflection radially by \sigma_\mathrm{bw}\approx \mathrm{bandwidth}\cdot\lambda/2d^2 (largest at high resolution), which also broadens the measured \Delta_\mathrm{Ewald} spread. Since prediction re-applies the bandwidth term per reflection (§8.2), this contribution is deconvolved from the estimate — R^2 = \langle\Delta_\mathrm{Ewald}^2\rangle - \langle\sigma_\mathrm{bw}^2\rangle — so that R is the intrinsic width and bandwidth is not double-counted. Still predictions use an excitation-error cutoff proportional to R.

11.2 Mosaicity from rotation data

For rotation data the mosaicity \sigma_M is estimated by maximum likelihood from the rocking offsets \tau of indexed spots, using the XDS reflection-fraction model R(\tau;\sigma_M/\zeta) (Kabsch 2010): each spot's exact Bragg angle is located near its frame, \zeta (the rotation-axis Lorentz component) is computed, and \sigma_M is chosen to maximize \sum_i \log R(\tau_i;\sigma_M/\zeta_i).

The \phi search window for the Bragg angle is set wider than the oscillation, so that reflections recorded at large rocking offset are included. These tail reflections carry most of the information about the mosaic width; a window limited to the oscillation range would truncate the \tau distribution and bias \sigma_M low.

The fit uses only the strongest 250 spots of an image, whatever the indexing spot budget (--max-spots) is. A spot is detected when I_\mathrm{full}R(\tau) clears the finder threshold, so selecting spots by intensity censors on R(\tau): a deeper list holds proportionally more large-\tau partially recorded spots and the fit widens with it. Left uncapped, \sigma_M therefore tracks the spot budget rather than the crystal — and since an over-wide mosaicity mis-states every partiality, the merge degrades sharply with it.

The estimated mosaicity feeds the rotation prediction (how many frames each reflection spans, §8.3) and the rotation partiality (§10.2). It is held fixed during scaling: in the per-image scale fit the mosaicity is degenerate with the scale G (both rescale the predicted intensity), so refining it there is unstable. A correct mosaicity matters because it controls both how much of each rocking curve is captured and the partiality used to form fulls (§10.6); too small a value truncates the captured curve and over-peaks the partiality, degrading the combined fulls.

---

12. Auxiliary statistics: ⟨I/σ(I)⟩ and Wilson plot

12.1 Per-shell ⟨I/σ(I)⟩

For monitoring integration quality, Jungfraujoch reports mean \langle I/\sigma(I)\rangle in a fixed number of resolution shells. Shelling is performed in 1/d^2 space (typical of crystallographic practice).

12.2 Wilson plot (B-factor proxy)

A Wilson-type analysis is computed by binning intensities by resolution and fitting: $ \langle I\rangle \propto \exp!\left(-\frac{B}{2}\frac{1}{d^2}\right), $ i.e. $ \log \langle I\rangle = \mathrm{const} - \frac{B}{2}\left(\frac{1}{d^2}\right). $ 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).

A dataset-wide Wilson B is also estimated over the merged reflections — restricted to the meaningful resolution range (skipping the low-resolution non-linear region below ~4 Å and shells past the signal limit \langle I/\sigma\rangle < 1, so it is insensitive to how far the merged data extend) — and written to the merged mmCIF as _reflns.B_iso_Wilson_estimate (and reported as WILSON_B=), the analogue of XDS's Wilson-line B. It is diagnostic only and is not fed back into scaling. It is the XDS convention, not the CCP4/Phenix one, and the two are not comparable. The slope here is fitted to \log\langle I\rangle directly; TRUNCATE, ctruncate and phenix.xtriage fit \log(\langle I\rangle/\Sigma), dividing out \Sigma=\sum_j f_j^2(s), the falloff of the atomic form factors for an assumed composition. Leaving \Sigma in the slope inflates B by roughly 2 to 8 Ų (measured across in-house merges; the arithmetic gives +6.9 Ų over 4.0-1.5 Å for a generic protein), and the fitted range accounts for more still: against ctruncate and phenix.xtriage on the same merged files this number runs 10 to 36 Ų high, in the same direction every time - though those two disagree with each other by 7 to 16 Ų, so there is a band rather than a right answer. Read it as a relative quantity, comparable between rugnux runs and against XDS, and do not compare it with a value quoted from a CCP4 or Phenix log. The same caveat applies to _reflns.B_iso_Wilson_estimate in the merged mmCIF, whose deposited values are conventionally the $\Sigma$-normalised kind. It is also where the flux the fixed integration disk clips lands: that loss is degenerate with an overall B, so the reported number carries an $r_1$-dependent contribution and is not a property of the crystal alone (§9.1). The per-image estimate (used for the live radiation-damage plot) is accepted only when the fit is well-correlated and physically plausible (0 < B < 200 Ų); on a bad frame (an indexing glitch, too few reflections) the Wilson line runs wildly steep, so an implausible B is reported as NaN rather than a spurious hundreds-of-Ų value.