diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp
index 98664cef7..401c8882d 100644
--- a/image_analysis/scale_merge/StillsPartialityRefine.cpp
+++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp
@@ -434,7 +434,7 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
     // intensities of thousands of crystals average out, and fitted as s^2 = gamma0^2 + gamma_e^2 d*^2 +
     // (bandwidth_sigma lambda d*^2 / 2)^2.
     const auto &indexing = experiment_.GetIndexingSettings();
-    struct Point { double x, de2, w; };
+    struct Point { double x, de2, w, var_w; };
     std::vector<Point> pts;
     for (const auto &o: outcomes) {
         if (!o.image_scale_g || !(*o.image_scale_g > 0.0))
@@ -445,11 +445,12 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
         for (const auto &r: o.reflections) {
             if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_))
                 continue;
-            const double w = r.I * r.prescaling_corr * r.qe_corr * r.flight_corr / *o.image_scale_g;
-            if (!std::isfinite(w) || !std::isfinite(r.dist_ewald) || !(r.d > 0.0f))
+            const double corr = r.prescaling_corr * r.qe_corr * r.flight_corr / *o.image_scale_g;
+            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});
+                           static_cast<double>(r.dist_ewald) * r.dist_ewald, w, sigma_w * sigma_w});
         }
     }
     constexpr size_t N_SHELLS = 8;
@@ -458,24 +459,32 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
     std::sort(pts.begin(), pts.end(), [](const Point &a, const Point &b) { return a.x < b.x; });
 
     Logger logger("StillsPartialityRefine");
-    std::array<double, N_SHELLS> sx{}, sy{};
+    // A shell's moment is no better than the summed intensity it is weighted by, so a shell is used only
+    // where that sum is measured to 10 % (sum I / sigma(sum I) >= 10): the shells beyond the diffraction
+    // limit, where the sum is noise and can even be negative, say nothing about the width.
+    constexpr double MIN_SHELL_SIGNAL = 10.0;
+    std::vector<double> sx, sy;
     for (size_t k = 0; k < N_SHELLS; ++k) {
-        double sw = 0.0, swde2 = 0.0, sde2 = 0.0, sumx = 0.0;
+        double sw = 0.0, svar = 0.0, swde2 = 0.0, sde2 = 0.0, sumx = 0.0;
         const size_t lo = k * pts.size() / N_SHELLS, hi = (k + 1) * pts.size() / N_SHELLS;
         for (size_t i = lo; i < hi; ++i) {
             sw += pts[i].w;
+            svar += pts[i].var_w;
             swde2 += pts[i].w * pts[i].de2;
             sde2 += pts[i].de2;
             sumx += pts[i].x;
         }
-        if (!(sw > 0.0))
-            return std::nullopt;
         const double n = static_cast<double>(hi - lo);
-        sx[k] = sumx / n;
-        sy[k] = UntruncatedVariance(swde2 / sw, sde2 / n);
-        logger.Debug("Stills rocking width at d = {:.2f} A: {:.2e} 1/A ({:.3f} deg)", 1.0 / std::sqrt(sx[k]),
-                     std::sqrt(sy[k]), std::sqrt(sy[k] / sx[k]) * kRadToDeg);
+        if (!(sw >= MIN_SHELL_SIGNAL * std::sqrt(svar)) || !(swde2 > 0.0))
+            continue;
+        const double x = sumx / n, y = UntruncatedVariance(swde2 / sw, sde2 / n);
+        sx.push_back(x);
+        sy.push_back(y);
+        logger.Debug("Stills rocking width at d = {:.2f} A: {:.2e} 1/A ({:.3f} deg), I/sigma {:.0f}", 1.0 / std::sqrt(x),
+                     std::sqrt(y), std::sqrt(y / x) * kRadToDeg, sw / std::sqrt(svar));
     }
+    if (sx.size() < 3)
+        return std::nullopt;
 
     // Least squares for s^2 = A + B d*^2 + C d*^4 over the shells. Every term is a variance, so a negative
     // one is unphysical: of the fits on each subset of the three terms, the best with none negative is kept.
@@ -489,7 +498,7 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
         const int m = static_cast<int>(terms.size());
         Eigen::MatrixXd M = Eigen::MatrixXd::Zero(m, m);
         Eigen::VectorXd v = Eigen::VectorXd::Zero(m);
-        for (size_t k = 0; k < N_SHELLS; ++k)
+        for (size_t k = 0; k < sx.size(); ++k)
             for (int i = 0; i < m; ++i) {
                 const double fi = std::pow(sx[k], terms[i]);
                 v(i) += fi * sy[k];
@@ -503,7 +512,7 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
         for (int i = 0; i < m; ++i)
             coef[terms[i]] = c(i);
         double residual = 0.0;
-        for (size_t k = 0; k < N_SHELLS; ++k) {
+        for (size_t k = 0; k < sx.size(); ++k) {
             const double r = sy[k] - coef[0] - coef[1] * sx[k] - coef[2] * sx[k] * sx[k];
             residual += r * r;
         }