mx/Jungfraujoch
1
1
Fork 0
Code Issues Pull Requests 1 Actions Packages Projects Releases 52 Wiki Activity

accurate sigma estimation on merging with partiality-weighting #47

Open
takaba_k wants to merge 1 commits from 2605-sigma-correction into main
pull from: 2605-sigma-correction
merge into: mx:main
Conversation 0 Commits 1 Files Changed 1
+7 -15
1 changed files with 7 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/ScaleAndMerge.cpp
+7 -15
View File
@@ -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
});
}
}