Jungfraujoch/rugnux/Rugnux.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "Rugnux.h"
#include "ModelValidation.h"
#include "SpotWidth.h"

#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <condition_variable>
#include <deque>
#include <functional>
#include <future>
#include <mutex>
#include <numeric>
#include <set>
#include <sstream>
#include <thread>

#include "../reader/JFJochHDF5Reader.h"
#include "../common/JFJochMath.h"
#include "../common/ParallelFor.h"
#include "../common/Logger.h"
#include "../common/AzimuthalIntegrationMapping.h"
#include "../common/AzimuthalIntegrationProfile.h"
#include "../common/CUDAWrapper.h"
#include "../common/JFJochException.h"
#include "../common/time_utc.h"
#include "../writer/FileWriter.h"
#include "../image_analysis/MXAnalysisWithoutFPGA.h"
#include "../image_analysis/beam_stop/ShadowFinder.h"
#include "../image_analysis/IndexAndRefine.h"
#include "../image_analysis/indexing/AnalyzeIndexing.h"
#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h"
#include "../image_analysis/geom_refinement/BeamCenterFromSpots.h"
#include "../image_analysis/geom_refinement/GeometryRefiner.h"
#include "../image_analysis/indexing/IndexerThreadPool.h"
#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"
#include "gemmi/twin.hpp"
#include "../image_analysis/azint/AzIntEngineCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h"
#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/RfreeFlags.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
#include "../image_analysis/scale_merge/ScalingResult.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/geom_refinement/PostRefine.h"
#include "../image_analysis/lattice_search/LatticeSearch.h"
#include "../image_analysis/scale_merge/AnisotropyAnalysis.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../image_analysis/scale_merge/HKLKey.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/bragg_integration/CalcISigma.h"
#include "../common/Definitions.h"
#include "../common/CorrelationCoefficient.h"
#include <array>
#include <map>

namespace {
    // A CUDA out-of-memory, or a host allocation that failed, says nothing about the image that
    // happened to be in flight and everything about the machine - the next image will hit it too.
    // Skipping it the way a genuinely undecodable image is skipped drops frames out of the middle of
    // the dataset and still exits 0, so the run looks complete while every merged number downstream
    // has quietly changed. Fail the run instead, and let the operator free the card and repeat it.
    bool IsFatalResourceError(const std::exception &e) {
        if (dynamic_cast<const std::bad_alloc *>(&e)) return true;
        const auto *jf = dynamic_cast<const JFJochException *>(&e);
        return jf != nullptr && (jf->Category() == JFJochExceptionCategory::GPUCUDAError
                                 || jf->Category() == JFJochExceptionCategory::MemAllocFailed);
    }

    // How precisely the beam centre has to be placed before the estimate is used at all. An estimate
    // that fails it is treated as no estimate: the spot symmetry falls through to the background and
    // the background to the header. A fit's own scatter will not serve as the gate - the ways either
    // estimator goes badly wrong all leave a small one - so the spot symmetry reports instead how
    // far its answer moves when the search is started elsewhere, which does separate them. On the 38
    // rotation regression crystals that stays below 0.35 px against errors of at most 0.74; over a
    // ring of injected geometries it rises past 1.0 px on 96 % of the answers that are 2 px or more
    // out, and on none of the 38 as they are.
    constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f;

    // The longest cell axis assumed where no cell is known yet. It has to be LARGE: a long axis makes
    // the tolerance below small, and the ceiling is the larger of it and the flat bound above, so a
    // large default is the one that changes nothing.
    constexpr float UNKNOWN_CELL_A_MAX_A = 200.0f;

    // How far the beam centre may be wrong before the FIRST PASS stops indexing. A centre error is
    // fixed in the LAB frame, so accumulating a sweep smears every reciprocal-lattice point around a
    // circle and the FFT amplitude at an axis of length a is multiplied by J0(2 pi delta p a/(D
    // lambda)); this is the displacement at which that has fallen to 0.70, and the first zero is
    // 2.09x further out. Over the in-house and non-SLS corpora it runs from 0.6 px to 9 px - a
    // fifteen-fold spread, which is why the flat bound above is only ever used as a FLOOR under it.
    //
    // It is printed, and it is deliberately not a relevance test. It says how wrong the header may
    // be; the estimator's sigma says how well the estimator knows its own answer, and gating one on
    // the other rejects a centre correct to 0.03 px because its sigma was 0.98. Nor is a move that is
    // small against it a move that can be ignored: measured on real data, a 0.12 px change of centre
    // - 0.03x of this - is the difference between the deposited cell and a halved axis.
    float BeamCenterNeed_pxl(const DiffractionExperiment &experiment) {
        const auto cell = experiment.GetUnitCell();
        const float a_max = cell ? std::max({cell->a, cell->b, cell->c}) : UNKNOWN_CELL_A_MAX_A;
        return 0.183f * experiment.GetDetectorDistance_mm() * experiment.GetWavelength_A()
               / (experiment.GetPixelSize_mm() * a_max);
    }

    // A detector-plane displacement resolved ACROSS the spindle - the only component the tolerance
    // above is about. Along the spindle the error translates the derotated cloud rigidly and the FFT
    // amplitude cannot see it at all; across it the peaks smear and the amplitude falls. (The lattice
    // fit that follows the FFT is not translation-invariant, so "along the spindle" is free for the
    // transform and not for the run - which is why this is reported and not used to decide anything.)
    // The detector direction wanted is the one whose lab image is perpendicular to the spindle n:
    // (ux e1 + uy e2).n = 0, so u is proportional to (-(e2.n), e1.n).
    std::optional<float> AcrossSpindle_pxl(const DiffractionExperiment &experiment, float dx, float dy) {
        const auto goniometer = experiment.GetGoniometer();
        if (!goniometer || !goniometer->IsScanning())
            return {};
        const auto geom = experiment.GetDiffractionGeometry();
        const float beam_x = geom.GetBeamX_pxl(), beam_y = geom.GetBeamY_pxl();
        const Coord n = goniometer->GetAxis().Normalize();
        const Coord origin = geom.LabCoord(beam_x, beam_y);
        const Coord e1 = geom.LabCoord(beam_x + 1.0f, beam_y) - origin;
        const Coord e2 = geom.LabCoord(beam_x, beam_y + 1.0f) - origin;
        const float ux = -(e2 * n), uy = e1 * n;
        const float length = std::hypot(ux, uy);
        if (length < 1e-3f)
            return {};
        return std::abs(dx * ux + dy * uy) / length;
    }

    // Images the projection is built from when the beam centre is wanted but the beam-stop pre-pass
    // is off; with it on, that pre-pass's image count is used and one projection serves both.
    constexpr int BEAM_CENTER_PROJECTION_IMAGES = 60;

    // Spots per pre-scan frame kept for the beam centre, the strongest first.
    constexpr size_t BEAM_CENTER_SPOTS_PER_IMAGE = 500;

    // Frames read for the spot symmetry alone when the first sample could not measure it. Both
    // estimators are counting statistics over PAIRS of frames - the Friedel one over frames half a
    // turn apart, the second crossing over any two frames a reflection is recorded on - so a sweep
    // whose frames carry few spots can only be measured by reading more of them. On the two
    // rotation regression crystals that decline, this is what turns the decline into an answer
    // 0.04 px from the post-refined centre; the sample the beam stop is built from is left alone.
    constexpr int BEAM_CENTER_SPARSE_IMAGES = 400;

    // Frames the spot symmetry reads, as a multiple of the beam stop's. Its two estimators are
    // counts of matched PAIRS and the pairs grow faster than the frames do, so this is the one knob
    // that buys independence from WHICH frames were drawn - and that dependence is the estimator's
    // largest remaining failure. Asking the same code for one frame more or less re-draws the whole
    // sample: over five draws of 38 sweeps, 60 frames leaves three answers more than 1 px out, none
    // of which reports an unusual sigma, and 120 leaves none in 186. It is also what makes the end
    // margin below free: at 60 frames excluding the ends moved three centres, at 120 it moves none.
    constexpr int BEAM_CENTER_IMAGE_FACTOR = 2;

    // Images left out at each end of the sweep. Shutter synchronisation disturbs the first and last
    // frame or two of a run, and a partly exposed frame is dark over the whole detector - which is
    // exactly the kind of error a projection cannot average away and a background comparison reads
    // as signal. The margin is capped at a tenth of the sweep so a short run still has a sample.
    constexpr int PRESCAN_END_MARGIN_IMAGES = 5;

    // Workers reading the pre-scan sample. Each owns a shard of the beam-stop projection so no two
    // threads touch the same accumulator, and a shard costs 20 bytes per pixel - 362 MB on a 16M
    // detector - so this is capped well below the worker count of the run proper. The accumulation
    // is memory-bound rather than compute-bound, so a handful of workers already saturates it.
    constexpr size_t PRESCAN_MAX_WORKERS = 8;

    // The spot width is measured on a GROWING share of the pre-scan sample: every eighth frame of
    // it, then every fourth, then every second, then all of it - an eighth, a quarter, a half and
    // the whole. Each stride is a multiple of the next, so a tier's sample really does contain the
    // one before it, and the measurement stops at the first tier whose answer has settled
    // (spot_width::WidthSettled, in the loop below). A crystal whose frames carry plenty of isolated
    // strong spots is characterised from a quarter of the sample; a sparse one still gets all of it.
    //
    // The tiers add PRESCAN_MAX_WORKERS-sized steps as well: on a 60 frame sample they are 8, 15, 30
    // and 60 frames, i.e. 8 + 7 + 15 + 30 read in the four passes, so no pass leaves most of the
    // workers idle waiting for the others to finish theirs. A ladder that starts smaller measures
    // fewer frames but costs more wall clock, because a pass of five frames occupies eight workers
    // for as long as a pass of eight does.
    //
    // Only the width work is tiered. Every frame of the sample is still READ once - the frames are
    // simply visited in tier order - and the beam-stop projection still sees all of them, so the
    // mask and the beam centre are what they were. What the tiers cut is the per-frame work the
    // width needs on top of the read: the host decompression, the preprocessing and the spot
    // finding, which together are two thirds of the pre-scan's cycles once the projection itself is
    // decoding on the GPU. The width measurement proper is under 1 % of them.
    //
    // Measured over the rotation battery: the radius chosen is the whole-sample answer on 38 of 38
    // crystals, from an average of 32 frames of the 60 instead of all of them.
    constexpr std::array<int, 4> WIDTH_TIER_STRIDE = {8, 4, 2, 1};

    // Pick up to requested_images ordinals spread evenly across [0, images_to_process) for the
    // first pass of two-pass rotation indexing.
    std::vector<int> select_equally_spaced_image_ordinals(int images_to_process, int requested_images) {
        std::vector<int> ret;
        if (images_to_process <= 0 || requested_images <= 0)
            return ret;

        const int n = std::min(images_to_process, requested_images);
        if (n == 1) {
            ret.push_back(0);
            return ret;
        }

        std::set<int> unique_ordinals;
        for (int i = 0; i < n; i++)
            unique_ordinals.insert(static_cast<int>(
                std::llround(static_cast<double>(i) * static_cast<double>(images_to_process - 1) /
                             static_cast<double>(n - 1))));

        ret.assign(unique_ordinals.begin(), unique_ordinals.end());
        return ret;
    }

    // A sample arranged as pairs of images half a turn apart: half the ordinals spread over the
    // part of the sweep that has a partner 180 deg later, and those partners. That is what the
    // Friedel mirror needs. It is NOT the plain selection rearranged - even on a full turn the two
    // sets are all but disjoint, sharing 2 of 60 ordinals on a 3600 image sweep - so the caller
    // keeps them apart and gives each consumer the one it needs. Empty when the sweep is shorter
    // than half a turn, which is the caller's signal that the spot symmetry cannot be measured here.
    std::vector<int> select_half_turn_paired_ordinals(int images_to_process, int requested_images,
                                                      float degrees_per_ordinal) {
        if (!(degrees_per_ordinal > 0.0f))
            return {};
        const int half_turn = static_cast<int>(std::lround(180.0 / degrees_per_ordinal));
        if (half_turn < 1 || images_to_process - half_turn < 1)
            return {};

        std::set<int> ordinals;
        for (const int ordinal : select_equally_spaced_image_ordinals(images_to_process - half_turn,
                                                                      requested_images / 2)) {
            ordinals.insert(ordinal);
            ordinals.insert(ordinal + half_turn);
        }
        return {ordinals.begin(), ordinals.end()};
    }

    // Apply a goniometer rotation SCALE: the stage turned k times the angle the file records. What is
    // wrong is the SWEEP, not where it began, so the per-frame increment and the per-frame oscillation
    // width both scale by k while the starting angle is left alone - the stage reached that position
    // before the sweep started. A constant offset in phi is in any case exactly degenerate with a
    // rotation of the crystal orientation about the spindle, which indexing refines away, so anchoring
    // the stretched sweep at image 0 costs nothing and keeps the first frame at the angle the file gives
    // it. Every angle in the pipeline comes from this object, so scaling it is the whole correction.
    GoniometerAxis ScaleRotation(const GoniometerAxis &g, float k) {
        GoniometerAxis scaled(g.GetName(), g.GetStart_deg(), g.GetIncrement_deg() * k,
                              g.GetAxis(), g.GetHelicalStep());
        if (const auto wedge = g.GetScreeningWedge())
            scaled.ScreeningWedge(*wedge * k);
        return scaled;
    }

    // Writing an image to the process file takes the global HDF5 mutex, which is the same mutex every
    // worker thread needs to fetch its next image. The write itself is short, but with a worker per
    // hardware thread they spend longer queueing for the lock than the writing takes: on a 16 Mpx
    // rotation dataset the per-image loop runs 7.7 s with 8 workers and 13.1 s with 48. Give the file
    // to one thread and let the workers post to it, and the lock has one taker again.
    //
    // A DataMessage does not own its pixels - it points into the reader's buffer - so the raw image is
    // parked in the queue beside it. Without that the worker frees the pixels on its next iteration and
    // the writer reads whatever landed there.
    class ProcessFileWriter {
    public:
        ProcessFileWriter(FileWriter &writer, size_t capacity)
            : writer_(writer), capacity_(capacity),
              thread_([this] { Run(); }) {}

        // Post one image. Blocks while the queue is full, which is what stops a run whose analysis
        // outpaces its writer from holding every image it has ever processed in memory.
        void Post(const DataMessage &msg, std::shared_ptr<JFJochReaderRawImage> img) {
            std::unique_lock lock(m_);
            space_.wait(lock, [this] { return queue_.size() < capacity_ || failed_; });
            if (failed_)
                return;   // the error is rethrown by Finish(); dropping the rest is deliberate
            queue_.push_back({msg, std::move(img)});
            work_.notify_one();
        }

        // Drain, join, and rethrow whatever the writer thread hit. Must be called before the
        // FileWriter is used for anything else - the end message, Finalize().
        void Finish() {
            {
                std::lock_guard lock(m_);
                done_ = true;
            }
            work_.notify_one();
            if (thread_.joinable())
                thread_.join();
            if (error_)
                std::rethrow_exception(error_);
        }

        ~ProcessFileWriter() {
            {
                std::lock_guard lock(m_);
                done_ = true;
                failed_ = true;   // an unwinding run should not wait for the backlog to be written
            }
            work_.notify_one();
            space_.notify_all();
            if (thread_.joinable())
                thread_.join();
        }

    private:
        struct Job {
            DataMessage msg;
            std::shared_ptr<JFJochReaderRawImage> img;
        };

        void Run() {
            while (true) {
                Job job;
                {
                    std::unique_lock lock(m_);
                    work_.wait(lock, [this] { return !queue_.empty() || done_; });
                    if (queue_.empty())
                        return;
                    job = std::move(queue_.front());
                    queue_.pop_front();
                }
                space_.notify_one();
                try {
                    writer_.Write(job.msg);
                } catch (...) {
                    std::lock_guard lock(m_);
                    if (!error_)
                        error_ = std::current_exception();
                    failed_ = true;
                    space_.notify_all();
                    return;
                }
            }
        }

        FileWriter &writer_;
        const size_t capacity_;
        std::mutex m_;
        std::condition_variable work_;
        std::condition_variable space_;
        std::deque<Job> queue_;
        bool done_ = false;
        bool failed_ = false;
        std::exception_ptr error_;
        std::thread thread_;
    };

}

// What the file says about the beam centre, before a frame is read. Nothing here measures anything -
// it reads the value that is there and says where it looks like it came from - and it is free, which
// is why it runs on every run whether or not a centre was asked for.
//
// A centre ON A MASKED PIXEL is the one of these worth acting on: a real beam does not sit on a dead
// pixel or in a module gap, so the value is wrong about something. Over the 29 non-SLS datasets whose
// outcome is on record (base rate of a bad outcome 48%) it fires on 4 and 3 of those 4 went wrong,
// against 2 false alarms in 42 in-house masters. What it does NOT promise is that correcting the
// centre repairs the run - on one of the three it does not.
//
// The other two are PROVENANCE, not prediction, and are said rather than warned about. A value equal
// to the geometric detector centre, or a whole number of pixels in both coordinates, was typed rather
// than measured - true of 29 of 42 in-house masters, and those are out by a median 3.5 px against
// 0.7 px for the rest - but a typed value is often close enough that as predictors of harm both sit
// at or below the base rate. They are still worth saying, because the spot-symmetry search reaches
// only about +-55 px from wherever it starts: a placeholder tens of pixels out is unreachable by
// construction, all five of its searches agree, and the small sigma they report says nothing is
// wrong. Only the background estimator recovers that case, and only a human reading this line knows
// to look.
// The lattice METRIC can carry more rotational symmetry than the group the intensities supported.
// That is not by itself a mistake - a pseudo-symmetric metric is ordinary, and only the intensities
// can say whether the extra rotations are real - but it is where the corpus's symmetry under-calls
// live, so the run says it out loud rather than leaving the user to notice. It decides nothing: no
// threshold, no promotion, no demotion.
//
// The metric symmetry is Le Page's, and its obliquity is the continuous measure of how nearly a
// two-fold holds; 1 degree is the middle of the band over which the answer is stable.
// Following Le Page (1982) J. Appl. Cryst. 15, 255-259
static void WarnIfMetricHostsMoreSymmetry(const std::optional<UnitCell> &cell,
                                          const gemmi::SpaceGroup &sg) {
    if (!cell)
        return;
    Logger logger("Rugnux");
    const gemmi::UnitCell gc(cell->a, cell->b, cell->c, cell->alpha, cell->beta, cell->gamma);
    const size_t metric_rotations =
        gemmi::find_lattice_symmetry(gc, sg.centring_type(), 1.0).sym_ops.size();
    const size_t group_rotations = sg.operations().sym_ops.size();
    if (metric_rotations <= group_rotations)
        return;
    // The proper rotations of each holohedry, which is what find_lattice_symmetry returns.
    const char *metric_system = "a higher";
    switch (metric_rotations) {
        case 2:  metric_system = "monoclinic";   break;
        case 4:  metric_system = "orthorhombic"; break;
        case 6:  metric_system = "trigonal";     break;
        case 8:  metric_system = "tetragonal";   break;
        case 12: metric_system = "hexagonal";    break;
        case 24: metric_system = "cubic";        break;
        default: break;
    }
    logger.Warning("The cell metric is {} - it admits {} rotations where {} has {}. The intensities "
                   "did not support the extra ones, so they are either absent or too weak to see; "
                   "this is a note, not a correction. Where the higher symmetry is real, merging in "
                   "it would raise multiplicity and completeness.",
                   metric_system, metric_rotations, sg.xhm(), group_rotations);
}

static void LogHeaderBeamCenter(const DiffractionExperiment &experiment, const PixelMask &pixel_mask) {
    Logger logger("Rugnux");
    const float beam_x = experiment.GetBeamX_pxl();
    const float beam_y = experiment.GetBeamY_pxl();
    const auto width = static_cast<int>(experiment.GetXPixelsNumConv());
    const auto height = static_cast<int>(experiment.GetYPixelsNumConv());

    if (std::hypot(beam_x - static_cast<float>(width) / 2.0f,
                   beam_y - static_cast<float>(height) / 2.0f) < 1.0f)
        logger.Info("Beam centre ({:.2f},{:.2f}) is the geometric centre of the detector - a value "
                    "written rather than measured", beam_x, beam_y);
    else if (beam_x == std::round(beam_x) && beam_y == std::round(beam_y))
        logger.Info("Beam centre ({:.0f},{:.0f}) is a whole number of pixels in both coordinates - a "
                    "value written rather than measured", beam_x, beam_y);

    const int x = static_cast<int>(std::lround(beam_x));
    const int y = static_cast<int>(std::lround(beam_y));
    if (x >= 0 && y >= 0 && x < width && y < height
        && pixel_mask.GetMask()[static_cast<size_t>(y) * width + x] != 0)
        logger.Warning("Beam centre ({:.2f},{:.2f}) lands on a masked pixel. A real beam does not, so "
                       "the value in the file is suspect; --estimate-beam-center measures it",
                       beam_x, beam_y);
}

Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment,
                             PixelMask pixel_mask, ProcessConfig config)
    : reader_(reader), experiment_(std::move(experiment)),
      pixel_mask_(std::move(pixel_mask)), config_(std::move(config)),
      user_fixed_sg_(experiment_.GetGemmiSpaceGroup()) {
    // Bit 9 describes where THIS run found the beam stop, so a mask read back from a file that
    // already carries one starts clear; the user mask (bit 8) is left as it was loaded.
    pixel_mask_.ClearBeamStopMask(experiment_);

    LogHeaderBeamCenter(experiment_, pixel_mask_);

    // A manually asserted stage calibration applies to everything, before anything reads an angle.
    if (config_.rotation_scale.has_value())
        if (const auto g = experiment_.GetGoniometer()) {
            Logger("Rugnux").Info("Goniometer rotation scale {:.5f} applied from the command line",
                                  *config_.rotation_scale);
            experiment_.Goniometer(ScaleRotation(*g, *config_.rotation_scale));
        }
}

// Report the detector geometry the way XDS states it, so a run can be compared line by line against
// the IDXREF.LP / CORRECT.LP of XDS on the same data. XDS is never given this geometry - the durin
// plugin hands it image data only (plugin_get_header returns dimensions, bytes per pixel, pixel size
// and frame count, nothing else) and XDS refines its own from XDS.INP - which is precisely why
// having ours in the same convention is worth printing.
//
// The two laboratory frames coincide: x along increasing detector column, y along increasing row,
// z along the beam. XDS places a pixel at
//     x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis)
// and DiffractionGeometry::LabCoord is that same expression with X_axis = poni_rot*(1,0,0) and
// Y_axis = poni_rot*(0,1,0). So the tilt goes out as the two axis vectors rather than as angles,
// which is also how XDS reports it after refinement. Taking the axes as differences of LabCoord
// keeps them in step with whatever the geometry currently is, tilt included.
//
// ORGX/ORGY are 1-based - XDS counts pixels from 1 - and are the foot of the perpendicular from the
// crystal rather than the direct beam, which is what our beam centre is too
// (docs/DETECTOR_GEOMETRY.md), so they map across without correction.
//
// ROTATION_AXIS is printed as stored. Its direction is right, the frames being shared; its SIGN has
// not been cross-checked against an XDS refinement, so treat a flip there as unconfirmed rather than
// as a discrepancy in the data.
static void LogXDSGeometry(const DiffractionExperiment &experiment) {
    Logger logger("Rugnux");

    const auto geom = experiment.GetDiffractionGeometry();
    const float beam_x = geom.GetBeamX_pxl();
    const float beam_y = geom.GetBeamY_pxl();
    const float pixel_mm = geom.GetPixelSize_mm();

    const Coord origin = geom.LabCoord(beam_x, beam_y);
    const Coord x_axis = ((geom.LabCoord(beam_x + 1.0f, beam_y) - origin) / pixel_mm).Normalize();
    const Coord y_axis = ((geom.LabCoord(beam_x, beam_y + 1.0f) - origin) / pixel_mm).Normalize();

    logger.Info("XDS convention: ORGX= {:.2f} ORGY= {:.2f} DETECTOR_DISTANCE= {:.4f}",
                beam_x + 1.0f, beam_y + 1.0f, geom.GetDetectorDistance_mm());
    logger.Info("XDS convention: DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}",
                x_axis.x, x_axis.y, x_axis.z);
    logger.Info("XDS convention: DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}",
                y_axis.x, y_axis.y, y_axis.z);
    logger.Info("XDS convention: INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= {:.6f} QX= QY= {:.6f}",
                geom.GetWavelength_A(), pixel_mm);
    if (const auto goniometer = experiment.GetGoniometer())
        logger.Info("XDS convention: ROTATION_AXIS= {:.6f} {:.6f} {:.6f}",
                    goniometer->GetAxis().x, goniometer->GetAxis().y, goniometer->GetAxis().z);
}

// The FFT recovers a basis vector only up to fft_max_unit_cell_A: FFTIndexer sizes its projected
// histogram from that value and the transform's last usable bin IS that length. So a reference cell
// longer than the bound cannot be represented by the search meant to find it - the run then reports
// "found no lattice", which is strictly worse than not giving the cell at all (measured on a 545 A
// axis: a partial triclinic solution without -C, nothing with it). Raise the bound to cover the cell,
// with the margin DIALS uses when a cell is known (max_cell = 1.3 * longest axis, see
// dials/algorithms/indexing/indexer.py find_max_cell).
//
// The floor at the short end is adjusted here too, for the same reason and against the same cell -
// see the second half of the body.
//
// ROTATION ONLY. Stills fire the FFT once per image across every worker and routinely run with a
// known cell, so widening the transform there would cost the whole serial run for nothing; the
// caller applies this inside its rotation branch.
static void RaiseFFTBoundForKnownCell(IndexingSettings &settings, const DiffractionExperiment &x) {
    const auto cell = x.GetUnitCell();
    if (!cell.has_value() || !cell->is_finite())
        return;

    constexpr float KNOWN_CELL_MARGIN = 1.3f;
    const float longest = std::max({cell->a, cell->b, cell->c});
    // Clamp rather than let the setter throw. A cell that needs more than the search can ever reach
    // is a real possibility - the long-axis rescue calls this with whatever its coarse pass recovered,
    // and that has been seen to be an implausible harmonic - and a bad candidate must not take the
    // whole run down. Searching as far as the indexer goes is the best available answer; if the cell
    // really is that long the pass simply fails to index, which the caller already handles.
    const float needed = std::min(KNOWN_CELL_MARGIN * longest,
                                  IndexingSettings::fft_max_unit_cell_limit_A);
    if (needed > settings.GetFFT_MaxUnitCell_A()) {
        settings.FFT_MaxUnitCell_A(needed);
        Logger("Rugnux").Info("FFT search bound raised to {:.0f} A to cover the given cell "
                              "(longest axis {:.1f} A)", needed, longest);
    }

    // The same argument at the other end of the search. Every candidate whose reduced cell has an
    // axis below fft_min_unit_cell_A is discarded, and the 10 A default is shorter than any protein
    // axis but longer than a small molecule's: a given cell with an 8 A axis was thrown away before
    // it could be compared against the reference it was supposed to match. Lower the floor to reach
    // the cell that was asked for, with the reciprocal of the margin used above.
    constexpr float KNOWN_CELL_FLOOR_MARGIN = 1.0f / KNOWN_CELL_MARGIN;
    const float shortest = std::min({cell->a, cell->b, cell->c});
    const float floor_A = std::max(KNOWN_CELL_FLOOR_MARGIN * shortest,
                                   IndexingSettings::fft_min_unit_cell_limit_A);
    if (floor_A < settings.GetFFT_MinUnitCell_A()) {
        settings.FFT_MinUnitCell_A(floor_A);
        Logger("Rugnux").Info("FFT shortest-axis floor lowered to {:.1f} A to cover the given cell "
                              "(shortest axis {:.1f} A)", floor_A, shortest);
    }
}

std::optional<CalibrationResult> Rugnux::RebinAndRefit(const CalibrationResult &first,
                                                       int start_image, int end_image,
                                                       int images_to_process,
                                                       const AzimuthalIntegrationMapping &mapping,
                                                       const std::vector<SpotToSave> &spots) {
    Logger logger("Rugnux");
    const DiffractionGeometry before = experiment_.GetDiffractionGeometry();

    // How far a geometry moves a ring ON THE DETECTOR, against the radial width of one profile bin.
    // Below one bin, re-binning would put every pixel back where it already is. A beam-centre shift
    // moves a ring by its own size; a distance change moves it in proportion to the ring's radius.
    const auto [radius_min, radius_max] = ProfileRadiusRange_pxl(mapping, before);
    if (!(radius_max > radius_min))
        return std::nullopt;
    const float bin_pxl = (radius_max - radius_min) / static_cast<float>(mapping.GetQBinCount());
    const auto ring_gap_pxl = [&](const DiffractionGeometry &a, const DiffractionGeometry &b) {
        return std::hypot(a.GetBeamX_pxl() - b.GetBeamX_pxl(), a.GetBeamY_pxl() - b.GetBeamY_pxl())
                + radius_max * std::abs(a.GetDetectorDistance_mm() - b.GetDetectorDistance_mm())
                  / std::max(1e-3f, b.GetDetectorDistance_mm());
    };
    const auto ring_shift_pxl = [&](const DiffractionGeometry &g) { return ring_gap_pxl(g, before); };

    // Which geometries are worth re-binning about. The fitted one is the obvious candidate, but it was
    // itself fitted on the profile now believed mis-binned, so where the beam centre is badly out it is
    // wrong and re-binning about it only digs deeper - measured, rms 3.75 -> 5.66 on an exposure 200 px
    // off. The spots' geometry is the one that does not degrade there, having never read the header. So
    // try both where they differ and keep whichever comes back better; they cost a read each, and a
    // calibration run is a handful of images.
    std::vector<DiffractionGeometry> candidates;
    const auto add = [&](const DiffractionGeometry &g) {
        if (!(ring_shift_pxl(g) > bin_pxl))
            return;   // this geometry bins the images the way they are already binned
        for (const auto &c : candidates)
            if (ring_gap_pxl(g, c) < bin_pxl)
                return;   // and this one bins them the way another candidate already does
        candidates.push_back(g);
    };
    add(first.geometry);
    if (first.spots_geometry)
        add(*first.spots_geometry);

    if (candidates.empty()) {
        logger.Info("Powder calibration: the fitted geometry moves a ring by {:.2f} px, under the "
                    "{:.2f} px a profile bin spans - re-binning would change nothing",
                    ring_shift_pxl(first.geometry), bin_pxl);
        return std::nullopt;
    }

    std::optional<CalibrationResult> best;
    for (const auto &geometry : candidates) {
        if (cancelled_)
            break;
        logger.Info("Powder calibration: re-integrating binned about beam ({:.1f},{:.1f}) at {:.3f} mm, "
                    "which moves a ring {:.2f} px against a {:.2f} px bin",
                    geometry.GetBeamX_pxl(), geometry.GetBeamY_pxl(),
                    geometry.GetDetectorDistance_mm(), ring_shift_pxl(geometry), bin_pxl);

        // The mapping is the only thing the new geometry is needed for. Preprocessing does not depend on
        // where the beam is, so it keeps experiment_ and nothing here mutates the run's own geometry.
        DiffractionExperiment rebinned_experiment(experiment_);
        rebinned_experiment.BeamX_pxl(geometry.GetBeamX_pxl()).BeamY_pxl(geometry.GetBeamY_pxl())
                           .DetectorDistance_mm(geometry.GetDetectorDistance_mm())
                           .PoniRot1_rad(geometry.GetPoniRot1_rad())
                           .PoniRot2_rad(geometry.GetPoniRot2_rad());
        const AzimuthalIntegrationMapping rebinned_mapping(rebinned_experiment, pixel_mask_);

        std::atomic<int> next_ordinal = 0;
        std::mutex merge_m;
        AzimuthalIntegrationProfile summed(rebinned_mapping);

        const int nthreads = std::max(1, config_.nthreads);
        std::vector<std::future<void>> workers;
        workers.reserve(nthreads);
        for (int t = 0; t < nthreads; ++t) {
            workers.emplace_back(std::async(std::launch::async, [&]() {
                std::vector<uint8_t> decompression_buffer;
                ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_);
                ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum());
                AzIntEngineCPU azint(rebinned_mapping);
                AzimuthalIntegrationProfile local(rebinned_mapping);

                while (!cancelled_) {
                    const int ordinal = next_ordinal.fetch_add(1);
                    const int image_idx = start_image + ordinal * config_.stride;
                    if (image_idx >= end_image) break;
                    try {
                        auto img = reader_.GetRawImage(image_idx);
                        if (!img) continue;
                        const uint8_t *image_ptr = img->image.GetUncompressedPtr(decompression_buffer);
                        preprocessor.Analyze(buffer, image_ptr, img->image.GetMode());
                        azint.Run(buffer, local);
                    } catch (const std::exception &e) {
                        if (IsFatalResourceError(e)) throw;
                        logger.Error("Failed to re-integrate image {}: {}", image_idx, e.what());
                    }
                }
                std::lock_guard lock(merge_m);
                summed += local;
            }));
        }
        for (auto &w : workers)
            w.get();
        if (cancelled_)
            return std::nullopt;

        // The profile is measured about this geometry, so this is also what the fit starts from and what
        // its ring points are recovered against.
        try {
            auto refit = CalibrateFromProfile(summed.GetResult(), rebinned_mapping, geometry,
                                              config_.calibrant_ring_q,
                                              config_.calibration_refine_tilt, spots);
            logger.Info("Powder calibration: that gives {} ring points at rms {:.3f} px",
                        refit.ring_points, refit.rms_radial_pxl);
            if (!best || (refit.ring_points * 2 >= best->ring_points
                          && refit.rms_radial_pxl < best->rms_radial_pxl))
                best = std::move(refit);
        } catch (const std::exception &e) {
            logger.Warning("Powder calibration: that profile could not be fitted ({})", e.what());
        }
    }

    if (!best)
        return std::nullopt;

    // Only if it is actually better. Re-binning about a geometry that was itself wrong makes the fit
    // worse, not better, and the first pass is then the answer to keep - it is measured on a profile
    // binned about a geometry now believed wrong, which is weak, but it is the one the run has.
    if (best->ring_points * 2 < first.ring_points || best->rms_radial_pxl >= first.rms_radial_pxl) {
        logger.Info("Powder calibration: the re-binned fit is no better than the first "
                    "({} points at rms {:.3f} against {} at {:.3f}) - keeping the first pass",
                    best->ring_points, best->rms_radial_pxl, first.ring_points, first.rms_radial_pxl);
        return std::nullopt;
    }
    // The header the answer is compared against is the run's own, not the pass before it: a reader wants
    // to know how far the calibration moved the FILE.
    best->header_distance_mm = first.header_distance_mm;
    return best;
}

void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) {
    Logger logger("Rugnux");

    // Two-pass rotation runs this twice. The shadow does not move, and re-detecting would find
    // nothing (its pixels are masked by now) and so clear the mask the first pass established;
    // the beam centre must not be re-read either, because by then the geometry is post-refined.
    const auto &current = pixel_mask_.GetMask();
    const bool have_shadow = std::any_of(current.begin(), current.end(),
                                         [](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; });
    const bool want_shadow = config_.detect_beam_stop.has_value() && !have_shadow;
    const bool want_beam_center = config_.estimate_beam_center && !beam_center_placed_;
    // The spot width is a property of the crystal, not of the geometry, so the second pass reuses the
    // radius the first one measured rather than re-reading the frames.
    //
    // It is READ OFF the geometry as it stands, though, and a beam centre this same pass is about to
    // correct leaves that geometry provisional: candidates are excluded within MIN_BEAM_DISTANCE_PX
    // of the centre and stratified into resolution bands about it, so a placeholder header centre
    // both measures a different set of spots and labels them at the wrong resolution - and
    // spot_width_measured_ would then stop the corrected pass from ever revisiting it. Where a
    // second pass follows, leave the width to that one, which starts from the corrected centre. A
    // run with no second pass measures it here: a provisional radius beats none.
    const bool second_pass_follows = want_beam_center && config_.rotation_indexing
                                     && config_.two_pass_rotation && experiment_.IsRotationIndexing();
    const bool want_width = config_.adaptive_integration_radius && !spot_width_measured_
                            && !second_pass_follows;
    if (!want_shadow && !want_beam_center && !want_width)
        return;

    // The two consumers want different frames and each gets its own set. The shadow is built from
    // the even spread over the sweep it has always been built from, so that asking for a beam
    // centre cannot move the mask; the spot symmetry needs frames paired half a turn apart, where
    // the sweep reaches that far. The two sets barely overlap, so the loop below walks their union
    // and every frame is read once.
    //
    // Both samples are taken from the middle of the sweep, both ends left out: an even spread over
    // [0, images_to_process) always picks the very first and the very last image, and those are the
    // ones a shutter out of step spoils - the last carries under 90 % of the sample's median counts
    // on 12 of 38 regression sweeps and is essentially unexposed on 5.
    const auto goniometer = experiment_.GetGoniometer();
    const int margin = std::min(PRESCAN_END_MARGIN_IMAGES, images_to_process / 10);
    auto shadow_sample = select_equally_spaced_image_ordinals(images_to_process - 2 * margin,
                                                              frame_count);
    auto spot_sample = want_beam_center && goniometer
                          ? select_half_turn_paired_ordinals(
                                    images_to_process - 2 * margin,
                                    BEAM_CENTER_IMAGE_FACTOR * frame_count,
                                    // Degrees per ORDINAL: the shifted goniometer's increment already
                                    // is that, and rotation data cannot be strided anyway.
                                    std::abs(goniometer->GetIncrement_deg()))
                          : std::vector<int>();
    for (int &ordinal : shadow_sample)
        ordinal += margin;
    for (int &ordinal : spot_sample)
        ordinal += margin;
    const bool want_spot_symmetry = !spot_sample.empty();
    const std::set<int> shadow_set(shadow_sample.begin(), shadow_sample.end());
    const std::set<int> spot_set(spot_sample.begin(), spot_sample.end());
    std::set<int> sample = shadow_set;
    sample.insert(spot_set.begin(), spot_set.end());
    if (sample.empty())
        return;

    // Say which of the two this is reading for. It runs before the first image of the run proper is
    // processed, so without this the job sits at its initial status for the whole pre-scan - which is
    // tens of seconds once the beam centre asks for its own frames on top of the projection.
    if (observer)
        observer->OnPhase(want_shadow ? (want_beam_center ? "Beam stop and beam centre" : "Beam stop")
                                      : (want_beam_center ? "Beam centre" : "Spot width"));

    ShadowFinder finder(experiment_, pixel_mask_);

    // The spot symmetry and the spot width are read in the same pass as the projection rather than in
    // a pre-pass of their own. The image is preprocessed a second time for them: the shadow is
    // accumulated from raw counts and the spot finder works on the converted image.
    const bool want_spots = want_spot_symmetry || want_width;
    std::unique_ptr<AzimuthalIntegrationMapping> prescan_mapping;
    if (want_spots)
        prescan_mapping = std::make_unique<AzimuthalIntegrationMapping>(experiment_, pixel_mask_);

    // Everything reading one frame needs that cannot be shared between threads: its own scratch for
    // each of the two decompressions, its own preprocessed image and its own spot finder.
    struct PreScanWorker {
        std::vector<uint8_t> shadow_buffer;
        std::vector<uint8_t> decompression_buffer;
        std::unique_ptr<ImagePreprocessorCPU> preprocessor;
        std::unique_ptr<ImagePreprocessorBuffer> preprocessed;
        std::unique_ptr<ImageSpotFinder> spot_finder;
    };
    const auto make_worker = [&] {
        PreScanWorker w;
        if (want_spots) {
            w.preprocessor = std::make_unique<ImagePreprocessorCPU>(experiment_, pixel_mask_);
            w.preprocessed = std::make_unique<ImagePreprocessorBuffer>(experiment_.GetPixelsNum());
            if (config_.spot_finding.adaptive_threshold)
                w.spot_finder = std::make_unique<AdaptiveSpotFinderCPU>(*prescan_mapping);
            else
                w.spot_finder = std::make_unique<ImageSpotFinderCPU>(experiment_.GetXPixelsNumConv(),
                                                                      experiment_.GetYPixelsNumConv());
        }
        return w;
    };

    // The spots of one image, and what the two consumers take off them. Taken as a step of its own
    // because the sparse-sweep pass below reads images the shadow does not. The frame they belong to
    // is stamped by the caller, which is what keeps the pool independent of the order the workers
    // finish in.
    const auto find_spots = [&](PreScanWorker &w, CompressedImage &image, int image_idx,
                                bool for_beam_center, std::vector<BeamCenterSpot> &out,
                                bool for_width, std::vector<spot_width::FluxCurve> &curves) {
        try {
            w.preprocessor->Analyze(*w.preprocessed,
                                    image.GetUncompressedPtr(w.decompression_buffer), image.GetMode());
        } catch (const std::exception &e) {
            if (IsFatalResourceError(e)) throw;
            logger.Warning("Pre-scan: failed to preprocess image {}: {}", image_idx, e.what());
            return;
        }
        auto spots = w.spot_finder->Run(*w.preprocessed, config_.spot_finding);
        // The width goes first: it wants the whole spot list, to know which spots are isolated.
        if (for_width)
            MeasureSpotFluxCurves(*w.preprocessed, static_cast<int>(experiment_.GetXPixelsNumConv()),
                                  static_cast<int>(experiment_.GetYPixelsNumConv()),
                                  experiment_.GetDiffractionGeometry(), spots, curves);
        if (!for_beam_center)
            return;
        // The strongest of a crowded frame: the symmetry is over-determined either way, and the
        // matching is quadratic in the spots of one frame.
        if (spots.size() > BEAM_CENTER_SPOTS_PER_IMAGE) {
            std::partial_sort(spots.begin(), spots.begin() + BEAM_CENTER_SPOTS_PER_IMAGE, spots.end(),
                              [](const DiffractionSpot &a, const DiffractionSpot &b) {
                                  return a.Count() > b.Count(); });
            spots.resize(BEAM_CENTER_SPOTS_PER_IMAGE);
        }
        for (const auto &spot : spots) {
            const Coord centroid = spot.RawCoord();
            out.push_back({centroid.x, centroid.y, static_cast<float>(spot.Count()), 0});
        }
    };

    std::vector<float> frame_angle_deg;
    std::vector<BeamCenterSpot> beam_center_spots;
    std::vector<spot_width::FluxCurve> width_curves;
    // Images the width ended up being measured on, for the log: the tiers below stop as soon as the
    // answer has settled, so this is a property of the crystal and worth reporting.
    size_t width_images = 0;
    // Frames that could actually be read. Counted here rather than asked of the shadow finder, which
    // is only fed when the shadow is wanted and would report zero for a pass that read every frame
    // for the beam centre alone.
    std::atomic<size_t> frames_read{0};

    // Read the sample on several workers. The reader serialises on the HDF5 lock, but the
    // decompression, the projection and the spot finding - which is all of the cost on a large
    // detector - run in parallel. Each worker accumulates into a shard of its own, so nothing is
    // locked while an image is added, and the per-frame results are stitched together in sample
    // order below so the beam centre sees the same input however the workers interleaved.
    {
        const std::vector<int> ordinals(sample.begin(), sample.end());
        std::vector<std::vector<BeamCenterSpot>> spots_of(ordinals.size());
        std::vector<std::vector<spot_width::FluxCurve>> curves_of(ordinals.size());
        // A frame joins the pool only if it could be read, as it did when this was a serial loop.
        std::vector<char> spot_read(ordinals.size(), 0);
        const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
                                                 std::min(PRESCAN_MAX_WORKERS, ordinals.size()));
        finder.SetShardCount(nworkers);

        // The order the frames are visited in, and where each width tier ends in it. The tiers
        // stride the PROJECTION's sample rather than the union read below, so asking for a beam
        // centre - which brings frames of its own - cannot move which frames the width is measured
        // on; the frames only the beam centre wants are appended to the last tier. Without a width
        // to measure there is one tier and the sample is read in its own order.
        std::vector<size_t> visit;
        std::vector<size_t> tier_end;
        visit.reserve(ordinals.size());
        {
            std::vector<size_t> width_pos;
            for (size_t i = 0; i < ordinals.size(); i++)
                if (shadow_set.contains(ordinals[i])) width_pos.push_back(i);
            std::vector<char> taken(ordinals.size(), 0);
            const size_t n_tier = want_width ? WIDTH_TIER_STRIDE.size() : 1;
            for (size_t tier = 0; tier < n_tier; tier++) {
                const size_t stride = want_width ? WIDTH_TIER_STRIDE[tier] : 1;
                for (size_t j = 0; j < width_pos.size(); j += stride)
                    if (!taken[width_pos[j]]) {
                        taken[width_pos[j]] = 1;
                        visit.push_back(width_pos[j]);
                    }
                tier_end.push_back(visit.size());
            }
            for (size_t i = 0; i < ordinals.size(); i++)
                if (!taken[i]) visit.push_back(i);
            tier_end.back() = visit.size();
        }

        // Kept across the tiers rather than rebuilt for each. A worker's preprocessed image is a
        // buffer the size of the detector, so each is still built on its own worker the first time
        // that worker runs and its pages are touched by the thread that reads them.
        std::vector<PreScanWorker> workers(nworkers);
        std::vector<char> worker_built(nworkers, 0);

        bool width_settled = !want_width;
        std::optional<float> previous_r80;
        size_t phase_begin = 0;
        for (const size_t phase_end : tier_end) {
            std::atomic<size_t> next{phase_begin};
            std::vector<std::future<void>> futures;
            futures.reserve(nworkers);
            for (size_t t = 0; t < nworkers; t++)
                futures.emplace_back(std::async(std::launch::async, [&, t] {
                    if (!worker_built[t]) {
                        workers[t] = make_worker();
                        worker_built[t] = 1;
                    }
                    PreScanWorker &w = workers[t];
                    for (size_t v = next.fetch_add(1); v < phase_end; v = next.fetch_add(1)) {
                        const size_t i = visit[v];
                        const int ordinal = ordinals[i];
                        const int image_idx = start_image + ordinal * config_.stride;
                        std::shared_ptr<JFJochReaderRawImage> img;
                        try {
                            img = reader_.GetRawImage(image_idx);
                        } catch (const std::exception &e) {
                            if (IsFatalResourceError(e)) throw;
                            logger.Warning("Pre-scan: failed to load image {}: {}", image_idx,
                                           e.what());
                            continue;
                        }
                        if (!img) continue;

                        DataMessage msg{};
                        msg.image = img->image;
                        msg.number = ordinal;
                        msg.original_number = image_idx;
                        frames_read.fetch_add(1, std::memory_order_relaxed);
                        // Gated on want_shadow, not on the frame set alone: AddImage projects the
                        // whole image, and with --detect-beam-stop=off (or a shadow already in the
                        // mask) that is the pass's largest per-frame cost paid for a result nobody
                        // reads. The set still decides WHICH frames feed it, so the mask is built
                        // from the frames it always was.
                        if (want_shadow && shadow_set.contains(ordinal))
                            finder.AddImage(msg, w.shadow_buffer, t);

                        // The width is measured on the projection's own even spread over the sweep,
                        // so asking for a beam centre cannot move it either.
                        const bool for_beam_center = spot_set.contains(ordinal);
                        const bool for_width = !width_settled && shadow_set.contains(ordinal);
                        if (for_beam_center) spot_read[i] = 1;
                        if (for_beam_center || for_width)
                            find_spots(w, msg.image, image_idx, for_beam_center, spots_of[i],
                                       for_width, curves_of[i]);
                    }
                }));
            for (auto &f : futures) f.get();
            phase_begin = phase_end;

            // What this tier makes of the width, pooled in sample order so the answer does not
            // depend on how the workers interleaved. A tier too sparse to measure settles nothing:
            // the test is against what the smaller sample actually said, not against the radius it
            // fell back to.
            if (width_settled) continue;
            width_images = static_cast<size_t>(
                    std::count_if(visit.begin(), visit.begin() + static_cast<std::ptrdiff_t>(phase_end),
                                  [&](size_t i) { return shadow_set.contains(ordinals[i]); }));
            width_curves.clear();
            for (const auto &c : curves_of)
                width_curves.insert(width_curves.end(), c.begin(), c.end());
            const auto r80 = spot_width::R80AtReference(width_curves);
            width_settled = r80 && previous_r80 && spot_width::WidthSettled(*r80, *previous_r80);
            previous_r80 = r80;
        }

        // Frame numbering follows the sample order, exactly as the serial read did.
        for (size_t i = 0; i < ordinals.size(); i++) {
            if (!spot_read[i]) continue;
            const int frame = static_cast<int>(frame_angle_deg.size());
            // The goniometer has been shifted so that a LOCAL ordinal maps to the angle of the
            // original image it came from (see RunPipeline), so the angle is asked for by ordinal,
            // as everywhere else.
            frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast<float>(ordinals[i])));
            for (auto &spot : spots_of[i]) {
                spot.frame = frame;
                beam_center_spots.push_back(spot);
            }
        }
    }

    // The integration radius from the recorded spot width. r1 is the aperture the integrator learns
    // the profile width over and a second moment over a disk of radius a saturates at a^2/4, so at the
    // shipped r1 = 4 the learned sigma cannot exceed 2 px and a wider spot is fitted with a profile
    // the model cannot represent. r2 follows r1, and r3 is taken so that the r2..r3 background ring
    // keeps the area it has at the shipped 4/6/13 - a ring that shrinks with the disk is what makes a
    // bare --integration-radius worse than the default it replaces.
    if (want_width) {
        spot_width_measured_ = true;
        const auto r80 = spot_width::R80AtReference(width_curves);
        if (!r80) {
            logger.Info("Spot width: not measurable on {} spots from {} images, keeping the "
                        "integration radius at r1={:.1f}", width_curves.size(), width_images,
                        experiment_.GetBraggIntegrationSettings().GetR1());
        } else {
            const float r1 = spot_width::R1ForWidth(*r80);
            const float r2 = r1 + 2.0f;
            const float r3 = std::sqrt(r2 * r2 + 133.0f);
            BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
            // Remembered only where the radius actually moved, so the guard in RunAllPasses has
            // something to give back and cannot announce a fallback that changes nothing.
            if (r1 != bis.GetR1()) bragg_before_adaptive_ = bis;
            bis.R1(r1).R2(r2).R3(r3);
            experiment_.ImportBraggIntegrationSettings(bis);
            logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots from {} images) => "
                        "integration radii r1={:.1f} r2={:.1f} r3={:.2f}", *r80, spot_width::D_REF_A,
                        width_curves.size(), width_images, r1, r2, r3);
        }
    }

    if (frames_read.load() == 0) {
        logger.Warning("Pre-scan: no image could be read. Skipping.");
        return;
    }

    if (want_shadow) {
        const auto shadow = finder.GetMask(config_.nthreads);
        const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u);
        pixel_mask_.LoadBeamStopMask(experiment_, shadow);
        // Not "this beamline has no beam stop". The search floods outwards from seeds within 4 px of
        // the ASSUMED beam centre, so a centre far enough out puts every seed in a module gap and the
        // flood finds nothing - which makes an empty shadow a statement about the centre. It is the
        // sharpest free predicate there is: over 29 non-SLS datasets it fires on 3, all 3 of which had
        // a bad outcome, with no false alarm in 39 in-house sweeps (which report 0.44-1.51 % masked).
        if (shadowed == 0)
            logger.Warning("Beam stop: no shadow found in {} images. The search starts from the beam "
                           "centre in the file, so this points at that centre rather than at the "
                           "experiment; --estimate-beam-center measures it", finder.GetFrameCount());
        else
            logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images",
                        shadowed, 100.0 * static_cast<double>(shadowed) / static_cast<double>(shadow.size()),
                        finder.GetFrameCount());
    }

    // The beam centre comes after, so that the shadow is in the mask by the time it is placed and
    // the holder arm - the largest azimuthal asymmetry on the detector - is out of the way of the
    // background the centre is read from.

    // What the scattered background makes the centre, measured on EVERY run and committed on none.
    // The fit reads the projection the beam-stop pre-pass has just built and no frame of its own, so
    // it costs a fraction of a second on top of a pass that has already run - which is what makes it
    // affordable to do always, and doing it always is the point. A header beam centre is the metadata
    // field the field jokes about, and nothing in any package says how wrong a given one is: on a run
    // that indexes and merges perfectly well this line still reports that the file claims a centre the
    // data put 73 px away. The run is not touched by it. What consumes the number is the second first
    // pass in RunPipeline, and --estimate-beam-center's fall-through below, which is the same fit.
    //
    // Only where the projection exists: with --detect-beam-stop=off nothing accumulates it, and
    // projecting every frame for this alone is the largest per-frame cost of the pre-scan. On the
    // second pass of a two-pass run the shadow is already in the mask, so this does not run again -
    // by then the centre has been post-refined and the file's value is no longer what is in it.
    if (want_shadow && config_.beam_center_check) {
        background_center_ = FindBeamCenterFromBackground(experiment_, pixel_mask_,
                                                          finder.GetMeanProjection());
        const float need = BeamCenterNeed_pxl(experiment_);
        if (!background_center_) {
            logger.Info("Beam centre check: the background is too flat to place the centre on this "
                        "run - the file's ({:.2f},{:.2f}) is neither confirmed nor contradicted",
                        experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
        } else {
            const float dx = background_center_->beam_x_pxl - experiment_.GetBeamX_pxl();
            const float dy = background_center_->beam_y_pxl - experiment_.GetBeamY_pxl();
            const float moved = std::hypot(dx, dy);
            const auto across = AcrossSpindle_pxl(experiment_, dx, dy);
            // Everything a reader needs to judge the header, on one line: what the file says, what the
            // background says, how far apart they are, how well the background knows its own answer,
            // and what this particular geometry asks of the centre. The last is what turns a bare pixel
            // count into a verdict - 3 px is nothing at 100 mm and a lost lattice at 500 mm.
            logger.Info("Beam centre check: the file says ({:.2f},{:.2f}), the scattered background "
                        "says ({:.2f},{:.2f}) +- {:.2f} px - a difference of {:.2f} px, against the "
                        "{:.2f} px this geometry asks the centre to be right to. Measured only; the "
                        "run keeps the centre in the file",
                        experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
                        background_center_->beam_x_pxl, background_center_->beam_y_pxl,
                        background_center_->sigma_pxl, moved, need);
            // Three sigma, because the fit is a PRECISION and not an accuracy: on a background with no
            // curvature there is nothing to separate shift from amplitude and it follows the noise in
            // g', confidently. A move it cannot resolve is not evidence about the header either way.
            if (moved <= 3.0f * background_center_->sigma_pxl)
                logger.Info("Beam centre check: that is under three times the fit's own sigma, so the "
                            "file's centre is as good as this measurement can tell - the second first "
                            "pass below indexes it anyway, because a fraction of a pixel is enough to "
                            "decide a cell and no sigma small enough to gate on would notice");
            else if (across && *across > need)
                // Info and not a warning: `need` is where the FFT amplitude at a 200 A axis has fallen
                // to 0.70, which is a floor under what matters and not a threshold for failure - it is
                // exceeded on two of three ordinary in-house sweeps whose lattice comes out right. What
                // is worth a warning is what the second pass below then finds, not this.
                logger.Info("Beam centre check: {:.2f} px of the difference is ACROSS the spindle, more "
                            "than the {:.2f} px this geometry absorbs, and a centre error is fixed in "
                            "the lab frame - it smears the accumulated reciprocal-space cloud, so the "
                            "first pass can lose peaks or take an axis harmonic. The second first pass "
                            "below indexes both centres and reports whether it did", *across, need);
            else
                logger.Info("Beam centre check: that is a real difference, and {} - the second first "
                            "pass below indexes both centres and compares what they give",
                            across ? "it is mostly along the spindle, where nothing announces itself"
                                   : "this is not a rotation sweep, so there is no spindle to resolve it on");
        }
    }

    if (!want_beam_center)
        return;
    beam_center_placed_ = true;

    // The spot symmetry first where the sweep supports it: it is exact geometry rather than a
    // property of the scattering, and it measures the sets the background cannot - a flat radial
    // profile at long wavelength leaves that fit nothing to bite on. The background answers
    // wherever the sweep does not reach half a turn, which is every still and every screening
    // wedge, so the two cover between them what neither covers alone.
    std::string source = "spot symmetry";
    std::optional<BeamCenterEstimate> estimate;
    SpindleEstimate spindle;
    // The bound an estimate has to meet. One-sided by construction: the flat constant was calibrated
    // on a 244 mm / 0.95 A / 0.15 mm geometry and is far too tight for a loose one, so where the
    // geometry asks less of the centre the bound is relaxed to what it asks and never below the
    // constant. Nothing it accepts today is lost.
    const float need = BeamCenterNeed_pxl(experiment_);
    const float ceiling = std::max(MAX_BEAM_CENTER_SIGMA_PXL, need);
    if (want_spot_symmetry)
        estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, beam_center_spots,
                                                  config_.fit_spindle ? &spindle : nullptr);

    // Where it does not come out - no answer at all, or one that moves when the search is started
    // elsewhere - read more of the sweep and ask again. Both estimators are counts of matched spot
    // PAIRS, the Friedel one over frames half a turn apart and the second crossing over any two
    // frames one reflection is recorded on, so a sweep whose frames carry few spots, or whose pairs
    // match rarely, is measurable only from more of them. More data first and a weaker estimator
    // second: this runs before the fall-through, and only what survives it falls through. The extra
    // frames go to the spot finder alone, so the beam-stop projection keeps the images it was
    // validated on.
    if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > ceiling)) {
        // What has already been read FOR THE SPOTS, which is not the whole sample: the frames the
        // shadow was built from carry none of them, and testing against the union would silently
        // skip the ones this pass exists to read.
        std::vector<int> extra;
        for (int ordinal : select_half_turn_paired_ordinals(
                     images_to_process - 2 * margin, BEAM_CENTER_SPARSE_IMAGES,
                     std::abs(goniometer->GetIncrement_deg()))) {
            ordinal += margin;
            if (!spot_set.contains(ordinal))
                extra.push_back(ordinal);
        }
        if (!extra.empty()) {
            logger.Info("Beam centre: the spot symmetry does not come out on {} frames "
                        "({} spots); reading {} more frames",
                        frame_angle_deg.size(), beam_center_spots.size(), extra.size());
            // Its own phase: this reads several times the frames the pass above did, and is the one
            // step of the pre-scan a user is likely to sit through wondering what is happening.
            if (observer)
                observer->OnPhase("Beam centre (reading more frames)");
            // Read on workers as above - this pass is several times the size of the first one - and
            // append in `extra` order so the pool does not depend on how they interleaved.
            std::vector<std::vector<BeamCenterSpot>> extra_spots(extra.size());
            std::vector<char> extra_read(extra.size(), 0);
            const size_t nworkers = std::min<size_t>(std::max<size_t>(config_.nthreads, 1),
                                                     std::min(PRESCAN_MAX_WORKERS, extra.size()));
            std::atomic<size_t> next{0};
            std::vector<std::future<void>> futures;
            futures.reserve(nworkers);
            for (size_t t = 0; t < nworkers; t++)
                futures.emplace_back(std::async(std::launch::async, [&] {
                    PreScanWorker w = make_worker();
                    for (size_t i = next.fetch_add(1); i < extra.size(); i = next.fetch_add(1)) {
                        const int image_idx = start_image + extra[i] * config_.stride;
                        std::shared_ptr<JFJochReaderRawImage> img;
                        try {
                            img = reader_.GetRawImage(image_idx);
                        } catch (const std::exception &e) {
                            if (IsFatalResourceError(e)) throw;
                            logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what());
                            continue;
                        }
                        if (img) {
                            extra_read[i] = 1;
                            // The extra frames go to the beam centre alone; the width has already
                            // been measured on the projection's own sample.
                            find_spots(w, img->image, image_idx, true, extra_spots[i],
                                       false, width_curves);
                        }
                    }
                }));
            for (auto &f : futures) f.get();

            for (size_t i = 0; i < extra.size(); i++) {
                if (!extra_read[i]) continue;
                const int frame = static_cast<int>(frame_angle_deg.size());
                frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast<float>(extra[i])));
                for (auto &spot : extra_spots[i]) {
                    spot.frame = frame;
                    beam_center_spots.push_back(spot);
                }
            }
            estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg,
                                                      beam_center_spots,
                                                      config_.fit_spindle ? &spindle : nullptr);
        }
    }
    if (config_.fit_spindle && std::isfinite(spindle.azimuth_rad))
        logger.Info("Spindle: azimuth about the beam {:.3f} mrad, tip towards it {:.3f} mrad "
                    "(Friedel vote {:.0f} -> {:.0f})",
                    1e3 * spindle.azimuth_rad, 1e3 * spindle.tip_rad,
                    spindle.vote_excess_nominal, spindle.vote_excess);
    if (!estimate || estimate->sigma_pxl > ceiling) {
        // Say so. Below about 220 deg of sweep the spot symmetry never clears this bound and the
        // background answers instead - measured within 0.23 px of the full-sweep centre at every span
        // from 120 to 220 deg - so the fall-through is what covers short sweeps, and until this line
        // existed nothing in the log said the spot arm had even been tried.
        if (want_spot_symmetry && estimate)
            logger.Info("Beam centre: the spot symmetry answers ({:.2f},{:.2f}) at sigma {:.2f} px, "
                        "over the {:.2f} px ceiling - falling through to the background",
                        estimate->beam_x_pxl, estimate->beam_y_pxl, estimate->sigma_pxl, ceiling);
        else if (want_spot_symmetry)
            logger.Info("Beam centre: the spot symmetry does not come out on this sweep - falling "
                        "through to the background");
        source = "background";
        // The same fit the check above has already done wherever the projection exists - same mask,
        // same geometry, same projection - so this asks for it again only where it was not run.
        estimate = background_center_
                       ? background_center_
                       : FindBeamCenterFromBackground(experiment_, pixel_mask_, finder.GetMeanProjection());
    }
    if (!estimate) {
        logger.Info("Beam centre: not measurable, keeping ({:.2f},{:.2f})",
                    experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
        return;
    }
    const float dx = estimate->beam_x_pxl - experiment_.GetBeamX_pxl();
    const float dy = estimate->beam_y_pxl - experiment_.GetBeamY_pxl();
    const float moved = std::hypot(dx, dy);
    const auto across = AcrossSpindle_pxl(experiment_, dx, dy);
    // Two conditions, and the second is the new one. A move under three times the estimator's own
    // sigma is not a measurement of anything, and a centre that is not moved cannot move the
    // two-pass loop off its fixed point either. On the 39 rotation regression crystals this adopts
    // 38 and keeps the header on one, a 0.74 px move at sigma 0.32.
    const bool commit = estimate->sigma_pxl <= ceiling && moved > 3.0f * estimate->sigma_pxl;
    const char *verdict = "COMMIT";
    if (!commit)
        verdict = estimate->sigma_pxl > ceiling ? "reject: sigma over the ceiling (kept header)"
                                                : "reject: the move is under 3 sigma (kept header)";
    logger.Info("Beam centre from {}: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved {:.2f} px, "
                "sigma {:.2f} px against a {:.2f} px ceiling => {}",
                source, experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
                estimate->beam_x_pxl, estimate->beam_y_pxl, moved, estimate->sigma_pxl,
                ceiling, verdict);
    if (across)
        logger.Info("Beam centre: {:.2f} px of that move is across the spindle, where this geometry "
                    "asks for {:.2f} px. That is a floor and not a tolerance - a move 0.03x of it has "
                    "decided a halved axis", *across, need);
    if (across && *across > need)
        logger.Warning("Beam centre: the file's centre is {:.2f} px out across the spindle, more than "
                       "the {:.2f} px this geometry absorbs - expect the first pass to lose peaks or "
                       "to take an axis harmonic", *across, need);
    if (commit)
        experiment_.BeamX_pxl(estimate->beam_x_pxl).BeamY_pxl(estimate->beam_y_pxl);
}

void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process,
                                  RugnuxObserver *observer) {
    Logger logger("Rugnux");

    // Rotation has its own two-pass; this is stills only.
    if (experiment_.IsRotationIndexing()) {
        logger.Warning("--refine-geometry is stills-only; rotation data uses two-pass indexing. Ignoring.");
        return;
    }
    // The bundle anchors a known cell; without one there is nothing to refine against.
    const auto cell = experiment_.GetUnitCell();
    if (!cell.has_value()) {
        logger.Warning("--refine-geometry needs a reference unit cell (-C / -S / reference MTZ). Skipping.");
        return;
    }

    const int refine_frames = std::max(1, config_.refine_geometry.value());
    if (observer)
        observer->OnPhase("Geometry refinement (first pass)");

    const auto dataset = reader_.GetDataset();

    // Index a spread sample large enough to yield ~refine_frames strong frames even at a low hit rate.
    const int sample_budget = std::min(images_to_process, std::max(refine_frames * 50, 8000));
    const std::vector<int> sample = select_equally_spaced_image_ordinals(images_to_process, sample_budget);

    // First-pass engines, built from the current (nominal) geometry: an index-only pass (no
    // integration) to obtain each frame's spots + assigned HKL + orientation.
    AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_);
    IndexerThreadPool pool(experiment_.GetIndexingSettings(), IndexerConstruction::OnFirstUse);
    IndexAndRefine indexer(experiment_, &pool, /*retain_outcomes=*/false);

    auto pass_settings = config_.spot_finding;
    pass_settings.enable = true;
    pass_settings.indexing = true;
    pass_settings.quick_integration = false;

    // A frame is worth bundling only if it indexed this many spots; the bundle then takes the
    // refine_frames strongest of those.
    constexpr int MIN_STRONG_SPOTS = 20;
    // Stop sampling once there are comfortably more strong frames than the bundle will use, instead of
    // draining the whole budget. The budget is sized for a low-hit-rate dataset (refine_frames * 50, at
    // least 8000), so on data that indexes well it meant re-indexing the ENTIRE run to keep 200 frames -
    // 99% of this pass was sampling. Four times the bundle leaves the "strongest N" selection a real
    // pool to choose from.
    const int strong_target = refine_frames * 4;

    // The sample is cut into a FIXED number of interleaved stripes - stripe s is sample positions
    // s, s+STRIPES, s+2*STRIPES, ... - and each stripe stops once it has contributed its share of the
    // strong frames. Two things follow, both of which the previous shared-cursor-with-a-shared-stop
    // arrangement got wrong. Every stripe is spread over the WHOLE run, so stopping early no longer
    // means fitting the geometry to the beginning of it. And which frames get examined depends only
    // on the data: it is not a race between the workers, and it does not change with -N, because a
    // stripe is processed identically whichever worker happens to claim it.
    constexpr int STRIPES = 32;
    const int stripe_strong_target = std::max(1, (strong_target + STRIPES - 1) / STRIPES);

    std::vector<GeomRefineFrame> frames;
    std::mutex frames_mutex;
    std::atomic<int> next_stripe = 0;
    std::atomic<int> examined = 0;

    auto worker = [&]() {
        pin_gpu();   // round-robin per worker thread; must precede engine construction
        MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, indexer,
                                       /*enable_fused_adaptive_gpu=*/true);
        AzimuthalIntegrationProfile profile(mapping);

        for (int stripe = next_stripe.fetch_add(1); stripe < STRIPES && !cancelled_;
             stripe = next_stripe.fetch_add(1)) {
        int stripe_strong = 0;
        for (int idx = stripe; idx < static_cast<int>(sample.size()) && !cancelled_
                               && stripe_strong < stripe_strong_target; idx += STRIPES) {
            examined.fetch_add(1, std::memory_order_relaxed);
            const int ordinal = sample[idx];
            const int image_idx = start_image + ordinal * config_.stride;

            std::shared_ptr<JFJochReaderRawImage> img;
            try {
                img = reader_.GetRawImage(image_idx);
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                logger.Warning("Geometry refinement: failed to load image {}: {}", image_idx, e.what());
                continue;
            }
            if (!img) continue;

            DataMessage msg{};
            msg.image = img->image;
            msg.number = ordinal;
            msg.original_number = image_idx;
            if (dataset->efficiency.size() > image_idx)
                msg.image_collection_efficiency = dataset->efficiency[image_idx];

            try {
                analysis.Analyze(msg, profile, pass_settings);
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                continue;
            }
            if (!msg.indexing_result.value_or(false) || !msg.indexing_lattice.has_value())
                continue;

            GeomRefineFrame f;
            f.lattice = *msg.indexing_lattice;
            f.ordinal = ordinal;
            for (const auto &s : msg.spots)
                if (s.indexed && s.lattice == 0)
                    f.spots.push_back(GeomRefineSpot{s.x, s.y,
                        static_cast<int32_t>(s.h), static_cast<int32_t>(s.k), static_cast<int32_t>(s.l)});
            if (f.spots.size() >= 6) {
                if (static_cast<int>(f.spots.size()) >= MIN_STRONG_SPOTS)
                    stripe_strong++;
                const std::unique_lock ul(frames_mutex);
                frames.push_back(std::move(f));
            }
        }
        }
    };

    std::vector<std::future<void>> futures;
    futures.reserve(config_.nthreads);
    for (int i = 0; i < config_.nthreads; ++i)
        futures.push_back(std::async(std::launch::async, worker));
    for (auto &fut : futures)
        fut.get();

    if (cancelled_)
        return;

    // Keep the strongest frames (most indexed spots) for the bundle - the true cell indexes many
    // spots per frame, and orientation diversity comes for free from independent serial stills.
    std::vector<GeomRefineFrame> strong;
    for (auto &f : frames)
        if (static_cast<int>(f.spots.size()) >= MIN_STRONG_SPOTS)
            strong.push_back(std::move(f));
    // frames is in the order the workers happened to finish, so the spot count alone does not order
    // it: without the ordinal tie-break, equally strong frames would swap places between runs and a
    // different bundle would be refined.
    std::sort(strong.begin(), strong.end(),
              [](const GeomRefineFrame &a, const GeomRefineFrame &b) {
                  if (a.spots.size() != b.spots.size()) return a.spots.size() > b.spots.size();
                  return a.ordinal < b.ordinal;
              });
    if (static_cast<int>(strong.size()) > refine_frames)
        strong.resize(refine_frames);

    logger.Info("Geometry refinement: indexed {} of {} sampled frames (budget {}), bundling {} strong "
                "frames (>= {} spots)",
                frames.size(), std::min<int>(examined.load(), static_cast<int>(sample.size())),
                sample.size(), strong.size(), MIN_STRONG_SPOTS);

    GeometryRefinerSettings gr_settings;
    gr_settings.crystal_system = experiment_.GetCrystalSystem();
    gr_settings.num_threads = config_.nthreads;

    const GeometryRefinerResult r = RefineGlobalGeometry(
        experiment_.GetDiffractionGeometry(), *cell, strong, gr_settings);
    if (!r.ok) {
        logger.Warning("Geometry refinement did not run (too few strong frames / did not converge); "
                       "keeping the input geometry");
        return;
    }

    logger.Info("Geometry refinement: beam ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, distance {:.4f} -> {:.4f} mm",
                experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), r.beam_x_px, r.beam_y_px,
                experiment_.GetDetectorDistance_mm(), r.distance_mm);
    logger.Info("Geometry refinement: cell a/b/c {:.3f}/{:.3f}/{:.3f} -> {:.3f}/{:.3f}/{:.3f} A "
                "(median residual {:.3f} px, {} frames, {} spots)",
                cell->a, cell->b, cell->c, r.cell.a, r.cell.b, r.cell.c,
                r.median_residual_px, r.frames_used, r.spots_used);

    experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm);
    experiment_.SetUnitCell(r.cell);
    // The beam centre now comes from indexed spots, which the background cannot better.
    beam_center_placed_ = true;
}

ProcessResult Rugnux::Run(RugnuxObserver *observer) {
    // Each pass times only itself, so the canonical result's processing_time_s is the last pass alone -
    // on the rotation two-pass that is under a third of what the run actually took, the pre-scan and the
    // first pass being invisible in it. Time the whole thing here, where every pass is inside.
    const auto run_start = std::chrono::steady_clock::now();
    ProcessResult result = RunAllPasses(observer);
    result.total_time_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - run_start).count();
    return result;
}

ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
    // Rotation two-pass geometry post-refinement: the first pass integrates at the header geometry and
    // post-refines the detector geometry (distance + beam from the observed spot positions, cell scale + axis
    // from phi_obs); the second pass re-indexes and re-integrates with the refined geometry. Only the second
    // (refined-geometry) pass writes merged files, under the plain "<prefix>_*" name. The first pass still runs
    // under a "<prefix>_01" prefix so that the _process.h5 it writes, when one is asked for, does not collide
    // with the second pass's - the ice-ring flags and geometry stored in each are the ones that pass used.
    if (config_.mode == ProcessMode::FullAnalysis && config_.rotation_postrefine_geometry
        && experiment_.IsRotationIndexing()) {
        Logger logger("Rugnux");
        const std::string base_prefix = config_.output_prefix;
        const auto gonio_snapshot = experiment_.GetGoniometer();
        const int64_t max_spot_count_snapshot = experiment_.GetDatasetSettings().GetMaxSpotCount();
        prepass_detector_geometry_.reset();
        prepass_rotation_scale_.reset();
        prepass_result_.reset();
        force_rotation_result_.reset();
        bragg_adaptive_.reset();

        logger.Info("Rotation two-pass geometry post-refinement: first pass (header geometry) -> {}_01_*",
                    base_prefix);
        // Keep an empty prefix empty (the "compute stats, persist nothing" mode): the "_01" suffix would
        // make it non-empty and RunPipeline would then write stray _01_* files despite no output wanted.
        config_.output_prefix = base_prefix.empty() ? base_prefix : base_prefix + "_01";
        auto pass1 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/true);
        pass1.pass_number = 1;
        pass1.pass_count = 2;
        if (cancelled_) { config_.output_prefix = base_prefix; return pass1; }
        // Undo the pre-pass goniometer shift (a sub-range or strided run moves start/increment inside
        // the pass), but KEEP an axis direction the pass changed. The rotation-axis sign rescue runs
        // inside the pass and adopts the opposite axis when the file's own indexes almost nothing;
        // restoring the snapshot wholesale threw that away, handing pass 2 the sign that failed and
        // then applying pass 1's fitted rotation scale to the axis it was not fitted on. It
        // self-healed - pass 2's rescue fires again and a wrong sign indexes almost no frames - but
        // only by repeating a whole first pass, and only while "a wrong sign always scores below
        // half" holds, which is an invariant nothing states or tests.
        if (gonio_snapshot) {
            GoniometerAxis restored = *gonio_snapshot;
            if (const auto after_pass1 = experiment_.GetGoniometer())
                restored.Axis(after_pass1->GetAxis());
            experiment_.Goniometer(restored);
        }
        // A measured spot budget is the pass's own, not the run's: give the second pass the same list to
        // measure from, so the two passes cannot ratchet each other down. The canonical pass then reads
        // the budget off the refined geometry, where the spots that do lie on the lattice actually do.
        experiment_.MaxSpotCount(max_spot_count_snapshot);

        // Apply the post-refined detector geometry for the second pass, keeping the header geometry so
        // the run can go back to it if the refined pass turns out worse (see the quality guard below).
        const std::array<float, 3> header_geometry{experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
                                                   experiment_.GetDetectorDistance_mm()};
        if (prepass_detector_geometry_) {
            const auto &g = *prepass_detector_geometry_;
            experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]);
        }
        // ... and the goniometer rotation scale, on the same measure-then-re-integrate footing: the angles
        // in the file are the commanded ones, so a stage that ran fast is a geometry error like any other.
        if (prepass_rotation_scale_ && gonio_snapshot) {
            experiment_.Goniometer(ScaleRotation(*gonio_snapshot, *prepass_rotation_scale_));
            // The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has
            // just stopped using, and the override can only ever raise the second pass's own estimate (it
            // takes the larger of the two). Carrying it over would hold the second pass at the rocking
            // width the uncorrected angles produced - the correction half-applied. Drop it and let the
            // second pass fit its own.
            prepass_mosaicity_.clear();
            logger.Info("Two-pass: goniometer rotation scale {:.5f} will drive the second integration pass "
                        "(oscillation {:.4f} -> {:.4f} deg per image)", *prepass_rotation_scale_,
                        gonio_snapshot->GetIncrement_deg(),
                        gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_);
        }
        // The integration radius the pre-scan measured from the crystal's own spot width, applied now
        // that the geometry is fitted. The pre-pass ran at the radius the run started with (see
        // bragg_adaptive_), so everything pass 2 consumes from it - the post-refined distance and beam,
        // the space group, the smoothed mosaicity - is what the fixed radius produces, and the widened
        // disk changes only the intensities of the canonical pass.
        if (bragg_adaptive_) {
            experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_);
            logger.Info("Two-pass: the canonical pass integrates at the measured spot width, "
                        "r1={:.1f} r2={:.1f} r3={:.2f}", bragg_adaptive_->GetR1(),
                        bragg_adaptive_->GetR2(), bragg_adaptive_->GetR3());
        }

        logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*",
                    base_prefix);
        config_.output_prefix = base_prefix;   // the refined pass is the canonical result (no _02 suffix)
        auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
        // The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever
        // result is returned - it is a result of the run, not of the pass that happened to fit it.
        pass2.post_refine = pass1.post_refine;
        pass2.pass_number = 2;
        pass2.pass_count = 2;

        // r2 = r1 + 2 is also the inner edge of the background ring, so widening the signal disk pushes
        // that ring out into the neighbours; a reflection whose ring is left with five or fewer clean
        // pixels has no background and is dropped whole. On a pattern dense enough in three dimensions
        // at once that costs a quarter of the observations - and NOTHING available before integrating
        // predicts it: the crystal it happens to has a mid-table predicted spot spacing, wider than
        // crystals that survive r1 = 12. What does measure it is the integrator counting its own drops.
        // The pass that integrates at the widened radius is now the canonical one, so the measurement
        // and the remedy are both here: re-run it at the radius that keeps the data, the same remedy the
        // two guards below use and, like them, a pass only on a crystal that was going to be wrong.
        if (!cancelled_ && bragg_adaptive_ && bkg_starved_fraction_
            && *bkg_starved_fraction_ > spot_width::BKG_STARVED_MAX_FRACTION) {
            logger.Info("Two-pass: at the widened integration radius r1={:.1f} the neighbouring "
                        "reflections left {:.2f}% of the predicted reflections without a background "
                        "ring (bound {:.2f}%), so this pass is re-run at r1={:.1f} r2={:.1f} r3={:.2f}",
                        bragg_adaptive_->GetR1(), 100.0 * *bkg_starved_fraction_,
                        100.0 * spot_width::BKG_STARVED_MAX_FRACTION,
                        bragg_before_adaptive_->GetR1(), bragg_before_adaptive_->GetR2(),
                        bragg_before_adaptive_->GetR3());
            const double starved = *bkg_starved_fraction_;
            const float wide_r1 = bragg_adaptive_->GetR1();
            const int passes = pass2.pass_count + 1;
            experiment_.ImportBraggIntegrationSettings(*bragg_before_adaptive_);
            bragg_adaptive_.reset();
            pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
            pass2.post_refine = pass1.post_refine;
            pass2.pass_number = passes;
            pass2.pass_count = passes;
            pass2.pass_decision = fmt::format(
                    "post-refined geometry adopted, re-run at the fixed integration radius: at r1={:.1f} "
                    "neighbouring reflections left {:.2f}% of the predicted reflections without a "
                    "background ring", wide_r1, 100.0 * starved);
        }

        // Supercell-collapse guard for the second pass. The de-novo lattice search can be BISTABLE on a
        // periodically-modulated dataset: the tiny post-refine geometry move tips it out of the true-cell
        // basin onto a spurious SUPERCELL (a large multiple of pass-1's cell volume), and the space group
        // reused from pass-1 is then stamped onto that wrong-shaped cell - catastrophic. Pass 1 (header
        // geometry, un-perturbed) found the correct lattice, so re-run the second pass with pass-1's FULL
        // indexing result forced (its refined lattice + orientation + rotation axis + geometry - not just the
        // cell, which loses the per-frame orientation and integrates worse): the re-run then matches de-novo
        // quality but at the refined geometry. Only fires when pass 2 is markedly LARGER than pass 1 - a
        // genuine de-novo demotion to a smaller primitive (the pseudo-symmetry case the second pass re-indexes
        // de novo to catch) has pass2 <= pass1 and is left untouched.
        if (!cancelled_ && prepass_result_ && pass1.consensus_cell && pass2.consensus_cell) {
            const double v1 = gemmi::UnitCell(*pass1.consensus_cell).volume;
            const double v2 = gemmi::UnitCell(*pass2.consensus_cell).volume;
            if (v1 > 0.0 && v2 > 1.5 * v1) {
                logger.Info("Two-pass: second-pass cell volume {:.0f} A^3 is {:.2f}x pass-1 ({:.0f} A^3) - a "
                            "spurious supercell; re-running the second pass with pass-1's result forced",
                            v2, v2 / v1, v1);
                force_rotation_result_ = *prepass_result_;
                // Pass 1's result carries pass 1's goniometer, and the per-image path takes its angles
                // from there - forcing it whole would put the uncorrected angles back and silently undo
                // the rotation scale for the re-run.
                if (prepass_rotation_scale_ && force_rotation_result_->axis)
                    force_rotation_result_->axis =
                            ScaleRotation(*force_rotation_result_->axis, *prepass_rotation_scale_);
                pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
                pass2.post_refine = pass1.post_refine;
                pass2.pass_number = 3;
                pass2.pass_count = 3;
                pass2.pass_decision = fmt::format(
                        "post-refined geometry adopted, re-run with pass-1's lattice forced: the de-novo "
                        "cell volume {:.0f} A^3 was {:.2f}x pass-1's ({:.0f} A^3), a spurious supercell",
                        v2, v2 / v1, v1);
                force_rotation_result_.reset();
            }
        }
        // Pass 2 is normally the better answer, which is why it is the canonical output - but it is not
        // guaranteed to be. Two ways it can be wrong: it merges more unique reflections than the cell it
        // settled on can hold (completeness above 100% is arithmetically impossible and means the cell is
        // wrong), or its CC1/2 collapses relative to pass 1. Both bounds are set where only a failure
        // reaches them, so a normal run - where pass 2 is a little better - keeps pass 2.
        //
        // The comparison is made on each pass's SEARCH merge, not on the final one. Each pass
        // determines its own space group, so the two final merges can be in different groups, and
        // completeness and CC1/2 in different groups are not the same measurement. The search merge is
        // P1 on both sides, over the full resolution range and before any correction surface is fitted
        // - the one merge the two passes make in the same terms.
        if (!cancelled_ && pass1.has_merge_statistics && pass2.has_merge_statistics) {
            const double compl1 = pass1.search_merge_completeness;
            const double compl2 = pass2.search_merge_completeness;
            const double cc1 = pass1.search_merge_cc_half;
            const double cc2 = pass2.search_merge_cc_half;
            constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5;   // rounding headroom, nothing more
            constexpr double MAX_CC_HALF_LOSS = 0.05;
            if (pass2.lattice_conflicts_with_fixed_sg
                || compl2 > MAX_CREDIBLE_COMPLETENESS || cc2 < cc1 - MAX_CC_HALF_LOSS) {
                logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass "
                               "(completeness {:.1f}% vs {:.1f}%, CC1/2 before corrections {:.3f} vs "
                               "{:.3f}) - going back to the header geometry. The refined geometry did "
                               "not help this crystal.",
                               compl2, compl1, cc2, cc1);
                // The refined pass has already written the canonical files, so re-run at the header
                // geometry to replace them - the same remedy the supercell collapse above uses, and it
                // only costs a pass on a crystal that was going to be wrong otherwise.
                experiment_.BeamX_pxl(header_geometry[0]).BeamY_pxl(header_geometry[1])
                           .DetectorDistance_mm(header_geometry[2]);
                if (prepass_rotation_scale_ && gonio_snapshot) {
                    // The header rotation angles with it - but not the axis DIRECTION, which is not
                    // part of the geometry being reverted and may have been corrected by the sign
                    // rescue (see the restore before pass 2).
                    GoniometerAxis restored = *gonio_snapshot;
                    if (const auto current = experiment_.GetGoniometer())
                        restored.Axis(current->GetAxis());
                    experiment_.Goniometer(restored);
                }
                config_.output_prefix = base_prefix;
                auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
                redo.post_refine = pass1.post_refine;
                redo.pass_number = pass2.pass_count + 1;
                redo.pass_count = pass2.pass_count + 1;
                redo.pass_decision = fmt::format(
                        "header geometry re-adopted: the post-refined pass was worse "
                        "(completeness {:.1f}% vs {:.1f}%, CC1/2 before corrections {:.3f} vs {:.3f})",
                        compl2, compl1, cc2, cc1);
                return redo;
            }
            if (pass2.pass_decision.empty())
                pass2.pass_decision = fmt::format(
                        "post-refined geometry adopted (completeness {:.1f}% vs {:.1f}%, "
                        "CC1/2 before corrections {:.3f} vs {:.3f})",
                        compl2, compl1, cc2, cc1);
        }
        if (pass2.pass_decision.empty())
            pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_)
                    ? "post-refined geometry adopted"
                    : "the post-refinement committed no geometry change, so this pass reproduces the first";
        return pass2;
    }
    return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
}

namespace {
    // Angle between the crystal's symmetry axes and the spindle. A rotation sweep never records the
    // reflections whose reciprocal vector lies close to the spindle - the blind cusp - and normally
    // symmetry fills it in from an equivalent elsewhere in reciprocal space. It cannot when a symmetry
    // axis IS the spindle: the cusp then maps onto itself, every reflection in it is equivalent only to
    // other reflections in it, and the hole stays empty however long the sweep runs. Returns the
    // smallest angle in degrees between any proper rotation axis of the space group and the goniometer
    // axis, with the order of the axis that achieves it.
    std::optional<std::pair<double, int>> ClosestSymmetryAxisToSpindle(const gemmi::SpaceGroup &sg,
                                                                      const CrystalLattice &lattice,
                                                                      const Coord &spindle) {
        const double spindle_len = spindle.Length();
        if (spindle_len < 1e-9)
            return std::nullopt;
        std::optional<std::pair<double, int>> best;
        for (const auto &op : sg.operations().derive_symmorphic().sym_ops) {
            if (op.rot == gemmi::Op::identity().rot)
                continue;
            // Order of the rotation, then project onto its invariant direction by summing its powers:
            // (1/n) sum_k W^k annihilates everything except the axis.
            const gemmi::Op w{op.rot, {0, 0, 0}, op.notation};
            int order = 0;
            gemmi::Op cur = gemmi::Op::identity();
            std::array<std::array<double, 3>, 3> sum{};
            for (int k = 0; k < 6; ++k) {
                for (int i = 0; i < 3; ++i)
                    for (int j = 0; j < 3; ++j)
                        sum[i][j] += static_cast<double>(cur.rot[i][j]) / gemmi::Op::DEN;
                cur = cur.combine(w).wrap();
                ++order;
                if (cur.rot == gemmi::Op::identity().rot)
                    break;
            }
            if (order < 2)
                continue;
            // Any non-degenerate column of the projector is the axis in fractional direct coordinates.
            for (int c = 0; c < 3; ++c) {
                const Coord axis = lattice.Vec0() * static_cast<float>(sum[0][c])
                                 + lattice.Vec1() * static_cast<float>(sum[1][c])
                                 + lattice.Vec2() * static_cast<float>(sum[2][c]);
                const double len = axis.Length();
                if (len < 1e-6 * lattice.Vec0().Length())
                    continue;
                const double cosang = std::abs((axis * spindle) / (len * spindle_len));
                const double ang = std::acos(std::min(1.0, cosang)) * 180.0 / PI;
                if (!best.has_value() || ang < best->first)
                    best = std::make_pair(ang, order);
                break;
            }
        }
        return best;
    }
}

ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) {
    Logger logger("Rugnux");
    ProcessResult result;

    // Each pass determines its own space group, so start every one of them from the group the run was
    // ASKED for. A pass that determined one leaves it on experiment_, and the pass after it - the
    // re-runs in Run(), or a second Run() on the same object - would otherwise inherit that answer and
    // skip the search that is the whole point of running again.
    experiment_.SetSpaceGroup(user_fixed_sg_);

    const auto dataset = reader_.GetDataset();
    if (!dataset)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "No experiment dataset found in the input file");

    if (config_.stride <= 0)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image stride must be positive");

    const auto total_images_in_file = static_cast<int>(reader_.GetNumberOfImages());
    int end_image = config_.end_image;
    if (end_image < 0 || end_image > total_images_in_file)
        end_image = total_images_in_file;
    const int start_image = config_.start_image;
    // Round UP: the worker loop below walks ordinals until start + ordinal * stride reaches end_image,
    // so a stride that does not divide the range evenly still processes a final partial step. Rounding
    // down here would size the per-image arrays one short of the ordinals that index them.
    const int images_to_process = (end_image - start_image + config_.stride - 1) / config_.stride;
    if (images_to_process <= 0) {
        logger.Warning("No images to process (start {}, end {}, stride {}, total {})",
                       start_image, end_image, config_.stride, total_images_in_file);
        return result;
    }

    const bool full = (config_.mode == ProcessMode::FullAnalysis);
    // Powder calibration runs the same per-image engine as the full analysis (spot finding lives inside
    // MXAnalysisWithoutFPGA) with indexing switched off in its settings. By rings it needs only the
    // azimuthal profile and the azint worker would do - but the spots are worth their half second even
    // then, because the circle through them places the beam centre without reference to the header at
    // all, which is the one hypothesis that survives a header the profile cannot correct.
    const bool calibration = (config_.mode == ProcessMode::Calibration);
    const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots;
    const bool per_image_analysis = full || calibration;
    const bool write_files = write_output && !config_.output_prefix.empty();

    // Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization,
    // space group, unit cell, ...) are configured on experiment_ by the caller.
    experiment_.BitDepthImage(32).PixelSigned(true);
    experiment_.Mode(DetectorMode::Standard);
    experiment_.OverwriteExistingFiles(true);
    // Offline processing: the output prefix is a local path the operator chose (rugnux -o or the
    // viewer's "next to the input file"), so use the trusted setter that permits an absolute path.
    experiment_.FilePrefixTrusted(config_.output_prefix.empty() ? "output" : config_.output_prefix);
    experiment_.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
    experiment_.ImagesPerTrigger(images_to_process);
    experiment_.NumTriggers(1);
    if (full)
        experiment_.Compression(CompressionAlgorithm::BSHUF_LZ4);

    // The pipeline indexes images 0..N-1 within this run; if we process a sub-range/strided
    // selection, shift the goniometer so local index i maps to the angle of original image
    // start+i*stride (keeping the per-image rotation wedge), otherwise rotation angles would be
    // wrong for any start_image != 0.
    if (const auto g = experiment_.GetGoniometer();
        g.has_value() && (start_image != 0 || config_.stride != 1)) {
        const float incr = g->GetIncrement_deg();
        GoniometerAxis shifted(g->GetName(),
                               g->GetStart_deg() + incr * static_cast<float>(start_image),
                               incr * static_cast<float>(config_.stride),
                               g->GetAxis(), g->GetHelicalStep());
        shifted.ScreeningWedge(g->GetScreeningWedge().value_or(incr));
        experiment_.Goniometer(shifted);
    }

    // Stills global geometry-refinement first pass: bundle-adjust the shared beam/distance/cell from a
    // sample of strongly-indexed frames and apply it to experiment_, so everything below (azimuthal
    // mapping, the main index+integrate pass, scaling) uses the refined geometry. Rotation is skipped
    // inside (it has its own two-pass).
    if (full && config_.refine_geometry.has_value())
        RefineStillsGeometry(start_image, end_image, images_to_process, observer);

    // After any geometry refinement, so the shadow is found about the beam centre actually used,
    // and before the azimuthal mapping and the output mask below, which both read pixel_mask_.
    if (config_.detect_beam_stop.has_value() || config_.estimate_beam_center
        || config_.adaptive_integration_radius)
        PreScan(start_image, images_to_process,
                config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES), observer);

    // The geometry pre-pass integrates at the radius the run started with; the width the pre-scan just
    // measured is applied to the canonical pass instead (RunAllPasses, between the two). See
    // bragg_adaptive_ for why - in short, the pre-pass keeps only the geometry it fits to the observed
    // reflection positions, and those positions move with the signal disk.
    if (geometry_prepass && bragg_before_adaptive_) {
        bragg_adaptive_ = experiment_.GetBraggIntegrationSettings();
        experiment_.ImportBraggIntegrationSettings(*bragg_before_adaptive_);
    }

    AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_);

    JFJochReceiverPlots plots;
    plots.Setup(experiment_, mapping);

    // Output file (NXmxIntegrated master that links back to the original images).
    StartMessage start_message;
    experiment_.FillMessage(start_message);
    start_message.arm_date = dataset->arm_date;
    start_message.az_int_bin_to_q = mapping.GetBinToQ();
    start_message.az_int_bin_to_two_theta = mapping.GetBinToTwoTheta();
    start_message.az_int_q_bin_count = mapping.GetQBinCount();
    start_message.az_int_phi_bin_count = mapping.GetAzimuthalBinCount();
    if (mapping.GetAzimuthalBinCount() > 1)
        start_message.az_int_bin_to_phi = mapping.GetBinToPhi();
    if (full) {
        start_message.rois = experiment_.ROI().ExportMetadata();
        start_message.max_spot_count = experiment_.GetMaxSpotCount();
    }
    start_message.master_suffix = "process";
    start_message.file_format = FileWriterFormat::NXmxIntegrated;
    start_message.write_master_file = true;
    start_message.write_images = false;
    // The process file LINKS to the original images rather than writing its own, which only an HDF5
    // input can offer: there is nothing for a virtual dataset to point at in a directory of CBFs. So
    // ask for the links only when the source is HDF5, and otherwise write the images into the file -
    // the per-image analysis is the point of it either way, and without pictures nothing downstream
    // could show the frame a reflection came from.
    if (auto *hdf5 = dynamic_cast<JFJochHDF5Reader *>(&reader_)) {
        // With the stride: image i of the process file is source image start_image + i * stride, so
        // the linked pictures line up with the per-image analysis written beside them.
        start_message.hdf5_source_data =
            hdf5->GetHDF5DataSource(start_image, images_to_process, config_.stride);

        // Everything describing the pixels must describe those linked files. experiment_ cannot: it
        // is pinned to signed 32-bit, the container the reader hands images out in, which would type
        // the virtual dataset int32 over unsigned 16- or 32-bit sources and silently convert every
        // value on read. Only consumers going through the virtual view ever saw it - our own reader
        // resolves the mapping and opens the source itself.
        const auto stored = hdf5->GetStoredPixelFormat();
        start_message.bit_depth_image = stored.bit_depth;
        // Unsigned only, for the reason FillMessage gives: DIALS remaps the top two codes whenever
        // the field is present, and on a signed image those land inside the trusted range.
        start_message.bit_depth_readout.reset();
        if (!stored.is_signed)
            start_message.bit_depth_readout = stored.bit_depth;
        start_message.pixel_signed = stored.is_signed;
        start_message.error_value = stored.is_signed
                                    ? -(int64_t(1) << (stored.bit_depth - 1))      // INTx_MIN
                                    : (int64_t(1) << stored.bit_depth) - 1;        // UINTx_MAX
        start_message.underload_value = stored.is_signed ? start_message.error_value.value() + 1 : 0;
    } else {
        start_message.write_images = true;   // no file to link to; the process file carries them
    }

    // The full-detector mask is 4 B/px and the ROI map 2 B/px - tens of megabytes each - and nothing
    // but the process file reads them. Copy them where a writer is actually built rather than here:
    // the default merging path writes no process file at all.
    auto fill_writer_maps = [&] {
        start_message.pixel_mask["default"] = pixel_mask_.GetMask(experiment_);
        if (full && !experiment_.ROI().empty())
            start_message.roi_map = experiment_.ExportROIMap();
    };

    std::unique_ptr<FileWriter> writer;
    std::unique_ptr<ProcessFileWriter> writer_queue;
    // Whether anything is in place to break the indexing ambiguity, for the warning further down: a
    // reference MTZ, or (serial stills) the reference computed from a model just below.
    bool ambiguity_reference = !config_.reference_data.empty();
    if (write_files && config_.write_process_h5) {
        fill_writer_maps();
        writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
                                              /*trusted_path=*/true);
        // Deep enough that a worker never waits for the writer in the normal case, shallow enough that
        // the backlog is bounded by the worker count rather than by the length of the run.
        writer_queue = std::make_unique<ProcessFileWriter>(*writer, 4 * std::max(1, config_.nthreads));
    }

    const char *mode_name = full ? "full analysis"
                                 : calibration ? "powder calibration" : "azimuthal integration";
    logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]",
                images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name);
    LogXDSGeometry(experiment_);
    if (observer)
        observer->OnPhase(full ? "Full analysis" : calibration ? "Powder calibration" : "Azimuthal integration");

    // Full-analysis shared engines.
    std::unique_ptr<IndexerThreadPool> indexer_pool;
    std::unique_ptr<IndexAndRefine> indexer;
    // The FFT search bound indexer_pool was built with. The long-axis rescue needs it to tell whether a
    // cell it recovered still fits, since a pool's transform is sized once, at construction.
    float indexing_bound_A = 0.0f;
    if (per_image_analysis) {
        auto indexing_settings = experiment_.GetIndexingSettings();
        // Two-pass rotation fires the FFT indexer only twice - the two first-pass schemes, run in
        // parallel (the main per-image loop then uses the forced lattice, no FFT). Size the pool to
        // those two so a rotation run doesn't pay to initialise (cuFFT plan + GPU alloc) indexers it
        // never uses. Stills fire the FFT once per image across all workers, so keep the default there.
        if (config_.rotation_indexing && config_.two_pass_rotation && experiment_.IsRotationIndexing()) {
            indexing_settings.IndexingThreads(2);
            // Two indexer threads leave the rest of the machine idle for the whole first pass, which is
            // a third of the run. Split each scheme's candidate-cell refinement over -N/2 threads so the
            // two together use -N. The split is by candidate cell, so the result does not depend on it.
            // Clamped to what IndexingSettings accepts - a big host would otherwise ask for more than
            // the setting's limit and throw before a single image is read.
            indexing_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
            RaiseFFTBoundForKnownCell(indexing_settings, experiment_);
        }
        indexing_bound_A = indexing_settings.GetFFT_MaxUnitCell_A();
        indexer_pool = std::make_unique<IndexerThreadPool>(indexing_settings, IndexerConstruction::OnFirstUse);
        indexer = std::make_unique<IndexAndRefine>(experiment_, indexer_pool.get());
        // With no per-image file to write, nothing reads the message's reflection list.
        if (!writer_queue)
            indexer->KeepReflectionsInMessage(false);
        // The reference that breaks the per-image indexing ambiguity: a reference MTZ where there is
        // one, otherwise intensities computed from the model, which carry the same information. Serial
        // stills need it here or not at all - each crystal is indexed in its own hand, so by the time
        // there is a merge for --model to fit, the hands have already been averaged together. Rotation
        // data have one lattice and one hand, so their ambiguity is settled at the merge instead.
        std::vector<MergedReflection> model_reference;
        if (!config_.reference_data.empty()) {
            indexer->ReferenceIntensities(config_.reference_data);
        } else if (!config_.model_path.empty() && !experiment_.IsRotationIndexing()
                   && experiment_.GetUnitCell().has_value() && experiment_.GetGemmiSpaceGroup().has_value()) {
            const auto d_min = experiment_.GetScalingSettings().GetHighResolutionLimit_A()
                .value_or(experiment_.GetDetectorMaxResolution_A());
            model_reference = ModelReferenceIntensities(
                config_.model_path, experiment_.GetUnitCell(),
                &experiment_.GetSpaceGroupOrP1(), d_min, logger);
            if (!model_reference.empty()) {
                indexer->ReferenceIntensities(model_reference);
                ambiguity_reference = true;
            }
        }
        // Second pass of the rotation two-pass: predict with the smoothed mosaicity fitted in the pre-pass.
        if (!geometry_prepass && !prepass_mosaicity_.empty())
            indexer->SetPredictionMosaicityOverride(prepass_mosaicity_);
        // The geometry pre-pass fits the detector geometry to spot POSITIONS and throws its intensities
        // away, so a fixed group's centring absences cost it half its events for nothing. A no-op when
        // no group is fixed - prediction is already in P there - so the de-novo path is untouched.
        if (geometry_prepass)
            indexer->PredictAllCentringNodes(true);
    }

    const auto start_time = std::chrono::steady_clock::now();

    // First pass of two-pass rotation indexing (full analysis only).
    if (full && force_rotation_result_.has_value()) {
        // Supercell-collapse fallback: force pass-1's WHOLE indexing result (lattice + refined orientation /
        // search metadata / axis), not just its lattice - a lattice-only force loses that metadata and
        // integrates markedly worse.
        indexer->ForceRotationIndexerResult(*force_rotation_result_);
        logger.Info("Rotation indexer result forced from pass 1 (lattice + geometry) - skipping first pass");
    } else if (full && config_.forced_rotation_lattice.has_value()) {
        indexer->ForceRotationIndexerLattice(*config_.forced_rotation_lattice);
        logger.Info("Rotation indexer lattice forced externally - skipping first pass");
    } else if (full && config_.rotation_indexing && config_.two_pass_rotation) {
        if (observer)
            observer->OnPhase("Rotation indexing (first pass)");

        // Mid-exposure goniometer angle for an image ordinal (matches the per-image path).
        auto rot_angle = [&](int ordinal) -> std::optional<float> {
            if (const auto g = experiment_.GetGoniometer())
                return g->GetAngle_deg(static_cast<float>(ordinal)) + g->GetWedge_deg() / 2.0f;
            return std::nullopt;
        };

        // Spots for a first-pass image ordinal, cached (a frame is reused across schemes and the
        // validation set; re-finding is expensive on the --redo-rotation-spots path). The workers below
        // hand their spots back and this thread puts them in, so the map itself is still only ever
        // touched from here and needs no locking.
        std::map<int, std::vector<SpotToSave>> spot_cache;
        // The spots of one frame, found with the caller's own analysis engine. A frame's spots are a
        // pure function of that frame and the settings - the engine carries nothing from one image to
        // the next, which is why the main image loop can give one to each of its workers - so it does
        // not matter which engine, or which thread, finds them.
        const auto find_spots = [&](MXAnalysisWithoutFPGA &analysis, AzimuthalIntegrationProfile &profile,
                                    int ordinal) {
            const int image_idx = start_image + ordinal * config_.stride;
            std::vector<SpotToSave> spots;
            try {
                if (auto img = reader_.GetRawImage(image_idx)) {
                    DataMessage m{};
                    m.number = ordinal;
                    m.original_number = image_idx;
                    auto first_pass = config_.spot_finding;
                    first_pass.indexing = false;
                    first_pass.quick_integration = false;
                    m.image = img->image;
                    if (dataset->efficiency.size() > image_idx)
                        m.image_collection_efficiency = dataset->efficiency[image_idx];
                    analysis.Analyze(m, profile, first_pass);
                    spots = std::move(m.spots);
                }
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                logger.Warning("First-pass spot read failed for image {}: {}", image_idx, e.what());
            }
            return spots;
        };

        // Find the spots of every ordinal in the list that is not cached yet, on several workers, and
        // put them in the cache.
        //
        // The cache ends up holding exactly what a serial loop would have put there: the results are
        // collected into a vector indexed by position and inserted afterwards, so neither the spots of
        // a frame nor the set of frames in the cache depends on how the workers interleaved.
        //
        // Worker t keeps its analysis engine for the whole first pass. Building one is a CUDA stream
        // and a full set of device buffers - measured at 0.13 s against 4 ms for a frame - and the
        // "wedge" scheme only discovers how far into the sweep it has to read as the frames arrive,
        // so this is called several times and rebuilding them each time cost a quarter of the pass.
        // The engine belongs to a card, not to a thread, so a later call only has to make that card
        // current again; workers take a card by index so an engine always meets the one it was
        // built on.
        const int spot_workers = get_gpu_count() > 0
                ? std::min(config_.nthreads, std::max(8, 4 * get_gpu_count()))
                : config_.nthreads;
        struct SpotEngine {
            std::unique_ptr<MXAnalysisWithoutFPGA> analysis;
            std::unique_ptr<AzimuthalIntegrationProfile> profile;
        };
        std::vector<SpotEngine> engines(std::max(spot_workers, 1));
        const auto use_gpu_of_worker = [](size_t t) {
            if (const int32_t gpus = get_gpu_count(); gpus > 0)
                set_gpu(static_cast<int32_t>(t) % gpus);
        };
        int prefetch_batches = 0;
        double prefetch_time_s = 0.0;
        const auto prefetch_spots = [&](const std::vector<int> &ordinals) {
            std::vector<int> wanted;
            for (const int ordinal : ordinals)
                if (!spot_cache.contains(ordinal))
                    wanted.push_back(ordinal);
            if (wanted.empty())
                return;
            const auto prefetch_start = std::chrono::steady_clock::now();
            ++prefetch_batches;

            std::vector<std::vector<SpotToSave>> found(wanted.size());
            std::atomic<size_t> next{0};
            std::vector<std::future<void>> futures;
            const size_t nworkers = std::min<size_t>(engines.size(), wanted.size());
            futures.reserve(nworkers);
            for (size_t t = 0; t < nworkers; t++)
                futures.emplace_back(std::async(std::launch::async, [&, t] {
                    use_gpu_of_worker(t);   // must precede engine construction
                    SpotEngine &e = engines[t];
                    if (!e.analysis) {
                        e.analysis = std::make_unique<MXAnalysisWithoutFPGA>(
                                experiment_, mapping, pixel_mask_, *indexer,
                                /*enable_fused_adaptive_gpu=*/true);
                        e.profile = std::make_unique<AzimuthalIntegrationProfile>(mapping);
                    }
                    for (size_t i = next.fetch_add(1); i < wanted.size(); i = next.fetch_add(1))
                        found[i] = find_spots(*e.analysis, *e.profile, wanted[i]);
                }));
            for (auto &f : futures)
                f.get();

            for (size_t i = 0; i < wanted.size(); i++)
                spot_cache.emplace(wanted[i], std::move(found[i]));
            prefetch_time_s += std::chrono::duration<double>(
                    std::chrono::steady_clock::now() - prefetch_start).count();
        };

        // Release the engines on workers of their own, so a stream and its device buffers are freed
        // with the card they were made on current.
        const auto release_spot_engines = [&] {
            std::vector<std::future<void>> futures;
            futures.reserve(engines.size());
            for (size_t t = 0; t < engines.size(); t++)
                futures.emplace_back(std::async(std::launch::async, [&, t] {
                    use_gpu_of_worker(t);
                    engines[t] = SpotEngine{};
                }));
            for (auto &f : futures)
                f.get();
        };

        // How many validation frames a candidate global lattice actually indexes per image, scored on
        // the *real* per-image path (force the candidate, then run the same indexing+refinement the
        // main loop uses). A cell that merely fits the accumulated first-pass cloud (a spurious
        // sub-lattice from an under-sampled scheme) indexes few frames, while the true cell indexes
        // many - this is the discriminator that lets the best scheme win.
        const std::vector<int> validation = select_equally_spaced_image_ordinals(
            images_to_process, std::min(images_to_process, 60));
        auto validation_settings = config_.spot_finding;
        validation_settings.indexing = true;
        validation_settings.quick_integration = false;
        auto count_indexed = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int {
            idx.ForceRotationIndexerResult(r);
            // Fill the spot cache first. It is a plain map filled on demand, and the loop below runs the
            // frames in parallel, where a lookup racing an insert is not something the map survives.
            // After the first candidate this costs nothing.
            prefetch_spots(validation);
            // Each frame's verdict is its own - the score is just how many of them index - so they run
            // together. Refining one frame's geometry is a few Ceres solves and no GPU, and it is the
            // same call the main image loop already makes from every one of its workers on this same
            // IndexAndRefine, writing only unit_cells[] under its own mutex. With the candidate forced,
            // GetLattice() returns it and the branch that would advance the indexer's own state is never
            // reached. The offline solver stops on an iteration count rather than a clock
            // (IndexAndRefine::RefineGeometryIfNeeded), so a loaded machine cannot change a verdict.
            std::atomic<int> count = 0;
            ParallelFor(static_cast<int>(validation.size()),
                        std::min<size_t>(validation.size(), config_.nthreads), [&](int i) {
                DataMessage m{};
                m.number = validation[i];
                m.spots = spot_cache.at(validation[i]);
                if (idx.IndexFrameOnly(m, validation_settings))
                    ++count;
            });
            return count.load();
        };

        // How deep into an image's intensity-ordered spot list this lattice is still being seen - the
        // measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the
        // spots rather than the frames: each frame's spots are tallied by rank and the run keeps the
        // ranks the tally supports (AddSpotBudgetEvidence / SpotBudgetFromEvidence). A budget deeper than
        // that admits detections that are not this crystal's reflections, and the frame gate then counts
        // them against the frame - which is how a larger budget can integrate FEWER images.
        auto measure_budget = [&](IndexAndRefine &idx, const RotationIndexerResult &r) -> int64_t {
            idx.ForceRotationIndexerResult(r);
            prefetch_spots(validation);
            std::mutex evidence_mutex;
            std::vector<int64_t> indexed(experiment_.GetMaxSpotCount(), 0);
            std::vector<int64_t> counted(indexed.size(), 0);
            const bool index_ice_rings = experiment_.GetIndexingSettings().GetIndexIceRings();
            ParallelFor(static_cast<int>(validation.size()),
                        std::min<size_t>(validation.size(), config_.nthreads), [&](int i) {
                DataMessage m{};
                m.number = validation[i];
                m.spots = spot_cache.at(validation[i]);
                if (!idx.IndexFrameOnly(m, validation_settings))
                    return;
                std::unique_lock ul(evidence_mutex);
                AddSpotBudgetEvidence(m.spots, index_ice_rings, indexed, counted);
            });
            return SpotBudgetFromEvidence(indexed, counted);
        };

        // Feed one first-pass scheme (a set of image ordinals) into its own rotation indexer, ready to
        // be indexed. Spots are pulled from the cache here; the FFT + refinement runs separately
        // (RunIndexing) so the two schemes' indexing can overlap.
        //
        // The feed itself stays in frame order and stops where it always did, but it keeps the cache a
        // chunk of frames ahead of itself so those frames' spots are found on several workers. A chunk
        // rather than the whole list because the "wedge" scheme stops on the accumulation cap, and how
        // far into the sweep that is only becomes known as the frames arrive.
        constexpr int SPOT_PREFETCH_CHUNK = 64;
        auto feed_scheme = [&](IndexerThreadPool &pool, const std::vector<int> &ordinals) {
            auto ri = std::make_unique<RotationIndexer>(experiment_, pool);
            for (size_t i = 0; i < ordinals.size(); i++) {
                if (cancelled_ || ri->AccumulationFull())
                    break;
                if (!spot_cache.contains(ordinals[i]))
                    prefetch_spots({ordinals.begin() + i,
                                    ordinals.begin() + std::min(ordinals.size(), i + SPOT_PREFETCH_CHUNK)});
                ri->ProcessImage(ordinals[i], spot_cache.at(ordinals[i]), rot_angle(ordinals[i]));
            }
            return ri;
        };

        // Sampling schemes, tried in order of preference. We index each and keep the one that indexes
        // the most frames; the schemes cover complementary failure modes and no single one is best for
        // every crystal:
        //  - "spread": images spread over the whole rotation - full angular range, coarse stride; the
        //    long-standing default, best for typical crystals.
        //  - "wedge": consecutive frames from the clean start - native stride keeps each reflection's
        //    rocking curve continuous across frames, so the FFT can trace a long axis whose fine
        //    reciprocal spacing the coarse spread cannot resolve (e.g. an F-cubic cell viewed down
        //    [111]); bounded by the accumulation cap to the least radiation-damaged early wedge.
        std::vector<int> all_consecutive(images_to_process);
        std::iota(all_consecutive.begin(), all_consecutive.end(), 0);
        const std::vector<int> spread_ordinals = select_equally_spaced_image_ordinals(
            images_to_process, config_.rotation_indexing_image_count);
        const std::vector<std::pair<std::string, std::vector<int>>> schemes = {
            {"spread", spread_ordinals},
            {"wedge",  all_consecutive},
        };

        // indexer_error: what the indexer threw, if it threw. An empty result with no error is a
        // verdict on the crystal; an empty result with one is a verdict on the machine.
        struct FirstPass { std::optional<RotationIndexerResult> result; int score = -1; std::string name; double vol = 0.0;
                           std::optional<std::string> indexer_error; };

        // Feed both schemes (single-threaded spot cache), run their RunIndexing() concurrently, then keep
        // the lattice that indexes the most validation frames on the real per-image path.
        auto pick_best = [&](IndexerThreadPool &pool, IndexAndRefine &idx) -> FirstPass {
            std::vector<std::unique_ptr<RotationIndexer>> ris;
            ris.reserve(schemes.size());
            for (const auto &[name, ordinals] : schemes) {
                if (cancelled_)
                    break;
                ris.push_back(feed_scheme(pool, ordinals));
            }
            std::vector<std::future<void>> index_futures;
            index_futures.reserve(ris.size());
            for (auto &ri : ris) {
                RotationIndexer *rp = ri.get();
                index_futures.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
            }
            for (auto &f : index_futures)
                f.get();

            FirstPass bp;
            for (size_t i = 0; i < ris.size(); i++) {
                if (cancelled_)
                    break;
                if (!bp.indexer_error)
                    bp.indexer_error = ris[i]->GetIndexerError();
                auto found = ris[i]->GetLattice();
                if (!found.has_value())
                    continue;
                RotationIndexerResult result = std::move(*found);
                int score = count_indexed(idx, result);
                const std::string &name = schemes[i].first;
                logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score,
                            static_cast<int>(validation.size()));

                // Metric symmetry that costs nearly every frame is not this crystal's symmetry. The
                // Bravais class is decided on the UNREFINED FFT candidate against a fixed angular
                // tolerance, so a lattice that is pseudo-symmetric to a few tenths of a degree gets
                // promoted a class too far; the constraint then snaps a real angle to the ideal one,
                // which throws nearly every reflection of every frame out of tolerance (measured: 2/60
                // frames for a monoclinic lattice promoted to C-orthorhombic, 39/60 for its primitive
                // cell). The indexer's own guard compares the two on the accumulated first-pass cloud,
                // where they differ by less than a factor 2 - far too little to act on.
                // The bar is a clear majority, not a margin: the 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 does index a few more (measured: 47->54 frames on a real
                // I-centred orthorhombic crystal, 49->58 on an F-cubic one - a 10% margin demotes both).
                // Only a constrained cell that fails outright while its unconstrained cell works is
                // evidence of a false promotion, so demand exactly that. It is the same "fails to index
                // half the frames" test the long-axis rescue below uses, and it leaves genuine symmetry
                // alone; the space group is in any case decided from the intensities later.
                const int majority = static_cast<int>(validation.size()) / 2;
                if (result.unconstrained && score < majority) {
                    RotationIndexerResult alt = *result.unconstrained;
                    const int alt_score = count_indexed(idx, alt);
                    if (alt_score > majority) {
                        logger.Info("Scheme '{}': {}-centred {} indexes {}/{} frames but its unconstrained "
                                    "cell indexes {}/{} - the metric symmetry is a false promotion, dropping it",
                                    name, result.search_result.centering,
                                    gemmi::crystal_system_str(result.search_result.system), score,
                                    static_cast<int>(validation.size()), alt_score,
                                    static_cast<int>(validation.size()));
                        result = std::move(alt);
                        score = alt_score;
                    }
                }

                // Compare PRIMITIVE volumes: two schemes can find the same lattice in different
                // settings, and a centred setting's cell is an exact integer multiple of its primitive
                // one - a rhombohedral lattice in hexagonal axes is exactly 3x its primitive
                // rhombohedral cell. Comparing the centred volumes makes the supercell test below fire
                // on that pair and "demote" a perfectly good setting to a threefold-smaller merge, which
                // is enough to change the space group the search then picks.
                const double vol = std::abs(
                        result.lattice.ToPrimitive(result.search_result.centering).CalcVolume());

                // A later scheme wins if it indexes clearly more frames (>10%).
                const bool clearly_more = static_cast<float>(score) > bp.score * 1.1f + 0.5f;

                // Integer-supercell tie-break. The validation-frame count saturates - a spurious axis
                // multiple (2x/3x...) indexes every frame its true cell does, so both schemes reach the
                // same frame total and the count alone cannot tell them apart; the default then keeps
                // whichever ran first. When the two schemes tie on frames but their cell volumes are
                // related by an integer factor >=2, the larger cell is that spurious supercell and the
                // smaller is the true reduced cell: take it, regardless of scheme order. A near-integer
                // ratio separates a real axis multiplication from a centering coincidence. Volumes are
                // primitive (see above), so a pure setting difference is a ratio of 1 and never fires.
                bool integer_subcell = false;
                double harmonic_ratio = 0.0;
                if (bp.result.has_value() && vol > 1.0 && bp.vol > 1.0) {
                    const bool tied = static_cast<float>(score) >= bp.score * 0.9f - 0.5f;
                    const double ratio = (vol < bp.vol) ? bp.vol / vol : vol / bp.vol;
                    const double nearest = std::round(ratio);
                    const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15;
                    integer_subcell = tied && !clearly_more && integer_multiple && vol < bp.vol;
                    // sqrt(3) as well as an integer: a hexagonal supercell that lengthens one axis by
                    // sqrt(3) is the other harmonic seen on real data and its volume ratio is not a
                    // whole number. This is a REPORT and not a decision - it is deliberately taken
                    // over the whole pair, `clearly_more` included, because the harmonic is often the
                    // one that indexes more frames (see below).
                    if (integer_multiple || std::abs(ratio - std::sqrt(3.0)) < 0.15)
                        harmonic_ratio = ratio;
                }

                // Two candidates whose primitive volumes differ by a small integer or by sqrt(3) are
                // one lattice and its axis harmonic, and what decides between them is the BEAM CENTRE.
                // A centre error along the spindle translates the derotated cloud rigidly, which the
                // FFT amplitude cannot see at all, so the peaks stay sharp and the indexer commits -
                // with confidence - to a sub-multiple. Measured on real data: a 0.12 px change of
                // centre, 0.03x of what the J0 law asks across the spindle, is the whole difference
                // between the deposited cell and a halved axis, and the halved one is the answer that
                // indexed MORE frames (100.00 % against 99.23 %). So the frame count cannot arbitrate
                // this pair, and until this line nothing in the run said what the pair meant.
                if (harmonic_ratio > 0.0)
                    logger.Warning("Scheme '{}' (primitive cell volume {:.0f} A^3) and '{}' ({:.0f} A^3) "
                                   "differ by {:.2f}x - one is an axis harmonic of the other. Which of "
                                   "the two is right is settled by the beam centre to a fraction of a "
                                   "pixel, and the harmonic is often the one that indexes MORE frames, "
                                   "so the frame counts above cannot decide it. Check the beam centre "
                                   "(--estimate-beam-center)", name, vol, bp.name, bp.vol, harmonic_ratio);

                if (!bp.result.has_value() || clearly_more || integer_subcell) {
                    if (integer_subcell)
                        logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at "
                                    "equal frame count - adopting the smaller true cell", name, vol,
                                    std::round(bp.vol / vol), bp.name, bp.vol);
                    bp.score = score;
                    bp.vol = vol;
                    bp.name = name;
                    bp.result = std::move(result);
                }
            }
            return bp;
        };

        // Find the spots of everything the first pass is certain to want in one go - both schemes'
        // opening frames and the validation set. A batch of workers costs an analysis engine each,
        // which is far more than a frame is to analyse, so one batch for the lot beats one per scheme.
        {
            std::vector<int> upfront = spread_ordinals;
            upfront.insert(upfront.end(), validation.begin(), validation.end());
            for (int i = 0; i < std::min(SPOT_PREFETCH_CHUNK, images_to_process); i++)
                upfront.push_back(i);
            std::sort(upfront.begin(), upfront.end());
            upfront.erase(std::unique(upfront.begin(), upfront.end()), upfront.end());
            prefetch_spots(upfront);
        }

        FirstPass best = pick_best(*indexer_pool, *indexer);

        // Rotation-axis SIGN rescue. The sign is a convention the input cannot settle: a PILATUS
        // miniCBF header carries only "Oscillation_axis omega" and no direction at all, and facilities
        // disagree about the sign of the NXmx vector (measured across ten of them: most want the
        // stored [-1,0,0], APS 19-ID and SPring-8 BL41XU want the opposite). Nothing downstream can
        // repair it either - the default spindle fit refines the axis DIRECTION by fractions of a
        // degree and cannot flip it, that being a discrete change.
        //
        // It is, however, decidable from the data: the wrong sign does not index at all. Measured on
        // one dataset by changing nothing but the sign: 100% of frames indexed one way, 0 of 60 the
        // other. So try the other sign and keep whichever indexes more frames - the same
        // validation-frame count the scheme choice above already uses, so no new metric and no new
        // threshold. Runs before the long-axis rescue: with the sign wrong every candidate lattice is
        // wrong, and the coarse pass below would spend itself on a mis-signed cloud.
        //
        // Only after a poor pass, so a correctly-signed file costs nothing; and spot finding is not
        // repeated, the cache being keyed by image and the spots not depending on the axis.
        if (!cancelled_ && best.score < 0.5 * static_cast<double>(validation.size())) {
            if (const auto gon = experiment_.GetGoniometer(); gon.has_value() && gon->IsScanning()) {
                GoniometerAxis flipped = *gon;
                flipped.Axis(-gon->GetAxis());
                experiment_.Goniometer(flipped);
                FirstPass alt = pick_best(*indexer_pool, *indexer);
                if (alt.result.has_value() && alt.score > best.score) {
                    logger.Info("Rotation axis sign: the file's axis indexes {}/{} validation frames, "
                                "the opposite sign {}/{} - adopting the opposite sign for this run",
                                best.score, static_cast<int>(validation.size()),
                                alt.score, static_cast<int>(validation.size()));
                    best = alt;
                } else {
                    experiment_.Goniometer(*gon);   // the file was right; put it back
                }
            }
        }

        // The background centre, tried as a SECOND FIRST PASS.
        //
        // The estimate itself is free - it is read off the projection --detect-beam-stop has already
        // built - and a first pass re-uses the spots it has already found, so indexing the run a
        // second time at the measured centre costs about what one rung of the ladder below costs.
        // That buys the comparison the on-failure ladder structurally cannot make: a centre wrong
        // ALONG the spindle does not fail. It keeps 96-100 % of frames indexed and quietly returns a
        // 2x, 3x or sqrt(3) axis harmonic, so no failure trigger ever fires on it. Indexing both
        // centres and comparing the two ANSWERS is what can see it.
        //
        // What must NOT arbitrate the two is the frame count. Acceptance is a fractional-Miller test,
        // so a cell twice as long has to place every spot twice as accurately to score the same:
        // measured on real data, two centres 0.12 px apart gave the deposited cell at 99.23 % and a
        // halved axis at 100.00 %, and the WRONG answer indexed better. A rule of the form "take the
        // centre that indexes more" would therefore pick wrong in exactly the case this comparison
        // exists for. So only what needs no arbiter is decided here:
        //
        //   * the file's centre does not index a majority -> nothing defends it, so try the measured
        //     one (and the axis sign with it, see below). A run that indexes nothing has nothing to
        //     lose; where the measured one is no better, the file's centre is put back;
        //   * the file's centre indexes and the two agree on the lattice -> keep the file's centre and
        //     say they agree. This is the overwhelming majority, and saying it is worth the pass: it
        //     is a free statement that the header is good enough for this crystal;
        //   * they index and disagree -> say so, with both cells and their primitive-volume ratio,
        //     and change nothing. Naming a disagreement is useful; resolving it needs an arbiter
        //     nobody has.
        //
        // Where the file's centre indexes nothing and the measured one does not either, the axis-sign
        // question is asked AGAIN at the measured centre before giving up. The two unknowns are
        // discrete and coupled, and the rescue above asked its question at a centre that may be 73 px
        // wrong: on the two depositions here whose header is the geometric centre of the detector,
        // each error hides the other, both runs sit at 0/60 until both are moved, and moving both
        // takes them to 60/60 on the deposited lattice. Asking the pair costs one more first pass on
        // a run that has already failed twice.
        if (!cancelled_ && config_.beam_center_check && background_center_) {
            // Consumed, not re-read: a two-pass rotation run comes back through here with the geometry
            // the first pass post-refined, and a pre-scan estimate of the FILE's centre says nothing
            // about that one. The pre-scan does not re-measure it either - by then the beam-stop shadow
            // is already in the mask, so the projection is not rebuilt.
            const BeamCenterEstimate measured = *background_center_;
            background_center_.reset();
            const float header_x = experiment_.GetBeamX_pxl(), header_y = experiment_.GetBeamY_pxl();
            const float measured_x = measured.beam_x_pxl;
            const float measured_y = measured.beam_y_pxl;
            const float moved = std::hypot(measured_x - header_x, measured_y - header_y);
            // No gate on how significant the move is. It is tempting to skip the second pass when the
            // two centres are within the fit's own noise - it would save it on a quarter of the runs
            // here - but that is deciding at a threshold before the evidence is in, and the evidence
            // is cheap. It would also gate away exactly the case worth asking about: measured on real
            // data, a 0.12 px change of centre is the whole difference between the deposited cell and
            // a halved axis, and 0.12 px is well inside the sigma of any fit that could gate it.
            {
                const int majority = static_cast<int>(validation.size()) / 2;
                experiment_.BeamX_pxl(measured_x).BeamY_pxl(measured_y);
                FirstPass alt = pick_best(*indexer_pool, *indexer);
                const bool header_indexes = best.result.has_value() && best.score > majority;
                bool measured_indexes = alt.result.has_value() && alt.score > majority;
                if (!header_indexes) {
                    // Nothing defends the file's centre: the data do not index at it, and the
                    // background says it is `moved` px away. So try the measured one, and if that
                    // does not index either, ask the AXIS-SIGN question again there. The two unknowns
                    // are discrete and coupled, and the sign rescue above asked its question at a
                    // centre that may itself be wrong - each error hides the other, and only moving
                    // both reveals anything.
                    bool joint = false;
                    const auto gonio = experiment_.GetGoniometer();
                    if (!measured_indexes && gonio.has_value() && gonio->IsScanning()) {
                        GoniometerAxis flipped = *gonio;
                        flipped.Axis(-gonio->GetAxis());
                        experiment_.Goniometer(flipped);
                        const FirstPass both = pick_best(*indexer_pool, *indexer);
                        if (both.result.has_value() && both.score > majority) {
                            logger.Warning("Beam centre from the background: the file's centre "
                                           "({:.2f},{:.2f}) indexes {}/{} validation frames and the "
                                           "measured ({:.2f},{:.2f}) {}/{}, but the measured centre "
                                           "WITH the opposite rotation-axis sign indexes {}/{}. The "
                                           "file has both wrong, and neither shows while the other "
                                           "stands; the run takes the measured centre and the opposite "
                                           "sign", header_x, header_y, best.score,
                                           static_cast<int>(validation.size()), measured_x, measured_y,
                                           alt.score, static_cast<int>(validation.size()), both.score,
                                           static_cast<int>(validation.size()));
                            alt = both;
                            measured_indexes = true;
                            joint = true;
                        } else {
                            experiment_.Goniometer(*gonio);
                        }
                    }
                    if (measured_indexes) {
                        if (!joint)
                            logger.Warning("Beam centre from the background: the file's ({:.2f},{:.2f}) "
                                           "indexes {}/{} validation frames and the measured "
                                           "({:.2f},{:.2f}) indexes {}/{} - the beam centre in this "
                                           "file is wrong, and the run adopts the measured one",
                                           header_x, header_y, best.score,
                                           static_cast<int>(validation.size()), measured_x, measured_y,
                                           alt.score, static_cast<int>(validation.size()));
                        best = alt;
                    } else {
                        // Nothing worked. Put the file's centre back, so the ladder below searches
                        // around the value the run started with rather than around a fit that has just
                        // been tested and found no better - and so a run that fails for some other
                        // reason fails at the geometry it was given.
                        experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y);
                        logger.Info("Beam centre check: neither the file's centre ({}/{} frames) nor "
                                    "the measured one {:.2f} px away ({}/{}) indexes this run - the "
                                    "file's centre is kept and this is not a beam-centre failure",
                                    best.score, static_cast<int>(validation.size()), moved, alt.score,
                                    static_cast<int>(validation.size()));
                    }
                } else {
                    experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y);
                    if (measured_indexes) {
                        // Primitive volumes, as everywhere else here: two passes 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 ratio = alt.vol < best.vol ? best.vol / alt.vol : alt.vol / best.vol;
                        const auto &hc = best.result->search_result.conventional.GetUnitCell();
                        const auto &mc = alt.result->search_result.conventional.GetUnitCell();
                        const bool same_class =
                                best.result->search_result.system == alt.result->search_result.system
                                && best.result->search_result.centering == alt.result->search_result.centering;
                        // The primitive volume is what the beam centre decides. A centre error along
                        // the spindle makes the FFT take an axis harmonic, which changes the volume by
                        // an integer factor or by sqrt(3); the Bravais class can differ for a reason
                        // that has nothing to do with the centre, so the two are reported separately.
                        if (ratio < 1.02 && same_class) {
                            logger.Info("Beam centre check: indexing at the measured centre "
                                        "({:.2f},{:.2f}) returns the same lattice as the file's - "
                                        "{}-centred {}, cells agreeing to {:.2f} % in volume. The {:.2f} px "
                                        "difference does not decide this crystal's cell",
                                        measured_x, measured_y, alt.result->search_result.centering,
                                        gemmi::crystal_system_str(alt.result->search_result.system),
                                        100.0 * (ratio - 1.0), moved);
                        } else if (ratio < 1.02) {
                            // Same cell, different metric symmetry read off it. That is a promotion
                            // decided on an angular tolerance against the unrefined FFT candidate, not
                            // a lattice the beam centre chose, so it is worth naming and is not the
                            // alarming case: the cell is not in question. Both cases in this corpus
                            // are of this kind, and in both the file's centre is the one that finds
                            // the symmetry.
                            logger.Info("Beam centre check: indexing at the measured centre "
                                        "({:.2f},{:.2f}) returns the same CELL as the file's - the "
                                        "primitive volumes agree to {:.2f} % - but a different metric "
                                        "symmetry: {}-centred {} from the file's centre against "
                                        "{}-centred {} from the measured one. The {:.2f} px difference "
                                        "does not decide this crystal's cell; the run keeps the file's "
                                        "centre and the symmetry it found", measured_x, measured_y,
                                        100.0 * (ratio - 1.0), best.result->search_result.centering,
                                        gemmi::crystal_system_str(best.result->search_result.system),
                                        alt.result->search_result.centering,
                                        gemmi::crystal_system_str(alt.result->search_result.system),
                                        moved);
                        } else {
                            const double nearest = std::round(ratio);
                            const bool harmonic = (nearest >= 2.0 && std::abs(ratio - nearest) < 0.15)
                                                  || std::abs(ratio - std::sqrt(3.0)) < 0.15;
                            logger.Warning("Beam centre check: the two centres DISAGREE about the "
                                           "lattice. The file's ({:.2f},{:.2f}) gives {}-centred {}, "
                                           "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} (primitive volume "
                                           "{:.0f} A^3, {}/{} frames); the measured ({:.2f},{:.2f}) gives "
                                           "{}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} "
                                           "(primitive volume {:.0f} A^3, {}/{} frames){}",
                                           header_x, header_y, best.result->search_result.centering,
                                           gemmi::crystal_system_str(best.result->search_result.system),
                                           hc.a, hc.b, hc.c, hc.alpha, hc.beta, hc.gamma, best.vol,
                                           best.score, static_cast<int>(validation.size()),
                                           measured_x, measured_y, alt.result->search_result.centering,
                                           gemmi::crystal_system_str(alt.result->search_result.system),
                                           mc.a, mc.b, mc.c, mc.alpha, mc.beta, mc.gamma, alt.vol,
                                           alt.score, static_cast<int>(validation.size()),
                                           harmonic ? fmt::format(" - the volumes differ by {:.2f}x, so one "
                                                                  "is an axis harmonic of the other and the "
                                                                  "beam centre is what settles which",
                                                                  ratio)
                                                    : "");
                            // Deliberately no decision. The frame counts above are the one statistic
                            // that must not choose: on the measured case the harmonic indexed MORE.
                            logger.Warning("Beam centre check: the run keeps the file's centre. The "
                                           "frame counts cannot arbitrate this - a cell twice as long "
                                           "must place every spot twice as accurately to score the "
                                           "same, so the harmonic is often the one that indexes more. "
                                           "Re-run with --estimate-beam-center to measure the centre "
                                           "from the spot symmetry, which is exact geometry");
                        }
                    } else {
                        logger.Info("Beam centre check: the measured centre ({:.2f},{:.2f}) indexes "
                                    "{}/{} validation frames against the file's {}/{} - the file's "
                                    "centre stands", measured_x, measured_y, alt.score,
                                    static_cast<int>(validation.size()), best.score,
                                    static_cast<int>(validation.size()));
                    }
                }
            }
        }

        // Beam centre as an INDEXING HYPOTHESIS. A beam-centre error is not repairable downstream: it
        // is fixed in the LAB frame, so accumulating a sweep smears every reciprocal-lattice point
        // around a circle and the FFT amplitude at an axis of length a is multiplied by
        // J0(2 pi delta p a/(D lambda)). Past the first zero the true axis is gone and its harmonic
        // wins, which is what a first-pass axis doubling usually is. It is decidable from the data on
        // exactly the count the scheme choice and the axis-sign rescue already use: the right centre
        // indexes and the wrong one does not.
        //
        // BOTH detector directions are searched, and that is the one thing to keep. The J0 law is
        // about the FFT AMPLITUDE, which is translation-invariant, so a centre error along the spindle
        // is free for the transform - but the step after it, fitting a lattice whose origin is the
        // beam, is not, and a rigidly shifted lattice is fitted by a finer one. Measured by injection
        // on two crystals: across the spindle the cell stays right and the indexed fraction collapses
        // (99 % to 25 %), while along it the run keeps 96-100 % indexed and quietly adopts a 2x, 3x or
        // sqrt(3) supercell. The direction that announces itself is not the dangerous one.
        //
        // Only after a poor pass, so a correct header costs nothing; and spot finding is not repeated,
        // the cache being keyed by image and the spot POSITIONS not depending on the centre.
        if (!cancelled_ && !beam_center_searched_ && config_.beam_center_search_pxl > 0
            && best.score < 0.5 * static_cast<double>(validation.size())) {
            beam_center_searched_ = true;
            const float beam_x = experiment_.GetBeamX_pxl(), beam_y = experiment_.GetBeamY_pxl();
            const int reach = config_.beam_center_search_pxl;
            logger.Info("Beam-centre hypothesis: {}/{} validation frames at ({:.2f},{:.2f}); trying "
                        "the centre +-{} px along each detector axis, a pixel at a time",
                        best.score, static_cast<int>(validation.size()), beam_x, beam_y, reach);
            // Adopt only a hypothesis that INDEXES - the same majority the guards above test - rather
            // than whichever rung scores best. Taking any improvement was tried and is wrong: on a
            // dataset where no centre works the ladder wandered to the far end of its range on 10/60
            // against 5/60 and reported a 30 px move, which is the count's noise floor.
            const int adopt_bar = static_cast<int>(validation.size()) / 2;
            bool adopted = false;
            float won_x = beam_x, won_y = beam_y;
            for (int k = 1; k <= reach && !cancelled_ && !adopted; k++) {
                for (const auto &[dx, dy]: {std::pair<int, int>{-k, 0}, {k, 0}, {0, -k}, {0, k}}) {
                    const float trial_x = beam_x + static_cast<float>(dx);
                    const float trial_y = beam_y + static_cast<float>(dy);
                    experiment_.BeamX_pxl(trial_x).BeamY_pxl(trial_y);
                    FirstPass alt = pick_best(*indexer_pool, *indexer);
                    logger.Info("Beam-centre hypothesis ({:.2f},{:.2f}): {}/{} frames", trial_x,
                                trial_y, alt.score, static_cast<int>(validation.size()));
                    if (alt.result.has_value() && alt.score > adopt_bar && alt.score > best.score) {
                        best = alt;
                        won_x = trial_x;
                        won_y = trial_y;
                        adopted = true;
                        break;
                    }
                }
            }
            experiment_.BeamX_pxl(won_x).BeamY_pxl(won_y);
            if (adopted)
                logger.Info("Beam centre from indexing: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved "
                            "{:.2f} px, {}/{} validation frames", beam_x, beam_y, won_x, won_y,
                            std::hypot(won_x - beam_x, won_y - beam_y), best.score,
                            static_cast<int>(validation.size()));
        }

        // Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced
        // axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let
        // a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the
        // many high-order reflections off along the fine axis, so the true cell scores worst on the raw
        // cloud). Recover the true metric with a second search whose FFT unit-cell bound is widened past
        // the default, then RE-INDEX at full resolution constrained by the cell that comes back, as a
        // reference (the -C path): the reference filter drops the collapsed/supercell candidates and
        // refines an accurate global lattice. Only runs after a poor standard pass, so well-indexing
        // crystals never pay for it and their result is untouched.
        auto max_axis = [](const RotationIndexerResult &r) {
            const auto uc = r.lattice.GetUnitCell();
            return std::max({uc.a, uc.b, uc.c});
        };
        if (!cancelled_ && best.result.has_value()
            && best.score < 0.5 * static_cast<double>(validation.size())) {
            auto coarse_settings = experiment_.GetIndexingSettings();
            // The resolution setting only sizes the FFT histogram - it does not filter which spots go
            // into it - so this is here to pay for the widened bound below, not to make the long axis
            // easier to pick out. Measured, coarsening does the opposite: the acceptance angle widens
            // (0.159 -> 0.579 deg) but the true axis falls from rank 1996 to 16361 of 16384 directions.
            coarse_settings.FFT_HighResolution_A(3.5f);
            // The widened bound is what the rescue turns on, and it costs nothing: the histogram length is
            // 4*pi*max_length/high_resolution, so coarsening the resolution to 3.5 A has already
            // SHRUNK the transform by 1.75x - spending that back on reach leaves the coarse pass at
            // 1.14x the standard one. Without it the rescue cannot recover what it exists to recover:
            // a cell whose long axis is past the default 500 A bound is unrepresentable, so the
            // "recovered" metric would still be a sub-multiple. Only reached after a poor pass, so a
            // well-indexing crystal never pays. Stills never reach here.
            constexpr float RESCUE_MAX_UNIT_CELL_A = 1000.0f;
            if (coarse_settings.GetFFT_MaxUnitCell_A() < RESCUE_MAX_UNIT_CELL_A)
                coarse_settings.FFT_MaxUnitCell_A(RESCUE_MAX_UNIT_CELL_A);
            coarse_settings.IndexingThreads(2);
            coarse_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
            IndexerThreadPool coarse_pool(coarse_settings, IndexerConstruction::OnFirstUse);
            // Coarse first pass: keep the recovered cell with the LONGEST axis directly. Its full-
            // resolution per-frame validation would be low (the coarse cell is metrically right but
            // imprecise), so do NOT score it here - it is only the reference for the constrained re-index.
            std::vector<std::unique_ptr<RotationIndexer>> cris;
            cris.reserve(schemes.size());
            for (const auto &[name, ordinals] : schemes) {
                if (cancelled_) break;
                cris.push_back(feed_scheme(coarse_pool, ordinals));
            }
            std::vector<std::future<void>> cfut;
            cfut.reserve(cris.size());
            for (auto &ri : cris) {
                RotationIndexer *rp = ri.get();
                cfut.push_back(std::async(std::launch::async, [rp] { rp->RunIndexing(); }));
            }
            for (auto &f : cfut) f.get();
            std::optional<RotationIndexerResult> coarse_ref;
            for (auto &ri : cris) {
                auto r = ri->GetLattice();
                if (r.has_value() && (!coarse_ref.has_value() || max_axis(*r) > max_axis(*coarse_ref)))
                    coarse_ref = r;
            }
            // A recovered axis at least 30% longer than the standard pass found signals a rescued metric.
            if (!cancelled_ && coarse_ref.has_value()
                && max_axis(*coarse_ref) > 1.3 * max_axis(*best.result)) {
                logger.Info("Long-axis rescue: coarse pass recovered a {:.0f} A axis (was {:.0f} A); "
                            "re-indexing constrained by that cell", max_axis(*coarse_ref), max_axis(*best.result));
                experiment_.SetUnitCell(coarse_ref->lattice.GetUnitCell());
                // The re-index must be able to REPRESENT the cell the coarse pass just recovered. The
                // standard pool was sized at the default bound, so re-using it here would search for a
                // long axis it cannot hold and hand back the sub-multiple the rescue exists to replace.
                // count_indexed forces the lattice and never touches the pool, so the existing
                // IndexAndRefine is reused; only the FFT needs the wider one.
                auto wide_settings = experiment_.GetIndexingSettings();
                wide_settings.IndexingThreads(2);
                wide_settings.RefineThreads(std::clamp(config_.nthreads / 2, 1, MAX_REFINE_THREADS));
                RaiseFFTBoundForKnownCell(wide_settings, experiment_);
                FirstPass constrained;
                if (wide_settings.GetFFT_MaxUnitCell_A() > indexing_bound_A) {
                    IndexerThreadPool wide_pool(wide_settings, IndexerConstruction::OnFirstUse);
                    constrained = pick_best(wide_pool, *indexer);
                } else {
                    constrained = pick_best(*indexer_pool, *indexer);
                }
                experiment_.SetUnitCell(std::nullopt);   // leave the space-group / cell determination de-novo
                if (constrained.result.has_value() && constrained.score > best.score)
                    best = constrained;
            }
        }
        release_spot_engines();   // no more frames to read; give the cards their memory back

        if (!cancelled_) {
            // The two ways the first pass can end without a usable lattice that are about the DATA.
            // They are one situation from the user's side - no crystal was found in these frames - so
            // both name the cause and then give the same remedies.
            constexpr const char *next_steps =
                    " Nothing was integrated. Check the beam centre (--beam-x / --beam-y), raise "
                    "--max-spots if the frames are well populated with spots, and use --mode azint to "
                    "write the per-image analysis for inspection.";
            if (!best.result.has_value()) {
                // An indexer that FAILED is a third outcome, and none of those remedies fit it: it
                // never examined the frames, so nothing about the data or the processing options is
                // implicated and changing them cannot help. Name the failure instead.
                if (best.indexer_error)
                    throw JFJochException(JFJochExceptionCategory::IndexerError,
                                          fmt::format("Rotation indexing could not run: {}. Nothing was "
                                                      "integrated. The indexer never examined the frames, "
                                                      "so this is a failure of the machine and not of the "
                                                      "data - free the resource it needs (a GPU shared with "
                                                      "another process is the usual cause) and run again.",
                                                      *best.indexer_error));
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      std::string("Two-pass rotation indexing found no lattice.") + next_steps);
            }
            // A lattice that indexes almost none of the validation frames is not this crystal's lattice.
            // Forcing it anyway integrates every image against it; the few partials that survive yield no
            // fulls in the 3D combine, and the run then dies inside the merge complaining about resolution -
            // far from the cause. Stop here instead, while the reason is still visible. The floor is one
            // frame in six, several times below the weakest crystal that still merges. It is checked AFTER
            // the long-axis rescue, so a metric that rescue recovers is never rejected on its pre-rescue
            // score.
            if (best.score < static_cast<int>(validation.size()) / 6) {
                // Name the cell and Bravais class that was rejected. The commonest cause is a metric
                // symmetry promoted one class too far - the constrained cell then misses every
                // reflection by the small angle the constraint snapped away - and without the cell in
                // the message there is nothing to see that from.
                const auto &c = best.result->search_result.conventional.GetUnitCell();
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      fmt::format("Two-pass rotation indexing found only a lattice that "
                                                  "indexes {}/{} validation frames - too few for it to be "
                                                  "this crystal's lattice. It was {}-centred {}, "
                                                  "{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}.{}",
                                                  best.score, static_cast<int>(validation.size()),
                                                  best.result->search_result.centering,
                                                  gemmi::crystal_system_str(best.result->search_result.system),
                                                  c.a, c.b, c.c, c.alpha, c.beta, c.gamma, next_steps));
            }
            indexer->ForceRotationIndexerResult(*best.result);
            logger.Info("First-pass spot finding: {} frames in {} batches, {:.2f} s",
                        spot_cache.size(), prefetch_batches, prefetch_time_s);
            logger.Info("Two-pass rotation indexing found lattice (scheme '{}': {}/{} validation frames)",
                        best.name, best.score, static_cast<int>(validation.size()));

            // The spots in hand were already cut to the budget in force, so the measurement can only
            // shorten a budget, never lengthen one - which is why RunAllPasses hands the second pass the
            // list the first one started from rather than the list the first one settled on.
            if (config_.measure_spot_budget) {
                const int64_t budget = measure_budget(*indexer, *best.result);
                if (budget > 0 && budget < experiment_.GetMaxSpotCount()) {
                    logger.Info("Spot budget: keeping the strongest {} spots per image "
                                "(measured on {} frames; was {})", budget,
                                static_cast<int>(validation.size()), experiment_.GetMaxSpotCount());
                    experiment_.MaxSpotCount(budget);
                }
            }

            // Second pass: compare the de-novo lattice with pass 1's HERE, before integrating every image
            // with it. A markedly larger cell is the bistable supercell collapse, which was already
            // rejected in Run() - but only after the whole pass had run, so it cost an entire extra pass
            // on a dataset that ended up on pass-1's lattice anyway. Take pass-1's result now instead:
            // the same remedy, one pass earlier.
            if (!geometry_prepass && prepass_result_.has_value()) {
                // Primitive volumes, like the scheme comparison above: a centred conventional cell is an
                // exact integer multiple of its primitive one, so two settings of the SAME lattice differ
                // by that factor and comparing the conventional cells reads a mere change of setting as a
                // supercell - which then forces pass 1's setting and, with it, its lower symmetry.
                const double v1 = std::abs(prepass_result_->lattice
                        .ToPrimitive(prepass_result_->search_result.centering).CalcVolume());
                const double v2 = std::abs(best.result->lattice
                        .ToPrimitive(best.result->search_result.centering).CalcVolume());
                const bool supercell = v1 > 1.0 && v2 > 1.5 * v1;
                // Dropping to the primitive sub-cell of pass-1's centred lattice is the same
                // disagreement seen from the other side, and the volumes above cannot see it: the two
                // have the SAME primitive metric. Integrating there never predicts the reflections the
                // centring implies, and the pass then merges in P1 (measured on a C-centred monoclinic
                // crystal). Only this direction - a centring pass 2 finds where pass 1 had none is the
                // refined geometry doing its job, and the intensity test in the search confirms it.
                const bool lost_centring = prepass_result_->search_result.centering != 'P'
                                        && best.result->search_result.centering == 'P';
                if (supercell || lost_centring) {
                    logger.Info("Two-pass: that lattice ({}-centred, volume {:.0f} A^3) disagrees with "
                                "pass 1 ({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead",
                                best.result->search_result.centering, v2,
                                prepass_result_->search_result.centering, v1);
                    indexer->ForceRotationIndexerResult(*prepass_result_);
                }
            }
        }
    }

    // Main per-image loop, spread over N worker threads pulling from a shared counter. HDF5 reads
    // are serialized by the global hdf5_mutex; the analysis runs in parallel.
    std::atomic<int> next_ordinal = 0;
    std::atomic<int> finished_count = 0;
    // Images are written at their own ordinal, so a frame that fails leaves a HOLE rather than
    // shifting everything after it up. The extent of the written datasets is therefore the highest
    // ordinal reached, not the number that succeeded - counting successes makes the file one image
    // short for every failure, and the image it loses is the LAST one, not the one that failed.
    std::atomic<int> max_written_ordinal = -1;
    auto note_written = [&max_written_ordinal](int ordinal) {
        int prev = max_written_ordinal.load(std::memory_order_relaxed);
        while (prev < ordinal
               && !max_written_ordinal.compare_exchange_weak(prev, ordinal, std::memory_order_relaxed)) {
        }
    };
    std::atomic<uint64_t> total_uncompressed_bytes = 0;

    // What the integrator did to the reflections it was handed, summed over the workers as each one
    // finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is
    // done once per worker at the end rather than per image.
    std::atomic<uint64_t> bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0,
                          bragg_fallback = 0;

    // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and
    // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of
    // a fixed budget instead of everything it found - a powder ring is over-determined either way.
    constexpr size_t CALIBRATION_MAX_SPOTS = 20000;
    const size_t calibration_spots_per_image =
            std::max<size_t>(1, CALIBRATION_MAX_SPOTS / static_cast<size_t>(images_to_process));
    std::mutex calibration_spot_m;
    std::vector<SpotToSave> calibration_spot_list;

    auto azint_worker = [&]() {
        std::vector<uint8_t> decompression_buffer;
        ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_);
        ImagePreprocessorBuffer buffer(experiment_.GetPixelsNum());
        AzIntEngineCPU azint(mapping);
        AzimuthalIntegrationProfile profile(mapping);

        while (!cancelled_) {
            const int ordinal = next_ordinal.fetch_add(1);
            const int image_idx = start_image + ordinal * config_.stride;
            if (image_idx >= end_image) break;

            std::shared_ptr<JFJochReaderRawImage> img;
            try {
                img = reader_.GetRawImage(image_idx);
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                logger.Error("Failed to load image {}: {}", image_idx, e.what());
                continue;
            }
            if (!img) continue;

            DataMessage msg{};
            msg.image = img->image;
            msg.number = ordinal;
            msg.original_number = image_idx;
            if (dataset->efficiency.size() > image_idx)
                msg.image_collection_efficiency = dataset->efficiency[image_idx];
            total_uncompressed_bytes += msg.image.GetUncompressedSize();

            const auto t0 = std::chrono::steady_clock::now();
            try {
                const uint8_t *image_ptr = msg.image.GetUncompressedPtr(decompression_buffer);
                preprocessor.Analyze(buffer, image_ptr, msg.image.GetMode());
                azint.Run(buffer, profile);
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                logger.Error("Error integrating image {}: {}", image_idx, e.what());
                continue;
            }
            msg.azint_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
            msg.processing_time_s = msg.azint_time_s;
            msg.az_int_profile = profile.GetResult();
            msg.az_int_profile_count = profile.GetPixelCount();
            msg.az_int_profile_std = profile.GetStd();
            msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings());
            msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(),
                                                         config_.spot_finding.ice_ring_width_Q_recipA);
            msg.run_number = experiment_.GetRunNumber();
            msg.run_name = experiment_.GetRunName();

            plots.Add(msg, profile);
            if (writer_queue) writer_queue->Post(msg, img);
            note_written(ordinal);
            if (observer) observer->OnImageProcessed(msg);
            const int done = finished_count.fetch_add(1) + 1;
            if (observer) observer->OnProgress(done, images_to_process);
        }
    };

    auto full_worker = [&]() {
        pin_gpu();   // round-robin per worker thread; must precede engine construction
        MXAnalysisWithoutFPGA analysis(experiment_, mapping, pixel_mask_, *indexer,
                                       /*enable_fused_adaptive_gpu=*/true);
        AzimuthalIntegrationProfile profile(mapping);

        while (!cancelled_) {
            const int ordinal = next_ordinal.fetch_add(1);
            const int image_idx = start_image + ordinal * config_.stride;
            if (image_idx >= end_image) break;

            std::shared_ptr<JFJochReaderRawImage> img;
            try {
                img = reader_.GetRawImage(image_idx);
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                logger.Error("Failed to load image {}: {}", image_idx, e.what());
                continue;
            }
            if (!img) continue;

            DataMessage msg{};
            msg.image = img->image;
            msg.number = ordinal;
            msg.original_number = image_idx;
            if (dataset->efficiency.size() > image_idx)
                msg.image_collection_efficiency = dataset->efficiency[image_idx];
            total_uncompressed_bytes += msg.image.GetUncompressedSize();

            const auto t0 = std::chrono::steady_clock::now();
            try {
                analysis.Analyze(msg, profile, config_.spot_finding);
            } catch (const std::exception &e) {
                if (IsFatalResourceError(e)) throw;
                logger.Error("Error analyzing image {}: {}", image_idx, e.what());
                continue;
            }
            msg.processing_time_s = std::chrono::duration<float>(std::chrono::steady_clock::now() - t0).count();
            msg.run_number = experiment_.GetRunNumber();
            msg.run_name = experiment_.GetRunName();

            if (calibration) {
                std::lock_guard lock(calibration_spot_m);
                for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i)
                    calibration_spot_list.push_back(msg.spots[i]);
            }

            plots.Add(msg, profile);
            if (writer_queue) writer_queue->Post(msg, img);
            note_written(ordinal);
            if (observer) observer->OnImageProcessed(msg);
            const int done = finished_count.fetch_add(1) + 1;
            if (observer) observer->OnProgress(done, images_to_process);
        }

        const auto bragg = analysis.BraggCounts();
        bragg_predicted += bragg.predicted;
        bragg_starved += bragg.bkg_starved;
        bragg_starved_nb += bragg.bkg_starved_by_neighbour;
        bragg_fallback += bragg.profile_fallback;
    };

    if (observer)
        observer->OnPhase("Processing images");

    std::function<void()> worker = per_image_analysis ? std::function<void()>(full_worker)
                                                      : std::function<void()>(azint_worker);
    // How many workers the loop actually wants. Every one of them submits its own kernels to a card,
    // and a card runs out of room to accept them long before it runs out of work: measured on two
    // GPUs, the loop's own time falls to sixteen workers and then rises again, so forty-eight is
    // slower than sixteen on a 16 Mpx set. Sixteen is where that measurement put the minimum, and a
    // later one on a single card agrees - eight leaves the loop waiting on the queue rather than the
    // card. The cap is per card because that is what the queue depth belongs to, and it only applies
    // when -N was left at its default - an explicit -N is a deliberate instruction and is obeyed,
    // which is what a previous attempt at this got wrong.
    // Only the full-analysis worker drives a card; the azimuthal one preprocesses and integrates on
    // the CPU and wants every thread it can have, so the cap must not reach it.
    const int gpus = get_gpu_count();
    const int image_workers = (per_image_analysis && config_.nthreads_auto && gpus > 0)
            ? std::min(config_.nthreads, std::max(16, 16 * gpus))
            : config_.nthreads;
    if (image_workers < config_.nthreads)
        logger.Info("Per-image loop: {} of {} threads ({} GPUs) - past what a card can take queue "
                    "work from, more workers make it slower, not faster; pass -N to override",
                    image_workers, config_.nthreads, gpus);
    std::vector<std::future<void> > futures;
    futures.reserve(image_workers);
    const auto image_loop_start = std::chrono::steady_clock::now();
    for (int i = 0; i < image_workers; ++i)
        futures.push_back(std::async(std::launch::async, worker));
    for (auto &f: futures)
        f.get();
    // Everything the workers posted has to be on disk before the end message is written and the file
    // finalized, and a write that failed - out of space, above all - has to surface here rather than
    // leave a truncated file behind.
    if (writer_queue)
        writer_queue->Finish();
    // Wall time of the per-image loop alone. The per-stage means below cover only what runs inside it,
    // so without this there is nothing to compare them against and no way to see how much of a run is
    // spent outside it - on the first-pass indexing and on scaling/merging.
    result.image_loop_threads = image_workers;
    result.image_loop_time_s = std::chrono::duration<double>(
            std::chrono::steady_clock::now() - image_loop_start).count();

    result.cancelled = cancelled_;
    result.images_processed = finished_count.load();

    // What the stencil cost, as the integrator measured it: the reflections dropped whole for a
    // background ring left below six clean pixels, and - of those - the ones the NEIGHBOURS took
    // rather than the detector. The second is what a wider r1 can change, and it is what the two-pass
    // guard reads; the first is reported beside it because the difference is the detector's own floor.
    if (bragg_predicted > 0) {
        const auto frac = [&](uint64_t n) {
            return 100.0 * static_cast<double>(n) / static_cast<double>(bragg_predicted);
        };
        bkg_starved_fraction_ = static_cast<double>(bragg_starved_nb) / static_cast<double>(bragg_predicted);
        logger.Info("Integration at r1={:.1f} r2={:.1f} r3={:.2f}: of {} predicted reflections, "
                    "{:.3f}% lost their background ring - {:.3f}% to neighbouring reflections and the "
                    "rest to the detector; {:.3f}% profile-fit fallbacks",
                    experiment_.GetBraggIntegrationSettings().GetR1(),
                    experiment_.GetBraggIntegrationSettings().GetR2(),
                    experiment_.GetBraggIntegrationSettings().GetR3(),
                    bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb),
                    frac(bragg_fallback));
    }

    // Every image failing is a total failure, not a run that produced nothing: it used to be reported
    // only as per-image log lines while the process still exited 0 with no output file.
    if (!cancelled_ && images_to_process > 0 && result.images_processed == 0)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "No image could be analyzed - see the per-image errors above");

    result.mean_processing_time = plots.GetMeanProcessingTime();
    result.indexing_rate = plots.GetIndexingRate();

    // Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in
    // plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and
    // writes the PONI file; nothing here is stored in the _process.h5.
    if (calibration && !cancelled_) {
        if (observer)
            observer->OnPhase("Powder calibration");
        if (calibration_spots) {
            logger.Info("Powder calibration from {} pooled spots", calibration_spot_list.size());
            result.calibration = CalibrateFromSpots(calibration_spot_list,
                                                    experiment_.GetDiffractionGeometry(),
                                                    config_.calibrant_ring_q,
                                                    config_.calibration_refine_tilt);
        } else {
            const auto *summed = plots.GetAzIntProfileObject();
            if (!summed)
                throw JFJochException(JFJochExceptionCategory::CalibrationError,
                                      "No azimuthal profile was accumulated");
            result.calibration = CalibrateFromProfile(summed->GetResult(), mapping,
                                                      experiment_.GetDiffractionGeometry(),
                                                      config_.calibrant_ring_q,
                                                      config_.calibration_refine_tilt,
                                                      calibration_spot_list);

            // Re-bin, and fit again on a profile that was binned about the answer.
            //
            // Everything above reads a profile built with the geometry the run STARTED from, and a
            // calibration is run precisely because that geometry is in doubt. Binning is not something
            // the fit can undo: which pixels landed in which (q, phi) bin was decided when the images
            // were read. Get the centre wrong and every ring is smeared across its own sectors - past a
            // hundred pixels the fit leaves rms 4.8 px even when it is started from the exact answer.
            // Get the distance wrong and the radial sampling is simply the wrong scale, which is why a
            // run recovering 110 mm from a 250 mm header ended at 1.1 px where the same data from a
            // right header gives 0.4.
            //
            // So integrate the images once more, binned with what was just fitted. A calibration run is
            // a handful of images and the fit is worth far more than the second read.
            const auto rebinned = RebinAndRefit(*result.calibration, start_image, end_image,
                                                images_to_process, mapping, calibration_spot_list);
            if (rebinned)
                result.calibration = rebinned;
        }
    }

    // End message (also written to the file).
    EndMessage end_msg;
    // The EXTENT of the per-image datasets (see note_written), not the success count: those two
    // differ by exactly the number of frames that failed, and the writer sizes its virtual datasets
    // from this.
    end_msg.max_image_number = static_cast<uint64_t>(max_written_ordinal.load() + 1);
    end_msg.images_collected_count = result.images_processed;
    end_msg.images_sent_to_write_count = result.images_processed;
    end_msg.end_date = time_UTC(std::chrono::system_clock::now());
    end_msg.run_number = experiment_.GetRunNumber();
    end_msg.run_name = experiment_.GetRunName();
    end_msg.bkg_estimate = plots.GetBkgEstimate();
    result.spot_resolution_estimate_A = plots.GetResolutionEstimate();
    end_msg.ice_ring_score = plots.GetIceRingScoreArray();
    end_msg.ice_ring_score_mean = plots.GetIceRingScore();
    end_msg.az_int_result["dataset"] = plots.GetAzIntProfile();
    end_msg.indexing_rate = result.indexing_rate;

    if (full && !cancelled_) {
        if (const auto rot = indexer->FinalizeRotationIndexing(); rot.has_value()) {
            end_msg.rotation_lattice = rot->lattice;
            // The refined geometry, not the nominal StartMessage values. Two things about these five
            // numbers that are not obvious from their names. The beam centre is the PONI - the foot of
            // the perpendicular from the sample - so with a non-zero tilt it is NOT where the direct
            // beam lands; that point is D*tan(rot)/pixel away and is GetDirectBeam_pxl(). And the tilt
            // is not written here only for the record: IndexAndRefine puts it on outcome.experiment,
            // which is the geometry Bragg prediction and integration then run at. It is fitted with the
            // DISTANCE FROZEN and nothing checks the covariance, so the split between PONI and tilt is
            // not reproducible between passes even though their sum - the direct beam - is. Neither
            // number is meaningful on its own: report or feed back the whole set, or none of it.
            end_msg.refined_beam_center_x = rot->geom.GetBeamX_pxl();
            end_msg.refined_beam_center_y = rot->geom.GetBeamY_pxl();
            end_msg.refined_poni_rot1 = rot->geom.GetPoniRot1_rad();
            end_msg.refined_poni_rot2 = rot->geom.GetPoniRot2_rad();
            end_msg.refined_poni_rot3 = rot->geom.GetPoniRot3_rad();
            // rot3 is not carried: a rotation about the beam is an exact null of this experiment
            // (measured, smallest eigenvalue -5e-16), so the fit never moves it and reporting it would
            // suggest it had been determined.
            result.refined_detector_tilt_deg = std::array<double, 2>{
                rot->geom.GetPoniRot1_rad() * 180.0 / PI,
                rot->geom.GetPoniRot2_rad() * 180.0 / PI};
            if (rot->axis)
                end_msg.refined_rotation_axis = rot->axis->GetAxis();
            end_msg.rotation_lattice_type = LatticeMessage{
                .centering = rot->search_result.centering,
                .niggli_class = rot->search_result.niggli_class,
                .crystal_system = rot->search_result.system
            };
            result.rotation_lattice_found = true;
        }
        result.consensus_cell = indexer->GetConsensusUnitCell();
        end_msg.unit_cell = result.consensus_cell;
    }

    // Scaling and merging (full analysis only).
    if (full && !cancelled_ && result.indexing_rate.has_value() && result.indexing_rate > 0
        && (config_.run_scaling || !config_.reference_data.empty())) {
        // Scaling/merging is a long post-pass; report each sub-step as a phase so the GUI progress
        // bar reflects what is happening instead of freezing on one label.
        auto phase = [&](const std::string &p) {
            if (observer) observer->OnPhase(p);
        };
        phase("Scaling and merging");

        // Ice-ring handling (--detect-ice-rings): flag reflections sitting on a hexagonal-ice powder
        // ring. Their integrated intensity is contaminated by the strong, variable ice background, so
        // they drag the per-image scale fit. Flag them here so scaling (the per-image G fit and the
        // fulls refit) skips them, while the combine, the merge and the statistics keep them - dropping
        // them outright guts low/mid-resolution completeness on crystals that merge them fine.
        // The eleven hexagonal-ice bands cover 16-26% of the unique reflections at typical resolutions
        // whether or not the crystal has any ice, so handling ice unconditionally taxes clean data for
        // nothing. Gate the flagging on the run's own measured ice strength (1 = no ice; the strongest
        // ring over the smooth radial background, measured on the peak-excluded profile so Bragg peaks
        // cannot pose as ice). The score is logged either way, so a run that gates off still reports why.
        // Ice is looked for on TWO channels, because it arrives in two forms. Fine-grained ice is a
        // smooth powder ring and shows up in the radial profile (ice_ring_score). Ice in large
        // crystallites arrives as discrete reflections instead, leaves the profile flat, and is only
        // visible as a pile-up of found SPOTS on the ring positions - measured against the same q width
        // of ice-free flanks beside them, so the fall-off of spot density with resolution cancels. On
        // the rotation battery the two channels barely overlap: the crystals with smooth ice read
        // 2.1-2.4 / ~1.0, the ones whose spots are mostly ice read ~1.1 / 3.8-14.6, and a clean crystal
        // reads 1.04 on both. Either channel alone is blind to half the cases.
        const auto ice_score = plots.GetIceRingScore();
        const auto ice_spot_ratio = plots.GetIceRingSpotRatio();
        const float ice_min_score = experiment_.GetScalingSettings().GetIceMinScore();
        const float ice_min_spot_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio();
        const bool ice_present = (!ice_score.has_value() || *ice_score >= ice_min_score)
                              || (ice_min_spot_ratio > 0.0f && ice_spot_ratio.has_value()
                                  && *ice_spot_ratio >= ice_min_spot_ratio);
        if (experiment_.IsDetectIceRings() && !ice_present) {
            logger.Info("Ice-ring handling: measured ice score {:.2f} < {:.2f} and spot ratio {:.2f} < "
                        "{:.2f}, no ice detected - ice-ring handling skipped entirely",
                        *ice_score, ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio);
        } else if (experiment_.IsDetectIceRings()) {
            const float ice_width = config_.spot_finding.ice_ring_width_Q_recipA;
            size_t total = 0, flagged = 0;
            for (auto &outcome : indexer->GetIntegrationOutcome()) {
                for (auto &r : outcome.reflections) {
                    ++total;
                    r.on_ice_ring = IsOnIceRing(r.d, ice_width);
                    if (r.on_ice_ring)
                        ++flagged;
                }
            }
            logger.Info("Ice-ring handling: ice score {:.2f} (gate {:.2f}), spot ratio {:.2f} (gate "
                        "{:.2f}); flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); "
                        "excluded from scaling, kept for merging",
                        ice_score.value_or(NAN), ice_min_score, ice_spot_ratio.value_or(NAN),
                        ice_min_spot_ratio, flagged, total, ice_width);
        }

        // Scale the images and merge. Factored so it can run twice: first in P1 to give the
        // space-group search a merged dataset, then again in the determined space group so the
        // scaling sees symmetry equivalents and the final statistics are in the right symmetry.
        // (With a user-fixed space group it simply runs once, already in that symmetry.)
        struct ScaleMergeResult {
            std::vector<MergedReflection> merged;
            MergeStatistics statistics;
            // CC1/2 as measured before the correction surfaces (rotation) - the quantity the two-pass
            // quality guard in Run() compares the passes on. Equal to the reported CC1/2 wherever no
            // surface was fitted, which is every merge on the stills path.
            double cc_half_before_corrections = NAN;
        };
        // The reference path computes each image's G once (per-image scaling against the
        // reference); the scaling loop below is skipped, so G is stable across the two passes.
        // Smoothing it more than once would compound the correction, so do it only on the first
        // pass. The no-reference path recomputes G from scratch each pass and re-smooths correctly.

        // Rotation self-scaling + 3D combine + merge is done by the dedicated RotationScaleMerge (a single
        // allocate-once engine that recomputes partiality from the fitted mosaicity, combines the per-frame
        // partials, scales the fulls and merges; the whole hot path runs on the GPU when one is present).
        // Ingested once here (after ice flagging) and reused across both space-group passes; a fixed
        // (forced) mosaicity is handled by the recompute. The decay + absorption correction surfaces run as
        // post-scale-fulls stages (on by default, see ScalingSettings::CorrectionSurfaces); external-
        // reference scaling, the stills -B (per-image B-factor) and wedge refinement are unsupported here.
        // Re-seat the integrated reflections into a different setting of the SAME lattice, for the two
        // arms below. The change of basis is read straight off the two lattices: P[i][j] =
        // conv_real[i] . indexed_reciprocal[j] is the coefficient of the indexed cell's j-th axis in the
        // new cell's i-th axis, and indices transform with the axes, so it is also the matrix that takes
        // hkl across. Rounding it to integers - and declining when the residual is large - keeps it exact
        // even though the refined cell is not exactly LatticeSearch's Niggli cell.
        //
        // The per-image reflections and lattices are already in the _process.h5 by now - they go to file
        // as each image is processed, long before the group is known - so re-seating them here leaves the
        // stored ones in the OLD setting beside a cell in the new one. Accumulate the change of basis
        // (composed, in case more than one fires) and hand it to the writer, so the file says which setting
        // its per-image data are in and a reader can bring them across.
        gemmi::Mat33 stored_to_final;   // identity until something below re-seats
        const auto record_reindex = [&](const gemmi::Mat33 &m) {
            stored_to_final = m.multiply(stored_to_final);
            std::array<int32_t, 9> v{};
            for (int i = 0; i < 3; ++i)
                for (int j = 0; j < 3; ++j)
                    v[3 * i + j] = static_cast<int32_t>(std::lround(stored_to_final.a[i][j]));
            end_msg.reindex_matrix = v;
        };
        const auto reindex_into = [&](const LatticeSearchResult &cand) {
            const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()};
            const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
                                 end_msg.rotation_lattice->Cstar()};
            gemmi::Mat33 reindex;
            double reindex_res = 0.0;
            for (int i = 0; i < 3; ++i)
                for (int j = 0; j < 3; ++j) {
                    const double v = cv[i] * rs[j];
                    reindex.a[i][j] = std::round(v);
                    reindex_res = std::max(reindex_res, std::fabs(v - reindex.a[i][j]));
                }
            if (reindex_res >= 0.1)
                return false;
            for (auto &io : indexer->GetIntegrationOutcome()) {
                io.latt = io.latt.Multiply(reindex);
                for (auto &r : io.reflections) {
                    const double h = r.h, k = r.k, l = r.l;
                    r.h = static_cast<int32_t>(std::lround(reindex.a[0][0]*h + reindex.a[0][1]*k + reindex.a[0][2]*l));
                    r.k = static_cast<int32_t>(std::lround(reindex.a[1][0]*h + reindex.a[1][1]*k + reindex.a[1][2]*l));
                    r.l = static_cast<int32_t>(std::lround(reindex.a[2][0]*h + reindex.a[2][1]*k + reindex.a[2][2]*l));
                }
            }
            result.consensus_cell = cand.conventional.GetUnitCell();
            end_msg.unit_cell = result.consensus_cell;
            end_msg.rotation_lattice = end_msg.rotation_lattice->Multiply(reindex);
            end_msg.rotation_lattice_type = LatticeMessage{ .centering = cand.centering,
                .niggli_class = end_msg.rotation_lattice_type->niggli_class, .crystal_system = cand.system };
            record_reindex(reindex);
            return true;
        };

        // A space group the USER fixed carries its own Bravais lattice, and that is the setting its
        // reflections have to be indexed in - but the indexer answers to the metric, not to the group, and
        // LatticeSearch hands back the MOST symmetric setting the metric supports. A tetragonal-P lattice
        // therefore never comes back C-centred orthorhombic even though the same lattice has that setting,
        // and -S C222 on one merges with the C absence rule applied in the wrong frame: measured, half the
        // observations thrown away (2.14 M -> 1.08 M) and ISa 28.3 -> 24.6. Ask the character table for the
        // group's own class instead and re-seat the reflections into it.
        //
        // Keyed on the CENTRING differing, which is the whole of the harm: a fixed group of lower symmetry
        // than the lattice but with the same centring (-S P21 on an orthorhombic-P lattice, say) merges
        // correctly where it stands, and re-seating it would only permute axes that are already right.
        // A group the RUN determined is left alone - the arm above and the centring check below are its
        // handling, and a de-novo run must behave exactly as it did.
        if (user_fixed_sg_ && end_msg.rotation_lattice.has_value() && end_msg.rotation_lattice_type.has_value()) {
            const gemmi::SpaceGroup &fixed_sg = *user_fixed_sg_;
            const char indexed_centering = end_msg.rotation_lattice_type->centering;
            if (indexed_centering != fixed_sg.centring_type()) {
                // A trigonal-P group sits on a hexagonal-P lattice; every other system names its own.
                const auto want_system = fixed_sg.crystal_system() == gemmi::CrystalSystem::Trigonal
                                         && fixed_sg.centring_type() == 'P'
                                       ? gemmi::CrystalSystem::Hexagonal : fixed_sg.crystal_system();
                // Search the PRIMITIVE cell. rotation_lattice is the conventional one and its centring
                // is carried beside it, but the character table is reached through a Niggli reduction
                // that takes whatever it is handed to be primitive - so a centred conventional cell
                // reduces as a lattice with the wrong point set. This arm is entered exactly when the
                // centrings differ, which includes every centred indexed lattice, so the distinction is
                // not academic: a C-centred orthorhombic conventional cell read as primitive is itself a
                // perfectly good oP metric, and -S P212121 on one would match, reindex by the identity
                // and merge in a group whose absence rule deletes half the reflections that exist.
                // The other lattice arithmetic in this file goes through ToPrimitive for the same reason.
                const auto cand = LatticeSearchForClass(
                        end_msg.rotation_lattice->ToPrimitive(indexed_centering), want_system,
                        fixed_sg.centring_type());
                if (cand && reindex_into(*cand)) {
                    const auto &uc = *result.consensus_cell;
                    logger.Info("Reindexed the {}-centred indexed lattice into the {}-centred setting the "
                                "fixed space group {} needs: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} "
                                "beta={:.2f} gamma={:.2f}", indexed_centering, cand->centering,
                                fixed_sg.xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
                }
            }
        }

        // The lattice a pass indexed has to carry the Bravais lattice of the group it merges in. Only a
        // group the USER fixed can fail that - one the pass determined for itself was determined FROM
        // this lattice - and the re-seating above has already had its chance to find a setting that
        // does. What comes out otherwise is not a merely suboptimal answer: the absence rule is applied
        // across a frame the reflections are not in, and the statistics stop being arithmetic (measured:
        // 173.5% complete on a triclinic-P lattice merged in C2, and an undefined R_meas on an F-centred
        // cubic one merged in a trigonal-P group).
        //
        // The two-pass second pass is flagged and not adopted rather than refused: its lattice comes from
        // a de-novo re-index at the post-refined geometry, and the first pass - whose lattice did carry
        // the group - is still there to go back to. A pass with nothing to fall back on refuses instead;
        // the group is the user's assertion, and quietly determining a different one would answer a
        // question that was not asked. Name the cell that WAS indexed so the user can act on it.
        if (user_fixed_sg_ && end_msg.rotation_lattice_type.has_value() && end_msg.unit_cell.has_value()) {
            const gemmi::SpaceGroup &fixed_sg = *user_fixed_sg_;
            if (end_msg.rotation_lattice_type->centering != fixed_sg.centring_type()) {
                const auto &uc = *end_msg.unit_cell;
                if (prepass_result_) {
                    logger.Warning("Two-pass: the second pass indexed a {}-centred lattice, but the fixed "
                                   "space group {} has a {}-centred lattice - the group cannot describe "
                                   "this lattice, so this pass will not be adopted",
                                   end_msg.rotation_lattice_type->centering, fixed_sg.xhm(),
                                   fixed_sg.centring_type());
                    result.lattice_conflicts_with_fixed_sg = true;
                } else {
                    throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, fmt::format(
                            "The space group {} was fixed for this run, and its lattice is {}-centred - but "
                            "this crystal indexes as a {}-centred {} lattice (a={:.3f} b={:.3f} c={:.3f} "
                            "alpha={:.2f} beta={:.2f} gamma={:.2f}), and no setting of that lattice carries the "
                            "group's. Merging in {} would apply its absence rule to reflections that are not in "
                            "its frame, so nothing it reported would describe this crystal. Re-run without a fixed "
                            "space group to have it determined from the data, or fix one whose lattice this "
                            "crystal has.",
                            fixed_sg.xhm(), fixed_sg.centring_type(), end_msg.rotation_lattice_type->centering,
                            gemmi::crystal_system_str(end_msg.rotation_lattice_type->crystal_system),
                            uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, fixed_sg.xhm()));
                }
            }
        }

        const auto &rot_ss = experiment_.GetScalingSettings();
        const bool is_rotation = experiment_.IsRotationIndexing();   // rotation indexing -> rotation scaling/merge
        std::optional<RotationScaleMerge> rsm;
        if (is_rotation) {
            if (experiment_.GetRefineRotationWedgeInScaling()
                || rot_ss.GetRotationWedgeForScaling().has_value())
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      "Rotation scaling/merging (RotationScaleMerge) does not support "
                                      "wedge refinement");
            // A reference MTZ is allowed for rotation: it fixes the space group / cell (on the CLI) and
            // resolves the indexing ambiguity (below), but is NOT used to scale - the rotation merge stays
            // self-consistent.
            rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
                        static_cast<int>(config_.scaling_iter),
                        config_.nthreads, logger, config_.observation_dump_path);
            rsm->Ingest();
        }

        auto scale_and_merge = [&](const std::string &label, bool for_search) -> ScaleMergeResult {
            if (rsm) {
                phase("Scale/combine/merge (" + label + ")");
                // The geometry pre-pass merges only to choose the space group and to give the quality
                // guard something to judge the second pass against; its reflections are never written,
                // so the parts of the merge that only fill in an output file are skipped there.
                auto r = rsm->Run(for_search, /*full_stats=*/!geometry_prepass,
                                  /*measure_cc_before_corrections=*/true);
                result.error_model_isa = r.isa;
                result.error_model_isa_asymptotic = r.isa_asymptotic;
                if (r.resolution_fit_A) result.resolution_fit_A = r.resolution_fit_A;
                result.error_model_a = r.error_model_a;
                result.error_model_b = r.error_model_b;
                return ScaleMergeResult{std::move(r.merged), std::move(r.statistics),
                                        r.cc_half_before_corrections};
            }
            // Stills (rotation goes through RotationScaleMerge above): self-scale each image against the
            // running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even
            // with an external reference: the reference is used only to break the per-image indexing
            // ambiguity and to report CCref / inherit the R-free set, NOT as a scale anchor (scaling each
            // image against a foreign dataset injects cross-dataset systematics and is a worse reference
            // than the data's own merge). One pass: the per-image scale G is the exact one-pass solution,
            // so iterating would only rebuild the reference from the freshly-scaled (noisy) data and
            // degrade weak stills (measured CC1/2 collapse at the default 3 iters).
            {
                phase("Scaling images (" + label + ")");
                auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome());
                indexer->ScaleAllImages(merge_result);
            }
            // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
            // orientation tilt against the running merge and recompute each reflection's partiality + scale
            // correction (no re-integration). Never for the P1 search pass.
            if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) {
                phase("Partiality post-refine (" + label + ")");
                StillsPartialityRefine refiner(experiment_);
                const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads);
                logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
                // The refine fits its OWN per-image scale, with no floor, and adopts it unconditionally
                // when the image has no prior CC - which is exactly the state the guard above leaves. So
                // the guard was protecting nothing: measured on a jet dataset, of 367 images it left
                // unscaled only 8 were still unscaled in the output, and 53 reached the merge at or below
                // a fiftieth of the run median, the worst at a 4525th. Whatever produced the final scale,
                // it is the final scale that has to be credible.
                ScaleOnTheFly::RejectCollapsedScales(indexer->GetIntegrationOutcome());
            }
            const std::vector<IntegrationOutcome> &merge_input = indexer->GetIntegrationOutcome();

            phase("Merging");
            MergeOnTheFly merge_engine(experiment_);
            // For the de-novo P1 pass (for_search) drop the ice-ring reflections from the merged
            // intensities and the error model, so the space-group search sees clean data; the final
            // in-symmetry merge keeps them so completeness is not lost.
            merge_engine.ExcludeIceRings(for_search);
            if (result.consensus_cell.has_value())
                merge_engine.ReferenceCell(*result.consensus_cell);
            // Merge-consistency filter: on floods with many spurious lattices (e.g. XFEL large cells)
            // most indexed crystals do not correlate with the true structure; --min-image-cc drops the
            // crystals whose per-image CC to the reference is below the limit, so the merge keeps only
            // the coherent (real) ones. Set before the error model so the model, the merge and the
            // reported statistics are all fitted over the same set of images. The P1 search pass keeps
            // every image: its job is to find the symmetry, not to produce final intensities.
            merge_engine.FilterByImageCC(!for_search
                                         && experiment_.GetScalingSettings().GetMinCCForImage() > 0.0);
            merge_engine.RefineErrorModel(merge_input);
            if (merge_engine.ErrorModelActive())
            {
                // XDS convention, so the pair reads directly against a CORRECT.LP; the stills ISa is
                // already the whole-range 1/b, there being no asymptote tier on this path.
                const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
                logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} chi2={:.2f}",
                            em.a, em.b, em.isa, merge_engine.ErrorModelChi2());
            }
            for (size_t i = 0; i < merge_input.size(); ++i)
                merge_engine.AddImage(merge_input[i], static_cast<int64_t>(i));

            ScaleMergeResult out;
            out.merged = merge_engine.ExportReflections();

            // Automatic high-resolution cutoff (post-merge): trim the written reflections + reported
            // shells to the CC1/2 fall-off. The merge, scaling and error model above ran over the full
            // range, and the _process.h5 is written from the per-image outcomes, so no data is lost. A
            // manual --scaling-high-resolution wins; the P1 search merge (for_search) is never cut so the
            // space-group search still sees the full range. (Rotation is cut inside RotationScaleMerge.)
            const auto &cut_ss = experiment_.GetScalingSettings();
            const std::optional<double> effective_d_min = ApplyResolutionCutoff(
                out.merged, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(),
                cut_ss.GetResolutionCCTarget(), for_search, logger, &result.resolution_fit_A);

            phase("Computing statistics");
            out.statistics = merge_engine.MergeStats(out.merged, merge_input, config_.reference_data, effective_d_min);
            {   // stills: no asymptote tier, and its ISa already IS the whole-range 1/b
                const auto em = ToXdsErrorModel(merge_engine.ErrorModelA(), merge_engine.ErrorModelB());
                result.error_model_isa = em.isa;
                result.error_model_a = em.a;
                result.error_model_b = em.b;
            }
            out.cc_half_before_corrections = out.statistics.overall.cc_half;   // no surfaces on this path
            logger.Info("Merge complete ({} unique reflections, {})", out.merged.size(), label);
            return out;
        };

        // First pass: P1 when searching, or directly the user-fixed space group.
        const auto initial_sg = experiment_.GetGemmiSpaceGroup();
        auto sm = scale_and_merge(initial_sg ? initial_sg->short_name() : "P1", !initial_sg.has_value());

        // The two-pass quality guard judges the second pass against the first on THIS merge, which both
        // passes make in the same terms - P1 (or the group the user fixed), full range, no correction
        // surfaces - where their final merges need not even be in the same group. See RunAllPasses.
        {
            const auto &o = sm.statistics.overall;
            result.search_merge_completeness = o.possible_unique_reflections > 0
                    ? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0;
            result.search_merge_cc_half = sm.cc_half_before_corrections;
        }

        // Rotation two-pass geometry pre-pass: this first scale/merge has now fitted a frame-order-smoothed
        // mosaicity and written it back onto the per-frame outcomes. Capture it for the second pass's Bragg
        // prediction (so prediction uses the smoothed value, not a fresh per-image estimate) and post-refine +
        // apply the detector geometry from this integration. That is all the pre-pass is for: everything
        // the second pass consumes from it is fitted here, and the space group is left to the second pass
        // (see the search gate below).
        if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) {
            const auto &outcomes = indexer->GetIntegrationOutcome();
            prepass_mosaicity_.assign(outcomes.size(), NAN);
            for (size_t o = 0; o < outcomes.size(); ++o)
                if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg;
            if (result.consensus_cell.has_value()) {
                if (const auto rot = indexer->FinalizeRotationIndexing(); rot && rot->axis) {
                    prepass_result_ = *rot;   // full result (lattice+geom+axis), for the supercell fallback
                    PostRefineSettings prs;
                    prs.crystal_system = rot->search_result.system;
                    prs.num_threads = static_cast<int>(config_.nthreads);
                    prs.refine_geometry = true;
                    const auto pr = PostRefineRotationGeometry(outcomes, *rot->axis,
                                        experiment_.GetDiffractionGeometry(), rot->lattice, prs, logger);
                    result.post_refine = pr;
                    if (pr.ok && pr.detector_refined) {
                        // Store, don't apply: this pass continues to write its own (header-geometry) output;
                        // Run() applies the refined geometry before the second pass.
                        prepass_detector_geometry_ = std::array<float, 3>{
                            static_cast<float>(pr.beam_x_after_px),
                            static_cast<float>(pr.beam_y_after_px),
                            static_cast<float>(pr.distance_after_mm)};
                        logger.Info("Two-pass: post-refined detector geometry (distance {:.3f} mm, "
                                    "beam {:.2f},{:.2f} px) will drive the second integration pass",
                                    pr.distance_after_mm, pr.beam_x_after_px, pr.beam_y_after_px);
                    } else {
                        logger.Info("Two-pass: geometry post-refine committed no detector change "
                                    "- the second pass reproduces the first");
                    }
                    // DECISION POINT for the goniometer rotation scale. It does not ride on pr.ok - the
                    // scale is its own cross-validated fit and the crystals that have a stage fault are
                    // exactly the ones whose cell and detector steps do NOT pass, because the angle error
                    // is what their residual is made of. A calibration fault is rare (36 of 37 rotation
                    // datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset
                    // would silently damage it, so the asymmetry is deliberate: committed only when the
                    // fit is both cross-validated and outside the tolerance. A manual --rotation-scale is
                    // already on the goniometer and is left alone.
                    if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value())
                        prepass_rotation_scale_ = static_cast<float>(pr.rotation_scale);
                }
            }
        }

        // For a de-novo search, gate the P1 scaling+merge (the pass that feeds the space-group search) at
        // <I/sigma> >= 1 over thin resolution shells: noise-dominated high-res shells otherwise mislead
        // the determination, where a single artifact (a zinger the box sum integrates whole, or a
        // weak-reflection profile-fit runaway) has an astronomical resolution-normalised E and wrecks the
        // intensity-correlation / second-moment statistics. The final in-symmetry merge keeps the full
        // range (only this P1 search pass is cut); a manual --scaling-high-resolution, if coarser, wins.
        // The geometry pre-pass does not determine the space group. Its integration is at the header
        // geometry, which is the thing the second pass exists to improve on, and the symmetry evidence
        // moves with the geometry: measured on a rotation crystal, 0.33 mm of distance and 0.65 px of
        // beam centre took the screw-axis evidence from -19 to +311 nats, where merging the same
        // integration in symmetry instead of P1 moved it 2.5. Deciding on the worse integration also
        // fails in one direction only - every statistic the search weighs is a ratio of the candidate
        // merge's scatter to the subgroup's, so poorly-integrated data inflate them and REFUSE real
        // symmetry (corpus-wide the misses run ten under-calls to one over-call). So the search runs on
        // the canonical pass, at the refined geometry, and the pre-pass skips it along with the
        // in-symmetry re-merge that only exists to serve it.
        const bool search_space_group = !geometry_prepass && !initial_sg.has_value();

        double d_min_search = 0.0;
        if (search_space_group) {
            std::vector<std::pair<float, float>> rs;   // (d, I/sigma) over the full P1 merge
            rs.reserve(sm.merged.size());
            for (const auto &m : sm.merged)
                if (std::isfinite(m.I) && std::isfinite(m.sigma) && m.sigma > 0.0f &&
                    std::isfinite(m.d) && m.d > 0.0f)
                    rs.emplace_back(m.d, m.I / m.sigma);
            if (rs.size() >= 400) {
                std::sort(rs.begin(), rs.end(),
                          [](const auto &a, const auto &b) { return a.first > b.first; });   // low -> high res
                const int bins = std::clamp(static_cast<int>(rs.size() / 100), 1, 40);
                const size_t per = (rs.size() + bins - 1) / static_cast<size_t>(bins);
                for (size_t b = 0; b * per < rs.size(); ++b) {
                    const size_t lo = b * per, hi = std::min(rs.size(), lo + per);
                    double sum = 0.0;
                    for (size_t j = lo; j < hi; ++j) sum += rs[j].second;
                    if (sum / static_cast<double>(hi - lo) < 1.0) {
                        // Cut the noise-dominated high-res shells. When even the lowest-res shell fails,
                        // keep that shell alone rather than abandoning the cut. The bound is absolute
                        // while the merged I/sigma it tests saturates at the merge's own ISa (times the
                        // square root of the multiplicity), so a merge whose ISa has collapsed towards 1
                        // can fail it in every shell - and that is exactly the merge the search must not
                        // be handed whole, since the noise which drove ISa down is what the cut exists to
                        // remove. Abandoning it feeds back: more noise into the search, a weaker operator
                        // statistic, and a lost symmetry.
                        if (b > 0)
                            d_min_search = rs[lo].first;
                        else if (per >= 400)
                            d_min_search = rs[per - 1].first;
                        break;
                    }
                }
            }
            const double d_manual = experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0);
            if (rsm && d_min_search > d_manual) {
                logger.Info("Space-group search: re-scaling the P1 pass at {:.2f} A (<I/sigma> >= 1)", d_min_search);
                rsm->SetDMinLimit(d_min_search);
                sm = scale_and_merge("P1, space-group search", true);
                rsm->SetDMinLimit(d_manual > 0.0 ? std::optional<double>(d_manual) : std::nullopt);
            }
        }

        std::ostringstream stats_text;
        // True when the point group below was chosen BY THE SEARCH (rather than given by the user) and is
        // above triclinic - i.e. a promotion was made, which is what makes a later "the Laue class is
        // holohedral so there is no twin law" conclusion circular.
        bool promoted_point_group = false;
        if (search_space_group) {
            SearchSpaceGroupOptions sg_opts;
            sg_opts.nthreads = static_cast<size_t>(std::max(1, config_.nthreads));
            sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel();
            sg_opts.d_min_limit_A = std::max<double>(
                d_min_search, experiment_.GetScalingSettings().GetHighResolutionLimit_A().value_or(0.0));
            // Constrain the search to subgroups of the lattice (metric) symmetry found by rotation
            // indexing. Centering is not constrained here - it is determined from the absences.
            if (end_msg.rotation_lattice_type.has_value())
                sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system;
            // The cell the reflections are indexed on, which is what lets the search offer a candidate
            // in a setting other than gemmi's reference one: a setting names the symmetry directions by
            // axis, so a candidate the cell cannot host is refused before it is scored.
            if (result.consensus_cell.has_value())
                sg_opts.cell = gemmi::UnitCell(*result.consensus_cell);
            sg_opts.enumerate_all_settings = true;
            sg_opts.enumerate_all_rotation_sets = true;
            // The error model this merge was scaled with, so the search's "genuinely present" cut can be
            // read on the counting sigma instead of on the merged one (SearchSpaceGroupOptions::merge_isa).
            // Set per merge, not once: the two arms below are scaled separately and their ISa differ by up
            // to a factor of six. 0 (b not measurable) leaves the cut as the caller set it.
            sg_opts.merge_isa = result.error_model_isa;
            auto sg_search = SearchSpaceGroup(sm.merged, sg_opts);

            // Second opinion from a merge that keeps only well-measured observations (see
            // RotationScaleMerge::search_min_zeta). Where the two disagree the all-observation merge
            // decides. The filter earns its place because a reflection that crosses the Ewald sphere
            // near-tangentially is measured worst, and a handful of bad measurements can make a real
            // operator look like a twin law; but it discards 40-80% of the observations, which starves
            // the operator correlations. Both failure directions are real - the filtered arm has been
            // seen to lose an operator the full merge confirms, and, at geometries away from the ones
            // this comparison was first checked on, to confirm one the full merge refuses. So the two
            // are split rather than ranked: whichever arm finds the LARGER point group decides which
            // rotations are real, because a promotion is refused by the most damning statistic it can
            // be shown and only one arm has to be starved for that to happen, while the ABSENCES are
            // always judged on the all-observation merge, which is the only one holding the weak
            // reflections a systematic absence is made of. Both searches keep their vetoes, so a merohedral twin
            // stays refused in each (verified on a real twin at the three frame ranges where it
            // over-promotes: refused in both arms every time).
            if (rsm && rsm->GetSearchMinZeta() > 0.0 && !sg_search.point_group_hm.empty()) {
                const double zeta = rsm->GetSearchMinZeta();
                rsm->SetSearchMinZeta(0.0);
                auto sm_all = scale_and_merge("P1, all observations", true);
                rsm->SetSearchMinZeta(zeta);
                sg_opts.merge_isa = result.error_model_isa;   // this arm's own error model
                const auto alt = SearchSpaceGroup(sm_all.merged, sg_opts);
                // The order of the point group each arm confirmed. Read from the search, not from the
                // space group it went on to pick: Stage B leaves best_space_group unset when no
                // candidate has enough absences to be eligible, which is common on the filtered arm
                // precisely because the filter removes the weak reflections the absences live in. An
                // arm that confirmed 422 would then have scored 0 and lost to a P2 arm - a demotion,
                // reported as "taking the higher symmetry".
                const auto order_of = [](const SearchSpaceGroupResult &r) { return r.point_group_order; };
                logger.Info("Space-group search: all-observation merge -> {} (order {}), "
                            "Lorentz-filtered -> {} (order {})",
                            alt.point_group_hm.empty() ? "?" : alt.point_group_hm, order_of(alt),
                            sg_search.point_group_hm.empty() ? "?" : sg_search.point_group_hm,
                            order_of(sg_search));
                // The screw/centering determination in particular comes from the merge with all the
                // observations, because systematic absences are decided by the WEAK reflections and the
                // filter throws most of them away. Preferring the filtered arm on a tie cost four
                // crystals their screw axes (P2(1) read as P2) for exactly that reason, with the point
                // group and every intensity statistic identical.
                // Checked while both arms are still intact - below, the all-observation arm is adopted
                // by value and the two become the same result.
                // Which point group, by its symmorphic representative rather than by name - gemmi calls
                // both P321 and P312 "32", so two arms that disagree about WHICH 2-folds are real read
                // as agreeing. By the representative's full symbol rather than its NUMBER, for the same
                // reason one level down: P 1 2 1 and P 1 1 2 are both number 3 and put their 2-fold on
                // different axes. Empty when no point group was chosen.
                const auto pg_of = [](const SearchSpaceGroupResult &r) {
                    return r.point_group_representative.has_value() ? r.point_group_representative->xhm()
                                                                    : std::string();
                };
                if (order_of(alt) > 0 && order_of(alt) == order_of(sg_search)
                    && pg_of(alt) != pg_of(sg_search)) {
                    // Same order, different symmetry: the two merges disagree about WHICH operators are
                    // real, and neither is higher, so there is nothing to prefer. Say so instead of
                    // picking silently - the two imply different molecular-replacement searches, and
                    // trying both is cheap next to processing the data again.
                    const std::string a = alt.best_space_group.has_value()
                                          ? alt.best_space_group->short_name() : alt.point_group_hm;
                    const std::string b = sg_search.best_space_group.has_value()
                                          ? sg_search.best_space_group->short_name() : sg_search.point_group_hm;
                    const std::string msg = fmt::format(
                        "Space group is AMBIGUOUS between {} and {} (point groups {} and {}, the same "
                        "order): the merge of all observations and the merge of only the well-measured "
                        "ones each support one of them, and neither is the higher symmetry, so the data "
                        "do not decide. Processing continues in {} (the all-observation choice) - try "
                        "BOTH in molecular replacement, or re-run with -S {} to force the other one.",
                        a, b, alt.point_group_hm, sg_search.point_group_hm, a, b);
                    logger.Warning("{}", msg);
                    stats_text << "  !! " << msg << "\n\n";
                    result.warnings.push_back(msg);
                }
                if (order_of(alt) > order_of(sg_search))
                    logger.Info("Space-group search: all-observation merge supports {} where the "
                                "Lorentz-filtered one supports {} - taking the higher symmetry",
                                alt.point_group_hm, sg_search.point_group_hm);
                if (order_of(alt) < order_of(sg_search)) {
                    // The filtered arm found symmetry the all-observation arm did not. Discarding it
                    // was costing real point groups: over the corpus the two arms disagree on eight
                    // crystals and the higher one is right on all eight, because the filter removes
                    // the near-tangential observations that make a genuine operator look like a twin,
                    // and a promotion is refused by the most damning statistic it can be shown.
                    //
                    // But the two halves of the answer do not come from the same place. The filtered
                    // merge is deliberately starved, and ABSENCES live in the weak reflections it
                    // throws away - adopting its result whole recovers the point group and loses the
                    // screw (a P6_5 crystal comes back P6). So the point group is taken from the arm
                    // that found it and the absences are re-judged here, on every observation. That
                    // is what fixed_point_group is for: Stage A still runs and still reports what it
                    // would have refused, only the choice is overridden.
                    sg_opts.fixed_point_group = sg_search.point_group_representative;
                    auto split = SearchSpaceGroup(sm_all.merged, sg_opts);
                    sg_opts.fixed_point_group.reset();
                    logger.Info("Space-group search: Lorentz-filtered merge supports {} where the "
                                "all-observation one supports {} - taking the higher symmetry, and "
                                "judging its absences on all observations -> {}",
                                sg_search.point_group_hm, alt.point_group_hm,
                                split.best_space_group.has_value() ? split.best_space_group->xhm()
                                                                   : split.point_group_hm);
                    sg_search = std::move(split);
                } else {
                    sg_search = alt;
                }
                // The all-observation merge is the one carried forward: it is what the
                // pre-promotion twinning numbers below are measured on, and the filtered merge is the
                // second opinion reported above.
                sm = std::move(sm_all);
            }

            // Twinning evidence measured BEFORE any promotion, on the P1/subgroup merge the search was
            // given. Reported alongside the post-adoption analysis: once a higher Laue class has been
            // adopted, its own twinning test can only say "no twin law exists within this class", which
            // is circular - the promotion is what a twin would have caused. These numbers are not.
            const auto pre_promotion_twinning = AnalyzeTwinning(sm.merged, nullptr);
            result.pre_promotion_twinning = pre_promotion_twinning;
            if (pre_promotion_twinning.twinning_suspected)
                logger.Warning("Twinning indicators BEFORE the space-group decision (subgroup merge): "
                               "<|L|> = {:.3f}, second moment = {:.3f} (untwinned 0.500 / 2.00) - if the "
                               "search promoted the point group, treat the promotion with suspicion",
                               pre_promotion_twinning.mean_abs_l, pre_promotion_twinning.second_moment);
            promoted_point_group = !sg_search.point_group_hm.empty() && sg_search.point_group_hm != "1";
            if (!sg_search.refused_point_group_hm.empty())
                logger.Warning("Higher symmetry {} was confirmed by the operator correlations but refused: "
                               "{}. Processing in the lower symmetry (the recoverable direction).",
                               sg_search.refused_point_group_hm, sg_search.refused_reason);

            // Miller-index reindex under a change of basis: (h,k,l)_conv = reindex * (h,k,l)_prim.
            const auto reindex_hkl = [](auto &r, const gemmi::Mat33 &m) {
                const double h = r.h, k = r.k, l = r.l;
                r.h = static_cast<int32_t>(std::lround(m.a[0][0] * h + m.a[0][1] * k + m.a[0][2] * l));
                r.k = static_cast<int32_t>(std::lround(m.a[1][0] * h + m.a[1][1] * k + m.a[1][2] * l));
                r.l = static_cast<int32_t>(std::lround(m.a[2][0] * h + m.a[2][1] * k + m.a[2][2] * l));
            };

            // Intensity-based centred-lattice test (XDS CORRECT / POINTLESS style). When the indexer
            // demoted to a primitive/triclinic cell but the metric is pseudo-symmetric, the point-group
            // search above tested no operators (its candidate set is bounded by the indexer's lattice).
            // A genuinely centred lattice and a merely pseudo-symmetric triclinic one are indistinguishable
            // from spot POSITIONS - the 2-fold is an INTENSITY property - so decide it here, after P1
            // integration: reindex the P1-merged data to the metric candidate's conventional setting and
            // re-run the search. The I(h) vs I(Rh) correlation confirms a real centred lattice (high CC ->
            // commit the centring) and rejects a pseudo-symmetric triclinic (low CC -> stays P1, no reindex).
            std::optional<gemmi::Mat33> commit_reindex;
            std::optional<UnitCell> commit_cell;
            std::optional<CrystalLattice> commit_lattice;
            std::optional<LatticeMessage> commit_lattice_type;
            // Only a demoted (triclinic) indexing result reaches here with a primitive cell - that is both
            // the case the point-group search could not test (a triclinic holohedry offers no operators)
            // and the case where LatticeSearch's Niggli reduction is valid (it assumes a primitive input).
            const bool searched_trivial = !sg_search.best_space_group || sg_search.best_space_group->number == 1;
            const bool demoted_triclinic = end_msg.rotation_lattice_type.has_value()
                && end_msg.rotation_lattice_type->crystal_system == gemmi::CrystalSystem::Triclinic;
            if (searched_trivial && demoted_triclinic && end_msg.rotation_lattice.has_value()) {
                // The metric only generates the candidate; the intensity CC gate below makes the decision
                // (a false pseudo-symmetry is rejected by the correlation).
                const auto cand = LatticeSearch(*end_msg.rotation_lattice);
                // Exact integer reindex from the indexed setting straight to the candidate conventional
                // setting: P[i][j] = conv_real[i] . rot_reciprocal[j]. Derived from the two lattices, so it
                // holds even if the refined indexed cell is not exactly LatticeSearch's Niggli cell (whose
                // primitive_reduced->conventional reindex would otherwise be applied in the wrong setting).
                const Coord cv[3] = {cand.conventional.Vec0(), cand.conventional.Vec1(), cand.conventional.Vec2()};
                const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
                                     end_msg.rotation_lattice->Cstar()};
                gemmi::Mat33 reindex;
                double reindex_res = 0.0;
                for (int i = 0; i < 3; ++i)
                    for (int j = 0; j < 3; ++j) {
                        const double v = cv[i] * rs[j];
                        reindex.a[i][j] = std::round(v);
                        reindex_res = std::max(reindex_res, std::fabs(v - reindex.a[i][j]));
                    }
                if (cand.system != gemmi::CrystalSystem::Triclinic && reindex_res < 0.1) {
                    std::vector<MergedReflection> merged_c = sm.merged;
                    for (auto &m : merged_c)
                        reindex_hkl(m, reindex);
                    SearchSpaceGroupOptions o2 = sg_opts;
                    o2.lattice_system = cand.system;
                    o2.cell = gemmi::UnitCell(cand.conventional.GetUnitCell());
                    const auto s2 = SearchSpaceGroup(merged_c, o2);
                    if (s2.best_space_group.has_value() && s2.best_space_group->number > 1) {
                        // The centring cannot be confirmed from absences here (integrated in the primitive
                        // cell, so the centring-absent reflections do not exist) - it is metric-determined.
                        // Among the point-group-equivalent candidates the intensities cannot separate
                        // (e.g. P2 / P21 / C2), pick the one matching the metric candidate's centring.
                        gemmi::SpaceGroup chosen = *s2.best_space_group;
                        if (chosen.centring_type() != cand.centering)
                            for (const auto &alt : s2.alternatives)
                                if (alt.centring_type() == cand.centering) { chosen = alt; break; }
                        // Commit only if the adopted centring matches the metric candidate whose conventional
                        // cell we reindex into - otherwise the space group and the written cell disagree.
                        if (chosen.centring_type() == cand.centering) {
                            sg_search = s2;
                            // Swap the metric-matching candidate into the answer and leave the one it
                            // displaced among the alternatives, so the reported "or ..." list stays the
                            // set of groups the data cannot separate, each named once.
                            if (chosen.number != s2.best_space_group->number) {
                                sg_search.alternatives.push_back(*s2.best_space_group);
                                std::erase_if(sg_search.alternatives, [&](const gemmi::SpaceGroup &a) {
                                    return a.number == chosen.number;
                                });
                            }
                            sg_search.best_space_group = chosen;
                            commit_reindex = reindex;
                            commit_cell = cand.conventional.GetUnitCell();
                            commit_lattice = end_msg.rotation_lattice->Multiply(reindex);
                            commit_lattice_type = LatticeMessage{ .centering = cand.centering,
                                .niggli_class = end_msg.rotation_lattice_type->niggli_class,
                                .crystal_system = cand.system };
                            // Say so at the point of decision. The point group here is confirmed
                            // from the intensities, but the centering is not and cannot be - the
                            // reflections it extinguishes were never measured. (The report warns
                            // about it too, from the candidate's centering_untested flag, which
                            // survives into the pass that is written.)
                            // Only for a CENTRED candidate: a primitive one extinguishes nothing, so
                            // there is no untested absence class to warn about. Before the enumeration
                            // was widened this branch could only ever adopt a centred group, so the
                            // condition was implicit.
                            if (cand.centering != 'P')
                                logger.Warning("The {} centering of {} comes from the lattice METRIC, "
                                               "not from the intensities: the crystal is indexed on "
                                               "the primitive sub-cell, so the reflections a "
                                               "{}-centred lattice extinguishes were never predicted "
                                               "or integrated. The point group ({}) IS confirmed, by "
                                               "the operator correlations.",
                                               cand.centering, chosen.short_name(), cand.centering,
                                               s2.point_group_hm);
                        }
                    }
                }
            }

            // A DOUBLED DESCRIPTION of a primitive lattice. gemmi states every setting as a change of
            // basis from its reference one, and where that change has determinant < 1 the setting names
            // a cell of 1/det times the reference volume. In the monoclinic family exactly two settings
            // have it - C 1 1 2 and C 1 1 2_1 - and the reason is physical: with unique axis c the
            // C-centring vector (1/2,1/2,0) has no component along c, lies wholly in the oblique a-b
            // plane and reduces away, which is why the genuine c-unique centred settings are A/B/I 1 1 2.
            // (The centred settings of P1, and the C/F settings of the tetragonal groups, are the same
            // thing.) Such a group is never a distinct Bravais centring; it is a primitive lattice
            // written on twice the cell.
            //
            // It wins because the candidate table is scored ON the doubled cell, where half the
            // reflection positions are not reciprocal-lattice points at all and are absent BY
            // CONSTRUCTION, so the centred candidate collects a centring evidence that is tautological.
            // Refusing to enumerate those settings was measured and is worse, not better: the
            // volume-preserving alternative in the same table keeps the by-construction absences as real
            // observations and merges no better than P1. What the centred candidate gets RIGHT is
            // discarding them.
            //
            // So take the win as a SIGNAL rather than as the answer. Halve the cell, put the lattice back
            // in the primitive setting of the group's own crystal system, carry the reflections across -
            // the by-construction absences map to non-integral indices and are dropped, which is exactly
            // what the centred candidate was doing, so nothing that was merging is lost - and re-run the
            // search there. The enumeration is untouched, so no candidate that JUDGES another moves and
            // no chi2 reference changes; the doubled setting still wins its own table, and only the cell
            // its answer is expressed on is corrected.
            if (!commit_reindex && sg_search.best_space_group.has_value()
                && sg_search.best_space_group->number > 1
                && end_msg.rotation_lattice.has_value() && end_msg.rotation_lattice_type.has_value()) {
                const gemmi::SpaceGroup won = *sg_search.best_space_group;
                constexpr double den3 = static_cast<double>(gemmi::Op::DEN) * gemmi::Op::DEN * gemmi::Op::DEN;
                const double det = won.basisop().det_rot() / den3;
                const int multiple = static_cast<int>(std::lround(1.0 / det));   // 2 for C 1 1 2_1, 4 for F 1
                if (det < 0.99 && multiple >= 2 && multiple <= 4) {
                    // Reduce, then classify. The character table is reached through a Niggli reduction
                    // that takes whatever it is handed to be primitive, so the centred conventional cell
                    // has to be brought to its primitive one FIRST - asked on the doubled cell the same
                    // call answers with the same volume, merely permuted. Ask for the group's own system
                    // in 'P': the metric of a reducible lattice supports the doubled centred description
                    // too (that is how it got here), and an unrestricted search would hand it straight
                    // back.
                    const auto cand = LatticeSearchForClass(
                            end_msg.rotation_lattice->ToPrimitive(won.centring_type()),
                            won.crystal_system(), 'P');
                    // Change of basis from the doubled conventional setting straight to the reduced one,
                    // P[i][j] = conv_real[i] . doubled_reciprocal[j], as everywhere else here - except
                    // that its entries are multiples of 1/multiple rather than integers, which is the
                    // whole point of it.
                    gemmi::Mat33 reduce;
                    double reduce_res = 1.0;
                    if (cand) {
                        const Coord cv[3] = {cand->conventional.Vec0(), cand->conventional.Vec1(),
                                             cand->conventional.Vec2()};
                        const Coord rs[3] = {end_msg.rotation_lattice->Astar(), end_msg.rotation_lattice->Bstar(),
                                             end_msg.rotation_lattice->Cstar()};
                        reduce_res = 0.0;
                        for (int i = 0; i < 3; ++i)
                            for (int j = 0; j < 3; ++j) {
                                const double v = cv[i] * rs[j];
                                reduce.a[i][j] = std::round(v * multiple) / multiple;
                                reduce_res = std::max(reduce_res, std::fabs(v - reduce.a[i][j]));
                            }
                    }
                    // hkl across the same change of basis. A reflection whose image is not integral is
                    // not a reciprocal-lattice point of the reduced cell - it is one of the positions the
                    // doubled description invented - and it is dropped rather than reindexed.
                    const auto reduce_hkl = [&reduce](auto &r) {
                        const double h = r.h, k = r.k, l = r.l;
                        const double n[3] = {reduce.a[0][0]*h + reduce.a[0][1]*k + reduce.a[0][2]*l,
                                             reduce.a[1][0]*h + reduce.a[1][1]*k + reduce.a[1][2]*l,
                                             reduce.a[2][0]*h + reduce.a[2][1]*k + reduce.a[2][2]*l};
                        for (const double x : n)
                            if (std::fabs(x - std::round(x)) > 1e-3)
                                return false;
                        r.h = static_cast<int32_t>(std::lround(n[0]));
                        r.k = static_cast<int32_t>(std::lround(n[1]));
                        r.l = static_cast<int32_t>(std::lround(n[2]));
                        return true;
                    };
                    if (!cand || reduce_res >= 0.05) {
                        logger.Warning("{} describes a cell of {}x the volume its reference setting does - "
                                       "the lattice is a doubled description of a primitive one - but it "
                                       "could not be reduced ({}). Keeping the doubled cell.",
                                       won.short_name(), multiple,
                                       cand ? "the change of basis is not exact" : "no primitive character fits");
                    } else {
                        std::vector<MergedReflection> merged_p;
                        merged_p.reserve(sm.merged.size());
                        for (auto m : sm.merged)
                            if (reduce_hkl(m))
                                merged_p.push_back(m);
                        SearchSpaceGroupOptions o2 = sg_opts;
                        o2.lattice_system = cand->system;
                        o2.cell = gemmi::UnitCell(cand->conventional.GetUnitCell());
                        const auto s2 = SearchSpaceGroup(merged_p, o2);
                        // Only ever a re-expression: the reduced cell has to support the same symmetry
                        // the doubled one found, in the primitive setting the reduction produced.
                        // Anything less and the doubled answer stands - it is the one that merges.
                        const bool holds = s2.best_space_group.has_value()
                                           && s2.point_group_order >= sg_search.point_group_order
                                           && s2.best_space_group->centring_type() == cand->centering;
                        if (!holds) {
                            logger.Warning("{} describes a cell of {}x the volume its reference setting "
                                           "does, but the reduced primitive cell does not carry the same "
                                           "symmetry ({}, order {} against {}). Keeping the doubled cell.",
                                           won.short_name(), multiple,
                                           s2.best_space_group ? s2.best_space_group->short_name() : "nothing",
                                           s2.point_group_order, sg_search.point_group_order);
                        } else {
                            for (auto &io : indexer->GetIntegrationOutcome()) {
                                io.latt = io.latt.Multiply(reduce);
                                std::vector<Reflection> kept;
                                kept.reserve(io.reflections.size());
                                for (auto r : io.reflections)
                                    if (reduce_hkl(r))
                                        kept.push_back(r);
                                io.reflections = std::move(kept);
                            }
                            const size_t dropped = sm.merged.size() - merged_p.size();
                            result.consensus_cell = cand->conventional.GetUnitCell();
                            end_msg.unit_cell = result.consensus_cell;
                            end_msg.rotation_lattice = end_msg.rotation_lattice->Multiply(reduce);
                            end_msg.rotation_lattice_type = LatticeMessage{ .centering = cand->centering,
                                .niggli_class = cand->niggli_class, .crystal_system = cand->system };
                            sg_search = s2;
                            if (rsm) {
                                rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
                                            static_cast<int>(config_.scaling_iter),
                                            config_.nthreads, logger, config_.observation_dump_path);
                                rsm->Ingest();
                            }
                            const auto &uc = *result.consensus_cell;
                            logger.Info("{} names a cell of {}x the volume of its reference setting: the "
                                        "lattice is a doubled description of a primitive one, and {} of "
                                        "{} merged reflections were positions the doubled cell invented. "
                                        "Reduced to a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} "
                                        "gamma={:.2f} and searched again there: {}",
                                        won.short_name(), multiple, dropped, sm.merged.size(),
                                        uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma,
                                        s2.best_space_group->short_name());
                            // The reduction halves the cell, so its change of basis is not integral and
                            // end_msg.reindex_matrix (int32) cannot carry it. The per-image reflections
                            // already written to the _process.h5 therefore stay in the doubled setting.
                            logger.Warning("The per-image reflections in the process file are in the "
                                           "doubled setting the run integrated in; the cell, the merged "
                                           "reflections and the space group are in the reduced one. The "
                                           "change of basis between them is not integral, so it is not "
                                           "recorded in the file.");
                        }
                    }
                }
            }

            // Adopt the determined space group and re-scale/merge in it, so scaling uses symmetry
            // equivalents and the statistics come out in the right symmetry. P1 stands when nothing
            // is determined - and it has to stand as a MERGE of its own, not as whatever `sm` happens
            // to hold: everything merged up to here was a search merge, which drops the ice rings and
            // (on the Lorentz-filtered arm) most of the observations with them, is never resolution-cut
            // and fits no correction surfaces, because its intensities exist only to be searched.
            //
            // The P1 fallback is a belt, not a brace. SearchSpaceGroup cannot hand back an unset group
            // for a merge that has any reflections in it: finding no symmetry at all leaves it on the
            // identity point group, whose representative IS P1, and it returns that; finding a point
            // group but no screw or centering leaves the primitive symmorphic candidate, which has no
            // predicted absences to be inconsistent with and is therefore always eligible; and an empty
            // merge never reaches the search, because both merge engines throw on the resolution range
            // before returning one. So the fallback buys nothing today - what it buys is that the
            // promise above is kept HERE, rather than resting on eligibility gates in another file that
            // a later change could tighten without anyone noticing what else was leaning on them.
            {
                const gemmi::SpaceGroup sg = sg_search.best_space_group.value_or(gemmi::get_spacegroup_p1());
                // xhm(), not short_name(): the search can now return a non-reference setting, and
                // short_name() cannot tell P 1 2 1 from P 1 1 2.
                logger.Info("Adopting space group {} (number {})", sg.xhm(), sg.number);
                // A reindex is committed only when the intensity re-test confirmed a higher symmetry in a
                // centred setting: bring the integrated reflections + cell into that setting and re-ingest
                // the rotation merge, so the final scaling/merging folds the equivalents correctly.
                if (commit_reindex) {
                    for (auto &io : indexer->GetIntegrationOutcome()) {
                        io.latt = io.latt.Multiply(*commit_reindex);
                        for (auto &r : io.reflections)
                            reindex_hkl(r, *commit_reindex);
                    }
                    result.consensus_cell = commit_cell;
                    end_msg.unit_cell = commit_cell;
                    // Keep the lattice metadata (hence the master-file UB matrix and indexed-lattice
                    // vectors) in the same conventional setting as the cell, reflections and space group.
                    end_msg.rotation_lattice = commit_lattice;
                    end_msg.rotation_lattice_type = commit_lattice_type;
                    record_reindex(*commit_reindex);
                    if (rsm) {
                        rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
                                    static_cast<int>(config_.scaling_iter),
                                    config_.nthreads, logger, config_.observation_dump_path);
                        rsm->Ingest();
                    }
                }
                experiment_.SetSpaceGroup(sg);
                end_msg.space_group_name = sg.xhm();
                end_msg.space_group_number = static_cast<uint64_t>(sg.number);
                result.space_group = sg;
                WarnIfMetricHostsMoreSymmetry(result.consensus_cell, sg);
                phase("Re-scaling in space group " + sg.short_name());
                sm = scale_and_merge(sg.short_name(), false);
            }
            result.space_group_search = sg_search;
        } else {
            // A space group was fixed by the user; surface it so the viewer/CLI can still show it.
            result.space_group = experiment_.GetGemmiSpaceGroup();
            // ...and record it, exactly as the search arm above does with the group it determined, or
            // the _process.h5 comes out with a cell but no group and re-merging it defaults to P1.
            if (const auto &sg = experiment_.GetGemmiSpaceGroup(); sg.has_value()) {
                end_msg.space_group_name = sg->xhm();
                end_msg.space_group_number = static_cast<uint64_t>(sg->number);
            }
        }

        // Reference-based indexing-ambiguity resolution (rotation). When a reference MTZ is supplied and
        // the crystal has an indexing ambiguity (merohedral: lattice symmetry higher than the Laue group),
        // pick the reindexing that best correlates with the reference intensities and re-merge in it. The
        // twin-law reindex is metric-preserving, so only the hkl labels change (the cell is unchanged).
        // Stills resolve the ambiguity per image with ReindexAmbiguityResolver, not here.
        if (rsm && !config_.reference_data.empty() && result.consensus_cell
            && experiment_.GetSpaceGroupNumber().has_value()) {
            const gemmi::SpaceGroup &sg = experiment_.GetSpaceGroupOrP1();
            const auto choice = ChooseReindex(
                sm.merged, *result.consensus_cell, sg,
                [&](const std::vector<MergedReflection> &m) {
                    return ReferenceIntensityCC(m, config_.reference_data, sg);
                });
            if (!choice.is_identity) {
                logger.Info("Reference: resolved indexing ambiguity by reindexing to match the reference "
                            "(CC {:.3f}, vs {:.3f} unchanged)", choice.score, choice.identity_score);
                for (auto &io : indexer->GetIntegrationOutcome())
                    for (auto &r : io.reflections) {
                        const gemmi::Op::Miller h = choice.op.apply_to_hkl({{r.h, r.k, r.l}});
                        r.h = h[0]; r.k = h[1]; r.l = h[2];
                    }
                rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell,
                            static_cast<int>(config_.scaling_iter),
                            config_.nthreads, logger, config_.observation_dump_path);
                rsm->Ingest();
                phase("Re-merging in the reference indexing");
                sm = scale_and_merge(experiment_.GetGemmiSpaceGroup()->short_name(), false);
            }
        }

        const auto &twin_sg_opt = experiment_.GetGemmiSpaceGroup();
        const gemmi::SpaceGroup *twin_sg = twin_sg_opt ? &*twin_sg_opt : nullptr;
        // Not on the geometry pre-pass: the analysis goes into that pass's statistics text and its
        // written reflections, and neither survives the run. The promotion flag below is a different
        // thing - it is what the SEARCH did, the second pass reads it, and it is set either way.
        if (!geometry_prepass) {
            result.twinning = AnalyzeTwinning(sm.merged, twin_sg);
            stats_text << TwinningAnalysisToText(result.twinning) << "\n";
        }
        // Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the
        // search itself made, so the text cannot claim "no twin law exists" on its own say-so.
        result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group;

        // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the
        // crystal was mounted, which the user can change), warned about when the two nearly coincide.
        if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) {
            if (const auto gonio = experiment_.GetGoniometer()) {
                const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice,
                                                                  gonio->GetAxis());
                if (closest.has_value()) {
                    // Practical bound, not a derived one: the blind cone's half-angle is the maximum
                    // Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg
                    // case shows the loss while a 16.2 deg one is 99.7% complete.
                    constexpr double WARN_DEG = 15.0;
                    const auto [angle, order] = *closest;
                    if (angle < WARN_DEG) {
                        const std::string msg = fmt::format(
                            "The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep "
                            "never records the reflections whose reciprocal vector lies within the Bragg angle "
                            "of the spindle, and symmetry normally supplies them from an equivalent elsewhere; "
                            "it cannot here, because that axis maps the blind region onto itself. Those "
                            "reflections stay missing however long the sweep runs. The loss is confined to "
                            "that cone rather than spread over the data, so overall completeness may still "
                            "look reasonable - check it near the spindle direction. A second sweep on a "
                            "different axis, or re-mounting, recovers them.",
                            order, angle);
                        logger.Warning("{}", msg);
                        stats_text << "  !! " << msg << "\n\n";
                        result.warnings.push_back(msg);
                    } else {
                        stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at "
                                   << std::fixed << std::setprecision(1) << angle << " deg\n\n";
                    }
                }
            }
        }

        // Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds
        // the Laue symmetry the crystal can be validly indexed in several hands related by twin-law
        // operators. For an obvious merohedral case (P3/P4/P6...) users expect this; but a PSEUDO-merohedral
        // metric (e.g. a C2 crystal whose beta makes it pseudo-F-orthorhombic) is easy to miss, so surface it. The
        // reindex operator is applied directly to (h,k,l), so its triplet reads as an h,k,l transform.
        if (result.consensus_cell && twin_sg != nullptr) {
            const auto twin_ops = ReindexAmbiguityOperators(*result.consensus_cell, *twin_sg, 2.0);
            if (!twin_ops.empty()) {
                std::string laws;
                for (const auto &op : twin_ops) {
                    std::string t = op.triplet();
                    std::replace(t.begin(), t.end(), 'x', 'h');
                    std::replace(t.begin(), t.end(), 'y', 'k');
                    std::replace(t.begin(), t.end(), 'z', 'l');
                    laws += (laws.empty() ? "" : " ; ") + t;
                }
                // A model resolves it too: on rotation data at the validation step below, by R-free;
                // on serial stills only through the per-image reference computed from it at integration
                // time, which needs the cell and the space group up front - so a model on its own does
                // not always suffice there, and `ambiguity_reference` says whether it did.
                const bool will_be_resolved = ambiguity_reference
                    || (!config_.model_path.empty() && experiment_.IsRotationIndexing());
                if (!will_be_resolved) {
                    const std::string msg = fmt::format(
                        "Indexing ambiguity: this cell / space group admits alternative indexing "
                        "(reindex operator(s): {}). {} are indexed in one hand at random, and rugnux can only "
                        "break this against an external reference. WITHOUT one the merge mixes the hands and "
                        "CC1/2 is degraded - supply a reference MTZ (-z) or a model ({}) to resolve it.",
                        laws, experiment_.IsRotationIndexing() ? "Lattices" : "Serial-stills crystals",
                        experiment_.IsRotationIndexing() ? "--model" : "--model, which needs -C and -S here");
                    logger.Warning("{}", msg);
                    stats_text << "  !! " << msg << "\n\n";
                    result.warnings.push_back(msg);
                } else if (config_.reference_data.empty()) {
                    // Said before the model has been read, so it says what will be tried, not what
                    // came of it - a model that turns out to be unreadable is reported where it is
                    // read, and claiming success here would contradict it.
                    logger.Info("Indexing ambiguity present (reindex operator(s): {}); the supplied "
                                "model will be used to resolve it.", laws);
                } else {
                    logger.Info("Indexing ambiguity present (reindex operator(s): {}); resolved "
                                "against the supplied reference.", laws);
                }
            }
        }

        // Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not
        // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log, none
        // of which the geometry pre-pass produces.
        if (!geometry_prepass) {
            const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged);
            sm.statistics.wilson_b = wilson.b;
            sm.statistics.wilson_b_correlation = wilson.correlation;
            if (std::isfinite(wilson.b) && wilson.b > 0.0)
                logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)",
                            wilson.b, wilson.correlation, wilson.n_shells);

            // Diffraction anisotropy. Report-only, like the twinning and Wilson analyses above: it
            // describes the data and corrects nothing. The tensor and its resolution signature come from
            // the merged intensities; the error bar the verdict is gated on has to come from the
            // unmerged observations, because a merge has exact Laue symmetry by construction and the
            // tensor directions the symmetry forbids - the only place a dataset measures its own
            // systematic error - are identically zero in it.
            if (result.consensus_cell) {
                AnisotropyRunInfo aniso_run;
                if (sm.statistics.sweep_quality.measured && sm.statistics.sweep_quality.sweep_deg > 0.0f)
                    aniso_run.observed_rotation_deg = sm.statistics.sweep_quality.sweep_deg;
                aniso_run.dose_term_in_scale_model =
                        experiment_.GetScalingSettings().GetCorrectionSurfaces();
                aniso_run.radiation_damage_relative_b = sm.statistics.radiation_damage_delta_b;
                sm.statistics.anisotropy = AnalyzeAnisotropy(
                        sm.merged,
                        ScaledObservations(indexer->GetIntegrationOutcome(),
                                           experiment_.IsRotationIndexing(), twin_sg),
                        *result.consensus_cell, twin_sg, aniso_run);
                stats_text << AnisotropyToText(sm.statistics.anisotropy) << "\n";
            }
        }

        stats_text << sm.statistics;
        result.merge_statistics_text = stats_text.str();
        result.has_merge_statistics = true;
        result.merge_statistics = sm.statistics;
        result.cc_half_before_corrections = sm.cc_half_before_corrections;

        // Per-image form of the sweep-quality ranges, for the _process.h5: one code per image, 0 where
        // the image is in no flagged range, plus the vocabulary the codes index. Only filled when the
        // diagnostic ran, so absent datasets mean "not looked for" rather than "all clean".
        if (sm.statistics.sweep_quality.measured) {
            end_msg.sweep_quality.assign(end_msg.max_image_number, 0);
            for (const auto &r : sm.statistics.sweep_quality.ranges)
                for (int64_t i = std::max<int64_t>(0, r.first_image);
                     i <= r.last_image && i < static_cast<int64_t>(end_msg.sweep_quality.size()); ++i)
                    end_msg.sweep_quality[i] = static_cast<uint8_t>(r.reason) + 1;
            for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
                end_msg.sweep_quality_reasons.emplace_back(
                        SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
        }
        {
            // Stride rather than take the head: the merged list is ordered by hkl, so the first N
            // reflections are one corner of reciprocal space and would not show the intensity range.
            constexpr size_t MAX_POINTS = 4000;
            const size_t n = sm.merged.size();
            const size_t step = std::max<size_t>(1, n / MAX_POINTS);
            result.merged_i_sigma.reserve(std::min(n, MAX_POINTS + 1));
            for (size_t i = 0; i < n; i += step)
                if (std::isfinite(sm.merged[i].I) && std::isfinite(sm.merged[i].sigma)
                    && sm.merged[i].sigma > 0.0f)
                    result.merged_i_sigma.emplace_back(sm.merged[i].I, sm.merged[i].sigma);
        }
        result.has_reference = !config_.reference_data.empty();

        // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the
        // per-hkl hash that the merge assigned), so every dataset flags the same free reflections.
        if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) {
            const size_t matched = ApplyReferenceFreeFlags(sm.merged, experiment_.GetSpaceGroupOrP1(),
                                                           config_.reference_data);
            logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)",
                        matched, sm.merged.size());
        }

        // Radiation-damage report (rotation): the per-image scale CC-to-merge and mosaicity across the
        // sweep (the per-image scaling the rotation merge already fits, binned by frame = dose). A
        // per-image CC that falls, and/or a mosaicity that rises, with frame number is the classic
        // radiation-damage signature - a data-quality-vs-dose read complementary to the fitted decay
        // correction. The full per-image table is written to <prefix>_scaling.txt for detail.
        if (experiment_.IsRotationIndexing()) {
            const auto &outs = indexer->GetIntegrationOutcome();
            const int nf = static_cast<int>(outs.size());
            constexpr int nb = 10;
            std::array<double, nb> cc_sum{}, mos_sum{};
            std::array<int, nb> cc_n{}, mos_n{};
            for (int f = 0; f < nf; ++f) {
                const int b = nf > 1 ? std::min(nb - 1, f * nb / nf) : 0;
                if (outs[f].image_scale_cc && std::isfinite(*outs[f].image_scale_cc)) {
                    cc_sum[b] += *outs[f].image_scale_cc; cc_n[b]++;
                }
                if (outs[f].mosaicity_deg && std::isfinite(*outs[f].mosaicity_deg)) {
                    mos_sum[b] += *outs[f].mosaicity_deg; mos_n[b]++;
                }
            }
            int tot = 0;
            for (int n : cc_n) tot += n;
            if (tot >= 20) {
                std::ostringstream os;
                os << fmt::format("Radiation-damage report (per-image scale over {} frames):\n", nf);
                os << "     dose        <CC to merge>   <mosaicity deg>\n";
                for (int b = 0; b < nb; ++b) {
                    if (cc_n[b] == 0 && mos_n[b] == 0) continue;
                    os << fmt::format("     {:3d}-{:3d}%        {:>7}         {:>7}\n", b * 10, (b + 1) * 10,
                                      cc_n[b] ? fmt::format("{:.3f}", cc_sum[b] / cc_n[b]) : std::string("-"),
                                      mos_n[b] ? fmt::format("{:.3f}", mos_sum[b] / mos_n[b]) : std::string("-"));
                }
                if (cc_n[0] > 0 && cc_n[nb - 1] > 0) {
                    const double cc0 = cc_sum[0] / cc_n[0], cc9 = cc_sum[nb - 1] / cc_n[nb - 1];
                    os << fmt::format("   => per-image CC to merge {:.3f} (first 10%) -> {:.3f} (last 10%){}",
                                      cc0, cc9,
                                      cc9 < cc0 - 0.05 ? "  (falling: radiation damage / crystal decay)" : "");
                }
                // Relative B-factor across the sweep (measured before any decay correction, against the
                // low-dose start): the resolution-dependent complement to the CC/mosaicity read above.
                // Radiation damage fades the high-resolution intensity, so only a POSITIVE change is dose.
                // A curve no straight line describes gets no first->last number at all, and a batch whose
                // data could not measure one prints "-" rather than a value.
                if (!sm.statistics.radiation_damage_b_batch.empty()) {
                    if (std::isfinite(sm.statistics.radiation_damage_delta_b))
                        os << fmt::format("\n   => relative B-factor change over run = {:+.2f} A^2 (first->last){}",
                                          sm.statistics.radiation_damage_delta_b,
                                          sm.statistics.radiation_damage_delta_b > 5.0
                                              ? "  (significant dose-dependent scaling: radiation damage)" : "");
                    else
                        os << "\n   => relative B-factor: the per-batch curve below is not a trend, so no "
                              "first->last number describes it - dose does not come back, so whatever moved "
                              "here was not dose; see the sweep-quality report";
                    os << fmt::format("\n      per-batch relative-B (A^2, {:.0f} deg/batch):",
                                      sm.statistics.radiation_damage_batch_deg);
                    for (float bb : sm.statistics.radiation_damage_b_batch)
                        os << (std::isfinite(bb) ? fmt::format(" {:.1f}", bb) : std::string(" -"));
                }
                logger.Info("{}", os.str());
                result.radiation_damage_text = os.str();
            }

            // Sweep-quality report: the stretches of the sweep over which the crystal delivered much
            // less than the rest of the run, and what each one looks like. Nothing is excluded because
            // of it - the frames still carry signal, and this is a message for the beamline, not a
            // filter. Frame numbers are processed-image ordinals, inclusive, as in <prefix>_image.dat.
            const auto &sq = sm.statistics.sweep_quality;
            if (sq.measured) {
                std::ostringstream os;
                os << fmt::format("Sweep quality over {:.0f} deg (per-image scale with the incident flux, "
                                  "which varied {:.2f}x, already divided out):\n",
                                  sq.sweep_deg, sq.flux_peak_to_trough);
                if (sq.ranges.empty()) {
                    os << "     no stretch of the sweep is materially worse than the run";
                } else {
                    os << "     frames              rotation  diagnosis        severity scale    CC indexed\n";
                    for (const auto &r : sq.ranges)
                        os << fmt::format("     {:<17s} {:6.1f} deg  {:<16s}    {:5.2f} {:5.2f} {:5.2f}   {:4.0f}%\n",
                                          fmt::format("{}-{}", r.first_image, r.last_image), r.rotation_deg,
                                          SweepQualityReasonText(r.reason), r.severity, r.mean_relative_scale,
                                          r.mean_relative_cc, 100.0 * r.indexed_fraction);
                    os << "   => severity is the fraction of the run's typical diffracting power missing "
                          "over the range";
                    if (sq.modulation_peak_to_trough >= 1.05f)
                        os << fmt::format("\n   => once-per-revolution modulation of the per-image scale: "
                                          "{:.1f}x peak to trough", sq.modulation_peak_to_trough);
                }
                logger.Info("{}", os.str());
            }
        }

        // Model validation runs BEFORE the reflection files are written, because it is what settles the
        // frame they are written in: the enantiomorph, which merged intensities cannot choose, and -
        // with no reference MTZ - the alternative indexing. Both are relabelings of the same
        // measurements, so the statistics computed above are unaffected; what changes is that the file,
        // the R-factors and the maps then all describe one indexing instead of two.
        if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass) {
            phase("Validating against model");
            const auto data_sg = experiment_.GetGemmiSpaceGroup();
            // With a reference MTZ the merohedral indexing was already resolved against it (rotation
            // merge / stills scaling), so trust that; only probe indexing by R-free when model-only.
            const auto validation =
                ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path,
                                     config_.output_prefix, logger,
                                     data_sg ? &*data_sg : nullptr,
                                     /*probe_indexing_ambiguity=*/config_.reference_data.empty());
            // A model that was asked for and could not be used has to say so where anyone will see
            // it. Without this the run ends successfully with no R-free, no maps and nothing in the
            // report - indistinguishable from a run that was never given --model at all.
            result.model_validation = validation;
            if (!validation.failure_reason.empty())
                result.warnings.push_back("Model validation did not run: " + validation.failure_reason);
            if (data_sg.has_value()) {
                const gemmi::SpaceGroup *adopted =
                    AdoptModelFrame(validation, sm.merged, *data_sg,
                                    experiment_.GetScalingSettings().GetMergeFriedel(), logger);
                if (adopted->number != data_sg->number) {
                    // The written reflections and the report that describes them, but deliberately not
                    // end_msg: the _process.h5 already holds the per-image reflections as they were
                    // integrated, and a file whose recorded group did not match its own reflections
                    // would mis-merge on a later --mode scale. An enantiomorphic pair merges identically
                    // anyway - same Laue class, same absences - so nothing is lost by leaving it.
                    experiment_.SetSpaceGroup(*adopted);
                    result.space_group = *adopted;
                }
                // The unmerged export below comes from the integrated observations, not from the merge,
                // so it follows the same relabeling - otherwise the two files would describe different
                // indexings of one dataset. The enantiomorph is a label, so only the indexing moves it.
                const gemmi::Op &op = validation.indexing_op;
                if (!(op == gemmi::Op::identity()))
                    for (auto &io : indexer->GetIntegrationOutcome())
                        for (auto &r : io.reflections) {
                            const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}});
                            r.h = h[0]; r.k = h[1]; r.l = h[2];
                        }
            }
        }

        // Not on the geometry pre-pass. Pass 1 exists to choose the space group and post-refine the
        // geometry; pass 2 remakes these files seconds later at the refined geometry, and that is the
        // answer anyone reads. Writing them twice costs 0.6 s of a 15 s run - a fifth of the merged
        // mmCIF on a large crystal - for a file that is superseded before the run ends. The pass-2
        // quality guard is unaffected: has_merge_statistics is set well above this, at the merge.
        if (result.consensus_cell && write_files && config_.write_merged && !geometry_prepass) {
            phase("Writing reflections");
        const ErrorModelReport em_report{
                result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?",
                result.error_model_isa_asymptotic > 0 ? fmt::format("{:.2f}", result.error_model_isa_asymptotic)
                                                   : std::string(),
                result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(),
                result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()};
            WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics,
                             em_report, result.twinning, config_.output_prefix,
                             static_cast<size_t>(std::max(1, config_.nthreads)));
            // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS
            // comparison. The offline self-scaling result is otherwise not exposed (process.h5's
            // per-image arrays are only filled on the online per-image path). Sourced from the
            // partials, which carry the first-pass per-image scale.
            ScalingResult(indexer->GetIntegrationOutcome()).SaveToFile(config_.output_prefix);

            // P1 cross-check dataset. The group the files above are written in was chosen by the
            // search, and if that choice is wrong nothing in them says so - every statistic was
            // computed in the group that was assumed, and the only way back is to process the images
            // again. Merging the same integration once more in P1 removes that: the user can re-merge
            // it in any subgroup, run MR or refinement in it, or hand it to POINTLESS. It is a proper
            // merge - correction surfaces fitted, ice rings and near-tangential observations kept -
            // not the deliberately degraded one the search itself ran on, and it is the same file
            // `rugnux --mode scale -S P1` makes from a _process.h5. Written after the merged output,
            // so the answer of the run is already on disk in the group that was determined.
            //
            // Written on EVERY de-novo run, including one whose search concluded P1 and where the
            // file therefore repeats the merged output. Whether a file exists must not depend on
            // what the pipeline decided: a script harvesting results would otherwise have to
            // reproduce the search's decision to know whether to expect it, and a missing file
            // would not distinguish "the run chose P1" from "the run failed".
            //
            // A user-fixed space group (-S) is the one case that writes nothing, and the reason is
            // not policy: with a group fixed, prediction rejects that group's centring absences and
            // those reflections are never integrated (IndexAndRefine.cpp, `.centering = ...`). A P1
            // merge built from such a run would be missing whole centring classes - a misleading P1
            // dataset, not merely a smaller one. Do not "fix" this gap by writing the file anyway.
            //
            // Rotation only, for now. The stills merge re-fits each image's scale and each
            // reflection's partiality onto the integration outcomes, and _unmerged.mtz is written
            // from those afterwards - so on stills this merge changes a file that is the run's own
            // output (measured: same length, different bytes). The fix is to run this below the
            // unmerged export rather than here, which needs the merge lambda hoisted out of this
            // block; deferred, as rotation is what the online pipeline processes.
            if (config_.write_p1_crosscheck && search_space_group && is_rotation) {
                const auto determined_sg = experiment_.GetSpaceGroupNumber();
                const gemmi::SpaceGroup *determined = gemmi::find_spacegroup_by_number(
                        static_cast<int>(determined_sg.value_or(1)));
                // scale_and_merge reports the error model of whatever it merged; this merge is not
                // the run's answer, so put the determined group's model back afterwards.
                const double em_isa = result.error_model_isa;
                const double em_isa_asymptotic = result.error_model_isa_asymptotic;
                const double em_a = result.error_model_a;
                const double em_b = result.error_model_b;
                // Both the merge and the MTZ read the group from the experiment, so it is set for
                // the whole of it and restored after.
                experiment_.SpaceGroupNumber(1);
                const auto p1 = scale_and_merge("P1 cross-check", false);
                const std::string path = config_.output_prefix + "_P1.mtz";
                WriteMtzReflections(p1.merged, *result.consensus_cell, experiment_, path);
                experiment_.SpaceGroupNumber(determined_sg);
                result.error_model_isa = em_isa;
                result.error_model_isa_asymptotic = em_isa_asymptotic;
                result.error_model_a = em_a;
                result.error_model_b = em_b;
                if (determined_sg.value_or(1) > 1)
                    logger.Info("P1 cross-check dataset written to {} ({} unique reflections): the "
                                "same observations merged in P1 instead of {}, so a wrong space group "
                                "can be recovered from without reprocessing. It is not the result of "
                                "this run.", path, p1.merged.size(),
                                determined ? determined->short_name() : "?");
                else
                    logger.Info("P1 cross-check dataset written to {} ({} unique reflections). This "
                                "run determined P1, so it repeats the merged output; it is written "
                                "anyway, so the file is there on every run.", path, p1.merged.size());
            }
        }
    }

    // Unmerged observations. After the merge, because that is where the space group is settled and
    // the file has to be written in it, but from the integrated observations rather than the merged
    // ones: the partiality and the per-image scale are left for the reading program, which fits a
    // scale model of its own. write_files, not write_output - without an output prefix there is
    // nowhere to put it, and the name would come out as a bare "_unmerged.mtz".
    if (full && !cancelled_ && write_files && !geometry_prepass && result.consensus_cell) {
        if (config_.export_unmerged) {
            if (observer) observer->OnPhase("Writing unmerged reflections");
            const std::string path = config_.output_prefix + "_unmerged.mtz";
            WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell,
                                        experiment_, true, path);
            logger.Info("Unmerged observations written to {}", path);
        }
        if (config_.export_unmerged_partials) {
            if (observer) observer->OnPhase("Writing unmerged reflections");
            const std::string path = config_.output_prefix + "_unmerged_partials.mtz";
            WriteUnmergedMtzReflections(indexer->GetIntegrationOutcome(), *result.consensus_cell,
                                        experiment_, false, path);
            logger.Info("Unmerged partials written to {}", path);
        }
    }

    // When merging is on, the _process.h5 is skipped because the merged reflections are the wanted
    // output and the per-image file is large. If nothing indexed there are no merged reflections, so
    // the run would finish successfully having written no file at all. Write one now: the per-image
    // messages went past while there was no writer, but the end message still carries the per-image
    // scalars (spot counts, background, resolution estimate) and the dataset metadata - which is what
    // anyone asking "why did nothing index" needs, and it is small precisely because nothing indexed.
    if (!writer && full && write_files && !cancelled_ && result.indexing_rate.value_or(0.0f) <= 0.0f) {
        logger.Warning("No image indexed, so there are no merged reflections to write - writing "
                       "{}_process.h5 with the per-image analysis instead", config_.output_prefix);
        fill_writer_maps();
        writer = std::make_unique<FileWriter>(start_message, /*check_overwrite_at_start=*/true,
                                              /*trusted_path=*/true);
    }

    if (writer) {
        writer->WriteHDF5(end_msg);
        writer->Finalize();
        result.written_master_path = config_.output_prefix + "_process.h5";
    }

    if (observer)
        observer->OnPhase(cancelled_ ? "Cancelled" : "Done");

    // Total wall time for the whole run, including scaling/merging/writing above (not just the
    // per-image pass), so the reported rate reflects everything rugnux did.
    result.processing_time_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - start_time).count();
    if (result.processing_time_s > 0.0) {
        result.frame_rate_hz = static_cast<double>(result.images_processed) / result.processing_time_s;
        result.throughput_MBs = static_cast<double>(total_uncompressed_bytes) / (result.processing_time_s * 1e6);
    }

    logger.Info("{} {} images in {:.2f} s ({:.2f} Hz)", cancelled_ ? "Cancelled after" : "Processed",
                result.images_processed, result.processing_time_s, result.frame_rate_hz);

    // The geometry this pass integrated at, for the report. Taken from experiment_ rather than from the
    // caller's copy: the second pass of the rotation two-pass runs at the post-refined geometry.
    result.used_beam_x_pxl = experiment_.GetBeamX_pxl();
    result.used_beam_y_pxl = experiment_.GetBeamY_pxl();
    result.used_distance_mm = experiment_.GetDetectorDistance_mm();
    {
        const auto &g = experiment_.GetDiffractionGeometry();
        constexpr double DEG = 180.0 / PI;
        result.used_detector_tilt_deg = {g.GetPoniRot1_rad() * DEG, g.GetPoniRot2_rad() * DEG,
                                         g.GetPoniRot3_rad() * DEG};
        const auto direct = g.GetDirectBeam_pxl();
        result.used_direct_beam_x_pxl = direct.first;
        result.used_direct_beam_y_pxl = direct.second;
    }
    return result;
}