From b440e3792bbfa2e0cd5cfd9b41e12de969311737 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Oct 2026 21:23:28 +0200 Subject: [PATCH] StillsPartialityRefine: measure the run width only on shells with signal Without a resolution limit the outer d*^2 shells of the width measurement lie beyond the diffraction limit, where the summed intensity is noise or even negative; one such shell made the whole measurement fail, and 5 of 11 still-battery sets (sacla237, pal86, pal90, lcls242_s, lcls242_t) fell back to the per-crystal window fit. A shell now enters the fit only where its summed intensity is measured to 10 % (sum I / sigma >= 10), and at least three shells are needed. Weighting the fit by the shells' signal instead was tried and rejected: it lets the strong low-resolution shells set the width alone, and on the MicroMAX pink set the dimple R-free went 0.2442 -> 0.2748. The matrix subsets, which are cut at their resolution, are unchanged within 0.0010 dimple R-free (SwissFEL 0.2453, -C/-S 0.2513, mono 0.3188, pink 0.2442, SACLA 243 0.2455, SACLA 237 0.2327, LCLS 222 0.2335, PAL 86 0.3628). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- .../scale_merge/StillsPartialityRefine.cpp | 37 ++++++++++++------- 1 file changed, 23 insertions(+), 14 deletions(-) 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::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 pts; for (const auto &o: outcomes) { if (!o.image_scale_g || !(*o.image_scale_g > 0.0)) @@ -445,11 +445,12 @@ std::optional 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(r.d) * r.d), - static_cast(r.dist_ewald) * r.dist_ewald, w}); + static_cast(r.dist_ewald) * r.dist_ewald, w, sigma_w * sigma_w}); } } constexpr size_t N_SHELLS = 8; @@ -458,24 +459,32 @@ std::optional StillsPartialityRefine::MeasureR std::sort(pts.begin(), pts.end(), [](const Point &a, const Point &b) { return a.x < b.x; }); Logger logger("StillsPartialityRefine"); - std::array 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 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(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::MeasureR const int m = static_cast(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::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; }