diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 04796300..b4637fad 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -5,6 +5,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * Bragg integration: the background ring's outer radius default changes from 10 px to **13 px**, which roughly doubles the pixels behind each background estimate; the signal disk is unchanged. +* Scaling: when too few reflections are strong enough to constrain the error model's systematic term **b**, it is now held at zero and **ISa is reported as unmeasured** rather than as the spurious value the fit would otherwise extrapolate. * Scaling: the partials of one rocking event now share a single exact-Bragg angle, instead of each solving its own from its frame's refined orientation, which had been feeding per-frame refinement noise into partiality in proportion to 1/zeta. * Bragg integration: signal pixels shared with a neighbouring reflection are now dropped from the profile fit (`--overlap off|reject|exclude`, **default `exclude`**), so a crowded reflection no longer reads its neighbour's flux as its own; `reject` instead discards a reflection whose cleanly-observed profile fraction is below `--overlap-minpk` (default 0.75). * rugnux: The background ring can be **elongated radially per reflection** for broadband data (`--integration-stencil `, default 0 = the fixed circular ring), by `k` times the beam's radial streak; the `r1` signal box stays circular. diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index e494c8af..f9e43a73 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1795,6 +1795,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool std::vector reject_median(n_groups, NAN); double error_model_a = 1.0, error_model_b = 0.0, error_model_chi2 = 0.0; bool error_model_active = false; + bool error_model_b_unmeasured = false; // b had no leverage; ISa is not reported // One leverage-corrected sample per usable full: its raw variance, its group's mean intensity, its // squared deviation from that mean - and the resolution it sits at, because the fit is re-run below // over the samples that survive the automatic resolution cutoff. @@ -1822,6 +1823,22 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool bI2.push_back(median_of(vI2)); bd2.push_back(median_of(vd2) / CHI2_1_MEDIAN); } + // `b` is identified ONLY by the spread of I^2/sigma^2 across the bins, and the bins hold equal + // COUNTS - so when fewer reflections are strong than one bin holds (1/16 of the pool), the top + // bin's median sits at an intensity where b cannot be measured, and the fit assigns it the + // bins' own noise-selection slope instead. Measured: a crystal with 2.2% of its fulls above + // I/sigma 2 reports b = 5.6, i.e. sigma -> 2*I at the strong end and ISa 0.50, while XDS gets + // 6.13 from the same images; synthetic data with b_true = 0 and 1.8% strong reproduces that to + // two digits, and recovers the truth as soon as the strong fraction passes 1/16. sigma^2 is + // then a*sigma^2 + (b*mean)^2 with a per-group constant, which caps merged |I/sigma| at + // sqrt(n)/b - reporting an impossible systematic error on data whose CC1/2 is 99%. + // + // So refuse to report what was not measured: if the strongest bin's own (I/sigma)^2 is small, + // fit `a` alone and leave b at zero. The threshold is not delicate - over the rotation battery + // the two crystals this fires on sit at 0.22 and 0.84 while the next is 31.7 and a healthy one + // is 342, so anything from 4 to 25 selects the same two and nothing else changes. + constexpr double B_LEVER_MIN = 4.0; + const bool b_measurable = bI2.back() > B_LEVER_MIN * bs2.back(); std::vector bd2_sorted = bd2; const double dev2_floor = std::max(1e-30, 1e-3 * median_of(bd2_sorted)); double Ass = 0, AsI = 0, AII = 0, Bs = 0, BI = 0; @@ -1833,10 +1850,26 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool Bs += wgt * s2 * d2; BI += wgt * I2 * d2; } const double det = Ass * AII - AsI * AsI; - if (det > 1e-10 * Ass * AII) { + if (!b_measurable) { + // One parameter, same weights: dev2 = a*sigma^2 alone. + if (Ass > 0.0) { + error_model_a = std::clamp(Bs / Ass, 0.25, 100.0); + error_model_b = 0.0; + error_model_b_unmeasured = true; + error_model_active = true; + std::vector chi2; + chi2.reserve(smp.size()); + for (const auto &s : smp) { + const double v = error_model_a * s.s2; + if (v > 0.0) chi2.push_back(s.dev2 / v); + } + error_model_chi2 = chi2.empty() ? 0.0 : median_of(chi2) / CHI2_1_MEDIAN; + } + } else if (det > 1e-10 * Ass * AII) { error_model_a = std::clamp((Bs * AII - BI * AsI) / det, 0.25, 100.0); const double b2 = std::max((Ass * BI - AsI * Bs) / det, 0.0); error_model_b = std::sqrt(b2); + error_model_b_unmeasured = false; error_model_active = true; std::vector chi2; chi2.reserve(smp.size()); @@ -2136,9 +2169,17 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // is what XDS's ISa means, and the strong-reflection asymptote above, which XDS has no // equivalent of and which can only ever be the more optimistic of the two. const auto em = ToXdsErrorModel(error_model_a, error_model_b); - logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} " - "(strong-reflection asymptote {:.1f}) chi2={:.2f}", - em.a, em.b, em.isa, result.isa_asymptotic, error_model_chi2); + if (error_model_b_unmeasured) + logger.Warning("Error model (XDS convention): a={:.3f}, b NOT MEASURABLE - fewer than one " + "intensity bin's worth of reflections are strong enough to constrain it, so " + "b is held at 0 and ISa is not reported. chi2={:.2f}. This says the data do " + "not reach far enough for a systematic error to be seen, not that there is " + "none; a resolution range matched to the signal would measure it", + em.a, error_model_chi2); + else + logger.Info("Error model (XDS convention): a={:.3f} b={:.3e} ISa={:.1f} " + "(strong-reflection asymptote {:.1f}) chi2={:.2f}", + em.a, em.b, em.isa, result.isa_asymptotic, error_model_chi2); } AssignRfreeFlags(result.merged, x.GetSpaceGroupNumber().value_or(1), rfree_fraction);