Jungfraujoch/image_analysis/geom_refinement/PostRefine.cpp

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

#include "../../common/ParallelFor.h"
#include "PostRefine.h"

#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <map>
#include <memory>
#include <numeric>
#include <string>

#include "../../common/JFJochMath.h"   // PI
#include "XtalResidual.h"              // the positional detector<->reciprocal residual, and the cell it is parameterised by
#include "LatticeReduction.h"
#include "ceres/ceres.h"
#include "ceres/rotation.h"

namespace {

// How far the joint fit may move the beam centre away from a value something else already
// believes. The bound is not about how large a real correction can be - it is the backstop for a fit
// corrupted by something every cross-validation fold shares (a second lattice, most often), which the
// relative "the held-out residual improved" gate cannot see. The absolute size of the move is what
// separates that from a genuine header correction.
//
// It is measured from whichever of the nominal centre and the run's own MEASUREMENT of the centre
// (PostRefineSettings::measured_beam_px) is nearer. Anchored on the nominal centre alone it caps the
// correction at exactly the header's own error, which is the opposite of its job: the header is most
// worth correcting when it is most wrong. Measured on a rotation crystal whose header centre is 22.1 px
// out - a value the same run's beam-centre check had already placed to +-0.30 px and then discarded -
// the old bound rejected a fit that improved the held-out positional residual nine-fold.
constexpr double BEAM_BOUND_PXL = 15.0;

// How far the joint fit may move a cell ANGLE. Only monoclinic (beta) and triclinic leave one free,
// and those are the two systems whose conditioning is weakest, so the angles need the same absolute
// backstop the lengths and the distance already have rather than being left to the solver's box. The
// largest genuine angle move anywhere on the corpus is 0.35 deg, on a monoclinic beta, measured
// against the same run's own second pass; one degree is about twice that, and still well inside the
// +-2.86 deg box the solver works in, so the box reaches everywhere the gate accepts.
constexpr double ANGLE_BOUND_DEG = 1.0;

// One integrated partial, flattened across all images. Kept as narrow as the sort and the event
// split allow: on a large cell this array is gigabytes, and the scatter and every level of the
// per-bucket sort move all of it. The goniometer angle is not stored - it is a function of the
// image number alone, and is rebuilt from it where it is needed.
struct Partial {
    int h, k, l;
    float img;
    float I, sigma;
    float obs_x, obs_y;  // observed spot centroid (pixels); NAN if the box sum found no centroid
};

// A rocking event: one reflection's intensity-weighted centroid over the frames it spans.
struct Event {
    double phi_obs;     // rad, intensity-weighted rocking centroid
    double weight;      // sqrt(sum I / sum sigma)
    int h, k, l;
};

// Distance-INDEPENDENT Ewald excitation residual with the whole crystal free. The observed rocking
// centroid phi_obs says where the reflection actually crossed the Ewald sphere, which fixes the
// ABSOLUTE size of the reciprocal lattice with the detector never entering - that is what breaks the
// distance <-> cell-scale degeneracy the spot positions alone leave open, and it is why the two can be
// fitted together. The parameter blocks are XtalResidual's crystal half (rotation axis, orientation,
// cell lengths, cell angles), so one problem can share them between the two residuals.
// On the Ewald sphere <=> |p|^2 + 2 p_z/lambda == 0.
struct JointExcitationResidual {
    JointExcitationResidual(double lambda, double angle_rad, double weight,
                            int h, int k, int l, gemmi::CrystalSystem symmetry)
        : inv_lambda(1.0 / lambda), angle_rad(angle_rad), weight(weight),
          h(h), k(k), l(l), symmetry(symmetry) {}
    template<typename T>
    bool operator()(const T *const axis, const T *const p0, const T *const p1, const T *const p2,
                    T *residual) const {
        Eigen::Matrix<T, 3, 1> bxc, cxa, axb;
        T invV;
        XtalResidual::ReciprocalBasis(p1, p2, symmetry, bxc, cxa, axb, invV);
        const Eigen::Matrix<T, 3, 1> unrot = (bxc * T(h) + cxa * T(k) + axb * T(l)) * invV;
        const T recip_unrot[3] = {unrot[0], unrot[1], unrot[2]};
        T p_ref[3];
        const AngleAxisRotator<T> rot_p0(p0);
        rot_p0.Rotate(recip_unrot, p_ref);
        const T aa[3] = {T(-angle_rad) * axis[0], T(-angle_rad) * axis[1], T(-angle_rad) * axis[2]};
        T p_lab[3];
        ceres::AngleAxisRotatePoint(aa, p_ref, p_lab);
        const T zeta = p_lab[0] * p_lab[0] + p_lab[1] * p_lab[1] + p_lab[2] * p_lab[2]
                       + T(2.0) * p_lab[2] * T(inv_lambda);
        residual[0] = T(weight) * zeta * T(0.5) / T(inv_lambda);
        return true;
    }
    const double inv_lambda, angle_rad, weight;
    const double h, k, l;
    const gemmi::CrystalSystem symmetry;
};

// GONIOMETER ROTATION SCALE k: the same Ewald excitation residual, but with the crystal and the axis
// DIRECTION already committed by the joint fit, so the single free quantity is how far the stage actually
// turned per unit of commanded angle. Two differences from that fit matter:
//   * the angle is measured from the CENTRE of the sweep, not from the goniometer's zero. The reference
//     orientation is the one rotation indexing fitted against the commanded angles, so it has already
//     absorbed the MEAN angle error; only the part that varies across the sweep is left to fit. Scaling the
//     absolute angle instead - which is what reading k off the length of the fitted axis vector does - asks
//     the fit to also produce a constant offset it has no parameter for, and the least-squares compromise
//     shrinks k towards 1 by var(phi) / (var(phi) + phi_centre^2): exactly a factor of four for the common
//     case of a sweep starting at zero.
//   * e_mid is the committed reciprocal vector already turned to the sweep centre, so nothing but k is
//     free.
struct RotationScaleResidual {
    RotationScaleResidual(double lambda, double dangle_rad, const double u[3], const double e_mid[3])
        : inv_lambda(1.0 / lambda), dangle_rad(dangle_rad),
          ux(u[0]), uy(u[1]), uz(u[2]), ex(e_mid[0]), ey(e_mid[1]), ez(e_mid[2]) {}
    template<typename T>
    bool operator()(const T *const k, T *residual) const {
        const T a = T(-dangle_rad) * k[0];
        const T aa[3] = {a * T(ux), a * T(uy), a * T(uz)};
        const T p_ref[3] = {T(ex), T(ey), T(ez)};
        T p_lab[3];
        ceres::AngleAxisRotatePoint(aa, p_ref, p_lab);
        const T zeta = p_lab[0] * p_lab[0] + p_lab[1] * p_lab[1] + p_lab[2] * p_lab[2]
                       + T(2.0) * p_lab[2] * T(inv_lambda);
        residual[0] = zeta * T(0.5) / T(inv_lambda);
        return true;
    }
    const double inv_lambda, dangle_rad, ux, uy, uz, ex, ey, ez;
};

}  // namespace

PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome> &outcomes,
                                            const GoniometerAxis &axis,
                                            const DiffractionGeometry &nominal_geom,
                                            const CrystalLattice &reference_latt,
                                            const PostRefineSettings &settings,
                                            Logger &logger) {
    PostRefineResult result;
    result.geom = nominal_geom;
    result.cell = reference_latt.GetUnitCell();
    result.distance_before_mm = nominal_geom.GetDetectorDistance_mm();
    result.distance_after_mm = nominal_geom.GetDetectorDistance_mm();
    try {
        const double wedge_half = axis.GetWedge_deg() / 2.0;
        const double lambda = nominal_geom.GetWavelength_A();
        const Coord ax = axis.GetAxis();

        // The goniometer angle is a function of the image number alone, so it is recomputed where it
        // is used rather than carried through the array below: eight bytes per partial cost more in
        // the fill, the scatter and every level of the sort than the multiply-add that rebuilds them.
        const auto angle_rad = [&](float img) { return (axis.GetAngle_deg(img) + wedge_half) * PI / 180.0; };

        // Count first, then fill. Growing one vector by push_back over tens of millions of
        // reflections copies the whole thing every time it doubles - several gigabytes of pure
        // copying - and the counts are cheap to take. Each outcome then owns a slice, so the fill
        // runs on all threads and lands in the order the serial loop produced. The h range comes out
        // of the same sweep: the sort below buckets by h and needs to know how many buckets that is,
        // and this pass already reads every reflection.
        const size_t nthreads = std::max(1, settings.num_threads);
        const int n_out = static_cast<int>(outcomes.size());
        std::vector<size_t> pts_offset(n_out + 1, 0);
        std::vector<int> h_lo_of(n_out), h_hi_of(n_out);
        ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
            for (int o = lo; o < hi; o++) {
                size_t keep = 0;
                int lmin = std::numeric_limits<int>::max(), lmax = std::numeric_limits<int>::min();
                for (const auto &r : outcomes[o].reflections)
                    if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
                        keep++;
                        lmin = std::min(lmin, r.h);
                        lmax = std::max(lmax, r.h);
                    }
                pts_offset[o + 1] = keep;
                h_lo_of[o] = lmin;
                h_hi_of[o] = lmax;
            }
        });
        int h_lo = std::numeric_limits<int>::max(), h_hi = std::numeric_limits<int>::min();
        for (int o = 0; o < n_out; o++) {
            pts_offset[o + 1] += pts_offset[o];
            h_lo = std::min(h_lo, h_lo_of[o]);
            h_hi = std::max(h_hi, h_hi_of[o]);
        }
        const size_t n_pts = pts_offset[n_out];
        const int H = (h_lo <= h_hi) ? (h_hi - h_lo + 1) : 1;

        // A vector of n partials VALUE-initialises them: on a large cell that is gigabytes of zeroing
        // on one thread, and it is that one thread which first touches every page - which on a
        // multi-socket machine leaves the whole array on its node, so every pass that follows runs at
        // one node's memory bandwidth. new[] leaves the partials untouched, so the parallel fill is
        // the first touch and each page lands on the node of the thread that filled it.
        std::unique_ptr<Partial[]> pts(new Partial[n_pts]);
        // The bucket histogram the sort needs is taken here rather than in a pass of its own, since
        // the fill already has h in hand. Its chunks are the outcome chunks ParallelChunks makes, so
        // the scatter below has to be split the same way.
        const int nt = static_cast<int>(std::clamp<size_t>(nthreads, 1, std::max(1, n_out)));
        const int chunk = (n_out + nt - 1) / nt;
        std::vector<std::vector<int32_t>> hist(nt, std::vector<int32_t>(H, 0));
        ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
            std::vector<int32_t> &h_count = hist[lo / chunk];
            for (int o = lo; o < hi; o++) {
                size_t at = pts_offset[o];
                for (const auto &r : outcomes[o].reflections) {
                    if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) continue;
                    const float ox = std::isfinite(r.observed_x) ? r.observed_x : NAN;
                    const float oy = std::isfinite(r.observed_y) ? r.observed_y : NAN;
                    pts[at++] = Partial{r.h, r.k, r.l, r.image_number, r.I, r.sigma, ox, oy};
                    h_count[r.h - h_lo]++;
                }
            }
        });
        logger.Info("Post-refine: {} partials gathered", n_pts);
        if (n_pts < static_cast<size_t>(settings.min_events)) return result;

        // Where each bucket starts, and the buckets largest first: the sort lays the array out this
        // way and the event split below walks the same buckets.
        std::vector<int32_t> bstart(H + 1, 0);
        std::vector<int> order(H);

        // Bucket by h, then sort the buckets. h is the leading key, so the sorted array is the
        // buckets laid end to end, and each bucket sorts on its own thread. Sorting the whole thing
        // in one pass moved every partial through every level of a comparison sort, on one thread,
        // over tens of millions of reflections.
        {
            // Not a total order: two partials of one reflection on one image still tie, as they did
            // before this was bucketed. What makes the result reproducible is the scatter below
            // rather than the comparator - the prefix lays each bucket out chunk by chunk and a
            // chunk is a contiguous span of the gathered order, so a bucket reaches std::sort in
            // global gather order whatever the thread count. Ties therefore resolve the same way on
            // every run and at every -N; they are simply not resolved by rank.
            const auto part_less = [](const Partial &a, const Partial &b) {
                if (a.h != b.h) return a.h < b.h;
                if (a.k != b.k) return a.k < b.k;
                if (a.l != b.l) return a.l < b.l;
                return a.img < b.img;
            };
            int32_t acc = 0;
            for (int b = 0; b < H; ++b) {
                bstart[b] = acc;
                for (int t = 0; t < nt; ++t) { const int32_t c = hist[t][b]; hist[t][b] = acc; acc += c; }
            }
            bstart[H] = acc;

            std::unique_ptr<Partial[]> sorted(new Partial[n_pts]);
            ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
                std::vector<int32_t> fill = hist[lo / chunk];
                for (size_t i = pts_offset[lo]; i < pts_offset[hi]; ++i)
                    sorted[fill[pts[i].h - h_lo]++] = pts[i];
            });
            std::iota(order.begin(), order.end(), 0);
            std::sort(order.begin(), order.end(),
                      [&](int a, int b) { return (bstart[a + 1] - bstart[a]) > (bstart[b + 1] - bstart[b]); });
            ParallelFor(H, nthreads, [&](int oi) {
                const int b = order[oi];
                std::sort(sorted.get() + bstart[b], sorted.get() + bstart[b + 1], part_less);
            });
            pts.swap(sorted);
        }

        // Split into rocking events (same raw hkl, adjacent frames). Only >=2-frame events carry an
        // unbiased phi_obs (a single-frame centroid is just the frame centre).
        const float max_frame_gap = RockingEventFrameGap(axis.GetWedge_deg());
        const auto run_end = [&](size_t i, size_t end) {
            size_t j = i + 1;
            while (j < end && pts[j].h == pts[i].h && pts[j].k == pts[i].k && pts[j].l == pts[i].l
                   && pts[j].img - pts[j - 1].img <= max_frame_gap)
                ++j;
            return j;
        };
        // The event the partials [i, j) make, or false where their intensities cannot place a
        // centroid. Counting the events and writing them both walk the buckets, and both build the
        // event this way.
        const auto make_event = [&](size_t i, size_t j, Event &out) {
            double sumI = 0, sumIphi = 0, sumSig = 0;
            for (size_t m = i; m < j; ++m) {
                const double Ipos = std::max(0.0, static_cast<double>(pts[m].I));
                sumI += Ipos; sumIphi += Ipos * angle_rad(pts[m].img); sumSig += pts[m].sigma;
            }
            if (!(sumI > 0.0 && sumSig > 0.0)) return false;
            out = Event{sumIphi / sumI, std::sqrt(sumI / sumSig), pts[i].h, pts[i].k, pts[i].l};
            return true;
        };
        // An event never crosses an h boundary - h is the leading sort key - so the buckets can be
        // walked independently, and laying their events out in bucket order gives exactly the order
        // the serial walk produced. Counting first also sizes the array in one go, in place of a
        // push_back that grew a gigabyte by doubling.
        std::vector<int32_t> ev_count(H, 0);
        std::vector<size_t> ev_frames(H, 0);
        ParallelFor(H, nthreads, [&](int oi) {
            const int b = order[oi];
            const size_t end = bstart[b + 1];
            int c = 0;
            Event ev;
            for (size_t i = bstart[b]; i < end; ) {
                const size_t j = run_end(i, end);
                if (j - i >= 2 && make_event(i, j, ev)) ++c;
                i = j;
            }
            ev_count[b] = c;
        });
        std::vector<int32_t> ev_start(H + 1, 0);
        for (int b = 0; b < H; ++b) ev_start[b + 1] = ev_start[b] + ev_count[b];
        const size_t n_events = ev_start[H];
        std::unique_ptr<Event[]> events(new Event[n_events]);
        ParallelFor(H, nthreads, [&](int oi) {
            const int b = order[oi];
            const size_t end = bstart[b + 1];
            int at = ev_start[b];
            size_t frames = 0;
            for (size_t i = bstart[b]; i < end; ) {
                const size_t j = run_end(i, end);
                if (j - i >= 2 && make_event(i, j, events[at])) { frames += j - i; ++at; }
                i = j;
            }
            ev_frames[b] = frames;
        });
        size_t event_frames = 0;
        for (int b = 0; b < H; ++b) event_frames += ev_frames[b];
        // Frames per event is the phi_obs sampling: near 2 the reflections barely rock, so the angle
        // this refinement is fitted to is under-determined. It is a geometry count, so unlike an
        // intensity-weighted width it cannot be inflated by noise.
        logger.Info("Post-refine: {} multi-frame rocking events ({:.1f} frames per event)", n_events,
                    n_events == 0 ? 0.0 : static_cast<double>(event_frames) / n_events);
        if (static_cast<int>(n_events) < settings.min_events) return result;


        // The rotation-scale fit further down is a single scalar whose whole point is how the residual
        // varies ALONG the sweep, so it keeps every event. The cap below ranks by I/sigma, and on the
        // crystals that have a stage fault the strong events sit in the middle of the sweep - the part
        // that still indexes - so a capped set would leave the ends unrepresented in exactly the fit that
        // has to see them. Both sets come out of the one array by selecting on indices instead: with
        // the weights in the same places nth_element takes the same decisions it would take on the
        // events themselves, so the selection is the same one in the same order and the whole list no
        // longer has to be duplicated to survive it.
        constexpr size_t MAX_EVENTS = 20000;
        std::vector<int32_t> selected(n_events);
        std::iota(selected.begin(), selected.end(), 0);
        if (selected.size() > MAX_EVENTS) {
            std::nth_element(selected.begin(), selected.begin() + MAX_EVENTS, selected.end(),
                             [&](int32_t a, int32_t b) { return events[a].weight > events[b].weight; });
            selected.resize(MAX_EVENTS);
        }

        // ---- GEOMETRY REFINEMENT: ONE JOINT fit of the crystal (orientation, cell, rotation axis) and
        // the detector (distance, beam centre) against BOTH residuals at once:
        //   * the positional detector<->reciprocal residual at each partial's observed spot, and
        //   * the distance-INDEPENDENT Ewald excitation residual at each rocking centroid phi_obs.
        // It replaces a two-step fit that scaled the whole cell by ONE scalar against phi_obs and then
        // read the distance off that scaled cell. The two were separated because the positional residual
        // is degenerate with the cell scale - which is true of the positions ALONE, and is exactly what
        // the excitation residual breaks, so the degeneracy that motivated the split is already resolved
        // inside the same problem. Splitting it cost accuracy twice over: the first pass frees the whole
        // lattice against a frozen distance, so the distortion it absorbs is ANISOTROPIC and no single
        // scale can undo it; and whatever bias is left in that scale goes straight into the distance,
        // which is only ever determined relative to the cell. Measured over three wavelengths of one
        // crystal, that scale's bias changed SIGN with the wavelength and the distance error followed it
        // with a slope of 1.3 and a correlation of 0.99, while the rocking centroids on their own fixed
        // the cell volume to 0.06 %.
        // Detector tilt is held fixed (gauge-coupled to the orientation on a single crystal).
        // Committed only if it lowers a HELD-OUT (deterministic split-half) residual and the move stays
        // within the bounds below - otherwise the geometry is left at nominal ("quit when things go wrong").
        if (settings.refine_geometry) {
            const gemmi::CrystalSystem sys =
                (settings.crystal_system == gemmi::CrystalSystem::Trigonal) ? gemmi::CrystalSystem::Hexagonal
                                                                            : settings.crystal_system;
            const double ax0[3] = {ax.x, ax.y, ax.z};
            const double lambda_l = lambda;
            const double rot3 = nominal_geom.GetPoniRot3_rad();
            // Same for every observation, so taken once here rather than per residual.
            const double cos_rot3 = std::cos(rot3), sin_rot3 = std::sin(rot3);
            const DetectorOrientation orientation = nominal_geom.GetOrientation();
            const double pixel_mm = nominal_geom.GetPixelSize_mm();
            double det_rot[2] = {nominal_geom.GetPoniRot1_rad(), nominal_geom.GetPoniRot2_rad()};
            const UnitCell r0 = reference_latt.GetUnitCell();
            const double beam_x0 = nominal_geom.GetBeamX_pxl(), beam_y0 = nominal_geom.GetBeamY_pxl();
            const double dist0 = nominal_geom.GetDetectorDistance_mm();

            // Deterministic split of the reflections into a fit half and a held-out half. Avalanche-mix the
            // hkl hash so the split bit is decorrelated from the LSB - a plain h+k+l parity collides with the
            // lattice centering condition (e.g. an I-centred lattice has h+k+l even for EVERY present
            // reflection, so a parity split would leave the validation half empty). A reflection's rocking
            // event and its spot positions carry the same hkl, so both residual families split together.
            auto is_val = [](int h, int k, int l) {
                unsigned u = static_cast<unsigned>(h) * 2654435761u + static_cast<unsigned>(k) * 2246822519u
                             + static_cast<unsigned>(l) * 3266489917u;
                u ^= u >> 15; u *= 2246822519u; u ^= u >> 13;
                return (u & 1u) != 0u;
            };
            enum Subset { FIT, VAL, ALL };
            auto in = [&](int h, int k, int l, Subset s) {
                return s == ALL || (is_val(h, k, l) == (s == VAL)); };

            // The crystal as XtalResidual's parameter blocks - orientation (angle-axis), cell lengths,
            // cell angles - in the per-system parameterisation that residual's B matrix reads. Seeded from
            // the lattice rotation indexing settled on, and free from here on: all of it, not one scale.
            double p0[3] = {0, 0, 0}, p1[3] = {0, 0, 0}, p2[3] = {0, 0, 0};
            double beta = r0.beta;
            switch (sys) {
                case gemmi::CrystalSystem::Tetragonal:
                    LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1);
                    p1[0] = (p1[0] + p1[1]) / 2.0; break;
                case gemmi::CrystalSystem::Cubic:
                    LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1);
                    p1[0] = (p1[0] + p1[1] + p1[2]) / 3.0; break;
                case gemmi::CrystalSystem::Hexagonal:
                    LatticeToRodriguesAndLengths_Hex(reference_latt, p0, p1); break;
                case gemmi::CrystalSystem::Monoclinic:
                    LatticeToRodriguesLengthsBeta_Mono(reference_latt, p0, p1, beta);
                    p2[0] = beta; break;
                case gemmi::CrystalSystem::Orthorhombic:
                    LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1); break;
                default:
                    LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1);
                    p2[0] = r0.alpha * PI / 180.0; p2[1] = r0.beta * PI / 180.0; p2[2] = r0.gamma * PI / 180.0; break;
            }
            const double p0_0[3] = {p0[0], p0[1], p0[2]};
            const double p1_0[3] = {p1[0], p1[1], p1[2]};
            const double p2_0[3] = {p2[0], p2[1], p2[2]};
            // Only monoclinic (beta) and triclinic read p2; for every other system it is a block the
            // residual never touches, so it is held constant rather than left with a zero Jacobian column.
            const bool p2_free = (sys == gemmi::CrystalSystem::Monoclinic
                                  || sys == gemmi::CrystalSystem::Triclinic);

            // The observed spot positions, one entry per partial that produced a centroid. Count first,
            // then fill, exactly as the partial gather above does and for the same reason: this walks the
            // same tens of millions of partials, and a pointer vector grown by push_back copies itself
            // every time it doubles. new[] rather than a sized vector so the array is not zeroed on one
            // thread before the parallel fill overwrites it. The fill lands in the order the serial loop
            // produced, so the selection below sees the same sequence it always did.
            const int n_obs_chunks = static_cast<int>(std::clamp<size_t>(nthreads, 1,
                                                                         std::max<size_t>(1, n_pts)));
            const size_t obs_chunk = (n_pts + n_obs_chunks - 1) / n_obs_chunks;
            std::vector<size_t> obs_offset(n_obs_chunks + 1, 0);
            const auto keep_obs = [&](size_t i) {
                return std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y);
            };
            ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
                size_t keep = 0;
                for (int i = lo; i < hi; ++i) if (keep_obs(i)) keep++;
                obs_offset[static_cast<size_t>(lo) / obs_chunk + 1] = keep;
            });
            for (int c = 0; c < n_obs_chunks; ++c) obs_offset[c + 1] += obs_offset[c];
            size_t n_obs = obs_offset[n_obs_chunks];
            std::unique_ptr<const Partial *[]> obs(new const Partial *[n_obs]);
            ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
                size_t at = obs_offset[static_cast<size_t>(lo) / obs_chunk];
                for (int i = lo; i < hi; ++i) if (keep_obs(i)) obs[at++] = &pts[i];
            });
            constexpr size_t MAX_OBS = 20000;
            if (n_obs > MAX_OBS) {
                std::nth_element(obs.get(), obs.get() + MAX_OBS, obs.get() + n_obs,
                                 [](const Partial *a, const Partial *b) {
                                     return a->I / std::max(1e-9, static_cast<double>(a->sigma))
                                            > b->I / std::max(1e-9, static_cast<double>(b->sigma)); });
                n_obs = MAX_OBS;
            }
            result.obs_used = static_cast<int>(n_obs);

            // ===== The joint fit =====
            // Cost of a whole geometry over one subset, the two residual families kept apart so the log
            // can say which of them moved. Each is a mean per residual VALUE, and the number they are
            // summed over below is the same weighting the solver itself applies.
            auto joint_cost = [&](Subset s, const double bm[2], const double ds[1], const double rv[3],
                                  const double q0[3], const double q1[3], const double q2[3],
                                  double &pos_out, double &exc_out) {
                double cp = 0.0, ce = 0.0;
                size_t np = 0, ne = 0;
                for (size_t oi = 0; oi < n_obs; ++oi) {
                    const Partial *pp = obs[oi];
                    if (!in(pp->h, pp->k, pp->l, s)) continue;
                    XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3,
                                   angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation);
                    double res[3] = {0, 0, 0};
                    r(bm, ds, det_rot, rv, q0, q1, q2, res);
                    cp += res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
                    np += 3;
                }
                for (const int32_t i : selected) {
                    const Event &ev = events[i];
                    if (!in(ev.h, ev.k, ev.l, s)) continue;
                    JointExcitationResidual r(lambda_l, ev.phi_obs, settings.excitation_weight,
                                              ev.h, ev.k, ev.l, sys);
                    double res = 0.0;
                    r(rv, q0, q1, q2, &res);
                    ce += res * res;
                    ++ne;
                }
                pos_out = np ? cp / np : 0.0;
                exc_out = ne ? ce / ne : 0.0;
                return (np + ne) ? (cp + ce) / static_cast<double>(np + ne) : 0.0;
            };

            // One problem, both residual families, every block seeded at nominal. The detector tilt is
            // declared and held constant so that a non-zero rot1/rot2 still acts on the observed side.
            auto solve_joint = [&](Subset s, double bm[2], double ds[1], double rv[3],
                                   double q0[3], double q1[3], double q2[3]) {
                bm[0] = beam_x0; bm[1] = beam_y0; ds[0] = dist0;
                for (int j = 0; j < 3; ++j) {
                    rv[j] = ax0[j]; q0[j] = p0_0[j]; q1[j] = p1_0[j]; q2[j] = p2_0[j];
                }
                ceres::Problem p;
                for (size_t oi = 0; oi < n_obs; ++oi) {
                    const Partial *pp = obs[oi];
                    if (!in(pp->h, pp->k, pp->l, s)) continue;
                    p.AddResidualBlock(new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
                        new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3,
                                         angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation)),
                        new ceres::CauchyLoss(0.02), bm, ds, det_rot, rv, q0, q1, q2);
                }
                for (const int32_t i : selected) {
                    const Event &ev = events[i];
                    if (!in(ev.h, ev.k, ev.l, s)) continue;
                    p.AddResidualBlock(new ceres::AutoDiffCostFunction<JointExcitationResidual, 1, 3, 3, 3, 3>(
                        new JointExcitationResidual(lambda_l, ev.phi_obs, settings.excitation_weight,
                                                    ev.h, ev.k, ev.l, sys)),
                        new ceres::CauchyLoss(0.02), rv, q0, q1, q2);
                }
                if (p.NumResidualBlocks() == 0) return false;
                p.SetParameterBlockConstant(det_rot);
                if (!p2_free) p.SetParameterBlockConstant(q2);
                p.SetParameterLowerBound(ds, 0, dist0 * 0.95); p.SetParameterUpperBound(ds, 0, dist0 * 1.05);
                // The box has to reach everywhere the gate below would accept, or the gate is never the
                // thing that decides: a fit pinned at a box face lands exactly ON the bound and is then
                // refused for being there.
                for (int j = 0; j < 2; ++j) {
                    const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : bm[j];
                    p.SetParameterLowerBound(bm, j, std::min(bm[j], m) - BEAM_BOUND_PXL);
                    p.SetParameterUpperBound(bm, j, std::max(bm[j], m) + BEAM_BOUND_PXL);
                }
                for (int j = 0; j < 3; ++j) {
                    p.SetParameterLowerBound(rv, j, ax0[j] - 0.05);
                    p.SetParameterUpperBound(rv, j, ax0[j] + 0.05);
                    p.SetParameterLowerBound(q1, j, 0.97 * p1_0[j]);
                    p.SetParameterUpperBound(q1, j, 1.03 * p1_0[j]);
                }
                if (p2_free)
                    for (int j = 0; j < 3; ++j) {
                        p.SetParameterLowerBound(q2, j, p2_0[j] - 0.05);
                        p.SetParameterUpperBound(q2, j, p2_0[j] + 0.05);
                    }
                ceres::Solver::Options o; o.linear_solver_type = ceres::DENSE_QR; o.max_num_iterations = 60;
                o.num_threads = std::max(1, settings.num_threads); o.logging_type = ceres::LoggingType::SILENT;
                ceres::Solver::Summary sum; ceres::Solve(o, &p, &sum);
                return sum.IsSolutionUsable();
            };

            double beam[2] = {beam_x0, beam_y0}, dist[1] = {dist0};
            double axv[3] = {ax0[0], ax0[1], ax0[2]};
            bool commit = false;
            if (n_obs >= static_cast<size_t>(settings.min_events)) {
                double pos_nom = 0.0, exc_nom = 0.0, pos_ref = 0.0, exc_ref = 0.0;
                const double cv_nom = joint_cost(VAL, beam, dist, axv, p0, p1, p2, pos_nom, exc_nom);
                double bm_f[2], ds_f[1], rv_f[3], q0_f[3], q1_f[3], q2_f[3];
                const bool convJ = solve_joint(FIT, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f);
                const double cv_ref = joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, pos_ref, exc_ref);

                // Commit only for a small, credible move: distance < 1 % (a calibrated header needs
                // < ~0.6 %), each cell length < 1 %, each free cell angle < 1 deg (ANGLE_BOUND_DEG), and
                // the beam inside the bound measured from whichever centre anything already believes is
                // nearer. A larger move is the red flag for an unreliable fit - typically a second
                // lattice whose spots bias every cross-validation fold identically, so the relative "it
                // improved" gate is blind to it. The absolute size of the move discriminates a genuine
                // header correction from that failure far better than the absolute residual, which real
                // marginal (noisy / iced) data shares with it. Names the test it failed, or nullptr.
                const auto out_of_bounds = [&](const double bm[2], const double ds[1],
                                               const double q1c[3], const double q2c[3]) -> const char * {
                    if (std::fabs(ds[0] - dist0) >= 0.01 * dist0)
                        return "the distance moved more than 1 %";
                    double len_shift = 0.0, ang_shift = 0.0;
                    for (int j = 0; j < 3; ++j) {
                        len_shift = std::max(len_shift,
                                             std::fabs(q1c[j] - p1_0[j]) / std::max(1e-9, p1_0[j]));
                        ang_shift = std::max(ang_shift, std::fabs(q2c[j] - p2_0[j]));
                    }
                    if (len_shift >= 0.01)
                        return "a cell length moved more than 1 %";
                    if (ang_shift >= ANGLE_BOUND_DEG * PI / 180.0)
                        return "a cell angle moved more than 1 deg";
                    const double from_nominal = std::hypot(bm[0] - beam_x0, bm[1] - beam_y0);
                    const double from_measured = settings.measured_beam_px
                            ? std::hypot(bm[0] - (*settings.measured_beam_px)[0],
                                         bm[1] - (*settings.measured_beam_px)[1])
                            : std::numeric_limits<double>::infinity();
                    if (std::min(from_nominal, from_measured) >= BEAM_BOUND_PXL)
                        return "the beam moved further than the bound from every centre anything believes";
                    return nullptr;
                };

                // The two residual families are asked separately as well as together. Pooled, the
                // positional values outnumber the excitation ones about three to one where both caps
                // saturate, and the excitation residual is the only one that identifies the cell SCALE -
                // so a pooled mean can improve on the strength of the positions alone while the one
                // quantity the cell is committed for has got worse. Neither may degrade.
                const char *refused =
                        !convJ                        ? "the fit did not converge"
                        : !(cv_ref < 0.98 * cv_nom)   ? "the held-out residual did not improve enough"
                        : !(pos_ref < pos_nom)        ? "the held-out positional residual did not improve"
                        : !(exc_ref < exc_nom)        ? "the held-out excitation residual did not improve"
                                                      : out_of_bounds(bm_f, ds_f, q1_f, q2_f);

                // What the half that was held out earned the right to fit is re-fitted on all of it, and
                // that second solve moves the geometry - so the bounds are asked again of what will
                // actually be committed, not only of the half that passed the gate. Refused here, the
                // run keeps its header geometry exactly as the other refusals leave it.
                double ds_log[1] = {ds_f[0]}, bm_log[2] = {bm_f[0], bm_f[1]};
                double q1_log[3], q2_log[3];
                for (int j = 0; j < 3; ++j) { q1_log[j] = q1_f[j]; q2_log[j] = q2_f[j]; }
                if (refused == nullptr) {
                    solve_joint(ALL, beam, dist, axv, p0, p1, p2);   // commit: re-fit on all data
                    ds_log[0] = dist[0]; bm_log[0] = beam[0]; bm_log[1] = beam[1];
                    for (int j = 0; j < 3; ++j) { q1_log[j] = p1[j]; q2_log[j] = p2[j]; }
                    refused = out_of_bounds(beam, dist, p1, p2);
                    if (refused != nullptr) {
                        beam[0] = beam_x0; beam[1] = beam_y0; dist[0] = dist0;
                        for (int j = 0; j < 3; ++j) {
                            axv[j] = ax0[j]; p0[j] = p0_0[j]; p1[j] = p1_0[j]; p2[j] = p2_0[j];
                        }
                    }
                }
                commit = (refused == nullptr);
                // Name which test refused it. Several different things reject here and the geometry that
                // comes out is the same in all of them, so a run that silently keeps its header geometry
                // says nothing about whether the fit was bad, the improvement too small, or the move too
                // large for the bound - which is the one case where the number worth reading is the one
                // that was thrown away.
                const std::string verdict = commit ? std::string("COMMIT")
                                                   : "reject (" + std::string(refused) + ")";
                // The cell the log names is the EFFECTIVE one - what the residual's B matrix builds
                // from the blocks - not the blocks themselves, whose unused components a high-symmetry
                // system leaves at whatever the seed happened to put there. The geometry beside the
                // verdict is the one the verdict is about: the all-data re-fit where it committed, the
                // half-data fit where it did not. The residual pair stays the split-half gate's, since
                // the fit that is committed has no held-out half of its own.
                double len_nom[3], ang_nom[3], len_fit[3], ang_fit[3];
                EffectiveCellFromParams(sys, p1_0, p2_0, len_nom, ang_nom);
                EffectiveCellFromParams(sys, q1_log, q2_log, len_fit, ang_fit);
                logger.Info("Post-refine GEOM (joint crystal + detector): dist {:.3f} -> {:.3f} mm, beam "
                            "({:.2f},{:.2f}) -> ({:.2f},{:.2f}), cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} "
                            "{:.2f} -> {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, held-out positional "
                            "{:.3e} -> {:.3e}, excitation {:.3e} -> {:.3e} => {}",
                            dist0, ds_log[0], beam_x0, beam_y0, bm_log[0], bm_log[1],
                            len_nom[0], len_nom[1], len_nom[2],
                            ang_nom[0] * 180.0 / PI, ang_nom[1] * 180.0 / PI, ang_nom[2] * 180.0 / PI,
                            len_fit[0], len_fit[1], len_fit[2],
                            ang_fit[0] * 180.0 / PI, ang_fit[1] * 180.0 / PI, ang_fit[2] * 180.0 / PI,
                            pos_nom, pos_ref, exc_nom, exc_ref, verdict);
            } else {
                logger.Info("Post-refine GEOM: only {} positional observations - the joint fit needs the "
                            "spot positions as well as the rocking angles, so nothing is refined", n_obs);
            }
            result.cell_refined = commit;
            result.detector_refined = commit;

            const double axlen = std::sqrt(axv[0]*axv[0] + axv[1]*axv[1] + axv[2]*axv[2]);
            const double axdev = std::acos(std::clamp((axv[0]*ax0[0]+axv[1]*ax0[1]+axv[2]*ax0[2])
                                                      / std::max(1e-9, axlen), -1.0, 1.0)) * 180.0 / PI;
            // The committed crystal, as a lattice again. Where nothing was committed this is the lattice
            // rotation indexing handed in, unchanged - not its round trip through the symmetry-constrained
            // parameterisation, which would move the cell for a fit that was refused.
            double eff_len[3], eff_ang[3];
            EffectiveCellFromParams(sys, p1, p2, eff_len, eff_ang);
            const CrystalLattice committed_latt = commit
                    ? AngleAxisAndCellToLattice(p0, eff_len, eff_ang[0], eff_ang[1], eff_ang[2])
                    : reference_latt;
            const Coord As = committed_latt.Astar(), Bs = committed_latt.Bstar(), Cs = committed_latt.Cstar();
            if (commit) {
                result.cell = UnitCell{static_cast<float>(eff_len[0]), static_cast<float>(eff_len[1]),
                                       static_cast<float>(eff_len[2]),
                                       static_cast<float>(eff_ang[0] * 180.0 / PI),
                                       static_cast<float>(eff_ang[1] * 180.0 / PI),
                                       static_cast<float>(eff_ang[2] * 180.0 / PI)};
                logger.Info("Post-refine GEOM: committed cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, "
                            "rotation axis moved {:.3f} deg (the cell is the fit's own; the second pass "
                            "re-indexes at the refined detector geometry)", result.cell.a, result.cell.b,
                            result.cell.c, result.cell.alpha, result.cell.beta, result.cell.gamma, axdev);
            }

            // ===== Goniometer rotation SCALE k, its own one-parameter fit on the same rocking events =====
            // The angles stored in the file are the COMMANDED ones, so a stage that turned k times as far
            // is invisible in the header. Nothing else here can represent it: the cell scale, the axis
            // direction, the distance and the beam are all orthogonal to a rotation MAGNITUDE error. Fitted
            // after step A so the cell scale and the axis direction are fixed at their committed values and
            // k is the only free quantity.
            const double u[3] = {axv[0] / axlen, axv[1] / axlen, axv[2] / axlen};
            double phi_c = 0.0, phi_lo = events[0].phi_obs, phi_hi = events[0].phi_obs;
            for (size_t e = 0; e < n_events; ++e) {
                phi_c += events[e].phi_obs;
                phi_lo = std::min(phi_lo, events[e].phi_obs);
                phi_hi = std::max(phi_hi, events[e].phi_obs);
            }
            phi_c /= static_cast<double>(n_events);
            const double sweep_deg = (phi_hi - phi_lo) * 180.0 / PI;
            // The committed reciprocal vector turned to the sweep centre. The
            // angle then enters the fit measured FROM that centre. A constant crystal missetting about the
            // spindle is k with a slope in phi, so measuring the angle from the goniometer's zero instead
            // lets a missetting leak into k with gain <phi>/<phi^2> - which depends only on where the sweep
            // happens to sit. On a short sweep starting near zero that gain is enormous: a 0.14 deg
            // missetting on a 10 deg wedge fakes 1.4 % of k. Referred to the sweep centre the leak is
            // identically zero at any width, and no parameter has to be added to get it.
            // The residual is closed-form in k, so this is a one-parameter minimisation rather than a
            // solver problem. A rotation preserves length, so |p_lab| = |e_mid| whatever k is, and only
            // the z component moves; Rodrigues gives it exactly:
            //
            //   r(k) = C + A cos(a k) - B sin(a k) = C + R cos(a k + psi)
            //   C = lambda |e|^2 / 2 + u_z (u.e),  A = e_z - u_z (u.e),  B = (u x e)_z,  a = phi_obs - phi_c
            //
            // which is the same function the residual functor computes, to the last bit. Handing 8 million
            // one-parameter residual blocks to Ceres instead cost tens of millions of allocations and a
            // dense factorisation per iteration, for a fit that a scan over a bounded interval settles.
            // Coefficients are computed in double and stored narrowed: their rounding perturbs the
            // minimiser by ~1e-10, and k is carried downstream as a float.
            struct ScaleTerm { float a, C, R, psi; };
            std::vector<ScaleTerm> terms(n_events);
            std::vector<int> fifth_of(n_events);
            const double aa_c[3] = {-phi_c * u[0], -phi_c * u[1], -phi_c * u[2]};
            ParallelChunks(static_cast<int>(n_events), nthreads, [&](int lo, int hi) {
                for (int e = lo; e < hi; ++e) {
                    const Coord ec = As * static_cast<float>(events[e].h)
                                     + Bs * static_cast<float>(events[e].k)
                                     + Cs * static_cast<float>(events[e].l);
                    const double p[3] = {ec.x, ec.y, ec.z};
                    double em[3];
                    ceres::AngleAxisRotatePoint(aa_c, p, em);
                    const double ue = u[0] * em[0] + u[1] * em[1] + u[2] * em[2];
                    const double e2 = em[0] * em[0] + em[1] * em[1] + em[2] * em[2];
                    const double C = 0.5 * lambda_l * e2 + u[2] * ue;
                    const double A = em[2] - u[2] * ue;
                    const double B = u[0] * em[1] - u[1] * em[0];
                    terms[e] = ScaleTerm{static_cast<float>(events[e].phi_obs - phi_c),
                                         static_cast<float>(C), static_cast<float>(std::hypot(A, B)),
                                         static_cast<float>(std::atan2(B, A))};
                    fifth_of[e] = std::clamp(static_cast<int>(
                            5.0 * (events[e].phi_obs - phi_lo) / std::max(1e-9, phi_hi - phi_lo)), 0, 4);
                }
            });

            // Robust-loss scale from the scatter the events actually have: it varies by more than a decade
            // between datasets, so any fixed constant is either inert or throws away real data. Taken once,
            // over every event, so the all-data fit and every jackknife fold share it.
            const auto residual_at = [&](const ScaleTerm &t, double k) {
                return static_cast<double>(t.C)
                     + static_cast<double>(t.R) * std::cos(static_cast<double>(t.a) * k + t.psi);
            };
            double rms = 0.0;
            for (const auto &t : terms) {
                const double r = residual_at(t, 1.0);
                rms += r * r;
            }
            rms = std::sqrt(rms / static_cast<double>(terms.size()));
            const double huber_delta = std::max(1e-12, 2.0 * rms);
            const double huber_d2 = huber_delta * huber_delta;

            // Ceres minimises half the sum of the loss applied to the SQUARED residual, so that is what is
            // reproduced here. One pass yields the five per-fifth partial sums, which serve the all-data
            // fit and all five leave-a-fifth-out folds together.
            // Each chunk folds into its own slot and the slots are summed in chunk order, so the
            // sum does not depend on which worker finishes first: the same events always add up in
            // the same sequence, and the fit is reproducible run to run.
            const int n_terms = static_cast<int>(terms.size());
            const int cost_nt = static_cast<int>(std::max<size_t>(1, std::min(nthreads,
                                                                             static_cast<size_t>(n_terms))));
            const int cost_chunk = (n_terms + cost_nt - 1) / cost_nt;
            // Several k are always wanted at once (the grid below asks for 101), and they all sweep the
            // same event list, so sweep it ONCE and evaluate every k on each event while it is still in
            // registers. The per-thread accumulator is one slot per (k, fifth) - 4 kB for the grid, small
            // enough to stay in L1 - against re-reading the whole term array once per k. Each (k, fifth)
            // still receives its events in the same order and the chunks are still summed in chunk order,
            // so the sums are the ones a k-at-a-time loop produced, bit for bit.
            const auto cost_grid = [&](const std::vector<double> &ks) {
                const int nk = static_cast<int>(ks.size());
                std::vector<std::vector<std::array<double, 5>>> per_chunk(
                        cost_nt, std::vector<std::array<double, 5>>(nk));
                ParallelChunks(n_terms, nthreads, [&](int lo, int hi) {
                    std::vector<std::array<double, 5>> acc(nk);
                    for (int e = lo; e < hi; ++e) {
                        const ScaleTerm &term = terms[e];
                        const int fifth = fifth_of[e];
                        for (int g = 0; g < nk; ++g) {
                            const double r = residual_at(term, ks[g]);
                            const double s2 = r * r;
                            acc[g][fifth] += (s2 <= huber_d2) ? s2
                                                              : (2.0 * huber_delta * std::sqrt(s2) - huber_d2);
                        }
                    }
                    per_chunk[lo / cost_chunk] = std::move(acc);
                });
                std::vector<std::array<double, 5>> total(nk);
                for (const auto &acc : per_chunk)
                    for (int g = 0; g < nk; ++g)
                        for (int j = 0; j < 5; ++j) total[g][j] += acc[g][j];
                return total;
            };
            const auto cost_by_fifth = [&](double k) { return cost_grid({k})[0]; };

            // Scan the interval Ceres was bounded to, then close in. No event's phase can move by more than
            // a fraction of a period over an interval this narrow, so the objective has no structure the
            // grid could step over; the refinement is only there to place the minimum precisely.
            constexpr int SCALE_GRID = 101;
            constexpr double SCALE_K_LO = 0.95, SCALE_K_HI = 1.05;
            std::vector<double> grid_k(SCALE_GRID);
            for (int g = 0; g < SCALE_GRID; ++g)
                grid_k[g] = SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1);
            const std::vector<std::array<double, 5>> grid = cost_grid(grid_k);
            auto solve_scale = [&](int drop_fifth) {
                const auto total = [&](const std::array<double, 5> &f) {
                    double t = 0.0;
                    for (int j = 0; j < 5; ++j)
                        if (j != drop_fifth) t += f[j];
                    return t;
                };
                int best = 0;
                for (int g = 1; g < SCALE_GRID; ++g)
                    if (total(grid[g]) < total(grid[best])) best = g;
                const double step = (SCALE_K_HI - SCALE_K_LO) / (SCALE_GRID - 1);
                double a = std::max(SCALE_K_LO, SCALE_K_LO + step * (best - 1));
                double b = std::min(SCALE_K_HI, SCALE_K_LO + step * (best + 1));
                // Golden section: the objective is smooth but its curvature jumps wherever an event
                // crosses the Huber knee, which a derivative method would have to cope with.
                constexpr double INV_PHI = 0.6180339887498949;
                double c = b - INV_PHI * (b - a), d = a + INV_PHI * (b - a);
                double fc = total(cost_by_fifth(c)), fd = total(cost_by_fifth(d));
                // The fit is narrowed to a float before it is applied (Rugnux.h prepass_rotation_scale_),
                // so bracketing it below that type's epsilon, 6e-8, only buys about ten more full
                // passes over the events for a digit that cannot survive being stored.
                while (b - a > 1e-7) {
                    if (fc < fd) { b = d; d = c; fd = fc; c = b - INV_PHI * (b - a); fc = total(cost_by_fifth(c)); }
                    else         { a = c; c = d; fc = fd; d = a + INV_PHI * (b - a); fd = total(cost_by_fifth(d)); }
                }
                return 0.5 * (a + b);
            };
            const double k_fit = solve_scale(-1);
            result.rotation_scale = k_fit;

            // ----- Whether to COMMIT it. A stage fault is rare - 36 of 37 rotation datasets sit at 1.0000
            // on a direct scan - and a 1 % angle correction applied to a healthy dataset would damage it
            // silently, so every test below has to pass.
            // Preconditions: below these the fit is reported but never acted on. Under ~30 deg of sweep k
            // entangles with the axis direction and 10-20 deg truncations of a perfect dataset wander by
            // +-0.6 %; a screening wedge must not trigger a correction.
            constexpr int MIN_SCALE_EVENTS = 5000;
            constexpr double MIN_SCALE_SWEEP_DEG = 30.0;
            // T1 significance: 0.5 % is 18 sigma on the between-dataset scatter of healthy stages
            // (robust sd 2.8e-4) and still 3.5x below the one measured fault.
            constexpr double ROTATION_SCALE_TOL = 0.005;
            // T2 relevance: the misorientation the error produces at each end of the sweep. A large k over
            // a short sweep moves nothing and is not worth correcting.
            constexpr double MIN_SCALE_END_ERROR_DEG = 0.5;
            // T3 uniformity: a stage error is a ramp present in EVERY part of the sweep, so dropping any
            // fifth of it must leave the same k. A second lattice that dominates ONE END of the sweep -
            // exactly what happens where the primary stops indexing - fakes a k indistinguishable from a
            // real fault on T1 and T2, and is the reason this test is not optional. It replaces the
            // hkl-hash split used elsewhere here, which cannot see it: both halves of that split sit at
            // the same angles, so anything structured in phi survives in both folds.
            constexpr double MIN_SCALE_JACKKNIFE_FRAC = 0.5;
            const double end_error_deg = std::fabs(k_fit - 1.0) * sweep_deg / 2.0;
            const bool enough_data = static_cast<int>(n_events) >= MIN_SCALE_EVENTS
                                     && sweep_deg >= MIN_SCALE_SWEEP_DEG;
            const bool big_enough = enough_data && std::fabs(k_fit - 1.0) >= ROTATION_SCALE_TOL
                                    && end_error_deg >= MIN_SCALE_END_ERROR_DEG;
            double jackknife = 1.0;
            if (big_enough)
                for (int f = 0; f < 5; ++f)
                    jackknife = std::min(jackknife, (solve_scale(f) - 1.0) / (k_fit - 1.0));
            result.rotation_scale_suspect = big_enough && jackknife >= MIN_SCALE_JACKKNIFE_FRAC;
            logger.Info("Post-refine rotation SCALE: k = {:.5f} over {:.0f} deg of sweep centred on {:.1f} "
                        "deg ({} events): end error {:.2f} deg, leave-a-fifth-out {:.2f} => {}",
                        k_fit, sweep_deg, phi_c * 180.0 / PI, n_events, end_error_deg, jackknife,
                        result.rotation_scale_suspect ? "COMMIT"
                        : !enough_data ? "report only (too little sweep or too few events)"
                        : "reject (kept the stored angles)");
            if (result.rotation_scale_suspect)
                logger.Warning("Goniometer rotation scale looks off by {:+.2f} % (fitted {:.5f}): the stage "
                               "appears to have turned {} than the angles stored in the file, which are the "
                               "COMMANDED values. This is a hardware calibration fault, not a data problem - "
                               "left uncorrected it inflates mosaicity, biases the cell and loses "
                               "high-resolution reflections",
                               100.0 * (k_fit - 1.0), k_fit, k_fit > 1.0 ? "further" : "less far");

            // Assemble the committed geometry.
            result.distance_after_mm = dist[0];
            result.beam_x_before_px = beam_x0; result.beam_x_after_px = beam[0];
            result.beam_y_before_px = beam_y0; result.beam_y_after_px = beam[1];
            result.events_used = static_cast<int>(selected.size());
            result.ok = commit;
            if (!result.ok)
                logger.Info("Post-refine GEOM: the joint fit did not pass cross-validation - geometry "
                            "left at nominal");
            return result;
        }

        return result;   // refine_geometry is the only supported mode; nothing refined otherwise
    } catch (...) {
        result.ok = false;
        return result;
    }
}