Stills partiality: the run's width table is balanced on the partials before the post-refinement

The rocking width is measured from moments of the excitation error through the prediction window, and
its tails are where that is least sure. On rotation sweeps processed as stills, reflections predicted at
under a third of their peak were 30% fainter than the model at low resolution and up to 50% brighter at
high resolution, against reflections near their peak - a bias in I/p that no variance term removes.

Before the post-refinement each shell of the width table is now rescaled by the factor at which the
reflections predicted below half their peak and those above it put the same scale on the merged
intensities (ratio of sums, searched by bisection). On synthetic data with a known Gaussian width the
factors come out 0.98-1.00. On the two rotation sweeps the low-partiality/high-partiality ratio per
shell goes from 0.68-0.93 and 0.65-1.49 to 0.98-1.07 and 0.92-1.14 (one outer shell of the first set
over-corrected to 0.75).

Quick panel against the previous commit (CCref last shell / overall): XFEL serial 18k frames
0.653/0.665 -> 0.657/0.685, 1.5k frames 0.252/0.543 -> 0.291/0.559; rotation sweeps as stills
0.689/0.464 -> 0.673/0.483, 0.644/0.274 -> 0.645/0.280, Zn site sweep peaks 20.6 -> 21.0; monochromatic
jet 0.458/0.396 -> 0.464/0.355 (its lowest shell 0.33 -> 0.30).

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

image_analysis/scale_merge/StillsPartialityRefine.cpp

@@ -8,6 +8,7 @@
 #include <atomic>
 #include <cmath>
 #include <future>
+#include <limits>
 #include <thread>
 #include <vector>
 
@@ -25,6 +26,9 @@ namespace {
     // Model error of the per-crystal scale fit, as a fraction of the predicted intensity - ScaleOnTheFly's.
     constexpr double SCALE_MODEL_ERROR = 0.3;
 
+    // The width scales BalanceRunWidth searches: from 1/BALANCE_RANGE to BALANCE_RANGE times the table.
+    constexpr double BALANCE_RANGE = 3.0;
+
     // The smallest priors on a crystal's radial scale and log width scale (RunPrior).
     constexpr double MIN_PRIOR_ETA = 1e-5;
     constexpr double MIN_PRIOR_LN_WIDTH = 0.01;
@@ -743,6 +747,85 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
     return width;
 }
 
+StillsPartialityRefine::RunWidth StillsPartialityRefine::BalanceRunWidth(
+        const std::vector<IntegrationOutcome> &outcomes, const std::map<HKLKey, double> &reference,
+        const std::optional<UnitCell> &reference_cell, const RunWidth &width) const {
+    // The width table is measured from moments of the excitation error through the prediction window, and
+    // the tails of the curve are where that is least sure: on rotation sweeps processed as stills, the
+    // reflections recorded at under a third of their peak came out 30% fainter than the model at low
+    // resolution and up to 50% brighter at high resolution, against those recorded near their peak. So
+    // each shell of the table is scaled by the factor s at which its reflections predicted below half
+    // their peak and those above it put the same scale on the merged intensities:
+    // sum I / sum p_s I_ref is the same for both, p_s the Gaussian of width s sigma.
+    if (width.table_dstar.empty())
+        return width;
+    const auto &indexing = experiment_.GetIndexingSettings();
+    struct Point { double u, x, ref; };   // eps / sigma, I / G (corrected), the reference intensity
+    std::vector<std::vector<Point>> shells(width.table_dstar.size());
+    const Width w{.gamma0 = width.gamma0, .gamma_e = width.gamma_e, .bw = width.bandwidth_sigma,
+                  .table_dstar = &width.table_dstar, .table_s = &width.table_s};
+    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_) || !std::isfinite(r.I))
+                continue;
+            const auto it = reference.find(hkl_key_generator_(r));
+            if (it == reference.end() || !(it->second > 0.0))
+                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>(width.radial_scale);
+            const double dstar = VecLen(q.x, q.y, q.z);
+            const double de = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - inv_lambda;
+            const double x = r.I * r.prescaling_corr * r.qe_corr * r.flight_corr / *o.image_scale_g;
+            if (!std::isfinite(x))
+                continue;
+            // the table shell nearest in d*
+            size_t k = std::upper_bound(width.table_dstar.begin(), width.table_dstar.end(), dstar)
+                       - width.table_dstar.begin();
+            if (k == width.table_dstar.size() || (k > 0 && dstar - width.table_dstar[k - 1] < width.table_dstar[k] - dstar))
+                k = k > 0 ? k - 1 : 0;
+            shells[k].push_back({de / std::sqrt(w.Sigma2(dstar, q.z)), x, it->second});
+        }
+    }
+    // log(sum I / sum p_s I_ref) over the reflections below half the table's peak, less the same over those
+    // above it; falls as s grows.
+    auto gap = [](const std::vector<Point> &pts, double s) {
+        const double half = std::sqrt(2.0 * std::log(2.0));   // |u| at half the peak
+        double x_lo = 0.0, m_lo = 0.0, x_hi = 0.0, m_hi = 0.0;
+        for (const auto &p: pts) {
+            const double ps = std::exp(-0.5 * p.u * p.u / (s * s));
+            (std::fabs(p.u) > half ? x_lo : x_hi) += p.x;
+            (std::fabs(p.u) > half ? m_lo : m_hi) += ps * p.ref;
+        }
+        if (!(x_lo > 0.0 && m_lo > 0.0 && x_hi > 0.0 && m_hi > 0.0))
+            return std::numeric_limits<double>::quiet_NaN();
+        return std::log(x_lo / m_lo) - std::log(x_hi / m_hi);
+    };
+    RunWidth balanced = width;
+    std::string table;
+    for (size_t k = 0; k < shells.size(); ++k) {
+        double lo = 1.0 / BALANCE_RANGE, hi = BALANCE_RANGE;
+        const double g_lo = gap(shells[k], lo), g_hi = gap(shells[k], hi);
+        if (shells[k].size() < 400 || !(g_lo > 0.0) || !(g_hi < 0.0))
+            continue;
+        for (int i = 0; i < 40; ++i) {
+            const double mid = std::sqrt(lo * hi);
+            (gap(shells[k], mid) > 0.0 ? lo : hi) = mid;
+        }
+        balanced.table_s[k] *= std::sqrt(lo * hi);
+        table += fmt::format(" {:.2f} A x{:.2f}", 1.0 / width.table_dstar[k], std::sqrt(lo * hi));
+    }
+    Logger("StillsPartialityRefine").Info("Stills rocking width balanced on the partials, against resolution:{}",
+                                          table.empty() ? " unchanged" : table);
+    return balanced;
+}
+
 StillsPartialityRefine::CrystalPrior StillsPartialityRefine::RunPrior(const std::vector<CrystalFit> &fits) const {
     // The spread of a parameter over the crystals: the variance of the estimates less the median variance
     // they carry from their own data, floored at a small value (a crystal is never held exactly).
@@ -805,8 +888,6 @@ double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, si
         run_width = MeasureRunWidth(outcomes, reference_cell, radial_scale);
     }
     const double radial_scale = run_width ? run_width->radial_scale : 1.0;
-    if (used_width)
-        *used_width = run_width;
 
     // The priors on a crystal's radial scale and width are the run's own spread of them: the first cycle
     // refines every crystal against loose priors, and the spread of those estimates less their own
@@ -821,6 +902,8 @@ double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, si
         std::map<HKLKey, double> reference;
         for (const auto &m: merged)
             reference[hkl_key_generator_(m)] = m.I;
+        if (outer == 0 && run_width)
+            run_width = BalanceRunWidth(outcomes, reference, reference_cell, *run_width);
 
         std::atomic<size_t> next{0};
         auto worker = [&]() {
@@ -848,6 +931,9 @@ double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, si
             prior = RunPrior(fits);
     }
 
+    if (used_width)
+        *used_width = run_width;
+
     double tilt_sum = 0.0;
     size_t tilt_n = 0;
     for (const auto &f: fits)

image_analysis/scale_merge/StillsPartialityRefine.h

@@ -109,6 +109,9 @@ private:
     CrystalFit RefineOne(IntegrationOutcome &outcome, const std::map<HKLKey, double> &reference,
                          const std::optional<RunWidth> &run_width, double radial_scale,
                          const CrystalPrior &prior) const;
+    // The run's width table, each shell rescaled so that its weak and strong partials agree on the scale.
+    RunWidth BalanceRunWidth(const std::vector<IntegrationOutcome> &outcomes, const std::map<HKLKey, double> &reference,
+                             const std::optional<UnitCell> &reference_cell, const RunWidth &width) const;
     // The run's spread of eta and ln width over the crystals, net of their own uncertainty.
     CrystalPrior RunPrior(const std::vector<CrystalFit> &fits) const;
 };