Stills partiality: merge only partials of p >= 0.5, adopt the refined model on every crystal, centre the rocking curve on the run's radial scale

Three changes to the stills post-refinement and merge, decided on 0-cycle REFMAC R-free against the
deposited models (diff = ours - deposited data, noise ~0.002):

- Stills merge only observations whose model partiality is at least 0.5 (RugnuxDefaultScalingSettings;
  rotation keeps 0.02, and --min-partiality still overrides). Diagnosis on per-observation dumps: the
  empirical rocking curve has a floor of 5-15 % of its peak out on the tail (mis-oriented crystals), and
  dividing such observations by a small model p makes them biased, not noisy. Variance terms in (1-p) or
  (1-p)/p recovered a third of what the cut does; a narrower (correct-looking angular) width or a
  Lorentzian form was worse. On one set in scale mode the optimum is flat over 0.3-0.6.
- The CC adoption gate in RefineOne is gone: it left 25-33 % of crystals at p = 1 next to the refined ones,
  and collapsed to 3 % on a run whose curve is off-centre.
- The run's radial scale of the reciprocal lattice (a wavelength or cell a little off, invisible to spot
  positions since they fix only lambda d*) is measured from where the curve is centred against d*^2
  (MeasureRadialScale, iterated three times) and applied to q in the width measurement and the
  post-refinement. One run had its curve centred at +2-3e-3 d*^2 (scale -1.6e-3 measured); the others
  measure below 7e-4.

--mode scale on stills now starts where the full pipeline starts: p = 1, the excitation error recomputed
from the stored lattice (the window was read as zero before, so no run width was ever measured), and
collapsed scales rejected after the refinement.

Still battery, the 12 non-float sets, change of diff vs the previous build: -0.0035 mean, 9 better, three
within +-0.0006 except one +0.0053 (an edge win becoming a draw); the off-centre run -0.005 here and -0.017
in an earlier run of the same algorithm (that set varies by ~0.01 between runs). Rotation p.hkl
byte-identical on three sets. Test StillsPartialityRefine_MeasureRadialScale added.

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

image_analysis/scale_merge/StillsPartialityRefine.cpp

@@ -253,7 +253,8 @@ StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x)
 
 double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome,
                                          const std::map<HKLKey, double> &reference,
-                                         const std::optional<RunWidth> &run_width) const {
+                                         const std::optional<RunWidth> &run_width,
+                                         double radial_scale) const {
     if (outcome.reflections.empty())
         return 0.0;
 
@@ -264,8 +265,8 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome,
     const double inv_lambda = 1.0 / outcome.geom.GetWavelength_A();
 
     auto base_q = [&](const Reflection &r) {
-        return Astar * static_cast<float>(r.h) + Bstar * static_cast<float>(r.k)
-               + Cstar * static_cast<float>(r.l);
+        return (Astar * static_cast<float>(r.h) + Bstar * static_cast<float>(r.k)
+                + Cstar * static_cast<float>(r.l)) * static_cast<float>(radial_scale);
     };
 
     // Collect the reflections that constrain the fit (accepted, non-ice, finite, present in the reference).
@@ -378,16 +379,6 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome,
             break;
     }
 
-    // Keep what the crystal came in with. This refinement is on by default, so a crystal the tilt model
-    // happens to suit WORSE than the fixed partiality it replaces must not be made worse by it - and
-    // whether it suits is only known once the corrections are written and the CC re-measured.
-    const auto cc_before = outcome.image_scale_cc;
-    const auto g_before = outcome.image_scale_g;
-    std::vector<std::pair<float, float>> before;   // partiality, image_scale_corr
-    before.reserve(outcome.reflections.size());
-    for (const auto &r: outcome.reflections)
-        before.emplace_back(r.partiality, r.image_scale_corr);
-
     // Write the refined partiality + scale correction onto every reflection of the crystal (not only the
     // fit subset), so the merge sees a consistent model. image_scale_corr = prescaling_corr * qe_corr * flight_corr /
     // (partiality * G), the same composition ScaleOnTheFly writes.
@@ -410,20 +401,6 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome,
     const auto [cc, cc_n] = ImageReferenceCC(outcome.reflections, reference, hkl_key_generator_,
                                              d_min_limit_, d_max_limit_, min_partiality_);
 
-    // Adopt the refined model only if it correlates with the reference at least as well as the model it
-    // replaces. Rejecting puts the crystal back exactly as it arrived, which is the same state a
-    // crystal with too few reflections to fit ends in. A CC that cannot be measured at all counts as
-    // worse, not as no opinion: ImageReferenceCC returns NaN when too few reflections clear the
-    // partiality cut, which is precisely what a refinement that collapsed the partialities produces -
-    // and adopting it would leave image_scale_cc NaN, on which --min-image-cc drops the image outright.
-    if (cc_before.has_value() && std::isfinite(*cc_before) && !(std::isfinite(cc) && cc >= *cc_before)) {
-        for (size_t i = 0; i < outcome.reflections.size(); ++i) {
-            outcome.reflections[i].partiality = before[i].first;
-            outcome.reflections[i].image_scale_corr = before[i].second;
-        }
-        outcome.image_scale_g = g_before;
-        return 0.0;   // nothing adopted, so no tilt to report
-    }
     outcome.image_scale_cc = cc;
     outcome.image_scale_cc_n = cc_n;
 
@@ -431,8 +408,67 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome,
     return tilt_deg;
 }
 
+double StillsPartialityRefine::MeasureRadialScale(const std::vector<IntegrationOutcome> &outcomes,
+                                                  const std::optional<UnitCell> &reference_cell,
+                                                  double radial_scale) const {
+    // A scale of the reciprocal lattice - a wavelength, or a cell, a little off - leaves every spot where
+    // it was (positions fix only lambda d*), but moves each reflection's excitation error by
+    // eps lambda d*^2 / 2 for q -> (1 + eps) q, so the run's rocking curve sits off zero by an amount
+    // growing as d*^2. The prediction window is symmetric in de, so the intensity-weighted mean de of a
+    // shell is that offset, and a line through the origin in lambda d*^2 / 2 gives eps. A window that cuts
+    // an off-centre curve pulls the mean toward zero, so the caller iterates from the scale found.
+    const auto &indexing = experiment_.GetIndexingSettings();
+    struct Point { double x, de, w; };
+    std::vector<Point> pts;
+    double max_x = 0.0;
+    for (const auto &o: outcomes) {
+        if (!o.image_scale_g || !(*o.image_scale_g > 0.0))
+            continue;
+        if (reference_cell && !o.latt.GetUnitCell().is_close(*reference_cell, indexing.GetUnitCellDistTolerance(),
+                                                               indexing.GetUnitCellAngleTolerance_deg()))
+            continue;
+        const Coord S0 = o.geom.GetScatteringVector();
+        const double inv_lambda = 1.0 / o.geom.GetWavelength_A();
+        for (const auto &r: o.reflections) {
+            if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_) || !(r.d > 0.0f))
+                continue;
+            const double w = r.I * r.prescaling_corr * r.qe_corr * r.flight_corr / *o.image_scale_g;
+            if (!std::isfinite(w))
+                continue;
+            const Coord q = (o.latt.Astar() * static_cast<float>(r.h) + o.latt.Bstar() * static_cast<float>(r.k)
+                             + o.latt.Cstar() * static_cast<float>(r.l)) * static_cast<float>(radial_scale);
+            const double x = 0.5 / (inv_lambda * r.d * r.d);
+            pts.push_back({x, VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - inv_lambda, w});
+            max_x = std::max(max_x, x);
+        }
+    }
+    constexpr int N_SHELLS = 8;
+    if (pts.size() < 100 * N_SHELLS)
+        return radial_scale;
+    std::array<double, N_SHELLS> sw{}, swde{}, swx{};
+    for (const auto &p: pts) {
+        const int b = std::min(N_SHELLS - 1, static_cast<int>(p.x / max_x * N_SHELLS));
+        sw[b] += p.w;
+        swde[b] += p.w * p.de;
+        swx[b] += p.w * p.x;
+    }
+    double sxy = 0.0, sxx = 0.0;
+    for (int b = 0; b < N_SHELLS; ++b) {
+        if (!(sw[b] > 0.0))
+            continue;
+        const double x = swx[b] / sw[b], y = swde[b] / sw[b];
+        sxy += sw[b] * x * y;
+        sxx += sw[b] * x * x;
+    }
+    if (!(sxx > 0.0) || !std::isfinite(sxy))
+        return radial_scale;
+    // the slope is eps, removed by scaling q by (1 - eps)
+    return radial_scale * (1.0 - sxy / sxx);
+}
+
 std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureRunWidth(
-        const std::vector<IntegrationOutcome> &outcomes, const std::optional<UnitCell> &reference_cell) const {
+        const std::vector<IntegrationOutcome> &outcomes, const std::optional<UnitCell> &reference_cell,
+        double radial_scale) const {
     // The rocking width measured on the crystals rather than on the prediction window. A lattice point
     // sits anywhere across the window with equal likelihood, while what it records is I_full * p(de), so
     // the intensity-weighted second moment of the excitation error, sum I de^2 / sum I, is the width s^2
@@ -470,8 +506,16 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
             const double w = r.I * corr, sigma_w = r.sigma * corr;
             if (!std::isfinite(w) || !std::isfinite(sigma_w) || !std::isfinite(r.dist_ewald) || !(r.d > 0.0f))
                 continue;
-            pts.push_back({1.0 / (static_cast<double>(r.d) * r.d),
-                           static_cast<double>(r.dist_ewald) * r.dist_ewald, w, sigma_w * sigma_w, window});
+            // The excitation error the window was cut in is the prediction's; the curve's is taken at the
+            // run's radial scale (MeasureRadialScale), recomputed from the lattice when that is not one.
+            double de = r.dist_ewald;
+            if (radial_scale != 1.0) {
+                const Coord q = (o.latt.Astar() * static_cast<float>(r.h) + o.latt.Bstar() * static_cast<float>(r.k)
+                                 + o.latt.Cstar() * static_cast<float>(r.l)) * static_cast<float>(radial_scale);
+                const Coord S0 = o.geom.GetScatteringVector();
+                de = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - 1.0 / o.geom.GetWavelength_A();
+            }
+            pts.push_back({1.0 / (static_cast<double>(r.d) * r.d), de * de, w, sigma_w * sigma_w, window});
         }
     }
     constexpr size_t N_SHELLS = 8;
@@ -580,7 +624,12 @@ double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, si
     nthreads = std::max<size_t>(1, nthreads);
 
     double last_mean_tilt = 0.0;
-    const auto run_width = MeasureRunWidth(outcomes, reference_cell);
+    double radial_scale = 1.0;
+    for (int i = 0; i < 3; ++i)
+        radial_scale = MeasureRadialScale(outcomes, reference_cell, radial_scale);
+    Logger("StillsPartialityRefine").Info("Stills rocking curve centre: reciprocal lattice scaled by {:+.2e}",
+                                          radial_scale - 1.0);
+    const auto run_width = MeasureRunWidth(outcomes, reference_cell, radial_scale);
 
     for (int outer = 0; outer < settings_.outer_iterations; ++outer) {
         // Reference full intensities from the current corrections.
@@ -596,7 +645,7 @@ double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, si
         auto worker = [&]() {
             size_t i = next.fetch_add(1);
             while (i < outcomes.size()) {
-                const double t = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width) : 0.0;
+                const double t = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width, radial_scale) : 0.0;
                 if (t > 0.0) {
                     double prev = tilt_sum.load();
                     while (!tilt_sum.compare_exchange_weak(prev, prev + t)) {}

image_analysis/scale_merge/StillsPartialityRefine.h

@@ -58,7 +58,14 @@ public:
         double gamma0, gamma_e, bandwidth_sigma;
     };
     std::optional<RunWidth> MeasureRunWidth(const std::vector<IntegrationOutcome> &outcomes,
-                                            const std::optional<UnitCell> &reference_cell) const;
+                                            const std::optional<UnitCell> &reference_cell,
+                                            double radial_scale = 1.0) const;
+
+    // The run's radial scale of the reciprocal lattice (a wavelength or cell scale the spot positions cannot
+    // see), measured from where the rocking curve is centred against d*^2: one step from radial_scale, so
+    // call it a few times. Returns radial_scale unchanged when the run has too few reflections.
+    double MeasureRadialScale(const std::vector<IntegrationOutcome> &outcomes,
+                              const std::optional<UnitCell> &reference_cell, double radial_scale) const;
 
 private:
     const DiffractionExperiment experiment_;
@@ -72,5 +79,5 @@ private:
     // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped).
     double RefineOne(IntegrationOutcome &outcome,
                      const std::map<HKLKey, double> &reference,
-                     const std::optional<RunWidth> &run_width) const;
+                     const std::optional<RunWidth> &run_width, double radial_scale) const;
 };

rugnux/RugnuxDefaults.cpp

@@ -52,6 +52,13 @@ ScalingSettings RugnuxDefaultScalingSettings(bool rotation) {
     // I/sigma floor and falls back to P1.
     s.MinCapturedFraction(rotation ? 0.7 : 0.0);
 
+    // A still's intensity is divided by the partiality its excitation error gives. Below about one half
+    // the error of that model (the run's curve is an average over crystals, and mis-oriented crystals
+    // leave a floor far out on its tail) dominates, and dividing by a small p turns it into a bias that
+    // no variance term weights away. So stills merge only the reflections the model calls at least half
+    // recorded; rotation sums its partials into fulls and keeps the permissive floor.
+    s.MinPartiality(rotation ? 0.02 : 0.5);
+
     // Capture-aware systematic sigma: down-weights the over-extrapolated under-captured fulls.
     s.CaptureUncertaintyCoeff(rotation ? 1.0 : 0.0);
 

rugnux/rugnux_cli.cpp

@@ -229,7 +229,7 @@ void print_usage() {
     std::cout << "   --report-resolution <dmin>[,<dmax>] Also report the merging statistics over this resolution range (A, either order; dmax defaults to the run's own low-resolution limit) as a second table beside the run's own - the REFRES_* keys of the report - binned from the same merged reflections, so a run can be read against another program's table at that program's range. Report-only: the cut, the scaling, the error model and every decision stay the run's own, and shells finer than the run's own limit are empty rather than filled from data the run did not keep" << std::endl;
     std::cout << "   --ice-min-score <num>           Ice-presence gate: measured ice score (1 = no ice) a run must reach before ANY ice handling is applied - the flagging and the exclusion from scaling (default: 1.5). The eleven hexagonal bands cover 16-26% of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. 0 = no gate (always handle ice)" << std::endl;
     std::cout << "   --ice-min-spot-ratio <num>      Second ice-presence channel: found spots on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so --ice-min-score alone is blind to it. Default 2.0; 0 disables this channel" << std::endl;
-    std::cout << "   --min-partiality <num>          Minimum partiality to accept reflection (default: 0.02)" << std::endl;
+    std::cout << "   --min-partiality <num>          Minimum partiality to accept reflection (default: 0.02 for rotation, 0.5 for stills)" << std::endl;
     std::cout << "   --capture-uncertainty <num>     Systematic sigma num*(1-captured fraction)*<I>: on rot3d fulls, the share of the rocking curve not captured; on stills, the share of each partial its partiality does not see (default: 1.0 for rot3d, 0 otherwise)" << std::endl;
     std::cout << "   --min-captured-fraction <num>   rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << std::endl;
     std::cout << "   --mosaicity <num>               Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
@@ -1040,7 +1040,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     bool export_unmerged = true;            // on by default; --no-export-unmerged turns it off
     bool export_unmerged_partials = false;  // --export-unmerged-partials: the same, one row per image
     std::string dump_observations;   // diagnostic: dump unmerged -P rot3d fulls to this path
-    double min_partiality = 0.02;
+    std::optional<double> min_partiality_arg;   // --min-partiality; default per mode (RugnuxDefaultScalingSettings)
     std::optional<double> min_captured_fraction_arg;   // explicit --min-captured-fraction; default depends on rotation
     std::optional<double> capture_uncertainty_arg;   // explicit --capture-uncertainty; default depends on rot3d
     std::optional<double> forced_mosaicity_arg;      // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed
@@ -1468,7 +1468,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             no_scaling_corrections = true;
             break;
         case OPT_MIN_PARTIALITY:
-            min_partiality = parse_double_arg(optarg, "--min-partiality", logger);
+            min_partiality_arg = parse_double_arg(optarg, "--min-partiality", logger);
             break;
         case OPT_CAPTURE_UNCERTAINTY:
             capture_uncertainty_arg = parse_double_arg(optarg, "--capture-uncertainty", logger);
@@ -2260,7 +2260,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 ParseReportResolution(*report_resolution_arg, scaling_settings, logger));
         // Rotation merges decide with Friedel mates merged, -A or not (see the Rugnux constructor).
         scaling_settings.MergeFriedel(!anomalous_mode || rot);
-        scaling_settings.MinPartiality(min_partiality);
+        scaling_settings.MinPartiality(min_partiality_arg.value_or(scaling_settings.GetMinPartiality()));
         scaling_settings.MinCapturedFraction(
             min_captured_fraction_arg.value_or(scaling_settings.GetMinCapturedFraction()));
         scaling_settings.CaptureUncertaintyCoeff(
@@ -2468,6 +2468,20 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             // against a foreign dataset injects cross-dataset systematics and is a worse reference than
             // the data's own merge. The per-image scale G is the exact one-pass solution, so one pass
             // (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills).
+            // The file holds each still as it was merged - with its post-refined partiality - and not the
+            // excitation errors its prediction made. Put every reflection back where the prediction left it,
+            // at p = 1 and with the excitation error of the stored lattice and geometry, so that the scaling
+            // and the post-refinement start from the state they start from in the full pipeline.
+            for (auto &o: reflections) {
+                const Coord S0 = o.geom.GetScatteringVector();
+                const double inv_lambda = 1.0 / o.geom.GetWavelength_A();
+                for (auto &r: o.reflections) {
+                    const Coord S = o.latt.Astar() * static_cast<float>(r.h) + o.latt.Bstar() * static_cast<float>(r.k)
+                                    + o.latt.Cstar() * static_cast<float>(r.l) + S0;
+                    r.dist_ewald = static_cast<float>(std::fabs(std::sqrt(S.x * S.x + S.y * S.y + S.z * S.z) - inv_lambda));
+                    r.partiality = 1.0f;
+                }
+            }
             ScaleOnTheFly(experiment, MergeAll(experiment, reflections, nthreads, experiment.GetUnitCell())).Scale(reflections, nthreads);
             // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
             // orientation tilt against the merge and recompute each reflection's partiality + scale correction
@@ -2476,6 +2490,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 StillsPartialityRefine refiner(experiment);
                 const double mean_tilt = refiner.Run(reflections, nthreads, experiment.GetUnitCell());
                 logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
+                // As in the full pipeline: the refine fits its own per-image scale, with no floor.
+                ScaleOnTheFly::RejectCollapsedScales(reflections);
             }
             MergeOnTheFly merge_engine(experiment);
             merge_engine.Threads(nthreads);
@@ -3273,7 +3289,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         scaling_settings.ReportResolutionRange_A(
             ParseReportResolution(*report_resolution_arg, scaling_settings, logger));
     scaling_settings.MergeFriedel(!anomalous_mode);
-    scaling_settings.MinPartiality(min_partiality);
+    scaling_settings.MinPartiality(min_partiality_arg.value_or(scaling_settings.GetMinPartiality()));
     // Drop edge-of-sweep truncated fulls (rocking curve captured < this fraction) from the rot3d combine.
     // Defaults ON (0.7) for rotation - removes the low-capture fulls that inflate low-res R-meas and
     // slightly bias accuracy; off for non-rot3d (no combine). 0.7 (rather than 0.5) also strips the

tests/MergeScaleTest.cpp

@@ -789,6 +789,78 @@ namespace {
     }
 }
 
+// A reciprocal lattice scaled by (1 + eps) - a wavelength or cell a little off, which the spot positions
+// cannot see - centres every crystal's rocking curve at eps lambda d*^2 / 2 instead of zero. The windows are
+// cut symmetrically about the stored lattice's sphere, as prediction cuts them; the run's radial scale has
+// to come back from where the recorded intensity sits within them.
+TEST_CASE("StillsPartialityRefine_MeasureRadialScale") {
+    constexpr float CELL = 50.0f, WAVELENGTH = 1.0f, D_MIN = 2.0f;
+    constexpr double EPS = 4e-3, GAMMA0 = 2e-4, WINDOW = 1e-3;
+    DiffractionExperiment x;
+    x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1"));
+    DiffractionGeometry geom;
+    geom.Wavelength_A(WAVELENGTH);
+    const Coord S0 = geom.GetScatteringVector();
+
+    std::mt19937 rng(13);
+    std::normal_distribution<double> gauss(0.0, 1.0);
+    std::uniform_real_distribution<double> unit(0.0, 1.0);
+    std::exponential_distribution<double> wilson(1.0);
+    std::vector<IntegrationOutcome> outcomes(100);
+    const int hmax = static_cast<int>(CELL / D_MIN);
+    for (auto &o: outcomes) {
+        const auto axes = TurnedAxes(UnitCell(CELL, CELL, CELL, 90, 90, 90));
+        // a different orientation per crystal: turn the turned axes about a random axis
+        const double ax = gauss(rng), ay = gauss(rng), az = gauss(rng), n = std::sqrt(ax * ax + ay * ay + az * az);
+        const double angle = 2.0 * M_PI * unit(rng), c = std::cos(angle), s = std::sin(angle);
+        auto turn = [&](const Coord &v) {
+            const double k[3] = {ax / n, ay / n, az / n};
+            const double dot = k[0] * v.x + k[1] * v.y + k[2] * v.z;
+            const double cr[3] = {k[1] * v.z - k[2] * v.y, k[2] * v.x - k[0] * v.z, k[0] * v.y - k[1] * v.x};
+            return Coord(static_cast<float>(v.x * c + cr[0] * s + k[0] * dot * (1 - c)),
+                         static_cast<float>(v.y * c + cr[1] * s + k[1] * dot * (1 - c)),
+                         static_cast<float>(v.z * c + cr[2] * s + k[2] * dot * (1 - c)));
+        };
+        o.geom = geom;
+        o.latt = CrystalLattice(turn(axes[0]), turn(axes[1]), turn(axes[2]));
+        o.image_scale_g = 0.5 + unit(rng);
+        for (int h = -hmax; h <= hmax; ++h)
+            for (int k = -hmax; k <= hmax; ++k)
+                for (int l = -hmax; l <= hmax; ++l) {
+                    const Coord q = o.latt.Astar() * h + o.latt.Bstar() * k + o.latt.Cstar() * l;
+                    const float dstar = q.Length();
+                    if (dstar == 0.0f || dstar > 1.0f / D_MIN)
+                        continue;
+                    const double de_stored = (q + S0).Length() - 1.0 / WAVELENGTH;
+                    if (std::fabs(de_stored) > WINDOW)
+                        continue;
+                    const double de_true = (q * static_cast<float>(1.0 + EPS) + S0).Length() - 1.0 / WAVELENGTH;
+                    Reflection r{};
+                    r.h = h; r.k = k; r.l = l;
+                    r.d = 1.0f / dstar;
+                    r.dist_ewald = static_cast<float>(std::fabs(de_stored));
+                    r.I = static_cast<float>(*o.image_scale_g * 1000.0 * wilson(rng)
+                                             * std::exp(-0.5 * de_true * de_true / (GAMMA0 * GAMMA0)));
+                    r.sigma = 1.0f;
+                    r.prescaling_corr = 1.0f;
+                    r.partiality = 1.0f;
+                    o.reflections.push_back(r);
+                }
+    }
+    const StillsPartialityRefine refine(x);
+    double scale = 1.0;
+    for (int i = 0; i < 3; ++i)
+        scale = refine.MeasureRadialScale(outcomes, std::nullopt, scale);
+    CHECK(scale - 1.0 == Catch::Approx(EPS).epsilon(0.15));
+
+    // and a lattice that is right stays put
+    for (auto &o: outcomes)
+        for (auto &r: o.reflections)
+            r.I = static_cast<float>(*o.image_scale_g * 1000.0);
+    CHECK(std::fabs(refine.MeasureRadialScale(outcomes, std::nullopt, 1.0) - 1.0) < 5e-4);
+}
+
+
 // A frame indexed with its c axis doubled: its reciprocal lattice holds every node of the true one plus
 // the planes in between. Relabelling keeps exactly the true nodes, each at the same scattering vector,
 // and gives the frame the true cell.