Stills scaling: a model error in the per-crystal scale fit, and the rocking width measured per window group

Two changes to the stills scale model; rotation is untouched (p.hkl md5 identical on two rotation sets).

The per-crystal scale fits (ScaleOnTheFly and the partiality post-refinement's G and tilt) measured
residuals in counting sigmas. The strong reflections of a still scatter about the model by tens of
percent (robust sd of ln(I/<I>) 0.45-0.8 by shell on one run), so with a Cauchy loss at 3 counting
sigma every strong reflection was an outlier and the scale was set by the weak ones: re-fitting each
frame's scale against the final merge moved it by a factor 0.45-1.5 (5-95%). The fit's variance is now
the expected one at the model, background + counting + (0.3 * model)^2. The coefficient is not
delicate: 0.1-1.0 fit alike in an offline replay; what matters is that the term is there.

The run's rocking width came from the intensity-weighted moment of the excitation error over all
frames pooled. Frames predict through windows of their own that differ by an order of magnitude, and
the widest windows dominated both moments: on a run whose windows are mostly narrower than the curve
at high resolution the measurement read the window at every resolution, 2.8e-4 1/A with no angular
term, several times too wide at low resolution (low-partiality observations there merged at 0.1-0.3
of their reflection's mean). The moment is now taken per group of frames of similar window (4 per
shell), each against its own window, and the shell takes the median of the groups whose curve is
narrower than their window. On that run the width becomes 6.7e-5 1/A + 0.018 deg angular; where the
pooled measurement was already right (a SASE-bandwidth run), it stays where it was.

0-cycle REFMAC R-free against the deposited models, 9 still battery sets (6 used while developing, 3
held out), change from rc179: -0.0100, -0.0095, -0.0057, -0.0055, -0.0051, -0.0045, -0.0030,
-0.0024, +0.0009 (held out: -0.0055, -0.0051, +0.0009). The width change alone accounts for most of
the largest gain; the scale-fit variance for 0.000-0.003 per set.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit ccc0e01f56827bbe0648ed7c5b8049912d80650f)

tests/MergeScaleTest.cpp

@@ -957,6 +957,49 @@ TEST_CASE("StillsPartialityRefine_MeasureRunWidth") {
     CHECK(width->bandwidth_sigma < 1e-3);
 }
 
+// Frames predicted through a much wider window whose intensities do not follow the curve (a wrong lattice
+// on the frame) must not decide the measured width: pooled with the others they hold most of the de^2, and
+// the run's width then read as their window. Measured per group of similar window, they are a group whose
+// curve looks wider than its window, which is left out.
+TEST_CASE("StillsPartialityRefine_MeasureRunWidth_WideWindowFramesIgnored") {
+    constexpr double GAMMA0 = 1.0e-4, GAMMA_E = 0.035 * M_PI / 180.0;
+    DiffractionExperiment x;
+    x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
+
+    std::mt19937 rng(7);
+    std::uniform_real_distribution<double> unit(0.0, 1.0);
+    std::exponential_distribution<double> wilson(1.0);
+    std::normal_distribution<double> gauss(0.0, 1.0);
+    std::vector<IntegrationOutcome> outcomes(400);
+    for (size_t f = 0; f < outcomes.size(); ++f) {
+        auto &o = outcomes[f];
+        const bool wide = f % 4 == 3;
+        o.image_scale_g = 0.5 + unit(rng);
+        for (int j = 0; j < 200; ++j) {
+            const double d = 1.5 + 8.0 * unit(rng);
+            const double dstar = 1.0 / d;
+            const double s = std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E * dstar * dstar);
+            const double window = (wide ? 20.0 : 4.0) * std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E / (1.5 * 1.5));
+            const double de = window * (2.0 * unit(rng) - 1.0);
+            const double p = wide ? 1.0 : std::exp(-0.5 * de * de / (s * s));
+            Reflection r{};
+            r.h = j; r.k = 1; r.l = 1;
+            r.d = static_cast<float>(d);
+            r.dist_ewald = static_cast<float>(std::fabs(de));
+            r.I = static_cast<float>(*o.image_scale_g * 1000.0 * wilson(rng) * p + 5.0 * gauss(rng));
+            r.sigma = 5.0f;
+            r.prescaling_corr = 1.0f;
+            r.partiality = 1.0f;
+            r.image_scale_corr = static_cast<float>(1.0 / *o.image_scale_g);
+            o.reflections.push_back(r);
+        }
+    }
+    const auto width = StillsPartialityRefine(x).MeasureRunWidth(outcomes, std::nullopt);
+    REQUIRE(width.has_value());
+    CHECK(width->gamma0 == Catch::Approx(GAMMA0).epsilon(0.3));
+    CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.15));
+}
+
 // Outlier rejection on the stills merge: one observation at ten times its reflection's intensity among
 // nineteen consistent ones is dropped, and nothing else is.
 TEST_CASE("MergeOnTheFly_RejectsAnOutlierInTheBand") {