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

Rotation merging: keep the systematic sigma floor when ISa is unmeasurable

Browse Source 
The cap that refuses to report an impossible ISa was zeroing the asymptotic b
itself, and that same value is the floor passed to SigmaWithSystematicFloor -
where zero means "no floor". So on the degenerate low-multiplicity fit the guard
is written for, instead of capping merged I/sigma at 100 it removed the cap
entirely. Report the asymptote as unmeasured, keep the fitted value for the
floor.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-30 10:50:30 +02:00
co-authored by Claude Opus 5
parent 386f10ad07
commit 609650b061
1 changed files with 9 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/RotationScaleMerge.cpp
+9 -5
View File
@@ -1599,13 +1599,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
}
// Guard a degenerate low-multiplicity fit: with too few symmetry equivalents both the (a, b) fit and
// the per-group scatter collapse toward zero, and 1/b then reports an impossibly high asymptotic
// I/sigma. Real macromolecular data does not exceed ISa ~50; past a generous cap treat the asymptote
// as unmeasured (result.isa left undetermined) rather than emit a spurious extreme.
// I/sigma. Real macromolecular data does not exceed ISa ~50; past a generous cap report the asymptote
// as unmeasured rather than emit a spurious extreme. Only the REPORT is dropped: the sigma floor
// below keeps using the fitted value, because a degenerate fit is precisely where the raw
// 1/sqrt(sum_w) is least trustworthy and leaving the merged sigma uncapped would be the opposite of
// what this guard is for.
constexpr double MIN_ASYMPTOTIC_B = 0.01; // ISa cap 100
if (error_model_b_asymptotic < MIN_ASYMPTOTIC_B) error_model_b_asymptotic = 0.0;
const double isa_reported = error_model_b_asymptotic >= MIN_ASYMPTOTIC_B
? 1.0 / error_model_b_asymptotic : 0.0;
if (error_model_active)
logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b,
error_model_b_asymptotic > 0 ? 1.0 / error_model_b_asymptotic : 0.0, error_model_chi2);
isa_reported, error_model_chi2);
auto corrected_sigma = [&](float I_corr, float sigma_corr, int g) -> float {
if (!error_model_active) return sigma_corr;
@@ -1662,7 +1666,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
}
// ---- Export merged reflections (+ resolution-shell R-free flags). ----
Result result;
result.isa = error_model_b_asymptotic > 0 ? 1.0 / error_model_b_asymptotic : 0.0;
result.isa = isa_reported;
std::vector<double> merged_I(n_groups, NAN);
float d_min = std::numeric_limits<float>::max(), d_max = 0.0f;
for (int g = 0; g < n_groups; ++g) {