From a6e242fec78b64e4b7a7c10d22cea4d3406afb3e Mon Sep 17 00:00:00 2001 From: takaba_k Date: Mon, 27 Apr 2026 14:25:24 +0200 Subject: [PATCH] accurate sigma estimation on merging with partiality-weighting --- image_analysis/scale_merge/ScaleAndMerge.cpp | 22 +++++++------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 2d813178..d5baccf0 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -241,6 +241,7 @@ namespace { int hkl_slot; double I_corr; double sigma_corr; + double weight; }; void scale(const ScaleMergeOptions &opt, @@ -408,18 +409,16 @@ namespace { struct HKLAccum { double sum_wI = 0.0; double sum_w = 0.0; - double sum_wI2 = 0.0; - int count = 0; + double sum_wsigma2 = 0.0; }; std::vector accum(nhkl); for (const auto &co: corr_obs) { - const double w = 1.0 / (co.sigma_corr * co.sigma_corr); + const double w = co.weight / (co.sigma_corr * co.sigma_corr); auto &a = accum[co.hkl_slot]; a.sum_wI += w * co.I_corr; a.sum_w += w; - a.sum_wI2 += w * co.I_corr * co.I_corr; - a.count++; + a.sum_wsigma2 += w * w * co.sigma_corr * co.sigma_corr; } for (int h = 0; h < nhkl; ++h) { @@ -428,15 +427,7 @@ namespace { continue; out.merged[h].I = a.sum_wI / a.sum_w; - const double sigma_stat = std::sqrt(1.0 / a.sum_w); - - if (a.count >= 2) { - const double var_internal = a.sum_wI2 - (a.sum_wI * a.sum_wI) / a.sum_w; - const double sigma_int = std::sqrt(var_internal / (a.sum_w * (a.count - 1))); - out.merged[h].sigma = std::max(sigma_stat, sigma_int); - } else { - out.merged[h].sigma = sigma_stat; - } + out.merged[h].sigma = std::sqrt(a.sum_wsigma2) / a.sum_w; } } @@ -643,7 +634,8 @@ namespace { corr_obs.push_back({ o.hkl_slot, static_cast(r.I) / correction, - o.sigma / correction + o.sigma / correction, + 1 / correction }); } } -- 2.52.0