RESULT-CHANGING: FFT indexing finds rows shorter than its floor; the short-axis pass is retired

The per-direction peak search starts at the 10 A floor, so a lattice row shorter than that was
reported at a multiple of itself past it and the cell built from such rows was an exact integer
supercell. A second first-pass hypothesis at a 5 A floor (the short-axis pass) won the true cell
back on frame counts, and could not reach an axis below 5 A at all.

On rotation data (a scanning axis) every shortlisted row reported within three times the floor is
now asked, by a direct sum over the accumulated spots, at each whole fraction of its length down
to 3 A; the shortest fraction that keeps half the row's amplitude replaces it. The cells of those
rows are reduced (widened: a small cell's in-plane combinations fill the standard window) and
refined beside the standard cells, and ranked with them, so where a cell and its supercell index
alike the smaller one is taken. Protein shortlists have no rows that short, so nothing changes for
them; stills and grid scans keep the search they had. The short-axis pass, its speculative pool and
its live warm-up are removed; the refined-length bound of the cell fit goes from 5 to 2.5 A.

The 40 % overlap test of the candidate refinement no longer deletes a cell for its overlap with an
integer multiple or sub-cell of it - one lattice contains the other, so that overlap is total by
construction. Such candidates are kept after every cell the test accepts, which keep their order,
so the first lattice and the extra crystals are what they were. The rotation indexer carries the
refined integer supercells of the cell it chooses, and the twin-composite test asks them as well as
the standing cell: once the search offers a twin domain's own cell, the coincidence lattice it used
to stand on is one of those.

Measured (battery subsets vs rc181-b2): a small-molecule crystal with an axis under 5 A, indexed as an
index-5 supercell now gives its true cell and group; the 13 sets the short-axis pass rescued keep
their verdicts; the sparse reticular twin keeps its domain and overlap.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UqXYdL4qBGYRiSXSGXfCLu

docs/CPU_DATA_ANALYSIS_IMAGE.md

@@ -250,7 +250,7 @@ within 2 %. The outcomes:
      given.
 
 **Both centres judged on the merge (rotation, two-pass).** Where case 3 or 5 carried both centres, the
-whole first pass is run again at the measured centre — with its own short-axis pass and
+whole first pass is run again at the measured centre — with its own
 post-refinement — and the two arms are compared on their search merges ($P1$, whole range, before the
 correction surfaces). Where case 2 held a hypothesis, this is done only if the file's arm merges
 inconsistently: $\mathrm{CC}_{1/2}$ below 0.9 in the lowest-resolution shell of its search merge, where

docs/CPU_DATA_ANALYSIS_INDEXING.md

@@ -115,6 +115,8 @@ 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.
 
+**Rows shorter than the floor (rotation).** The per-direction peak search starts at $L_\mathrm{min}$, so a lattice row whose repeat is shorter — common in small molecules, absent in proteins — is reported at a multiple of itself past the floor, and a cell built from such rows is an exact integer supercell of the true one. The spots can tell a row from its multiple where the histogram cannot: summed at the row's own length every spot is in phase whichever of the two it is, but at a whole fraction $L/k$ of it only the spots of a lattice whose row that fraction is stay in phase. So on a rotation run every shortlisted row reported within $3L_\mathrm{min}$ is asked, by a direct sum over the accumulated spots, at each $L/k$ between 3 Å (the minimum cell of DIALS' fft3d) and the floor, and the shortest fraction that keeps at least half the row's amplitude replaces it. The cells of these rows are reduced and refined **beside** the standard ones, not instead of them, and ranked with them (§5.4): where a cell and its supercell index alike, the smaller is kept. Stills keep the standard search.
+
 ### 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.
@@ -134,6 +136,8 @@ Candidate bases are refined against observed spots using an iterative inlier‑f
 
 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.
 
+A candidate dropped only for its overlap with a lattice that **contains** it or is contained by it — an integer supercell or sub-cell, the rows of the larger cell integer combinations of the smaller's — is kept nonetheless, after every cell the test accepts: that overlap is total by construction and says nothing about which of the two is the crystal's. The rotation indexer refines these too and may choose one (a sub-cell that indexes as well, §5.2); the integer supercells of the cell it chooses are carried with it, so that the twin-composite test of the first pass, which asks whether a cell is the coincidence lattice of two twin domains, can also be asked of them when the search has offered a domain's own cell. Kept cells are never taken as extra crystals: those must share the first lattice's metric.
+
 An optional reference unit cell (if supplied) restricts acceptance to cells within a relative distance tolerance in edge lengths (permutation-invariant).
 
 ### 5.5 Spindle alignment: the part of the blind cone symmetry cannot repair
@@ -243,7 +247,7 @@ The rotation first pass refines the **twelve** best candidate lattices rather th
 
 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 first part of the sweep (the part the geometry pre-pass integrates), 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
+A further hypothesis is weighed at the same point, 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 — and no
 count of spots can settle the pair either: the larger lattice contains the smaller one, so it fits every
@@ -253,13 +257,8 @@ the smaller cell, the reflections an index-2 multiple adds are one of the superc
 classes, read against the six others as their null (see the
 [report](RUGNUX_REPORT.md#index-2-superstructure)), and the doubled-cell hypothesis merges the larger
 cell and adopts it where that class grows to a quarter of the lattice's intensity or holds it at every
-resolution. An index-3 or index-4 pair is kept at its smaller cell and is not tested. 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.
+resolution. An index-3 or index-4 pair is kept at its smaller cell and is not tested. A row shorter than the 10 Å
+floor of §5.2 needs no hypothesis of its own: the search offers its cell directly (§5.2).
 
 **A near tie is decided on every frame - from the start offline, late on a live run.** A run on a
 stored file has the whole sweep from the start, so its first pass decides each near tie that later

docs/RUGNUX_CLI.md

@@ -247,9 +247,9 @@ rotation explicitly and pick the pass or lattice.
 
 The FFT search looks for cell axes between 10 Å (`--fft-min-unit-cell`) and a default **longest
 axis of 500 Å**, which has no flag of its own — a reference cell (`-C`) moves both bounds to cover
-the cell it names. A de-novo rotation run also tries a second first-pass hypothesis with the short
-end lowered to 5 Å, so a small-molecule cell below the 10 Å floor is indexed on its true axes
-rather than as a supercell of them; the standard pass's answer stands unless that evidence takes it.
+the cell it names. On rotation data a lattice row shorter than the floor, down to 3 Å, is found under
+the multiple of it the search reports there, so a small-molecule cell below the 10 Å floor is
+indexed on its true axes rather than as a supercell of them.
 The long end is extended the same way: wherever the detector separates the lattice rows of an axis
 longer than 500 Å (distance × wavelength / (pixel × spot width)), a coarse search to that length,
 up to 1200 Å, offers its lattices as a further hypothesis, kept only where it puts more of the
@@ -265,7 +265,7 @@ For very long axes the direction grid's angular resolution binds well below 500
 | `-X, --indexing-algorithm <txt>` | `Auto` (default) \| `FFT` \| `cuFFT` \| `FFTW` \| `FFBIDX` \| `None`, case-insensitive. `FFT` is the de-novo FFT search, on the GPU (cuFFT) where there is one and `--no-gpu` is not given, otherwise on the CPU (FFTW); `cuFFT` and `FFTW` name one implementation each; `cuFFT` needs a GPU and is refused without one. `FFBIDX` (known cell, `-C`) runs on the GPU where there is one and otherwise on the CPU. `Auto` takes `FFBIDX` on stills with a known cell, and `FFT` otherwise; rotation data always use the FFT search |
 | `-F` | Short for `-X FFT`, kept for old command lines |
 | `-C, --unit-cell <cell>` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it. A candidate that is an integer supercell of the cell given (a twin's coincidence lattice, a row taken at a multiple of its length) is stepped down to it |
-| `--fft-min-unit-cell <num>` | Shortest cell axis the FFT search accepts, Å (default: 10). A candidate with a shorter axis is discarded — but a de-novo rotation run also runs a **second first-pass hypothesis with the floor lowered to 5 Å by default**, adopted when the standing cell turns out to be an integer supercell of what it finds, so a small-molecule cell is indexed without this flag. `-C` lowers the floor on its own to cover the cell given (and the second hypothesis then stays out of the way). On stills, where neither applies, a crystal below the floor cannot be indexed unless this is lowered |
+| `--fft-min-unit-cell <num>` | Shortest cell axis the FFT search accepts, Å (default: 10). A candidate with a shorter axis is discarded — but on rotation data a row shorter than the floor, down to 3 Å, is recognised under the multiple of it the search reports past the floor, so a small-molecule cell is indexed without this flag. `-C` lowers the floor on its own to cover the cell given. On stills a crystal below the floor cannot be indexed unless this is lowered |
 | `--min-indexed-spots <num>` | Spots a frame must have on the lattice before it counts as indexed (default: 6, minimum 4). It sets the reported indexing rate and the count the rotation first pass ranks candidate lattices by; whether the run has a lattice at all is decided on the pooled spots instead ([§4.1](CPU_DATA_ANALYSIS_INDEXING.md#41-indexed-spot-decision-inlier-test)), and integration is not gated by it |
 | `-S, --space-group <num\|symbol>` | Space group number (`96`) or Hermann-Mauguin symbol (`P43212`) — for indexing and scaling. On serial stills whose group admits two indexings, it also lets the crystals be put on one from the data ([the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)) |
 | `-r, --geom-refine <txt>` | Geometry refinement: `none` \| `orientation` \| `beam_and_lattice` (default) \| `flex` (try all three per image, keep whichever indexes the most spots; alias `multi`) |

docs/RUGNUX_OVERVIEW.md

@@ -30,8 +30,8 @@ the spots against each image's own per-resolution-ring noise
 FFT search looks for periodicity over thousands of directions; candidate cells are
 Niggli-reduced, classified by Bravais lattice, refined both constrained and triclinic, and decided
 on how many validation frames each actually indexes
-([§4–§7](CPU_DATA_ANALYSIS_INDEXING.md)). A de-novo run also tries a second hypothesis with the
-shortest accepted axis lowered from 10 to 5 Å, for small-molecule cells. The first pass is indexed a
+([§4–§7](CPU_DATA_ANALYSIS_INDEXING.md)). A lattice row shorter than the search's 10 Å floor, as
+small-molecule cells have, is recognised under the multiple of it the search reports, down to 3 Å. The first pass is indexed a
 second time at the beam centre the pre-scan measured; the file's centre is kept unless it indexes
 nothing, gives an axis harmonic of the measured centre's lattice, or loses to it when both first passes
 are merged ([§1.4](CPU_DATA_ANALYSIS_IMAGE.md)). A failed pass triggers the discrete rescues — the

docs/RUGNUX_TUTORIAL.md

@@ -214,10 +214,9 @@ rugnux -o xtal sweep_master.h5
 
 What differs from a protein is handled by the run itself:
 
-- **Short cell axes.** The FFT search's floor is 10 Å, but a de-novo rotation run also tries a second
-  first-pass hypothesis with the floor at 5 Å and takes it where the standing cell is an integer
-  supercell of what it finds, so a cell with short axes is indexed on its true axes. An axis shorter
-  than 5 Å needs `--fft-min-unit-cell` lowered, or the cell given with `-C`.
+- **Short cell axes.** The FFT search's floor is 10 Å, but on a rotation run a row shorter than that,
+  down to 3 Å, is recognised under the multiple of it the search reports, so a cell with short axes is
+  indexed on its true axes.
 - **Wide and split spots.** Spots far from the beam, at the high energies small-molecule data are often
   taken at, grow several times wider than the integration disk, and a crystal of slightly misaligned
   domains records each reflection as two or more spots. The pre-scan measures how wide the spots and

image_analysis/geom_refinement/XtalOptimizer.h

@@ -18,7 +18,9 @@ struct XtalOptimizerData {
     gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic;
     int64_t min_spots = 8;
 
-    float min_length_A = 5.0;
+    // Below the shortest row the FFT search reaches (SHORT_ROW_MIN_A), so that it bounds a fit that
+    // runs away rather than a cell the search found.
+    float min_length_A = 2.5;
     float max_length_A = 500.0;
     float min_angle_deg = 60.0f;
     float max_angle_deg = 120.0f;

image_analysis/indexing/FFTIndexer.cpp

@@ -96,7 +96,8 @@ void Sort(Coord &A, Coord &B, Coord &C) {
         std::swap(A, B);
 }
 
-std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &results, bool widen) const {
+std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &results, bool widen,
+                                                      float min_length) const {
     if (results.size() < 3)
         return {};
 
@@ -120,7 +121,7 @@ std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &
                     CrystalLattice reduced = raw.NiggliReduce(); // Reduce cell
                     const auto uc = reduced.GetUnitCell();
 
-                    if (uc.a < min_length_A || uc.b < min_length_A || uc.c < min_length_A)
+                    if (uc.a < min_length || uc.b < min_length || uc.c < min_length)
                         continue;
 
                     float alpha = uc.alpha, beta = uc.beta, gamma = uc.gamma;
@@ -163,7 +164,7 @@ std::vector<CrystalLattice> FFTIndexer::ReduceResults(const std::vector<Coord> &
                 CrystalLattice reduced = CrystalLattice(A, B, C).NiggliReduce();
                 const auto uc = reduced.GetUnitCell();
 
-                if (uc.a < min_length_A || uc.b < min_length_A || uc.c < min_length_A)
+                if (uc.a < min_length || uc.b < min_length || uc.c < min_length)
                     continue;
                 if (uc.alpha < min_angle_deg || uc.alpha > max_angle_deg ||
                     uc.beta  < min_angle_deg || uc.beta  > max_angle_deg ||
@@ -365,11 +366,35 @@ std::vector<Coord> FFTIndexer::FilterFFTResults(size_t max_vectors,
 // strongly. Where the partner is an integer multiple of the row (within a tenth), the row itself
 // explains it as a harmonic and the question is an axis-multiple one that is settled elsewhere, so
 // such a row is left as it is.
+namespace {
+    // The transform at one length, summed over the spots' projections on the row's direction, as the
+    // largest over +-`window_A` about `centre_A` in quarter steps; with the length it was found at.
+    std::pair<double, double> RowAmplitude(const std::vector<float> &proj, double centre_A, double window_A) {
+        double best = 0.0, best_len = centre_A;
+        for (int k = -4; k <= 4; k++) {
+            const double len = centre_A + 0.25 * k * window_A;
+            double re = 0.0, im = 0.0;
+            for (const float p : proj) {
+                const double phase = 2.0 * PI * len * p;
+                re += std::cos(phase);
+                im += std::sin(phase);
+            }
+            const double a = std::hypot(re, im);
+            if (a > best) {
+                best = a;
+                best_len = len;
+            }
+        }
+        return {best, best_len};
+    }
+
+    // The sum is a test between lengths, not a measurement, so a sample of the spots decides it as well
+    // as all of them: the coherent part grows as N and the incoherent one as sqrt(N).
+    constexpr size_t ROW_TEST_MAX_SPOTS = 4096;
+}
+
 void FFTIndexer::DeAlias(const std::vector<Coord> &coord, size_t nspots, std::vector<Coord> &rows) const {
-    // The sum is a test between two lengths, not a measurement, so a sample of the spots decides it as
-    // well as all of them: the coherent part grows as N and the incoherent one as sqrt(N).
-    constexpr size_t MAX_SPOTS = 4096;
-    const size_t stride = std::max<size_t>(1, nspots / MAX_SPOTS);
+    const size_t stride = std::max<size_t>(1, nspots / ROW_TEST_MAX_SPOTS);
     // The FFT gives a row's length to one bin, 2*max_length_A/histogram_size; the partner inherits that
     // error, and the spots' phases stay coherent only over a fraction of it at high resolution. Search
     // both lengths over +-1 bin in quarter-bin steps.
@@ -377,20 +402,7 @@ void FFTIndexer::DeAlias(const std::vector<Coord> &coord, size_t nspots, std::ve
     std::vector<float> proj;
     proj.reserve(nspots / stride + 1);
 
-    const auto amplitude = [&](double centre_A) {
-        double best = 0.0;
-        for (int k = -4; k <= 4; k++) {
-            const double len = centre_A + 0.25 * k * bin_A;
-            double re = 0.0, im = 0.0;
-            for (const float p : proj) {
-                const double phase = 2.0 * PI * len * p;
-                re += std::cos(phase);
-                im += std::sin(phase);
-            }
-            best = std::max(best, std::hypot(re, im));
-        }
-        return best;
-    };
+    const auto amplitude = [&](double centre_A) { return RowAmplitude(proj, centre_A, bin_A).first; };
 
     for (auto &row : rows) {
         const double len = row.Length();
@@ -414,6 +426,56 @@ void FFTIndexer::DeAlias(const std::vector<Coord> &coord, size_t nspots, std::ve
     }
 }
 
+// The per-direction peak search skips every length up to min_length_A, so a lattice row shorter than
+// that is reported at a multiple of itself past it - an exact harmonic - and the cell built from such
+// rows is an exact integer supercell of the crystal's. Small molecules have rows like that; no
+// protein does. The spots can tell the two apart where the histogram cannot: summed at the row's
+// own length every spot is in phase, whichever of the two the row is, but at a whole fraction of
+// that length only the spots of a lattice whose row it is are, and the others cancel. So each row
+// reported within a few multiples of the floor is asked at every fraction of its length between
+// SHORT_ROW_MIN_A and the floor, and takes the shortest one that keeps at least half its amplitude.
+// The transform of a lattice comb has all its harmonics at one height, less what positional scatter
+// takes from the higher ones, so the true row is the strongest of them, not merely half as strong.
+bool FFTIndexer::Fundamentals(const std::vector<Coord> &coord, size_t nspots, std::vector<Coord> &rows) const {
+    // A harmonic past the first few multiples past the floor is damped below the ones before it, so
+    // the peak search does not report a short row there.
+    constexpr float MAX_ROW_OVER_FLOOR = 3.0f;
+    constexpr double MIN_AMPLITUDE_RATIO = 0.5;
+
+    const size_t stride = std::max<size_t>(1, nspots / ROW_TEST_MAX_SPOTS);
+    const double bin_A = 2.0 * static_cast<double>(max_length_A) / static_cast<double>(histogram_size);
+    std::vector<float> proj;
+    proj.reserve(nspots / stride + 1);
+
+    bool changed = false;
+    for (auto &row : rows) {
+        const double len = row.Length();
+        if (!(len > min_length_A) || len > MAX_ROW_OVER_FLOOR * min_length_A)
+            continue;
+        const Coord dir = row.Normalize();
+        proj.clear();
+        for (size_t i = 0; i < nspots; i += stride)
+            proj.push_back(dir * coord[i]);
+        const auto [a_row, len_row] = RowAmplitude(proj, len, bin_A);
+        double fundamental = 0.0;
+        for (int k = 2; len_row / k >= SHORT_ROW_MIN_A; k++) {
+            if (len_row / k >= min_length_A)
+                continue;
+            // The row's length is known to a bin, its k-th fraction to a k-th of one.
+            const auto [a, l] = RowAmplitude(proj, len_row / k, bin_A / k);
+            if (a >= MIN_AMPLITUDE_RATIO * a_row)
+                fundamental = l;
+        }
+        if (fundamental > 0.0) {
+            Logger("FFTIndexer").Debug("FFT row of {:.2f} A is a multiple of a {:.2f} A row below the {:.1f} A floor",
+                                       len, fundamental, min_length_A);
+            row = dir * static_cast<float>(fundamental);
+            changed = true;
+        }
+    }
+    return changed;
+}
+
 float FFTIndexer::IndexedFraction(const CrystalLattice &latt,
                                   const std::vector<Coord> &coord, size_t nspots) const {
     if (nspots == 0)
@@ -437,8 +499,13 @@ float FFTIndexer::IndexedFraction(const CrystalLattice &latt,
 }
 
 std::vector<CrystalLattice> FFTIndexer::ReduceAndRefine(const std::vector<Coord> &coord, size_t nspots,
-                                                        const std::vector<Coord> &filtered, bool widen) {
-    const auto r = ReduceResults(filtered, widen);
+                                                        const std::vector<Coord> &filtered, bool widen,
+                                                        const std::vector<Coord> &short_rows) {
+    auto r = ReduceResults(filtered, widen, min_length_A);
+    // Widened: the short rows a small cell has many combinations of within one plane, enough to fill
+    // the standard reduction's window and leave the third axis out of every triple it forms.
+    for (auto &l : ReduceResults(short_rows, true, SHORT_ROW_MIN_A))
+        r.push_back(std::move(l));
 
     Eigen::MatrixX3<float> oCell(r.size() * 3u, 3u);
     Eigen::VectorX<float> scores(r.size());
@@ -465,7 +532,7 @@ std::vector<CrystalLattice> FFTIndexer::ReduceAndRefine(const std::vector<Coord>
         .dist_tolerance_vs_reference = dist_tolerance_vs_reference,
         .step_down_to_reference = step_down_to_reference,
         .reference_unit_cell = reference_unit_cell,
-        .min_length_A = min_length_A,
+        .min_length_A = short_rows.empty() ? min_length_A : SHORT_ROW_MIN_A,
         // A row DeAlias unfolded is longer than max_length_A by construction, and shorter than twice it.
         .max_length_A = std::any_of(filtered.begin(), filtered.end(),
                                     [&](const Coord &v) { return v.Length() > max_length_A; })
@@ -623,7 +690,25 @@ std::vector<CrystalLattice> FFTIndexer::RunInternal(const std::vector<Coord> &co
 
     auto rows = shortlist;
     DeAlias(coord, nspots, rows);
-    auto lattices = ReduceAndRefine(coord, nspots, rows, false);
+    // The cells of the short rows are offered beside the standard ones, not instead of them: the
+    // refinement ranks them all on the spots they index, and where the two index alike it keeps the
+    // smaller cell.
+    std::vector<Coord> short_rows;
+    if (search_short_rows) {
+        auto replaced = rows;
+        if (Fundamentals(coord, nspots, replaced)) {
+            // Shortest first, as ReduceResults combines the first few; keys computed once (see
+            // FilterFFTResults).
+            std::vector<std::pair<float, size_t>> by_length;
+            for (size_t i = 0; i < replaced.size(); i++)
+                by_length.emplace_back(replaced[i].Length(), i);
+            std::stable_sort(by_length.begin(), by_length.end(),
+                             [](const auto &a, const auto &b) { return a.first < b.first; });
+            for (const auto &[len, i] : by_length)
+                short_rows.push_back(replaced[i]);
+        }
+    }
+    auto lattices = ReduceAndRefine(coord, nspots, rows, false, short_rows);
 
     // If the best cell indexes few of the (un-refined) accumulated spots, the true cell may be large/
     // elongated with a long axis beyond the standard triple window (a superstructure, or a satellite-

image_analysis/indexing/FFTIndexer.h

@@ -13,6 +13,10 @@
 
 #define FFT_MAX_SPOTS (64*1024)
 
+// The shortest lattice row the short-row search looks for (FFTIndexer::Fundamentals) - the minimum
+// cell DIALS' fft3d searches down to.
+constexpr float SHORT_ROW_MIN_A = 3.0f;
+
 class FFTIndexer : public Indexer {
 
 protected:
@@ -37,7 +41,7 @@ protected:
     // widen=false: combine only the shortest few filtered vectors (fast, correct for compact cells).
     // widen=true: anchor the two short axes but let the third range over all vectors (reaches the long
     // axis of a large/elongated cell). Used only as a fallback when the standard reduction indexes poorly.
-    std::vector<CrystalLattice> ReduceResults(const std::vector<Coord> &results, bool widen) const;
+    std::vector<CrystalLattice> ReduceResults(const std::vector<Coord> &results, bool widen, float min_length) const;
     // `magnitudes`, when given, receives the FFT peak prominence of each returned row, in the
     // same order. The spindle severity needs it to tell a crystal's real rows from the short
     // spurious ones a long cell produces; nothing else does, so it is optional.
@@ -45,8 +49,14 @@ protected:
                                         std::vector<float> *magnitudes = nullptr) const;
     // Replace each row whose FFT length is an alias of a longer one by that longer row (see the definition).
     void DeAlias(const std::vector<Coord> &coord, size_t nspots, std::vector<Coord> &rows) const;
+    // Replace each row reported past min_length_A at a multiple of a shorter row by that row (see the
+    // definition); true where any was replaced.
+    bool Fundamentals(const std::vector<Coord> &coord, size_t nspots, std::vector<Coord> &rows) const;
+    // `short_rows`, when not empty, are reduced as well (widened, down to SHORT_ROW_MIN_A) and their cells
+    // refined and ranked together with those of `filtered`.
     std::vector<CrystalLattice> ReduceAndRefine(const std::vector<Coord> &coord, size_t nspots,
-                                                const std::vector<Coord> &filtered, bool widen);
+                                                const std::vector<Coord> &filtered, bool widen,
+                                                const std::vector<Coord> &short_rows = {});
     float IndexedFraction(const CrystalLattice &latt, const std::vector<Coord> &coord, size_t nspots) const;
     void SetupDirectionVectors();
     // Replace the search directions with a small cap around `axis` (same count, so the transform and

image_analysis/indexing/Indexer.cpp

@@ -8,6 +8,7 @@ void Indexer::Setup(const DiffractionExperiment& experiment) {
     indexing_tolerance = experiment.GetIndexingSettings().GetTolerance();
     dist_tolerance_vs_reference = experiment.GetIndexingSettings().GetUnitCellDistTolerance();
     step_down_to_reference = experiment.GetIndexingSettings().GetRotationIndexing();
+    search_short_rows = experiment.IsRotationIndexing();
     viable_cell_min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots();
     index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings();
     wavelength_A = experiment.GetWavelength_A();

image_analysis/indexing/Indexer.h

@@ -40,6 +40,10 @@ protected:
     // crystal form among several in the stream, and a sub-lattice of a frame's supercell matches a
     // NEIGHBOURING form within the tolerance as readily as the given one.
     bool step_down_to_reference = false;
+    // Rotation data only (a scanning axis, so not a grid scan either): rows shorter than the FFT floor
+    // are looked for under the rows the transform reports past it (FFTIndexer::Fundamentals). Stills
+    // keep the search they had; nothing has been measured on them.
+    bool search_short_rows = false;
     bool index_ice_rings = true;
 
     DiffractionGeometry geom;

image_analysis/indexing/PostIndexingRefinement.cpp

@@ -245,6 +245,17 @@ namespace {
     }
 }
 
+// One lattice contains the other: the rows of the larger cell are integer combinations of the rows
+// of the smaller, and the volumes differ by an integer factor of at least two.
+bool IntegerRelated(const Eigen::Matrix3f &a, const Eigen::Matrix3f &b) {
+    const bool a_larger = std::abs(a.determinant()) > std::abs(b.determinant());
+    const Eigen::Matrix3f &large = a_larger ? a : b, &small = a_larger ? b : a;
+    if (std::abs(large.determinant()) < 1.5f * std::abs(small.determinant()))
+        return false;
+    const Eigen::Matrix3f m = large * small.inverse();
+    return (m - m.array().round().matrix()).cwiseAbs().maxCoeff() < 0.1f;
+}
+
 std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
                                    size_t nspots,
                                    Eigen::MatrixX3<float> &oCell,
@@ -393,24 +404,38 @@ std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
               });
 
     std::vector<RefinedCandidate> accepted;
+    // Candidates the overlap test drops only for their overlap with a cell they are an integer multiple
+    // or sub-cell of: one of those two lattices contains the other, so that overlap is total by
+    // construction and says nothing about which of them is the crystal's. They are kept, in rank order,
+    // after every cell the test accepts - which come first and in the order they always had - for the
+    // caller to decide on.
+    std::vector<RefinedCandidate> multiples;
 
     for (const auto &candidate: candidates) {
-        int64_t overlap = 0;
+        int64_t overlap = 0, overlap_unrelated = 0;
 
         // Check all already selected lattices and see how many spots are already indexed for the candidate
         // If the overlap is more than 40% of indexed spots - we assume the lattice doesn't bring anything new
         for (const auto &selected: accepted) {
+            int64_t n = 0;
             for (size_t i = 0; i < candidate.indexed_mask.size(); ++i) {
                 if (candidate.indexed_mask[i] && selected.indexed_mask[i])
-                    overlap++;
+                    n++;
             }
+            overlap += n;
+            if (!IntegerRelated(candidate.cell, selected.cell))
+                overlap_unrelated += n;
         }
 
-        if (overlap < static_cast<int64_t>(REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD
-                                           * static_cast<float>(candidate.indexed_spot_count))) {
+        const auto bound = static_cast<int64_t>(REFINE_CANDIDATE_OVERLAP_RATIO_THRESHOLD
+                                                * static_cast<float>(candidate.indexed_spot_count));
+        if (overlap < bound)
             accepted.emplace_back(candidate);
-        }
+        else if (overlap_unrelated < bound)
+            multiples.emplace_back(candidate);
     }
+    for (auto &m : multiples)
+        accepted.push_back(std::move(m));
 
     ret.reserve(accepted.size());
 

image_analysis/indexing/PostIndexingRefinement.h

@@ -34,6 +34,8 @@ struct RefineParameters {
     unsigned refine_threads = 1;
 };
 
+bool IntegerRelated(const Eigen::Matrix3f &a, const Eigen::Matrix3f &b);
+
 std::vector<CrystalLattice> Refine(const std::vector<Coord> &in_spots,
                                           size_t nspots,
                                           Eigen::MatrixX3<float> &oCell,

image_analysis/rotation_indexer/RotationIndexer.cpp

@@ -8,6 +8,7 @@
 #include "../indexing/FFTIndexer.h"
 #include "../lattice_search/LatticeSearch.h"
 #include "../indexing/MultiLatticeSearch.h"
+#include "../indexing/PostIndexingRefinement.h"
 
 #include <future>
 #include <mutex>
@@ -203,6 +204,8 @@ namespace {
         lattice(kb, b.indexed_lattice);
         for (const auto &l : a.extra_lattices) lattice(ka, l);
         for (const auto &l : b.extra_lattices) lattice(kb, l);
+        for (const auto &m : a.multiples) lattice(ka, m->lattice);
+        for (const auto &m : b.multiples) lattice(kb, m->lattice);
         return ka == kb && a.tilt_walk.has_value() == b.tilt_walk.has_value()
                && (a.unconstrained != nullptr) == (b.unconstrained != nullptr);
     }
@@ -236,7 +239,8 @@ std::vector<double> RotationIndexer::InputKey() const {
 }
 
 RotationIndexer::IndexingOutcome RotationIndexer::Outcome() const {
-    return {axis_, updated_geom_, search_result_, extra_lattices_, unconstrained_, tilt_walk_, indexed_lattice};
+    return {axis_, updated_geom_, search_result_, extra_lattices_, unconstrained_, tilt_walk_, indexed_lattice,
+            multiples_};
 }
 
 void RotationIndexer::RunIndexing() {
@@ -262,6 +266,7 @@ void RotationIndexer::RunIndexing() {
         unconstrained_ = kept->unconstrained;
         tilt_walk_ = kept->tilt_walk;
         indexed_lattice = kept->indexed_lattice;
+        multiples_ = kept->multiples;
         return;
     }
 
@@ -298,6 +303,7 @@ void RotationIndexer::RunIndexingLocked() {
     if (indexer_result.error)
         std::rethrow_exception(indexer_result.error);
     tilt_walk_.reset();
+    multiples_.clear();
 
     if (!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) {
         DiffractionExperiment experiment_copy(experiment);
@@ -440,6 +446,7 @@ void RotationIndexer::RunIndexingLocked() {
         XtalOptimizerData best_data;
         LatticeSearchResult best_sr;
         std::shared_ptr<RotationIndexerResult> best_alt;
+        std::vector<std::pair<XtalOptimizerData, LatticeSearchResult>> solved;
         for (size_t ci = 0; ci < n_try; ci++) {
             if (!work[ci].viable)
                 continue;
@@ -488,6 +495,7 @@ void RotationIndexer::RunIndexingLocked() {
 
             if (!ok)
                 continue;
+            solved.push_back({data, sr});
 
             // Prefer the indexer's earlier (primary) candidate; adopt a later one only if it indexes
             // clearly more AND indexes reasonably well in absolute terms. The absolute floor stops a
@@ -733,6 +741,21 @@ void RotationIndexer::RunIndexingLocked() {
             updated_geom_ = best_data.geom;
             axis_ = best_data.axis;
             unconstrained_ = std::move(best_alt);
+
+            const auto rows = [](const CrystalLattice &l) {
+                Eigen::Matrix3f m;
+                m << l.Vec0().x, l.Vec0().y, l.Vec0().z, l.Vec1().x, l.Vec1().y, l.Vec1().z,
+                     l.Vec2().x, l.Vec2().y, l.Vec2().z;
+                return m;
+            };
+            const CrystalLattice best_prim = best_data.latt.ToPrimitive(best_sr.centering);
+            for (const auto &[d, s] : solved) {
+                const CrystalLattice prim = d.latt.ToPrimitive(s.centering);
+                if (std::abs(prim.CalcVolume()) > 1.5f * std::abs(best_prim.CalcVolume())
+                    && IntegerRelated(rows(prim), rows(best_prim)))
+                    multiples_.push_back(std::make_shared<RotationIndexerResult>(RotationIndexerResult{
+                        .lattice = d.latt, .search_result = s, .geom = d.geom, .axis = d.axis}));
+            }
         }
 
         // Extra (twin) lattices: MultiLatticeSearch derives each rotation by relating the FFT's primary
@@ -889,6 +912,7 @@ std::optional<RotationIndexerResult> RotationIndexer::GetLattice() const {
         .unconstrained = unconstrained_,
         .tilt_walk = tilt_walk_,
         .twin_overlap = twin_overlap_,
+        .multiples = multiples_,
     };
 }
 
@@ -902,6 +926,7 @@ void RotationIndexer::ForceResult(const RotationIndexerResult &result) {
     unconstrained_ = result.unconstrained;
     tilt_walk_ = result.tilt_walk;
     twin_overlap_ = result.twin_overlap;
+    multiples_ = result.multiples;
 }
 
 bool RotationIndexer::AccumulationFull() const {
@@ -912,6 +937,7 @@ bool RotationIndexer::AccumulationFull() const {
 void RotationIndexer::ForceLattice(const CrystalLattice &lattice) {
     indexed_lattice = lattice;
     twin_overlap_.reset();
+    multiples_.clear();
     const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1();
     search_result_ = LatticeSearchResult{
         .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one

image_analysis/rotation_indexer/RotationIndexer.h

@@ -51,6 +51,11 @@ struct RotationIndexerResult {
         float fraction = 0.0f;   // the other domain's share of the crystal (the twin fraction)
     };
     std::optional<TwinOverlap> twin_overlap;
+    // The candidates the FFT offered beside the chosen cell that are integer supercells of it, each as
+    // refined in the selection. Not chosen - their cell indexes the cloud no better than this one - but
+    // carried, because some questions are asked of the larger lattice (the twin-composite hypothesis:
+    // the coincidence lattice of two domains is a supercell of either domain's cell).
+    std::vector<std::shared_ptr<RotationIndexerResult>> multiples;
 };
 
 // What a memo verification throws where a reused result differs from its recomputation - its own
@@ -76,6 +81,7 @@ public:
         std::shared_ptr<RotationIndexerResult> unconstrained;
         std::optional<RotationIndexerResult::TiltWalk> tilt_walk;
         std::optional<CrystalLattice> indexed_lattice;
+        std::vector<std::shared_ptr<RotationIndexerResult>> multiples;
     };
 private:
     mutable std::mutex m;
@@ -106,6 +112,7 @@ private:
     std::shared_ptr<RotationIndexerResult> unconstrained_;
     std::optional<RotationIndexerResult::TiltWalk> tilt_walk_;
     std::optional<RotationIndexerResult::TwinOverlap> twin_overlap_;   // only ever forced, see ForceResult
+    std::vector<std::shared_ptr<RotationIndexerResult>> multiples_;
 
     IndexerThreadPool &indexer_;
 

rugnux/Rugnux.h

@@ -342,8 +342,7 @@ double ValidationEvidenceLog10Tail(const ValidationSpotEvidence &e);
 // The significance the run's no-crystal test (RejectChanceLattice) asks of the lattice the first
 // pass ended on - Z ~ 6 instead of the 0.001 the comparisons between candidate lattices use. The
 // binomial null holds for a lattice chosen independently of the validation frames, and this one is
-// not: it is the winner of some fifty indexing attempts (schemes, ladder rungs, beam centres, the
-// short-axis pass), each fitting the orientation, cell and geometry to an accumulated cloud that
+// not: it is the winner of some fifty indexing attempts (schemes, ladder rungs, beam centres), each fitting the orientation, cell and geometry to an accumulated cloud that
 // already holds the long-lived reflections it then re-finds on the validation frames, while the
 // displaced arrangements get none of that selection. That makes the 0.001 null optimistic by about
 // the number of attempts, ~1.7 decades, before the fit's own freedom is counted. A real crystal's

rugnux/RugnuxFirstPass.cpp

@@ -479,7 +479,6 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) {
         return false;
 
     run.RequireLattice(best);
-    best = run.ShortAxisHypothesis(std::move(best));
     best = run.SubLatticeHypothesis(std::move(best));   // and the twin-composite hypothesis
     if (!p.geometry_prepass && prepass_result_.has_value())
         best = run.Pass1Fallback(std::move(best));
@@ -1435,14 +1434,6 @@ void Rugnux::FirstPassRun::PrefetchCertainFrames() {
     PrefetchSpots(upfront);
 }
 
-IndexingSettings Rugnux::FirstPassRun::ShortAxisSettings(const DiffractionExperiment &x) const {
-    auto settings = x.GetIndexingSettings();
-    settings.FFT_MinUnitCell_A(IndexingSettings::fft_min_unit_cell_limit_A);
-    settings.IndexingThreads(rx.FirstPassIndexerThreads());
-    settings.RefineThreads(std::clamp(rx.config_.nthreads / 2, 1, MAX_REFINE_THREADS));
-    return settings;
-}
-
 // The beam-centre check's pass at the measured centre, started beside the file centre's (see
 // centre_spec).
 void Rugnux::FirstPassRun::StartCentreSpeculation() {
@@ -1507,31 +1498,8 @@ void Rugnux::FirstPassRun::StartCentreSpeculation() {
     });
 }
 
-// The short-axis pass, started beside the standard indexing (see short_spec).
-void Rugnux::FirstPassRun::StartShortAxisSpeculation() {
-    if (rx.cancelled_ || rx.experiment_.GetUnitCell().has_value())
-        return;
-    auto &sp = short_spec.emplace();
-    sp.key = rx.ExperimentKey();
-    sp.mapping = spot_mapping;
-    sp.spots_per_image = first_pass_spots_per_image;
-    sp.experiment = std::make_unique<DiffractionExperiment>(rx.experiment_);
-    sp.pool = std::make_unique<IndexerThreadPool>(ShortAxisSettings(*sp.experiment),
-                                                  IndexerConstruction::OnFirstUse, rx.gpu_work_beside_);
-    for (const auto &[name, ordinals] : schemes)
-        sp.sris.push_back(FeedScheme(*sp.pool, ordinals, sp.experiment.get()));
-    for (auto &ri : sp.sris) {
-        RotationIndexer *rp = ri.get();
-        sp.futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
-    }
-}
-
 auto Rugnux::FirstPassRun::IndexSchemes() -> FirstPass {
-    // The short axis first: the measured centre's pass waits for the beam-centre capture.
-    return PickBest(*indexer_pool, *indexer, [&] {
-        StartShortAxisSpeculation();
-        StartCentreSpeculation();
-    });
+    return PickBest(*indexer_pool, *indexer, [&] { StartCentreSpeculation(); });
 }
 
 // Twinning by reticular merohedry: two domains of one lattice whose nodes coincide on a sub-lattice
@@ -1562,81 +1530,89 @@ auto Rugnux::FirstPassRun::IndexSchemes() -> FirstPass {
 auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass {
     if (rx.cancelled_ || !best.result.has_value() || best.score <= 0)
         return best;
-    const CrystalLattice prim = best.result->lattice
-            .ToPrimitive(best.result->search_result.centering).NiggliReduce();
-    std::vector<size_t> kept;
-    const auto hkl = IndexedMillerIndices(prim, ValidationCloud(*best.result), index_tol, &kept);
-    if (hkl.empty())
-        return best;
-    // The spots' intensities in the cloud's order (ValidationCloud), for the twin fraction.
-    std::vector<float> cloud_intensity, intensity;
-    for (const int ordinal : validation)
-        if (const auto it = spot_cache.find(ordinal); it != spot_cache.end())
-            for (const auto &sp : it->second)
-                cloud_intensity.push_back(sp.intensity);
-    for (const size_t i : kept)
-        if (i < cloud_intensity.size())
-            intensity.push_back(cloud_intensity[i]);
-    const PooledEvidence at_standing = MeasurePooled(*indexer, *best.result);
-    const double chance_share = static_cast<double>(at_standing.by_chance)
-                                / static_cast<double>(std::max<int64_t>(1, at_standing.on_lattice));
-    for (int n = 2; n <= TWIN_COMPOSITE_MAX_INDEX; n++) {
-        const auto tc = TwinCompositeAt(prim, hkl, n, chance_share, intensity);
-        if (!tc)
+    // The standing cell is asked, and then each integer supercell of it the indexer carried: where the
+    // search offers a domain's own cell, that cell indexes as well as the coincidence lattice and stands,
+    // and the coincidence lattice is one of its multiples.
+    std::vector<RotationIndexerResult> composites{*best.result};
+    for (const auto &m : best.result->multiples)
+        composites.push_back(*m);
+    for (const auto &composite : composites) {
+        const CrystalLattice prim = composite.lattice
+                .ToPrimitive(composite.search_result.centering).NiggliReduce();
+        const double vol = std::abs(prim.CalcVolume());
+        std::vector<size_t> kept;
+        const auto hkl = IndexedMillerIndices(prim, ValidationCloud(composite), index_tol, &kept);
+        if (hkl.empty())
             continue;
-        logger.Info("Twin-composite hypothesis, index {}: of {} validation spots on the standing cell "
-                    "({:.0f} A^3), the congruent sub-lattices ({},{},{}) and ({},{},{}) hold {} and {} alone, {} "
-                    "together and {} on neither ({:.1f} by chance) - per node {:.1f} on the emptiest of the other "
-                    "classes against {:.1f} on neither => {}", n, hkl.size(), best.vol,
-                    tc->f1[0], tc->f1[1], tc->f1[2], tc->f2[0], tc->f2[1], tc->f2[2], tc->only1, tc->only2,
-                    tc->both, tc->neither, tc->neither_expected, tc->other_density, tc->neither_density,
-                    tc->composite ? "two domains on their coincidence lattice" : "one lattice");
-        if (!tc->composite)
-            continue;
-        RotationIndexerResult cand = *best.result;
-        cand.search_result = LatticeSearch(tc->domain);
-        cand.lattice = cand.search_result.conventional;
-        cand.unconstrained.reset();
-        cand.extra_lattices.clear();
-        // The overlap, the reflections on both sub-lattices, holds both domains' intensity: the
-        // stronger domain's reflections h on the weaker's sub-lattice are those with h . t divisible
-        // by n, t the weaker sub-lattice's form taken as a vector of the standing cell.
-        cand.twin_overlap = RotationIndexerResult::TwinOverlap{
-            .t = tc->f2[0] * prim.Vec0() + tc->f2[1] * prim.Vec1() + tc->f2[2] * prim.Vec2(),
-            .n = n,
-            .rotation = tc->rotation,
-            .fraction = static_cast<float>(tc->fraction)};
-        const int score = CountIndexed(*indexer, cand);
-        // The composite is the evidence that the domain's cell is this crystal's; what is left to ask
-        // of the cell is that it is a lattice at all - that the validation spots it holds beat what a
-        // wrong spindle angle puts on it. Not a frame count: on a sparse pattern with two domains in
-        // the beam a frame holds a handful of each domain's spots, and a bar on frames is a bar on
-        // how sparse the pattern is.
-        // And the domain must index a majority of the frames or as many as the standing cell, the bar the
-        // sub-lattice hypothesis asks.
-        const int majority = static_cast<int>(validation.size()) / 2;
-        const bool adopt = (score > majority || score >= best.score)
-                           && ValidationEvidenceBeatsChance(MeasurePooled(*indexer, cand));
-        const auto uc = cand.lattice.GetUnitCell();
-        logger.Info("Twin-composite hypothesis: the stronger domain's cell, {}-centred {} {:.2f} {:.2f} "
-                    "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the {:.0f}% of its "
-                    "reflections the other domain overlaps, h.({:.2f},{:.2f},{:.2f}) divisible by {}, is left out "
-                    "of integration{}; twin fraction {:.2f} from the intensity of the spots on one domain alone)",
-                    cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system),
-                    uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score,
-                    static_cast<int>(validation.size()), adopt ? "adopting it" : "keeping the standing cell",
-                    100.0 * tc->p_both * n, cand.twin_overlap->t * cand.lattice.Astar(),
-                    cand.twin_overlap->t * cand.lattice.Bstar(), cand.twin_overlap->t * cand.lattice.Cstar(), n,
-                    tc->rotation ? fmt::format(", but for the axial ones whose twin partner is the reflection itself "
-                                               "or its Friedel mate - the twin law a rotation of {:.1f} deg",
-                                               std::acos(std::clamp((tc->rotation->Column(0).x + tc->rotation->Column(1).y
-                                                                     + tc->rotation->Column(2).z - 1.0f) / 2.0f, -1.0f, 1.0f))
-                                               * 180.0f / static_cast<float>(PI))
-                                 : std::string(), tc->fraction);
-        if (adopt) {
-            best = FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name),
-                             std::abs(tc->domain.CalcVolume())};
-            break;
+        // The spots' intensities in the cloud's order (ValidationCloud), for the twin fraction.
+        std::vector<float> cloud_intensity, intensity;
+        for (const int ordinal : validation)
+            if (const auto it = spot_cache.find(ordinal); it != spot_cache.end())
+                for (const auto &sp : it->second)
+                    cloud_intensity.push_back(sp.intensity);
+        for (const size_t i : kept)
+            if (i < cloud_intensity.size())
+                intensity.push_back(cloud_intensity[i]);
+        const PooledEvidence at_standing = MeasurePooled(*indexer, composite);
+        const double chance_share = static_cast<double>(at_standing.by_chance)
+                                    / static_cast<double>(std::max<int64_t>(1, at_standing.on_lattice));
+        for (int n = 2; n <= TWIN_COMPOSITE_MAX_INDEX; n++) {
+            const auto tc = TwinCompositeAt(prim, hkl, n, chance_share, intensity);
+            if (!tc)
+                continue;
+            logger.Info("Twin-composite hypothesis, index {}: of {} validation spots on the cell "
+                        "({:.0f} A^3), the congruent sub-lattices ({},{},{}) and ({},{},{}) hold {} and {} alone, {} "
+                        "together and {} on neither ({:.1f} by chance) - per node {:.1f} on the emptiest of the other "
+                        "classes against {:.1f} on neither => {}", n, hkl.size(), vol,
+                        tc->f1[0], tc->f1[1], tc->f1[2], tc->f2[0], tc->f2[1], tc->f2[2], tc->only1, tc->only2,
+                        tc->both, tc->neither, tc->neither_expected, tc->other_density, tc->neither_density,
+                        tc->composite ? "two domains on their coincidence lattice" : "one lattice");
+            if (!tc->composite)
+                continue;
+            RotationIndexerResult cand = composite;
+            cand.multiples.clear();
+            cand.search_result = LatticeSearch(tc->domain);
+            cand.lattice = cand.search_result.conventional;
+            cand.unconstrained.reset();
+            cand.extra_lattices.clear();
+            // The overlap, the reflections on both sub-lattices, holds both domains' intensity: the
+            // stronger domain's reflections h on the weaker's sub-lattice are those with h . t divisible
+            // by n, t the weaker sub-lattice's form taken as a vector of the standing cell.
+            cand.twin_overlap = RotationIndexerResult::TwinOverlap{
+                .t = tc->f2[0] * prim.Vec0() + tc->f2[1] * prim.Vec1() + tc->f2[2] * prim.Vec2(),
+                .n = n,
+                .rotation = tc->rotation,
+                .fraction = static_cast<float>(tc->fraction)};
+            const int score = CountIndexed(*indexer, cand);
+            // The composite is the evidence that the domain's cell is this crystal's; what is left to ask
+            // of the cell is that it is a lattice at all - that the validation spots it holds beat what a
+            // wrong spindle angle puts on it. Not a frame count: on a sparse pattern with two domains in
+            // the beam a frame holds a handful of each domain's spots, and a bar on frames is a bar on
+            // how sparse the pattern is.
+            // And the domain must index a majority of the frames or as many as the standing cell, the bar the
+            // sub-lattice hypothesis asks.
+            const int majority = static_cast<int>(validation.size()) / 2;
+            const bool adopt = (score > majority || score >= best.score)
+                               && ValidationEvidenceBeatsChance(MeasurePooled(*indexer, cand));
+            const auto uc = cand.lattice.GetUnitCell();
+            logger.Info("Twin-composite hypothesis: the stronger domain's cell, {}-centred {} {:.2f} {:.2f} "
+                        "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the {:.0f}% of its "
+                        "reflections the other domain overlaps, h.({:.2f},{:.2f},{:.2f}) divisible by {}, is left out "
+                        "of integration{}; twin fraction {:.2f} from the intensity of the spots on one domain alone)",
+                        cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system),
+                        uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score,
+                        static_cast<int>(validation.size()), adopt ? "adopting it" : "keeping the standing cell",
+                        100.0 * tc->p_both * n, cand.twin_overlap->t * cand.lattice.Astar(),
+                        cand.twin_overlap->t * cand.lattice.Bstar(), cand.twin_overlap->t * cand.lattice.Cstar(), n,
+                        tc->rotation ? fmt::format(", but for the axial ones whose twin partner is the reflection itself "
+                                                   "or its Friedel mate - the twin law a rotation of {:.1f} deg",
+                                                   std::acos(std::clamp((tc->rotation->Column(0).x + tc->rotation->Column(1).y
+                                                                         + tc->rotation->Column(2).z - 1.0f) / 2.0f, -1.0f, 1.0f))
+                                                   * 180.0f / static_cast<float>(PI))
+                                     : std::string(), tc->fraction);
+            if (adopt)
+                return FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name),
+                                 std::abs(tc->domain.CalcVolume())};
         }
     }
     return best;

rugnux/RugnuxFirstPass.h

@@ -73,7 +73,6 @@ struct Rugnux::FirstPassRun {
     void ReleaseSpotEngines();
     // Throws where the pass ends with no lattice at all.
     void RequireLattice(const FirstPass &best) const;
-    FirstPass ShortAxisHypothesis(FirstPass best);
     FirstPass SubLatticeHypothesis(FirstPass best);
     FirstPass TwinCompositeHypothesis(FirstPass best);
     FirstPass Pass1Fallback(FirstPass best);
@@ -143,9 +142,7 @@ struct Rugnux::FirstPassRun {
     FirstPass PickBest(IndexerThreadPool &pool, IndexAndRefine &idx,
                        const std::function<void()> &alongside = {});
     void ApplyState(const FirstPassState &st);
-    IndexingSettings ShortAxisSettings(const DiffractionExperiment &x) const;
     void StartCentreSpeculation();
-    void StartShortAxisSpeculation();
 
     // --- Near ties ---
     //
@@ -307,21 +304,6 @@ struct Rugnux::FirstPassRun {
     const std::optional<float> first_pass_file_d_min;
     FirstPassState committed;
 
-    // The short-axis pass below indexes the same schemes at a lower floor, after the rescues in
-    // between - and on nearly every run they change nothing it reads. So it is started here, on a
-    // copy of the experiment, beside the standard indexing it would otherwise wait for, and taken
-    // at its place only if what it read is still what the run has there (ExperimentKey, the spot
-    // list's mapping and budget); otherwise it is run again there as before.
-    struct ShortAxisSpeculation {
-        std::vector<double> key;
-        const AzimuthalIntegrationMapping *mapping = nullptr;
-        size_t spots_per_image = 0;
-        std::unique_ptr<DiffractionExperiment> experiment;
-        std::unique_ptr<IndexerThreadPool> pool;
-        std::vector<std::unique_ptr<RotationIndexer>> sris;
-        std::vector<std::future<void>> futures;   // last, so it is waited for before the rest goes
-    };
-    std::optional<ShortAxisSpeculation> short_spec;
     // The beam-centre check below runs a second first pass at the centre measured from the
     // background, and it reads the first one's answer only after that pass is over. So the spots at
     // the measured centre are found, and both schemes fed and indexed, here beside the file centre's

rugnux/RugnuxFirstPassRescues.cpp

@@ -795,166 +795,6 @@ auto Rugnux::FirstPassRun::AddEndWedge(FirstPass bp, IndexerThreadPool &pool, In
     return ScoreSchemes(ris, idx, &end_wedge, std::move(bp));
 }
 
-// Short-axis hypothesis - the small-molecule case.
-//
-// LAST, after the axis-sign and beam-centre rescues and the long-axis hypothesis, deliberately. Those all
-// run the standard floor and compare their result against `best` on volume ratios; a cell
-// adopted here is a k-fold sub-cell, so running before them makes every one of those
-// comparisons read an integer ratio and report an axis harmonic that is not one (measured:
-// the beam-centre check announced a "3.00x harmonic ... the beam centre is what settles
-// which" on a pair that differs only by the FLOOR). Here the comparison is like with like,
-// and this pass sees whatever geometry those rescues settled on.
-//
-// fft_min_unit_cell_A (10 A) is a protein constant: shorter than any protein axis and longer
-// than a small molecule's. It is applied to the FFT's PER-DIRECTION PEAK SEARCH
-// (FFTIndexerCPU.cpp / FFTIndexerGPU.cu skip every histogram bin below it) and only then to
-// the assembled cell (FFTIndexer::ReduceAndRefine, PostIndexingRefinement). The first of
-// those is what does the damage: the transform reports ONE peak per direction, so a
-// lattice row whose real-space repeat is shorter than the floor is reported at its first
-// HARMONIC above it, and the cell built out of those harmonics is an exact integer
-// sublattice of the true one - a real-space supercell. Measured on chemical-crystallography
-// data: index 2, 3 and 6 exactly, and the one dataset whose true reduced cell clears 10 A is
-// the one that comes out right. Because the floor is applied to the REDUCED cell it rejects
-// the lattice and not merely a basis choice, so nothing downstream can re-describe its way
-// back.
-//
-// Lowering the floor by default is not available: it only ever ADDS candidates, and the
-// Bravais classification the candidate pool feeds is under-determined, so it moves crystals
-// whose cell it never touched (measured: 2 regressions and 1 improvement over 20 protein
-// datasets, all three in the metric symmetry rather than the cell). So the low floor runs as
-// a SECOND hypothesis and the standard pass's answer stands unless the evidence takes it.
-//
-// What the evidence can be, before integration: a sublattice's per-spot residuals are
-// dominated pointwise by the finer lattice it contains, so every margin statistic - squares,
-// counts, trimmed sums - is biased toward the LARGER cell, and a spurious supercell indexes
-// every frame its true cell does. Only the SATURATION is informative. Hence the same shape
-// the scheme ratchet in PickBest already uses: an integer primitive-volume ratio plus a tie
-// on validation frames, or an outright win, which runs against the bias and so is real. The
-// statistic that genuinely decides an axis multiple is the occupancy of the reflection class
-// the multiple adds - AbsenceEvidence in SearchSpaceGroup - and that needs integrated
-// intensities, so it cannot run here.
-// Not when a cell was given. RaiseFFTBoundForKnownCell has already lowered the floor for
-// the standard pass against that cell, so this pass has nothing left to add - and it takes
-// its settings from rx.experiment_, which never carried the RAISED max bound, so on a cell
-// longer than the default search it could only ever return a harmonic of what was asked
-// for, tie on frames the way a sub-multiple always does, and be adopted over it.
-auto Rugnux::FirstPassRun::ShortAxisHypothesis(FirstPass best) -> FirstPass {
-    if (!rx.cancelled_ && !rx.experiment_.GetUnitCell().has_value()) {
-        std::unique_ptr<IndexerThreadPool> short_pool;
-        std::vector<std::unique_ptr<RotationIndexer>> sris;
-        if (short_spec && short_spec->key == rx.ExperimentKey() && short_spec->mapping == spot_mapping
-            && short_spec->spots_per_image == first_pass_spots_per_image) {
-            for (auto &f : short_spec->futures)
-                f.get();
-            short_pool = std::move(short_spec->pool);
-            sris = std::move(short_spec->sris);
-        } else {
-            short_spec.reset();
-            short_pool = std::make_unique<IndexerThreadPool>(ShortAxisSettings(rx.experiment_),
-                                                             IndexerConstruction::OnFirstUse, rx.gpu_work_beside_);
-            sris = FeedSchemes(*short_pool);
-            RunIndexing(sris);
-        }
-
-        for (size_t i = 0; i < sris.size() && !rx.cancelled_; i++) {
-            auto found = sris[i]->GetLattice();
-            if (!found.has_value())
-                continue;
-            RotationIndexerResult alt = std::move(*found);
-            const std::string &name = schemes[i].first;
-            // Primitive volumes, for the reason PickBest gives: two searches can find the same
-            // lattice in different settings, and a centred setting's cell is an exact integer
-            // multiple of its primitive one.
-            const double vol = std::abs(
-                    alt.lattice.ToPrimitive(alt.search_result.centering).CalcVolume());
-            if (!(vol > 1.0) || !(best.vol > 1.0))
-                continue;
-            const double ratio = best.vol / vol;
-            const double nearest = std::round(ratio);
-            // The gate is here rather than after the scoring because scoring is the expensive
-            // half - CountIndexed runs the real per-image path over the validation set, seconds
-            // on a protein - while the volume ratio is free. On every protein measured the low
-            // floor returns the SAME lattice (ratio 1.00-1.01), so this costs them one extra FFT
-            // and nothing else.
-            if (!(vol < best.vol) || nearest < 2.0 || std::abs(ratio - nearest) >= 0.15) {
-                logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': primitive volume "
-                            "{:.0f} A^3 against the standing cell's {:.0f} A^3, ratio {:.2f} - "
-                            "not an integer sub-cell, keeping the standing cell",
-                            IndexingSettings::fft_min_unit_cell_limit_A, name, vol, best.vol,
-                            ratio);
-                continue;
-            }
-
-            // The premise of this pass, checked rather than asserted. The pass exists because a
-            // lattice row shorter than the standard floor is reported at its first HARMONIC
-            // above it, so a cell it may rightfully adopt has a reduced axis BELOW that floor -
-            // and the floor is applied to the reduced cell, which is the one to measure. Where
-            // the alternative has no such axis the floor cannot be what the standard pass
-            // tripped over, and an integer volume ratio is then an ordinary axis multiple: the
-            // question the axis-harmonic arbiter settles just above on the occupancy of the
-            // class the multiple adds. That is intensity-class evidence about the very
-            // reflections in dispute, where the tie below is a frame count this comment already
-            // says a sublattice cannot lose - so where the two disagree the arbiter is right and
-            // this block must not overrule it. Measured on the corpus: the four small-molecule
-            // sets the pass is for all have a reduced axis under 9 A, and the two protein sets
-            // it halved have 35 and 47 A.
-            const auto alt_reduced = alt.lattice.ToPrimitive(alt.search_result.centering)
-                                                .NiggliReduce().GetUnitCell();
-            const float alt_shortest = std::min({alt_reduced.a, alt_reduced.b, alt_reduced.c});
-            const float standard_floor = rx.experiment_.GetIndexingSettings().GetFFT_MinUnitCell_A();
-            if (alt_shortest >= standard_floor) {
-                logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': a {:.0f}x sub-cell "
-                            "(vol {:.0f} A^3), but its shortest reduced axis is {:.1f} A, above "
-                            "the {:.1f} A floor the standard pass already searched - the floor "
-                            "is not why that cell was missed, so this is an ordinary axis "
-                            "multiple and not this pass's to decide; keeping the standing cell",
-                            IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, vol,
-                            alt_shortest, standard_floor);
-                continue;
-            }
-
-            const int alt_score = CountIndexed(*indexer, alt);
-            NearTieFrames("short-axis hypothesis", "sub-cell", alt_score, "standing", best.score);
-            // A pass that indexes nothing is not tied with a pass that indexes nothing. Without
-            // this, 0 against 0 passes the tie band and the volume ratio alone decides, which is
-            // the one thing this block must never do.
-            if (alt_score <= 0 || best.score <= 0) {
-                logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': {}/{} validation "
-                            "frames against {} - nothing indexes, keeping the standing cell",
-                            IndexingSettings::fft_min_unit_cell_limit_A, name, alt_score,
-                            static_cast<int>(validation.size()), best.score);
-                continue;
-            }
-            const bool tied = static_cast<float>(alt_score) >= best.score * 0.9f - 0.5f;
-            const bool clearly_more = static_cast<float>(alt_score) > best.score * 1.1f + 0.5f;
-            if (!tied && !clearly_more) {
-                logger.Info("Short-axis pass (floor {:.0f} A), scheme '{}': {:.0f}x sub-cell "
-                            "(vol {:.0f} A^3) indexes {}/{} validation frames against {} - "
-                            "keeping the standing cell",
-                            IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, vol,
-                            alt_score, static_cast<int>(validation.size()), best.score);
-                continue;
-            }
-            const auto uc = alt.lattice.GetUnitCell();
-            logger.Warning("Short-axis pass (floor {:.0f} A), scheme '{}': the standing "
-                           "cell is a {:.0f}x supercell. Its {:.0f} A^3 primitive volume against "
-                           "this pass's {:.0f} A^3, at {}/{} validation frames against {} - "
-                           "adopting {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}. The shortest axis "
-                           "is below the {:.0f} A default FFT floor, which is why the standard "
-                           "pass could only build a multiple of it",
-                           IndexingSettings::fft_min_unit_cell_limit_A, name, nearest, best.vol,
-                           vol, alt_score, static_cast<int>(validation.size()), best.score,
-                           uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma,
-                           rx.experiment_.GetIndexingSettings().GetFFT_MinUnitCell_A());
-            best.score = alt_score;
-            best.vol = vol;
-            best.name = name + " (short-axis)";
-            best.result = std::move(alt);
-        }
-    }
-    return best;
-}
-
 // Sub-lattice hypothesis: the standing cell's seven index-2 sub-lattices, in its own
 // orientation and geometry, each scored on the validation frames.
 //

rugnux/RugnuxLive.cpp

@@ -109,11 +109,9 @@ RugnuxLiveSetup MakeRugnuxStillsPrePassSetup(const JFJochReaderDataset &dataset,
 void WarmUpRugnuxLiveGPU(const RugnuxLiveSetup &setup) {
 #ifdef JFJOCH_USE_CUDA
     // As Rugnux sets them for its first-pass pool (the measured centre's indexing and pass 2 use it as
-    // well) and its short-axis pool. Should they part, the run builds its own, as without this.
+    // well). Should they part, the run builds its own, as without this.
     IndexingSettings first_pass = setup.experiment.GetIndexingSettings();
     first_pass.RefineThreads(std::clamp(setup.config.nthreads / 2, 1, MAX_REFINE_THREADS));
-    IndexingSettings short_axis = first_pass;
-    short_axis.FFT_MinUnitCell_A(IndexingSettings::fft_min_unit_cell_limit_A);
 
     // Each pool has two workers, on consecutive cards.
     const int gpus = get_gpu_count();
@@ -123,7 +121,6 @@ void WarmUpRugnuxLiveGPU(const RugnuxLiveSetup &setup) {
         cards.push_back(std::async(std::launch::async, [&, i] {
             ScopedGpuPin pin(i);
             KeepFFTIndexers(first_pass, per_card);
-            KeepFFTIndexers(short_axis, per_card);
             // Sizing the capture's work area makes its two plans (CudaFFTPlan::WorkBytes).
             BeamCenterFFTGPU::DeviceMemoryNeeded(setup.experiment.GetXPixelsNumConv(),
                                                  setup.experiment.GetYPixelsNumConv());

rugnux/RugnuxLive.h

@@ -26,7 +26,7 @@ RugnuxLiveSetup MakeRugnuxLiveSetup(const JFJochReaderDataset &dataset, int nthr
 RugnuxLiveSetup MakeRugnuxStillsPrePassSetup(const JFJochReaderDataset &dataset, int nthreads);
 
 // Builds on every card, at the same time, the GPU FFT indexers a live run of this setup is known to
-// need - its first pass's and its short-axis pass's (the settings Rugnux gives them) - and the cuFFT
+// need - its first pass's (the settings Rugnux gives it) - and the cuFFT
 // plans of its beam-centre capture, and with them the card's CUDA context and cuFFT's one-time
 // initialisation. Meant for the time between the start message and the first image, while the cards
 // are idle, under an FFTIndexerReuse and an FFTPlanReuse that outlive it: they are kept for the

rugnux/RugnuxPasses.cpp

@@ -405,9 +405,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
         // again at the other centre and keep whichever arm MERGES better, each judged on its own search
         // merge: P1, whole range, before any correction surface, which is the merge the two passes are
         // already compared in at the end of this function. Because the alternative arm is a whole first
-        // pass it gets the short-axis pass and the geometry post-refinement of its own, which the
-        // question needs - on a crystal with an axis under the FFT floor the check runs before that
-        // pass and compares supercells.
+        // pass it gets the geometry post-refinement of its own, which the question needs.
         //
         // Adopting a different metric moves the run off the geometry it was given, so it has to be
         // decisively better rather than merely ahead.
@@ -537,9 +535,8 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
             } else
                 alt = run_arm(false);
             // The question is asked again of what the two arms ENDED on. The check opened it on the
-            // lattices the standard-floor search returned, before the short-axis pass and the geometry
-            // post-refinement; either can replace that lattice, and on a crystal with an axis under the
-            // FFT floor the check compared two supercells the short-axis pass then threw away. Where
+            // lattices the first indexing returned, before the rescues and the geometry post-refinement;
+            // either can replace that lattice. Where
             // both arms now hold the same lattice in the same class, the centre did not decide the
             // symmetry after all: the two merges are two integrations of ONE hypothesis, their CC1/2
             // difference is the noise between them, and it is not asked. Same lattice means what it

tests/IndexingUnitTest.cpp

@@ -350,8 +350,9 @@ TEST_CASE("FFTIndexer_CentredReferenceCell","[Indexing]") {
     const auto from_primitive = index_with(primitive_uc);
     const auto from_conventional = index_with(conventional);
 
-    REQUIRE(from_primitive.size() == 1);
-    REQUIRE(from_conventional.size() == 1);
+    // The first cell is the one asked about; an integer multiple of it may follow (Refine keeps those).
+    REQUIRE(!from_primitive.empty());
+    REQUIRE(!from_conventional.empty());
 
     for (const auto &found : {from_primitive[0], from_conventional[0]}) {
         const auto uc_out = found.GetUnitCell();
@@ -747,6 +748,91 @@ namespace {
     }
 }
 
+TEST_CASE("FFTIndexer_ShortRowsBelowTheFloor","[Indexing]") {
+    // Two axes below the 10 A floor, one of them below 5 A. Rotation indexing finds the rows under the
+    // multiples the peak search reports past the floor and offers their cell; without it every cell
+    // offered has its axes past the floor, so it can only be a supercell.
+    const UnitCell uc(4.7, 7.8, 21.0, 90, 90, 90);
+    const CrystalLattice cl(uc);
+
+    // With a little positional scatter, as measured spots have: on an exact comb every harmonic of a row
+    // is as strong as the row itself, and nothing then tells the two apart.
+    std::mt19937 rng(7);
+    std::normal_distribution<float> scatter(0.0f, 0.005f);
+    std::vector<Coord> vec;
+    for (int h = -8; h <= 8; h++)
+        for (int k = -14; k <= 14; k++)
+            for (int l = -30; l <= 30; l++) {
+                const Coord s = h * cl.Astar() + k * cl.Bstar() + l * cl.Cstar();
+                if ((h != 0 || k != 0 || l != 0) && s.Length() < 1.0f / 0.8f)
+                    vec.push_back(s + Coord(scatter(rng), scatter(rng), scatter(rng)));
+            }
+
+    for (const auto algorithm : FFTAlgorithms()) {
+        for (const bool rotation : {true, false}) {
+            DiffractionExperiment experiment;
+            IndexingSettings settings;
+            settings.Algorithm(algorithm).RotationIndexing(rotation);
+            experiment.ImportIndexingSettings(settings);
+            experiment.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(1, 0, 0), std::nullopt));
+            std::unique_ptr<Indexer> indexer = CreateIndexer(experiment);
+            REQUIRE(indexer);
+
+            const auto result = indexer->Run(vec);
+            REQUIRE(!result.lattice.empty());
+            bool true_cell = false;
+            float shortest = 1e9f;
+            for (const auto &l : result.lattice) {
+                const auto u = l.GetUnitCell();
+                std::array<float, 3> out = {u.a, u.b, u.c};
+                std::sort(out.begin(), out.end());
+                shortest = std::min(shortest, out[0]);
+                true_cell |= std::abs(out[0] - 4.7f) < 0.05f && std::abs(out[1] - 7.8f) < 0.09f
+                             && std::abs(out[2] - 21.0f) < 0.2f;
+            }
+            if (rotation)
+                CHECK(true_cell);
+            else
+                CHECK(shortest >= 10.0f);
+        }
+    }
+}
+
+TEST_CASE("PostIndexingRefinement_IntegerMultipleSurvivesTheOverlapTest","[Indexing]") {
+    // A cell and its supercell index the same spots; the overlap test must not delete either.
+    const UnitCell uc(30.0, 40.0, 50.0, 90, 90, 90);
+    const CrystalLattice cl(uc);
+
+    std::vector<Coord> spots;
+    for (int h = -6; h <= 6; h++)
+        for (int k = -8; k <= 8; k++)
+            for (int l = -10; l <= 10; l++)
+                if (h != 0 || k != 0 || l != 0)
+                    spots.push_back(h * cl.Astar() + k * cl.Bstar() + l * cl.Cstar());
+
+    Eigen::MatrixX3<float> cells(6, 3);
+    const Coord a = cl.Vec0(), b = cl.Vec1(), c = cl.Vec2(), a2 = a * 2.0f;
+    cells << a.x, a.y, a.z, b.x, b.y, b.z, c.x, c.y, c.z,
+             a2.x, a2.y, a2.z, b.x, b.y, b.z, c.x, c.y, c.z;
+    Eigen::VectorX<float> scores(2);
+    scores << 0.2f, 0.2f;
+
+    RefineParameters p{
+        .viable_cell_min_spots = 6,
+        .dist_tolerance_vs_reference = 0.05f,
+        .min_length_A = 10.0f,
+        .max_length_A = 500.0f,
+        .min_angle_deg = 30.0f,
+        .max_angle_deg = 150.0f,
+        .indexing_tolerance = 0.1f
+    };
+    const auto out = Refine(spots, spots.size(), cells, scores, p);
+    REQUIRE(out.size() == 2);
+    // The two index alike, so the smaller cell ranks first.
+    CHECK(std::abs(out[0].CalcVolume()) == Catch::Approx(cl.CalcVolume()).epsilon(0.01));
+    CHECK(std::abs(out[1].CalcVolume()) == Catch::Approx(2.0 * cl.CalcVolume()).epsilon(0.01));
+}
+
 TEST_CASE("FFTIndexer_NoiseFrame", "[Indexing]") {
     // The frame that works the indexer hardest is not a clean crystal but junk. A real lattice
     // prunes the FFT shortlist to four or five distinct directions; noise leaves dozens, which is
@@ -797,7 +883,7 @@ TEST_CASE("FFTIndexer_CoplanarNoiseFrameLeavesTheIndexerUsable", "[Indexing]") {
             CHECK(IndexedFraction(noise, lattice, 0.1f) < 0.05f);
 
         const auto result = indexer->Run(lattice_spots);
-        REQUIRE(result.lattice.size() == 1);
+        REQUIRE(!result.lattice.empty());
         const auto found = result.lattice[0].GetUnitCell();
         std::array<float, 3> lengths = {static_cast<float>(found.a), static_cast<float>(found.b),
                                         static_cast<float>(found.c)};