From c140b3567da00c8e9634783167202a625b54a93c Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Sat, 18 Jul 2026 17:31:45 +0200
Subject: [PATCH] rugnux: restrain the gauge-weak beam centre toward the header
 in geometry refinement

Single-axis rotation has a gauge degeneracy: rotating the whole experiment
about the spindle leaves every spot position unchanged, so the beam-centre
component parallel to the spindle is a null / under-determined direction.
Refining it freely lets it wander a few pixels and absorb centroid
systematics into a wrong beam that the co-refined orientation keeps
position-consistent.

Add a soft header prior on that one component in the shared XtalOptimizer:
residual = w*(beam[parallel] - header), with w set so the restraint behaves
like a sigma_px-pixel prior competing with the positional residuals. The
gauge direction has ~zero data sensitivity, so the prior pins it near the
LaB6-monitored header, while a real, well-supported drift can still overcome
it (the beam does drift). One restraint covers all three refinement sites at
once: the noisy accumulated-spot primary, the per-frame filter, and
post-refine.

Validated against independent XDS processing: neutral on well-determined
data (the gauge direction has no effect there), and a real low-resolution
accuracy gain where the free refinement had otherwise wandered off the
header - the merged intensities move back into agreement with XDS while the
internal CC1/2, which is blind to this, barely changes. No space-group
determinations change.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 .../geom_refinement/XtalOptimizer.cpp         | 38 +++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp
index acf4071f..45774851 100644
--- a/image_analysis/geom_refinement/XtalOptimizer.cpp
+++ b/image_analysis/geom_refinement/XtalOptimizer.cpp
@@ -10,6 +10,21 @@
 #include "ceres/rotation.h"
 #include "LatticeReduction.h"
 
+// Soft header prior on ONE beam-centre component (the spindle-parallel, gauge-weak one). Residual = w*(b - b0);
+// the caller sets w so the prior behaves like a sigma-pixel restraint that competes with the (unit-weight)
+// positional residuals - strong enough to pin the gauge direction, negligible in the well-constrained one.
+struct BeamComponentPrior {
+    BeamComponentPrior(int component, double b0, double weight)
+        : component(component), b0(b0), weight(weight) {}
+    template<typename T>
+    bool operator()(const T *const beam, T *residual) const {
+        residual[0] = T(weight) * (beam[component] - T(b0));
+        return true;
+    }
+    int component;
+    double b0, weight;
+};
+
 struct XtalResidualRotationOnlyPrecomp {
     XtalResidualRotationOnlyPrecomp(const Coord &recip_obs,
                                     const CrystalLattice &latt,
@@ -199,6 +214,29 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
 
         if (!data.refine_beam_center)
             problem.SetParameterBlockConstant(beam);
+        else if (data.axis) {
+            // Gauge handling (single-axis rotation): rotating the whole experiment about the spindle leaves every
+            // spot position unchanged, so the beam-centre component PARALLEL to the spindle is a null/gauge-weak
+            // direction. Refining it freely lets it wander (~+3 px) and absorb centroid systematics into a wrong
+            // beam that the co-refined orientation keeps position-consistent. Rather than freeze it (the beam
+            // does drift - it is only LaB6-monitored to ~a few px), RESTRAIN it toward the header with a soft
+            // prior: the gauge direction has ~zero data sensitivity so the prior pins it near the header, while a
+            // real, well-supported drift can still overcome it. The spindle is along a detector axis in standard
+            // geometry, so restrain the dominant of X / Y.
+            const Coord spindle = data.axis->GetAxis();
+            const int parallel = (std::fabs(spindle.x) >= std::fabs(spindle.y)) ? 0 : 1;
+            // Weight so the prior is a sigma_px-pixel restraint that competes with the positional residuals.
+            // k = d|recip|/d(beam_px) ~ pixel/(distance*lambda) [A^-1/px]; scaling by sqrt(#residuals) makes the
+            // prior's curvature ~ (1/9) of the well-constrained-data curvature at sigma_px=3, i.e. data wins the
+            // perpendicular direction, the prior wins the gauge one.
+            constexpr double sigma_px = 3.0;
+            const double k = data.geom.GetPixelSize_mm() / (distance_mm * data.geom.GetWavelength_A());
+            const double w = k * std::sqrt(static_cast<double>(problem.NumResidualBlocks())) / sigma_px;
+            problem.AddResidualBlock(
+                new ceres::AutoDiffCostFunction<BeamComponentPrior, 1, 2>(
+                    new BeamComponentPrior(parallel, beam[parallel], w)),
+                nullptr, beam);
+        }
 
         if (!data.refine_detector_angles) {
             problem.SetParameterBlockConstant(detector_rot);