accurate sigma estimation on merging with partiality-weighting

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<HKLAccum> 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<double>(r.I) / correction,
-                o.sigma / correction
+                o.sigma / correction,
+                1 / correction
             });
         }
     }