rugnux: rotation detector distance as the fixed point of the integration; beam bound on the direct beam

Where pass 1's post-refinement leaves the distance at the header but its free fit asks for a move of a
step (1 %) or more, the distance is now asked for as a fixed point: pre-passes at distances along the free
fit's move (the direct beam held where the free fit put it), the move each one asks for against the
distance it was integrated at, and the zero of the least-squares line through them - taken only where the
distances tried bracket it (interpolated, never extrapolated). That distance is integrated once more and
adopted only where it realises a lower held-out residual than the header pass, beyond its noise, or the
validation frames prefer it - the tests a round of the geometry walk passes. Otherwise the header stands.

Why: a fit is measured on reflections integrated at the geometry it starts from and is pulled towards it.
On a detector swung out on a 2theta arm with a header 4 % short, the fit asked for 104.4 mm from data
integrated at 103 and 108.6 mm from data integrated at 109; spot-finder centroids in place of the
integrated ones moved neither. Each round of the existing walk closes a third of the gap and the walk
stopped 3 mm short. Measured (CPU, from the headers, no -C):
- two 2theta = 35 deg sweeps of one cubic crystal: 103 -> 108.0 / 107.5 mm, cell 74.4 -> 77.8 A,
  ISa 7.2 -> 13.4 and 4.5 -> 25.9, the second's space group back from P212121 to the cubic group;
- a 2theta = 45 deg sweep of an orthorhombic crystal: 183 -> 189.9 mm, cell 4.5 % -> 0.4 % short,
  ISa 4.2 -> 8.2, R_meas 32 -> 18 %;
- low-resolution sweeps at 2theta = 10 deg ask for +5.5 mm from every distance between 183 and 205 mm
  and bracket nothing: header kept, output bit-identical;
- eight controls whose header is right (free fit within 0.3 %): output bit-identical.

The beam the post-refinement fits is the PONI; on a swung detector it moves with the distance at a fixed
direct beam (7 px per mm at 35 deg), so the 15 px bound on it capped the distance near 105 mm. The bound,
and the comparison with a measured centre (which measures the direct beam), now apply to the direct beam;
on an untilted detector nothing changes.

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

image_analysis/geom_refinement/PostRefine.cpp

@@ -53,7 +53,7 @@ constexpr double ANGLE_BOUND_DEG = 1.0;
 // solve reaches at most this far, and a fit that wants more takes another step from there. A move of
 // more than one step is then ratified by re-indexing at it, which is what actually tells a second
 // lattice from a wrong header (see PostRefineResult::large_move).
-constexpr double STEP_BOUND = 0.01;
+constexpr double STEP_BOUND = POSTREFINE_STEP_BOUND;
 
 // ... and the walk is bounded, so a fit that never settles still terminates. Five steps is the +-5 %
 // box the solver was always given on the distance, i.e. as far as this fit has ever been trusted.
@@ -424,6 +424,20 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
             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();
+            // The beam the fit refines is the PONI, the foot of the normal from the sample. On a detector
+            // swung out on a 2theta arm the direct beam lies a distance-proportional way from it, so at a
+            // fixed direct beam the PONI moves with the distance (7 px per mm at 35 deg on a 0.1 mm pixel),
+            // and a bound on the PONI caps the distance instead of the beam. What the bound is about is
+            // where the beam hits the detector, which is also what a measured centre measures, so it is
+            // applied there: direct beam = PONI + beam_per_mm * distance. Zero on an untilted detector.
+            double beam_per_mm[2] = {0.0, 0.0};
+            {
+                const auto [dbx, dby] = nominal_geom.GetDirectBeam_pxl();
+                if (std::isfinite(dbx) && std::isfinite(dby)) {
+                    beam_per_mm[0] = (dbx - beam_x0) / dist0;
+                    beam_per_mm[1] = (dby - beam_y0) / dist0;
+                }
+            }
 
             // 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
@@ -670,11 +684,15 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
                 // the bound and is then refused for being there.
                 p.SetParameterLowerBound(ds, 0, wide ? dist0 * 0.95 : ds_seed * (1.0 - STEP_BOUND));
                 p.SetParameterUpperBound(ds, 0, wide ? dist0 * 1.05 : ds_seed * (1.0 + STEP_BOUND));
+                // The PONI box reaches every PONI whose direct beam the gate below can accept, at any
+                // distance the walk can reach (MAX_STEPS steps).
                 for (int j = 0; j < 2; ++j) {
                     const double b0 = (j == 0) ? beam_x0 : beam_y0;
-                    const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : b0;
-                    p.SetParameterLowerBound(bm, j, std::min(b0, m) - BEAM_BOUND_PXL);
-                    p.SetParameterUpperBound(bm, j, std::max(b0, m) + BEAM_BOUND_PXL);
+                    const double m = settings.measured_beam_px
+                                     ? (*settings.measured_beam_px)[j] - beam_per_mm[j] * dist0 : b0;
+                    const double reach = std::fabs(beam_per_mm[j]) * dist0 * MAX_STEPS * STEP_BOUND;
+                    p.SetParameterLowerBound(bm, j, std::min(b0, m) - BEAM_BOUND_PXL - reach);
+                    p.SetParameterUpperBound(bm, j, std::max(b0, m) + BEAM_BOUND_PXL + reach);
                 }
                 for (int j = 0; j < 3; ++j) {
                     p.SetParameterLowerBound(rv, j, ax0[j] - 0.05);
@@ -770,6 +788,11 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
                 bool convJ = solve_joint(FIT, /*wide=*/true, /*fix_distance=*/false, bm_f, ds_f, rv_f,
                                          q0_f, q1_f, q2_f, distance_corr);
                 const bool held_converged = held_fit.valid() && held_fit.get();
+                if (convJ) {
+                    result.distance_free_mm = ds_f[0];
+                    result.beam_x_free_px = bm_f[0];
+                    result.beam_y_free_px = bm_f[1];
+                }
                 double exc_ref_se = NAN;
                 double cv_ref = convJ ? joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f,
                                                    pos_ref, exc_ref, nullptr, &exc_ref_se)
@@ -876,16 +899,19 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
                 // step is ratified by re-indexing at where it arrived (result.large_move), which is the
                 // discriminator the size of the move was standing in for. Names the test it failed, or
                 // nullptr.
-                const auto out_of_bounds = [&](const double bm[2], const double q2c[3]) -> const char * {
+                const auto out_of_bounds = [&](const double bm[2], const double ds[1],
+                                               const double q2c[3]) -> const char * {
                     double ang_shift = 0.0;
                     for (int j = 0; j < 3; ++j)
                         ang_shift = std::max(ang_shift, std::fabs(q2c[j] - p2_0[j]));
                     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 db[2] = {bm[0] + beam_per_mm[0] * ds[0], bm[1] + beam_per_mm[1] * ds[0]};
+                    const double from_nominal = std::hypot(db[0] - beam_x0 - beam_per_mm[0] * dist0,
+                                                           db[1] - beam_y0 - beam_per_mm[1] * dist0);
                     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::hypot(db[0] - (*settings.measured_beam_px)[0],
+                                         db[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";
@@ -949,7 +975,7 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
                                 ? "the held-out positional residual did not fall by more than its noise"
                         : !(change.exc_n < 2 || change.exc <= change.exc_se)
                                 ? "the held-out excitation residual rose by more than its noise"
-                                                      : out_of_bounds(bm_f, q2_f);
+                                                      : out_of_bounds(bm_f, ds_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
@@ -971,7 +997,7 @@ PostRefineResult PostRefineRotationGeometry(PostRefineObservations observations,
                     solve_joint(ALL, !walked, hold_distance, 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, p2);
+                    refused = out_of_bounds(beam, dist, p2);
                     if (refused != nullptr) {
                         beam[0] = beam_x0; beam[1] = beam_y0; dist[0] = dist0;
                         for (int j = 0; j < 3; ++j) {

image_analysis/geom_refinement/PostRefine.h

@@ -31,6 +31,10 @@
 // Where there are more reflections than the fit's caps, the strongest by I/sigma are the ones kept;
 // inside the fit every reflection carries the same weight.
 
+// How far one step of the joint fit may move the distance or a cell length, relative (see STEP_BOUND in
+// PostRefine.cpp, where it is explained).
+constexpr double POSTREFINE_STEP_BOUND = 0.01;
+
 struct PostRefineResult {
     bool ok = false;
     DiffractionGeometry geom;                 // refined (distance / beam left as configured)
@@ -53,6 +57,11 @@ struct PostRefineResult {
     // ever - because freeing it did not lower the held-out excitation residual by more than that
     // residual's noise (see the joint fit in PostRefine.cpp). The report names it.
     bool distance_held = false;
+    // Where the fit with the distance free went, held or not (NaN where it did not converge): the
+    // detector move these data ask for at the geometry they were integrated at.
+    double distance_free_mm = std::numeric_limits<double>::quiet_NaN();
+    double beam_x_free_px = std::numeric_limits<double>::quiet_NaN();
+    double beam_y_free_px = std::numeric_limits<double>::quiet_NaN();
     // The held-out residual at the geometry the fit STARTED from - what the integration this fit was
     // handed actually realises there, before the fit moves anything - and its standard error over the
     // held-out values. Two passes at two geometries are compared on this, not on the fit's own

rugnux/Rugnux.h

@@ -1820,7 +1820,14 @@ class Rugnux {
         std::optional<GoniometerAxis> gonio_snapshot;
         std::array<float, 5> header_geometry{};
         RotationScaleOutcome scale_walk;   // where the walk has ended before they start
+        int distance_passes = 0;           // the pre-passes the fixed-point distance ran (FixedPointDistance)
     };
+    // The detector distance as the fixed point of the integration, where pass 1 held the header but its
+    // free fit asked for a move of a step or more (see RunAllPasses). Returns the geometry to start the
+    // refined passes from, or nothing where the header stands; `passes` counts the pre-passes it ran.
+    std::optional<std::array<float, 5>> FixedPointDistance(RugnuxObserver *observer, const ProcessResult &pass1,
+                                                           const DiffractionExperiment &first_pass_start,
+                                                           int &passes);
     // The refined passes of the rotation two-pass, from the second pass on. Empty where they ran beside
     // the rotation-scale walk and it overturned them.
     // After the refined passes: where the supercell probe's class stands clear, integrate and merge the

rugnux/RugnuxPasses.cpp

@@ -627,12 +627,21 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
                                                    experiment_.GetDetectorDistance_mm(),
                                                    experiment_.GetPoniRot1_rad(),
                                                    experiment_.GetPoniRot2_rad()};
+        int distance_passes = 0;
+        if (const auto g = FixedPointDistance(observer, pass1, experiment_before_first_pass, distance_passes)) {
+            prepass_detector_geometry_ = g;
+            pass1.post_refine->distance_after_mm = (*g)[2];
+            pass1.post_refine->beam_x_after_px = (*g)[0];
+            pass1.post_refine->beam_y_after_px = (*g)[1];
+            pass1.post_refine->distance_held = false;
+        }
+        if (cancelled_) return pass1;
         if (prepass_detector_geometry_) {
             const auto &g = *prepass_detector_geometry_;
             experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2])
                        .PoniRot1_rad(g[3]).PoniRot2_rad(g[4]);
         }
-        RefinedPassStart start{gonio_snapshot, header_geometry, {}};
+        RefinedPassStart start{gonio_snapshot, header_geometry, {}, distance_passes};
 
         // The goniometer rotation scale: the angles in the file are the commanded ones, so a stage that
         // turned at the wrong rate is a geometry error like any other - but one pass 1's fit reads only
@@ -713,6 +722,138 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
     return RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
 }
 
+// The detector distance as the fixed point of the integration. Pass 1 integrated at the header, and its
+// post-refinement left the distance there - the rocking angles did not ask for the move the free fit
+// wanted, or the held-out residual did not pay for it - but a fit is measured on reflections integrated at the geometry it starts from, and those
+// carry that geometry with them: which reflections are predicted, and on which frames their rocking
+// centroids are taken. Measured on a detector swung out on a 2theta arm with a header 4 % short, the fit
+// asked for 104.4 mm from data integrated at 103 and for 108.6 mm from data integrated at 109, each a
+// third of the way to where the other started, so a walk that moves to what each round asks for crawls
+// and stalls short; spot-finder centroids in place of the integrated ones moved neither number.
+//
+// So the distance is asked for as what it is, the point where the data integrated there ask for no
+// move: a few pre-passes at distances along the free fit's direction (the direct beam held where the free
+// fit put it), the move each one asks for against the distance it was integrated at, and the zero of the
+// straight line through them. That distance is integrated once more and adopted only for what it
+// REALISES against the header pass: a lower held-out residual beyond its noise, or validation spots the
+// lattice takes better - the two tests a round of the geometry walk has to pass. A header that is right
+// asks for less than a step and never gets here.
+std::optional<std::array<float, 5>> Rugnux::FixedPointDistance(RugnuxObserver *observer, const ProcessResult &pass1,
+                                                               const DiffractionExperiment &first_pass_start,
+                                                               int &passes) {
+    if (cancelled_ || config_.fast_live || !pass1.post_refine)
+        return std::nullopt;
+    const PostRefineResult &held = *pass1.post_refine;
+    const double d0 = experiment_.GetDetectorDistance_mm();
+    if (prepass_detector_geometry_ && std::fabs((*prepass_detector_geometry_)[2] - d0) > 1e-3 * d0)
+        return std::nullopt;   // pass 1 moved the distance itself; the walk takes it from there
+    const double move = held.distance_free_mm - d0;
+    if (!std::isfinite(move) || std::fabs(move) < POSTREFINE_STEP_BOUND * d0)
+        return std::nullopt;
+    Logger logger("Rugnux");
+    if (observer)
+        observer->OnStep("Distance check");
+
+    // The direct beam the free fit put the beam at, held there at every distance tried: on a swung
+    // detector the PONI moves with the distance at a fixed direct beam (zero per mm on an untilted one).
+    double per_mm[2] = {0.0, 0.0};
+    const auto [dbx, dby] = experiment_.GetDiffractionGeometry().GetDirectBeam_pxl();
+    if (std::isfinite(dbx) && std::isfinite(dby)) {
+        per_mm[0] = (dbx - experiment_.GetBeamX_pxl()) / d0;
+        per_mm[1] = (dby - experiment_.GetBeamY_pxl()) / d0;
+    }
+    const double direct_beam[2] = {held.beam_x_free_px + per_mm[0] * held.distance_free_mm,
+                                   held.beam_y_free_px + per_mm[1] * held.distance_free_mm};
+    const auto geometry_at = [&](double d) {
+        return std::array<float, 5>{static_cast<float>(direct_beam[0] - per_mm[0] * d),
+                                    static_cast<float>(direct_beam[1] - per_mm[1] * d), static_cast<float>(d),
+                                    experiment_.GetPoniRot1_rad(), experiment_.GetPoniRot2_rad()};
+    };
+    // A pre-pass at a distance, on a copy of the run as pass 1 started it: it leaves nothing on the run
+    // but its first-pass memo.
+    const auto prepass_at = [&](double d) {
+        Rugnux probe(*this);
+        probe.experiment_ = first_pass_start;
+        if (bragg_adaptive_)
+            probe.experiment_.ImportBraggIntegrationSettings(*bragg_adaptive_);
+        const auto g = geometry_at(d);
+        probe.experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]);
+        probe.prepass_mosaicity_.clear();
+        probe.prepass_detector_geometry_.reset();
+        probe.prepass_result_.reset();
+        ProcessResult r;
+        try {
+            r = probe.RunPipeline(observer, /*write_output=*/false, /*geometry_prepass=*/true);
+        } catch (const std::exception &e) {
+            if (IsFatalResourceError(e)) throw;
+            logger.Info("Distance check: the pre-pass at {:.3f} mm did not complete ({})", d, e.what());
+        }
+        for (const auto &m : probe.first_pass_memo_)
+            AddFirstPassMemo(m);
+        ++passes;
+        return r;
+    };
+    const auto asked = [](const ProcessResult &r) {
+        if (!r.post_refine) return std::numeric_limits<double>::quiet_NaN();
+        return std::isfinite(r.post_refine->distance_free_mm) ? r.post_refine->distance_free_mm
+                                                              : r.post_refine->distance_after_mm;
+    };
+
+    // Distances along the free fit's move, one move apart, until the move asked for changes sign.
+    std::vector<std::pair<double, double>> points{{d0, held.distance_free_mm}};
+    std::string text = fmt::format("{:.2f} -> {:.2f}", d0, held.distance_free_mm);
+    for (int k = 1; k <= 4 && !cancelled_; ++k) {
+        const double d = d0 + k * move;
+        const double a = asked(prepass_at(d));
+        text += fmt::format(", {:.2f} -> {:.2f}", d, a);
+        if (!std::isfinite(a))
+            break;
+        points.emplace_back(d, a);
+        if ((a - d) * move <= 0.0)
+            break;
+    }
+    if (cancelled_)
+        return std::nullopt;
+    // The straight line through (distance, asked - distance), least squares, and its zero.
+    double sx = 0, sy = 0, sxx = 0, sxy = 0;
+    for (const auto &[d, a] : points) {
+        sx += d; sy += a - d; sxx += d * d; sxy += d * (a - d);
+    }
+    const double n = static_cast<double>(points.size());
+    const double slope = (n * sxy - sx * sy) / (n * sxx - sx * sx);
+    const double intercept = (sy - slope * sx) / n;
+    // A fixed point is taken only where the distances tried bracket it - the move asked for changed
+    // sign - so it is interpolated, never extrapolated. Where every distance asks for more of the same
+    // move the data do not pin the distance at all (measured: a crystal at 3 A on a detector 10 deg out
+    // asked for 5.5 mm more from every distance between 183 and 205 mm), and the header stands.
+    const double last_move = points.back().second - points.back().first;
+    if (points.size() < 2 || !(slope < 0.0) || last_move * move > 0.0) {
+        logger.Info("Distance check: the moves asked for at each distance integrated ({}) do not bracket a "
+                    "fixed point - the header stands", text);
+        if (observer)
+            observer->OnNote("the header distance stands");
+        return std::nullopt;
+    }
+    const double d_star = std::clamp(-intercept / slope, std::min(d0, points.back().first),
+                                     std::max(d0, points.back().first));
+    const ProcessResult at_star = prepass_at(d_star);
+    const bool residual_fell = at_star.post_refine && HeldOutResidualFell(held, *at_star.post_refine);
+    const bool evidence_prefers = ValidationEvidencePrefers(pass1.validation_evidence, at_star.validation_evidence);
+    const bool adopt = residual_fell || evidence_prefers;
+    logger.Info("Distance check: integrated at -> asked for (mm): {}; fixed point {:.3f} mm (asked minus integrated "
+                "falls {:.2f} mm per mm); there the held-out residual is {:.3e} against the header's {:.3e} and "
+                "the validation spots on the lattice {} against {} - {}", text, d_star, -slope,
+                at_star.post_refine ? at_star.post_refine->held_out_before : NAN, held.held_out_before,
+                at_star.validation_evidence.on_lattice, pass1.validation_evidence.on_lattice,
+                adopt ? "adopted" : "the header stands");
+    if (observer)
+        observer->OnNote(adopt ? fmt::format("detector distance {:.2f} mm, where the data ask for no move", d_star)
+                               : std::string("the header distance stands"));
+    if (!adopt)
+        return std::nullopt;
+    return geometry_at(d_star);
+}
+
 // The rotation-scale walk of RunAllPasses, run beside the refined passes on a copy of the run (`check`
 // is what the two share): it starts once their first pass has indexed, calls them off (output_cancel)
 // before a refit, and says at its end whether it overturns them.
@@ -918,7 +1059,7 @@ std::optional<ProcessResult> Rugnux::RefinedPasses(RugnuxObserver *observer, Pro
     RotationScaleOutcome scale_walk = start.scale_walk;
     // The probe passes below and the lattice arms after them exist ONLY to measure - see there - so
     // the count of passes the run made has to include them even though they write nothing.
-    int arm_passes = 0;
+    int arm_passes = start.distance_passes;
 
     // ... and apply the adopted scale, on the same measure-then-re-integrate footing as the geometry.
     if (prepass_rotation_scale_ && gonio_snapshot) {