Jungfraujoch/docs/CPU_DATA_ANALYSIS_INDEXING.md

Data analysis: indexing and geometry (§4–§7)

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

:local:
:depth: 2

4. Indexing overview

Indexing maps observed reciprocal-space vectors \mathbf{s}_i to a lattice such that: $ \mathbf{s}_i \approx h_i\mathbf{a}^* + k_i\mathbf{b}^* + l_i\mathbf{c}^*, $ with integer (h_i,k_i,l_i).

Jungfraujoch supports two complementary indexing strategies:

1. FFT-based indexing (Rossmann-type): does not require an a priori unit cell; suitable for unknown samples.
2. Fast-feedback indexing (TORO-like): requires an approximate unit cell; optimized for speed and feedback.

Both feed into a common robust refinement/selection stage which maximizes the number of inliers under an indexing tolerance, and which can return more than one lattice per image (multi-lattice indexing; see §5.4).

4.1 Indexed-spot decision (inlier test)

Given a trial lattice with direct basis vectors \mathbf{a},\mathbf{b},\mathbf{c} (used here as reciprocal-space dot-test vectors), fractional indices are estimated by: $ h_f = \mathbf{s}\cdot\mathbf{a},\quad k_f = \mathbf{s}\cdot\mathbf{b},\quad l_f = \mathbf{s}\cdot\mathbf{c}. $ Let (h,k,l)=(\mathrm{round}(h_f),\mathrm{round}(k_f),\mathrm{round}(l_f)) and define the fractional residual: $ \delta^2 = (h_f-h)^2 + (k_f-k)^2 + (l_f-l)^2. $ A spot is indexed if \delta^2 < \tau^2, where \tau is the configured tolerance.

For indexed spots, the reciprocal lattice point \mathbf{p} = h\mathbf{a}^*+k\mathbf{b}^*+l\mathbf{c}^* is used to compute \Delta_\mathrm{Ewald}(\mathbf{p}) (stored as a diagnostic and later used in profile-radius estimation).

A frame is taken to be this crystal's when at least a fraction g = 0.20 of its in-resolution, non-ice spots index. On rotation data that decision is what admits the frame to integration, so its denominator matters: every spot handed to it that is not a reflection of this crystal argues against the frame.

4.2 The spot budget

Only the strongest --max-spots spots of an image are kept (FilterSpotsByCount), and that budget therefore sets the denominator above. Detections are not all reflections — background structure, unlisted ice and detector artefacts are found too — so a budget deeper than an image's reflections makes the test above a measurement of the background rather than of the crystal, and a larger budget can integrate fewer images.

rugnux measures the budget instead of fixing it. With the sweep's lattice in hand, the first pass tallies the spots of a sample of frames by their rank in the intensity-ordered list: how many images carried a spot at that rank, and on how many of them it indexed. Weighting each indexed spot by 1-g and each unindexed one by -g — the same weighing the frame test applies to the list as a whole — the running sum over ranks

E(N) = n_\mathrm{indexed}(N) - g\, n_\mathrm{counted}(N)

rises exactly while the spots at that depth lie on the lattice more often than g, and falls after. The budget is \arg\max_N E(N). Its meaning is "as deep into the list as the image is still showing reflections of this crystal": deeper spots cannot help the frame test and can only push a frame towards rejection. On crystals whose spot lists are reflections all the way down the maximum is at the end of the list and the budget is unchanged.

The peak has to be one. Under the null — the spots lie on the lattice at the same rate at every depth — E is a driftless random walk in the counted spots, with per-spot variance g(1-g), and the maximum of such a walk is positive whatever the data; an \arg\max taken on its own would shorten every dataset, including one with nothing to shorten. What the budget acts on is the fall from the peak to the end of the list, E(N^*) - E(L), which is the maximum of the same walk read backwards; by the reflection principle its null law is P(\mathrm{fall} > z\sqrt{g(1-g)T}) = 2(1-\Phi(z)) over T counted spots in all, so the search over ranks is already accounted for and no further multiple-comparison correction applies. The budget is taken only where the fall clears that bar at z = 3.29, one false shortening in a thousand measurements; otherwise the whole list is kept.

---

5. FFT indexing (unknown unit cell)

FFT indexing follows a classical approach: detect dominant periodicities by projecting reciprocal-space points onto many directions and Fourier transforming the resulting 1D histograms.

5.1 Directional projections and histograms

Choose a set of unit vectors \{\mathbf{u}_d\} on a half-sphere (a near-uniform distribution generated via a golden-angle construction). For each direction d, form a histogram in the scalar projection: $ t_{id} = \left|\mathbf{u}_d\cdot \mathbf{s}_i\right|. $

Bin width is chosen approximately as: $ \Delta t \approx \frac{1}{2 L_\mathrm{max}}, $ where L_\mathrm{max} is the maximum expected real-space unit-cell edge (Å). The histogram extent is tied to the maximum q used (set by a high-resolution cutoff for indexing).

5.2 FFT peak picking and candidate vectors

For each direction, the FFT magnitude spectrum is computed; peaks correspond to periodicities along \mathbf{u}_d. Each direction yields a candidate real-space length L chosen not by raw magnitude but by maximum prominence above a running-mean local background (subtracting the broad low-frequency envelope that otherwise dominates on weak or pink-beam frames), subject to L\ge L_\mathrm{min}.

The running-mean background window keeps a constant width and is slid inward at the ends of the spectrum rather than truncated there, so a peak within half a window of either end — which is where the longest cells sit — is judged against as much background as any other. Both window bounds stay monotonically non-decreasing in the bin index, so the GPU kernel's running sum is still valid.

The longest basis vector the transform can return is fft_max_unit_cell_A, since the histogram is sized from it and its last usable bin is that length; the shortest is fft_min_unit_cell_A (rugnux --fft-min-unit-cell, default 10 Å), below which a candidate is discarded. The defaults are unchanged (500 Å and 10 Å), but the accepted range for the maximum now reaches 1200 Å, and a reference cell given with -C moves both bounds on its own — up to reach a long axis, down to admit a small-molecule cell — since a cell the search cannot represent cannot be found by it.

Candidate vectors are \mathbf{v}_d = L_d\,\mathbf{u}_d.

A collinearity filter removes nearly parallel vectors (e.g. within 5°) and attempts to resolve harmonic ambiguity: shorter “fundamental” vectors may be preferred over longer harmonics if their peak magnitude is sufficiently strong relative to the dominant peak.

5.3 Lattice reduction and cell candidates

Triples of candidate vectors are combined to form candidate bases (\mathbf{A},\mathbf{B},\mathbf{C}), each reduced to its Niggli-reduced cell (Gruber-vector reduction) before comparison, and filtered by allowed length and angle ranges. Two passes are run: a standard pass forms shortest-vector triples from the ~30 strongest filtered directions; if the best cell then indexes fewer than half the spots, a widened fallback anchors the two shortest axes and lets the third range over up to ~60 candidate vectors (deduplicated by Niggli cell), catching large, elongated or superstructure cells the first pass misses.

A triple whose three vectors are coplanar is rejected before refinement. The length and angle filters cannot see it — any flat combination satisfies them — and a cell that flat has a metric determinant small enough for float to get its sign wrong, after which the guard against a negative argument to the square root places \mathbf{c} in the $\mathbf{a}$-\mathbf{b} plane, the reciprocal volume diverges and the solver reports a not-a-number Jacobian. The test is the volume fraction |V|/(|\mathbf{a}||\mathbf{b}||\mathbf{c}|), which must reach 0.02 — about 1.1° off flat, well below the flattest genuine candidate observed and far above where float loses the sign — and it is applied both where triples are produced and at the optimizer's entry points.

A shortlist confined to one plane cannot close a cell at all, and the row it is missing is the plane normal. That is detected from the eigenvalue ratio of the shortlisted directions' scatter matrix, and one further transform is then spent with the same direction count inside a narrow cap about the normal. More directions do not substitute for it: at the exact true direction a very long axis can still rank far below the shortlist cut, so for this rescue the obstacle is the ranking rather than the sampling, and a denser grid costs several times the device memory for the same answer.

Sampling has a limit of its own, and it binds well before the 1200 Å the accepted range for fft_max_unit_cell_A admits (§5.2). A direction off a real-space axis of length a by an angle \theta smears each projected lattice plane by about \theta/d_\mathrm{min} in the projection, so the planes (spacing 1/a) stay resolved only while \theta \lesssim d_\mathrm{min}/(2a). The shipped grid of 16384 directions puts the nearest one within about 0.6° of any axis, which satisfies that bound only up to roughly 120–150 Å at typical indexing resolutions; a longer axis is not refused but returned as a plausible sub-cell or harmonic. Raising the maximum cell alone therefore does not extend the reach — the direction grid has to resolve the axis before the histogram can represent it.

5.4 Robust refinement and best-cell selection

Candidate bases are refined against observed spots using an iterative inlier‑focused least‑squares procedure (trimmed/contracting threshold). Candidates are then ranked:

1. more indexed spots wins — unless two candidates index within ~10 % of each other, in which case
2. the smaller-volume cell is preferred (when the volumes differ by more than ~5 %), avoiding a doubled supercell, then
3. the smaller refinement score, then the spot count again.

Selection is not limited to a single lattice: after the best cell is accepted, further lattices are added as separate crystals provided fewer than ~40 % of their indexed spots overlap an already-accepted lattice (up to two extra by default), so split or multi-lattice crystals are indexed rather than discarded.

An optional reference unit cell (if supplied) restricts acceptance to cells within a relative distance tolerance in edge lengths (permutation-invariant).

---

6. Bravais lattice / centering inference (“lattice search”)

If the space group is supplied by the user, its lattice constraints are assumed for refinement and subsequent processing.

If not, Jungfraujoch attempts to infer the most plausible Bravais lattice type from the metric tensor after Niggli reduction:

1. Niggli reduction is performed to obtain a reduced cell in G^6 representation (Gruber vector).
2. The reduced cell is compared against a list of Niggli classes corresponding to Bravais lattices and centerings.
3. The highest-symmetry class that matches within tolerances is selected (relative metric tolerance and angular tolerance). The list is walked in order of decreasing symmetry and the first class that fits is taken, so a class that only just fits can pre-empt a lower-symmetry one that fits exactly.

The output includes:

• a conventional cell,
• crystal system (triclinic, monoclinic, …),
• centering symbol (one of P, C, I, F, R; the A/B variants are not emitted here — they are handled only later as prediction absences, §8.4).

This stage provides centering information used for systematic absences in prediction (§8.4) and for reporting.

A metric symmetry has to earn itself. The class is chosen from the unrefined candidate against a fixed angular tolerance (3°), so a lattice that is pseudo-symmetric to a few tenths of a degree is promoted a class too far — and the constraint then snaps a real angle to the ideal one, which throws nearly every reflection of every frame outside tolerance. Measured on a monoclinic crystal pseudo-C-orthorhombic to 0.42°, the promoted cell indexed 2 of 60 frames where its own primitive cell indexed 39: the same lattice, \mathbf{b}_{oC}=-(\mathbf{a}+2\mathbf{c}), at exactly twice the volume. A more accurate candidate is more likely to be promoted, not less: a run escapes the promotion only when the raw candidate misses the tolerance window.

The rotation first pass refines the twelve best candidate lattices rather than four. The pre-refinement indexed fraction is an unreliable ranking, so a correct cell can sit below several degenerate ones and never be refined at all.

For rotation data the first pass therefore refines the constrained cell and an unconstrained (triclinic) one from the same spots — which it finds itself, over a sample spread across the sweep, rather than reading what the acquisition wrote — and settles the two on how many of a fixed set of validation frames each actually indexes. The bar is a clear majority, not a margin. An unconstrained refinement holds no cell parameter fixed, so it can only index at least as many frames as the constrained one, and on genuine symmetry it indexes a few more — a percentage margin therefore demotes real lattices (measured: a genuine $I$-centred orthorhombic to P1). Only a constrained cell that fails outright while its unconstrained cell works distinguishes a false promotion. The intensities settle the space group later regardless (§13).

Two further hypotheses are weighed at the same point, both by default. Where two first-pass answers have primitive volumes in a small integer ratio (2–4×) and tie on the validation frames, the frame count has saturated — a spurious axis multiple indexes every frame its true sub-cell does — so the tie is settled at the granularity that does not saturate: which cell accounts for more of the validation frames' spots. That comparison leans toward the smaller cell by construction (indexing is a fractional-Miller test, so multiplying an axis by n multiplies that axis's residual by n), and only a reflection class the larger cell genuinely adds — a real superstructure's satellite rows — can pay for the loss; the occupancy of that added class is computed and logged beside the decision, deliberately without a threshold, but the spot count is what decides. Separately, the 10 Å floor of §5.2 is applied to the FFT's per-direction peak search, so a lattice row whose true repeat is below it is reported at its first harmonic and the cell assembled from those harmonics is an exact integer supercell of the true one; a second first-pass hypothesis with the floor lowered to 5 Å therefore also runs (rotation only, and not when a cell was given — -C already lowers the floor to cover it), and its answer is adopted only when it is an integer sub-cell of the standing one and ties or beats it on the validation frames. Everywhere else the standard pass's answer stands.

Note. In ambiguous or special cases, forcing space group to P1 (no symmetry assumptions) is recommended.

---

7. Geometry and lattice refinement

Refinement adjusts experimental geometry and crystal parameters to minimize discrepancies between observed spot reciprocal vectors and those predicted by a lattice model with integer indices.

7.1 Parameterization

The refinement jointly optimizes, depending on mode and constraints:

• beam center (x_\mathrm{beam}, y_\mathrm{beam}),
• detector distance D,
• detector tilt angles (two-angle model; third rotation often held at 0),
• rotation axis direction (for rotation datasets),
• crystal orientation (a global rotation),
• unit-cell parameters, with constraints determined by inferred crystal system.

The detector distance is not refined against one crystal's spots at all: the positional residual leaves it degenerate with the cell scale, so it is fitted separately - by the rotation post-refinement, which holds the cell at the value its own cell/axis step settled on, and by the stills --refine-geometry bundle. Per image, the beam centre and the crystal orientation are refined, and the unit cell as well for stills. The first-pass rotation indexing refines the detector tilt and the rotation-axis direction too, against the spots accumulated across the sweep; everywhere else both are held fixed, because on a single crystal a tilt is absorbed almost exactly by the crystal orientation. A lighter orientation-only mode refines just the crystal orientation, for stills whose geometry is already trusted. It carries a weak small-rotation prior penalising the whole angle-axis vector (all three components, at a low weight); what it is there for is the poorly-determined out-of-plane component, which is the one the data barely constrain.

For higher symmetries, constraints are enforced, e.g.

• cubic: a=b=c,\ \alpha=\beta=\gamma=90^\circ,
• tetragonal: a=b,
• hexagonal: a=b,\ \gamma=120^\circ,
• monoclinic (unique axis b): \alpha=\gamma=90^\circ, \beta refined.

7.2 Residuals and objective

For each indexed spot assigned integer (h,k,l), compute:

• observed reciprocal vector \mathbf{s}_\mathrm{obs} from its detector position and current geometry,
• predicted reciprocal vector \mathbf{s}_\mathrm{pred}(h,k,l;\ \text{lattice params}).

Residual is: $ \mathbf{r} = \mathbf{s}\mathrm{obs} - \mathbf{s}\mathrm{pred}. $

A non-linear least squares solver minimizes \sum \|\mathbf{r}\|^2 over all selected inlier spots.

7.3 Rotation datasets: bringing observations to a common reference frame

For oscillation/rotation data, each image corresponds to a rotation angle \phi about an axis \mathbf{m}_2. Observed reciprocal vectors are rotated “back to start” so that all images are refined in a single reference crystal frame: $ \mathbf{s}\mathrm{obs,ref} = R(\phi),\mathbf{s}\mathrm{obs}, $ where R(\phi) is the rotation by +\phi about the goniometer axis as stored in the file. The sign is a convention and it is load-bearing: rotating the observations forward by +\phi means the crystal itself turns by -\phi about that stored axis, i.e. R(\phi) is the inverse of the crystal's own rotation from the reference orientation to frame \phi. The same convention is why the unmerged-MTZ batch headers and the XDS geometry echo carry the axis negated relative to the input file (rugnux ▸ the unmerged export) — a reimplementation that takes R(\phi) as the crystal rotation must use R(-\phi) here instead. The angle \phi is taken at the centre of each frame's oscillation (the frame angle plus half the oscillation width).

7.4 Multi-stage tightening of inlier tolerance

Refinement is performed in stages with decreasing acceptance tolerance for including reflections (three stages, indexing tolerance 0.3\to0.2\to0.1), which stabilizes convergence when starting from imperfect indexing and approximate geometry.

The loose first stage necessarily admits some spots that are not reflections of this lattice — the fraction of randomly placed spots inside a fractional-Miller tolerance t is \tfrac{4}{3}\pi t^3, i.e. 11 % at t=0.3 — and an unweighted fit lets them pull the orientation. Each residual is therefore weighted by how strong its spot is for its resolution: the frame's spots are cut into equal-count resolution shells and each intensity is divided by its shell median, mapped to w^2=r/(1+r). The shell normalisation is what makes this safe — genuine high-resolution spots are legitimately weaker and carry the cell and distance information, so an un-normalised intensity weight would suppress exactly the spots the fit needs. The weight is a property of the spot and never of the current residual, so it does not depend on how far the geometry is from convergence.

7.5 Rotation geometry post-refinement (two-pass)

The refinement above (§7.2) runs per image against that image's spots. For rotation data an additional post-refinement (on by default; --rotation-no-postrefine disables it) improves the detector distance, beam centre and crystal cell/axis using all frames at once, then re-integrates:

1. Pass 1 integrates, scales and merges at the header geometry.

2. From pass-1's integrated reflections, the geometry is refined over all frames (Ceres, robust loss) in two separate steps rather than one joint fit:

   • Step A: crystal cell scale + goniometer-axis direction, from the observed rotation angles (a distance-independent excitation residual).
   • Step B: shared detector distance + beam centre, from the observed spot positions, with the cell held at step A — so the positional residual is no longer degenerate with the cell scale.

   Each step is cross-validated on a deterministic split of the reflections (an avalanche-mixed hkl hash, not a frame split and not an h+k+l parity, which would collide with a centering condition and leave the held-out half empty): fitted on one half, committed only if it lowers the held-out residual, otherwise left at nominal. The solver bounds the move — distance within ±5 %, beam centre within ±15 px — and detector tilt is held fixed, being gauge-coupled to the crystal orientation on a single crystal.

3. Pass 2 re-indexes de novo and re-integrates at the committed geometry. Only the detector distance and beam centre carry over: the refined cell and axis are used to make step B well-posed, but pass 2 re-indexes from scratch, so they are not propagated.

The space group is determined after pass 2, on the geometry the run refined, and pass 1 does not search at all: a decision taken on the worse of the two passes and then carried forward is a constraint on the better one, and had to be reconciled with what pass 2 later found. The guard that chooses which pass is written compares each pass's first merge — P1 on both sides, full resolution range, before the correction surfaces — which both passes produce anyway, so it never compares statistics computed in two different space groups. One index-time veto remains and is keyed to pass 1's lattice rather than its group: a centred pass-1 lattice against a primitive pass-2 one.

Only pass 2 is written, as the canonical <prefix>_* output. Pass 1's merge exists to give the guard something to judge pass 2 against, so it stops short of the parts of the merge that only fill in a file — the correction surfaces, the twinning and radiation-damage analyses, the R-free flags and the amplitudes — and writes no merged files of its own.

Goniometer rotation scale (report only). A stage that turns further than it was commanded to leaves no trace in the file, because the stored \omega values are the commanded ones; the excess then presents as the crystal drifting, in this program and in others. Step A already measures it without a new degree of freedom: its residual rotates by -\phi\,\mathbf{u} with \mathbf{u} an unnormalised 3-vector, so |\mathbf{u}| is the factor by which the stage actually turned, and normalising the axis throws it away. It is reported, and warned about beyond 0.5 %, under the same cross-validation that gates the cell move — a fold that merely soaked up noise cannot raise the flag. It is a detector, not a calibration: nothing corrects the data, and it under-reads the true magnitude, because the fit only sees reflections that indexed at the nominal angle and per-frame orientation refinement has already absorbed part of the error.

7.6 Detector geometry from powder rings

Everything above fits the geometry to Bragg data, where the beam centre is the weakest parameter: it is gauge-coupled to the crystal orientation, which is why §7.5 restrains it toward the header value and commits only a sub-1 % move. A powder ring has no orientation to be coupled to. Where it falls on the detector depends on the geometry and on nothing else, which makes a calibrant — LaB₆, silver behenate, CeO₂, silicon — or even ice an independent constraint on exactly the quantity Bragg data cannot pin.

The ring positions are matched to the observed rings and the geometry is refined (Ceres, five parameters: beam centre, distance, and the two detector tilts) so that the |s| predicted at each observed ring point matches the ring it belongs to. The two tilts can be held fixed (rugnux --no-refine-tilt, the viewer's Refine detector tilt tick box), leaving a three-parameter fit: a tilt a downstream program cannot express is better left out of the fit than refined and then dropped, since the centre and the distance of a tilted fit have already absorbed it.

Calibrants. LaB₆, silver behenate, CeO₂ and silicon are held as unit cells and their rings enumerated from them. Ice is held as the hexagonal-ice ring positions of §3.3 instead — measured to 1.522 Å, calculated below it — because hexagonal ice is P6_3/mmc with oxygen on 4f and enumerating hkl from its cell would emit rings the oxygen sublattice extinguishes. A calibrant is therefore a list of ring q values throughout, not a cell.

What a ring can and cannot determine. A ring is a conic centred on the beam, so a wrong centre makes its apparent radius oscillate once per turn, r(\phi)=R+\delta_x\cos\phi+\delta_y\sin\phi, with the same amplitude on every ring. A detector tilt \beta produces a \cos\phi term too — not the \cos2\phi one might expect — but one that grows as the ring's radius squared, r(\phi)=R+(R^2/F)(\beta_x\cos\phi+\beta_y\sin\phi); the true \cos2\phi term is O(R^3\beta^2/F^2), hundredths of a pixel. The two are therefore separated by how the amplitude scales with radius, which needs at least two rings — on a single ring they are exactly degenerate. None of this uses the calibrant's $d$-spacings, so the centre is determined without assuming anything about the standard.

The distance is different: it follows from r=F\tan2\theta with \sin\theta=\lambda/2d, so a fractional error in the lattice constant passes straight into it, and the $\lambda$–F pair is separated only by the curvature of \tan(2\arcsin(\lambda/2d)) across the rings — \partial\ln r/\partial\ln F=1 at every ring against \partial\ln r/\partial\ln\lambda=4\tan\theta/\sin4\theta, which runs from about 1.05 at low angle to 1.43 at high. That lever collapses as the detector moves back and the rings crowd into small 2\theta, so distance is a short-distance measurement and the wavelength is better calibrated by other means.

Reading the rings. The ring points come from one of two measurements, both accumulated over every processed image rather than one. The default reads the azimuthally-binned profile (§2) summed over the run: for each ring and each azimuthal sector, the radial peak is fitted against a locally interpolated background and the measured (q,\phi) mapped back through the current geometry to the pixel it came from. The alternative pools the spot lists, which samples each arc wherever the spot finder's threshold happens to bite. The accumulated profile is the same size however many images went into it; the pooled spot list is capped, each image contributing an equal share.

A plain radial profile — one azimuthal sector — has averaged the ring over every direction and carries no centre at all, so the profile route requires at least four sectors and uses 32 by default. Sixteen to thirty-two are enough; beyond that the limit is the ring's own texture, not counting statistics.

The extraction window around a ring is capped at half the gap to its neighbour, because the background under a peak is taken from the ends of that window: hexagonal ice has a triplet of rings (1.947, 1.916 and 1.882 Å) whose neighbours sit only 0.05–0.06 Å⁻¹ apart in q = 2\pi/d, which a fixed window merges into a single peak. Where only one ring is in reach the two tilts are held at their input values rather than fitted, since on a single ring they are degenerate with the centre (above) and the fit would otherwise trade the centre away for them.