StillsPartialityRefine: the partiality width is the run's, measured on the crystals, not the prediction window's

The stills partiality width came from a per-crystal least-squares fit of de^2 against d*^2 over the
integrated reflections. Those fill the prediction window evenly, so the fit returned the window's own
width (W^2/3, constant in the linear Ewald distance, gamma_e ~ 0): at low resolution several times the
crystal's rocking width, which mis-corrected every partial there.

The width is now measured once per run from what the reflections RECORD: the intensity-weighted second
moment of the excitation error, sum I de^2 / sum I, which is the width of a Gaussian partiality whatever
the window, in 8 shells of d*^2 over all crystals. Each shell is corrected for the window cutting the
curve's tails (a curve wider than the window is taken as wide as the window), and the shells are fitted
as s^2 = gamma0^2 + (gamma_e d*)^2 + (bw lambda d*^2 / 2)^2 with every term non-negative. The per-crystal
fit stays as the fallback when the run is too small to measure (< 800 reflections). The prediction window
is unchanged (2x the per-frame profile radius).

Measured: on the SwissFEL lysozyme stills the width comes out ~0.035 deg angular at every resolution,
the ~0.03-0.045 deg that the empirical rocking profile of the observations gives; the window-moment
model was 3-4x wider at low resolution.

dimple R-free (deposited free set, matrix subsets/cuts, PARTIALITY_MATRIX.md), shipped -> this:
  SwissFEL run0237 4k     0.2692 -> 0.2454
  same, -C/-S given       0.2679 -> 0.2513
  MicroMAX mono 4k        0.3341 -> 0.3189
  MicroMAX pink 4k        0.2415 -> 0.2442
  SACLA 243 subset        0.2576 -> 0.2454
  SACLA 237 subset        0.2342 -> 0.2331
  LCLS 222 subset         0.2344 -> 0.2330
  PAL 86 subset           0.3499 -> 0.3629  (binned CCD, stencil-limited)
For comparison --bandwidth 0.003 by hand: 0.2616 / - / 0.3108 / 0.2377 / 0.2460 / 0.2270 / 0.2404 / 0.3469.

Stills only: rotation never runs StillsPartialityRefine.

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

tests/MergeScaleTest.cpp

@@ -720,6 +720,7 @@ TEST_CASE("StillsPartialityRefine_RecoversTheTilt") {
                     Reflection r{};
                     r.h = h; r.k = k; r.l = l;
                     r.d = 1.0f / dstar;
+                    r.dist_ewald = std::fabs((q + S0).Length() - 1.0f / WAVELENGTH);
                     r.I = static_cast<float>(G * p_true * full_I(h, k, l));
                     r.sigma = static_cast<float>(1.0 + 0.01 * r.I);
                     r.prescaling_corr = 1.0f;
@@ -745,7 +746,7 @@ TEST_CASE("StillsPartialityRefine_RecoversTheTilt") {
     const double after = mean_abs_error();
     CHECK(mean_tilt > 0.015);
     CHECK(mean_tilt < 0.06);
-    CHECK(after < 0.6 * before);   // measured 0.54: the prior holds the larger tilts back
+    CHECK(after < 0.6 * before);   // measured 0.41 (0.54 with the per-crystal window width): the prior holds the larger tilts back
 }
 
 namespace {
@@ -914,3 +915,44 @@ TEST_CASE("MergeOnTheFly_RsplitKeepsEveryFrameWhole") {
             CHECK(sh.r_split == Catch::Approx(expected).epsilon(1e-5));
     CHECK_FALSE(SerialShellTable(stats).empty());
 }
+
+// The run's rocking width is measured on the crystals, not on the prediction window: reflections spread
+// evenly across a window four times wider than the curve, each recording I_full * p(de), give back the
+// curve's gamma0 and angular width, while the plain moment of de would return the window's.
+TEST_CASE("StillsPartialityRefine_MeasureRunWidth") {
+    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(5);
+    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(200);
+    for (auto &o: outcomes) {
+        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 = 4.0 * std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E / (1.5 * 1.5));
+            const double de = window * (2.0 * unit(rng) - 1.0);
+            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) * std::exp(-0.5 * de * de / (s * s))
+                                     + 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.1));
+    CHECK(width->bandwidth_sigma < 1e-3);
+}