post-refine: bound the beam move on the measurement, not on the header

PostRefine step B refuses a fitted beam centre more than 15 px from the geometry the run started
from. Two things were wrong with that and they compounded.

The 15 px was BOTH the Ceres search box and the acceptance threshold, so a fit that wanted to move
further was pinned at the box face, landed at exactly hypot >= 15.0, and was then refused for being
there. The gate never saw the fit it was judging. Measured on one crystal at three detector
distances, the fit came back at 15.128, 15.133 and 15.130 px - the box corner, three times.

And it was measured from the header, which is the value most worth correcting exactly when it is
most wrong. The run already knows better: the default beam-centre check fits the centre from the
isotropy of the scattered background on every rotation run, reported it, and then discarded it.

The bound is now measured from whichever of the nominal centre and that measurement is nearer -
accept a move that is small relative to something independent, rather than small relative to the
header alone - and the search box is widened to match, so the gate decides rather than the box. The
measurement is used only where it was precise enough to be adopted as a centre in its own right; a
run without one is bounded exactly as before. Step B also now says WHICH of its four tests refused a
fit, since all four leave the same geometry behind.

Measured over 113 rotation datasets against the same corpus, one binary per arm, with a control arm
that reproduced the stored battery on 112 of 112 reports:

  the bound fired and refused a real correction   1   6ukf
  the bound fired and caught a runaway            0
  the bound never fired                         112

On 6ukf the header centre is 22.1 px from what the background fit measures. The old fit, pinned at
the box, improved the held-out positional residual 9-fold and was refused; the new fit lands 0.46 px
from the measurement and improves it 418-fold. ISa 4.42 -> 6.71, R_meas 18.5 -> 12.8 %, CC1/2 0.987
-> 0.995, d_min 1.119 -> 0.973 A; at matched resolution, 1.68 A goes I/sigma 8.9 -> 15.1 and R_meas
28.8 -> 13.4 %. The refined cell moves from 0.4 % off the deposited one to under 0.1 %.

It survives perturbation, and inverts it: over +-0.5 mm of detector distance the corrected run is
flat (ISa 6.62 / 6.71 / 6.94) while the uncorrected one is bistable (4.02 / 4.42 / 6.90). The change
removes an instability rather than exploiting one.

No space group anywhere in the corpus moves, and no other dataset's report changes at all - 112 of
113 are identical. The distance half of the gate is untouched: it is separately argued, and it is
what legitimately refuses the one dataset whose distance wants to move 1.05 %.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N

docs/CHANGELOG.md

@@ -18,6 +18,7 @@
 * The `rugnux` report gives the detector tilt and the direct beam beside the beam centre, and carries the tilt in the `dataset_settings` block.
 * The `rugnux` report gives the detector tilt rotation indexing measured, so it can be checked against a powder calibration.
 * `--no-refine-tilt` holds the detector tilt at the value in the file, instead of zeroing it, when the calibration starts from the spots.
+* A post-refined beam centre is judged against the measurement the run makes for itself, so a badly-placed header no longer caps how far it may be corrected.
 * The `jfjoch_viewer` grid scan view draws the cells in the proportion of the scan steps, so the map has the shape of the scanned area.
 
 ### 1.0.0-rc.165

image_analysis/geom_refinement/PostRefine.cpp

@@ -20,6 +20,20 @@
 
 namespace {
 
+// How far post-refinement's step B 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;
+
 // 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
@@ -693,8 +707,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
                 // above bakes those five blocks in rather than declaring them and freezing them, so there
                 // is nothing left to hold constant here.
                 p.SetParameterLowerBound(dist, 0, dist0 * 0.95); p.SetParameterUpperBound(dist, 0, dist0 * 1.05);
-                for (int j = 0; j < 2; ++j) { p.SetParameterLowerBound(beam, j, beam[j] - 15.0);
-                                              p.SetParameterUpperBound(beam, j, beam[j] + 15.0); }
+                // 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] : beam[j];
+                    p.SetParameterLowerBound(beam, j, std::min(beam[j], m) - BEAM_BOUND_PXL);
+                    p.SetParameterUpperBound(beam, j, std::max(beam[j], m) + BEAM_BOUND_PXL);
+                }
                 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);
@@ -716,15 +736,34 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
                 // degeneracy) far off. 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 the multi-lattice case.
+                const double from_nominal = std::hypot(beam_fit[0] - beam_x0, beam_fit[1] - beam_y0);
+                const double from_measured = settings.measured_beam_px
+                        ? std::hypot(beam_fit[0] - (*settings.measured_beam_px)[0],
+                                     beam_fit[1] - (*settings.measured_beam_px)[1])
+                        : std::numeric_limits<double>::infinity();
                 const bool in_bounds = std::fabs(dist_fit - dist0) < 0.01 * dist0
-                                       && std::hypot(beam_fit[0] - beam_x0, beam_fit[1] - beam_y0) < 15.0;
+                                       && std::min(from_nominal, from_measured) < BEAM_BOUND_PXL;
                 result.detector_refined = convB && cvB_ref < 0.98 * cvB_nom && in_bounds;
                 if (result.detector_refined) { double bo[2]; solve_detector(ALL, bo, dist); beam[0] = bo[0]; beam[1] = bo[1]; }
+                // Name which test refused it. Three different things reject here and the geometry that
+                // comes out is the same in all three, 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 char *verdict =
+                        result.detector_refined      ? "COMMIT"
+                        : !convB                     ? "reject (the fit did not converge; kept nominal detector)"
+                        : !(cvB_ref < 0.98 * cvB_nom) ? "reject (the held-out residual did not improve enough; "
+                                                        "kept nominal detector)"
+                        : std::fabs(dist_fit - dist0) >= 0.01 * dist0
+                                                     ? "reject (the distance moved more than 1 %; kept nominal detector)"
+                                                     : "reject (the beam moved further than the bound from every "
+                                                       "centre anything believes; kept nominal detector)";
                 logger.Info("Post-refine GEOM step B (distance/beam): dist {:.3f} -> {:.3f} mm, beam "
                             "({:.2f},{:.2f}) -> ({:.2f},{:.2f}), held-out pos {:.3e} -> {:.3e} => {}",
                             dist0, result.detector_refined ? dist : dist0, beam_x0, beam_y0,
                             result.detector_refined ? beam[0] : beam_x0, result.detector_refined ? beam[1] : beam_y0,
-                            cvB_nom, cvB_ref, result.detector_refined ? "COMMIT" : "reject (kept nominal detector)");
+                            cvB_nom, cvB_ref, verdict);
             } else {
                 logger.Info("Post-refine GEOM step B: only {} positional observations - skipped", n_obs);
             }

image_analysis/geom_refinement/PostRefine.h

@@ -3,6 +3,8 @@
 
 #pragma once
 
+#include <array>
+#include <optional>
 #include <vector>
 
 #include "../../common/DiffractionGeometry.h"
@@ -52,6 +54,11 @@ struct PostRefineSettings {
     double excitation_weight = 1.0;   // weight of the phi/excitation residual vs the positional one
     int min_events = 50;
     int num_threads = 1;
+    // An independent measurement of the beam centre in pixels, where the run has one (rugnux's pre-scan
+    // fit to the isotropy of the scattered background, which runs on every rotation run). Step B bounds
+    // how far it may move the beam from whichever of this and the nominal centre is NEARER; see
+    // BEAM_BOUND_PXL in PostRefine.cpp for why the nominal centre alone is not enough to bound it.
+    std::optional<std::array<double, 2>> measured_beam_px;
 };
 
 // nominal_geom / reference_latt: the current detector geometry and the phi=0 reference lattice (orientation

rugnux/Rugnux.cpp

@@ -1046,6 +1046,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
     if (want_shadow && config_.beam_center_check) {
         background_center_ = FindBeamCenterFromBackground(experiment_, pixel_mask_,
                                                           finder.GetMeanProjection());
+        measured_beam_center_ = background_center_;
         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 "
@@ -3512,6 +3513,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
                     prs.crystal_system = rot->search_result.system;
                     prs.num_threads = static_cast<int>(config_.nthreads);
                     prs.refine_geometry = true;
+                    // The centre the pre-scan measured, where the fit was precise enough to be used as
+                    // a centre at all (the same ceiling --estimate-beam-center adopts one on). It does
+                    // not move the geometry - it only tells the bound on step B's beam move where else
+                    // the centre is independently believed to be.
+                    if (measured_beam_center_
+                        && measured_beam_center_->sigma_pxl <= std::max(MAX_BEAM_CENTER_SIGMA_PXL,
+                                                                        BeamCenterNeed_pxl(experiment_)))
+                        prs.measured_beam_px = std::array<double, 2>{
+                            static_cast<double>(measured_beam_center_->beam_x_pxl),
+                            static_cast<double>(measured_beam_center_->beam_y_pxl)};
                     const auto pr = PostRefineRotationGeometry(outcomes, *rot->axis,
                                         experiment_.GetDiffractionGeometry(), rot->lattice, prs, logger);
                     result.post_refine = pr;

rugnux/Rugnux.h

@@ -380,6 +380,10 @@ class Rugnux {
     // the second first pass, and clears it, so the two-pass loop does not re-ask on a post-refined
     // geometry the pre-scan estimate says nothing about.
     std::optional<BeamCenterEstimate> background_center_;
+    // The same pre-scan measurement, kept for the whole run instead of being consumed by the first
+    // pass. Post-refinement measures how far it may move the beam from the nearer of the nominal
+    // centre and this one, so a header that is far out does not cap the correction at its own error.
+    std::optional<BeamCenterEstimate> measured_beam_center_;
 
     // The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the
     // two-pass rotation run measures it once and both passes integrate at the same radius.