MergeOnTheFly: capture term and multiplicative outlier band for the stills merge (both opt-in)

--capture-uncertainty now reaches the stills merge: each partial carries (c (1 - p) <I>)^2 in its
variance, the rotation combine's capture term per observation, in the merge weights and in the
error-model samples alike. --reject-outliers on stills now uses the rotation merge's band
(OutlierBand.h): exp(+-n s) about the reflection mean for the systematic part, s the weighted median
|ln(I / median)| of the merge's own observations, instead of a symmetric +-n sigma band that reaches
below zero when b is large (stills b 1.4-1.7). A median needs three observations.

Neither is switched on for stills: measured with the run-measured partiality width before its
window-truncation correction (dimple R-free, SwissFEL 4k / SACLA 243 / MicroMAX mono / MicroMAX pink):
no options 0.2567 / 0.2497 / 0.3279 / 0.2477, --capture-uncertainty 1 0.2510 / 0.2470 / 0.3261 /
0.2527, --reject-outliers 6 0.2619 / 0.2487 / 0.3317 / 0.2630 - mixed, and the photon-counting pink
set loses with both. Raising --min-partiality to 0.1 or 0.2 lost 0.008-0.019 on SwissFEL and mono.
Defaults unchanged, so stills output is unchanged; rotation does not use MergeOnTheFly.

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

tests/MergeScaleTest.cpp

@@ -956,3 +956,75 @@ TEST_CASE("StillsPartialityRefine_MeasureRunWidth") {
     CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.1));
     CHECK(width->bandwidth_sigma < 1e-3);
 }
+
+// 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") {
+    DiffractionExperiment x;
+    x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
+    ScalingSettings s = x.GetScalingSettings();
+    s.OutlierRejectNsigma(6.0);
+    x.ImportScalingSettings(s);
+
+    std::mt19937 rng(3);
+    std::normal_distribution<float> gauss(0.0f, 1.0f);
+    std::vector<IntegrationOutcome> outcomes(20);
+    for (size_t i = 0; i < outcomes.size(); ++i) {
+        for (int h = 1; h <= 50; ++h) {
+            Reflection r{};
+            r.h = h; r.k = 2; r.l = 3;
+            r.d = 2.0f;
+            const float I_true = 100.0f * static_cast<float>(h);
+            r.I = I_true * (1.0f + 0.1f * gauss(rng)) + 3.0f * gauss(rng);
+            if (i == 7 && h == 25)
+                r.I = 10.0f * I_true;
+            r.sigma = 3.0f;
+            r.prescaling_corr = 1.0f;
+            r.partiality = 1.0f;
+            r.image_scale_corr = 1.0f;
+            outcomes[i].reflections.push_back(r);
+        }
+    }
+    MergeOnTheFly merge(x);
+    merge.RefineErrorModel(outcomes);
+    REQUIRE(merge.ErrorModelActive());
+    merge.AddImages(outcomes);
+    CHECK(merge.RejectedCount() == 1);
+}
+
+// The capture term widens the sigma of a partial by (c (1 - p) <I>): a partial merges at a lower weight
+// than a full of the same intensity, and not at all differently with c = 0.
+TEST_CASE("MergeOnTheFly_CaptureTermDownweightsPartials") {
+    auto merged_I = [](double c) {
+        DiffractionExperiment x;
+        x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
+        ScalingSettings s = x.GetScalingSettings();
+        s.CaptureUncertaintyCoeff(c);
+        x.ImportScalingSettings(s);
+        std::mt19937 rng(9);
+        std::normal_distribution<float> gauss(0.0f, 1.0f);
+        std::vector<IntegrationOutcome> outcomes(40);
+        for (size_t i = 0; i < outcomes.size(); ++i)
+            for (int h = 1; h <= 40; ++h) {
+                Reflection r{};
+                r.h = h; r.k = 1; r.l = 1;
+                r.d = 2.0f;
+                // Every other image records the reflection as a 30 % partial that reads 20 % high.
+                const bool partial = i % 2 == 1;
+                r.partiality = partial ? 0.3f : 1.0f;
+                r.image_scale_corr = 1.0f / r.partiality;
+                r.I = r.partiality * 100.0f * static_cast<float>(h) * (partial ? 1.2f : 1.0f) * (1.0f + 0.05f * gauss(rng));
+                r.sigma = 2.0f;
+                r.prescaling_corr = 1.0f;
+                outcomes[i].reflections.push_back(r);
+            }
+        MergeOnTheFly merge(x);
+        merge.RefineErrorModel(outcomes);
+        merge.AddImages(outcomes);
+        const auto m = merge.ExportReflections();
+        return m.at(10).I / (100.0f * static_cast<float>(m.at(10).h));
+    };
+    const float without = merged_I(0.0), with = merged_I(1.0);
+    CHECK(with < without);
+    CHECK(with < 1.08f);
+}