diff --git a/image_analysis/bragg_integration/CalcISigma.cpp b/image_analysis/bragg_integration/CalcISigma.cpp index 27bf1ff36..9ed47147b 100644 --- a/image_analysis/bragg_integration/CalcISigma.cpp +++ b/image_analysis/bragg_integration/CalcISigma.cpp @@ -9,11 +9,11 @@ #include "Regression.h" #include "../../common/ResolutionShells.h" -// Upper bound on a physically plausible macromolecular isotropic B-factor (A^2). A fit that lands -// outside (0, WILSON_B_MAX) comes from a bad frame / bad dataset (too few or mis-indexed reflections, -// a degenerate resolution range) rather than real Debye-Waller falloff, so it is rejected as -// indeterminate rather than reported. Radiation-damaged / low-resolution data rarely exceeds ~120 A^2; -// 200 leaves generous head-room while still excluding the hundreds-of-A^2 garbage. +// Upper bound on a physically plausible per-image isotropic B-factor (A^2). A single frame's fit that +// lands outside (0, WILSON_B_MAX) comes from a bad frame (too few or mis-indexed reflections, a +// degenerate resolution range) rather than real Debye-Waller falloff, so it is rejected as +// indeterminate rather than reported. The dataset-wide estimate below has no such cap: a crystal +// diffracting to 4-7 A has a true Wilson B of 300 A^2 and more. static constexpr float WILSON_B_MAX = 200.0f; // What each of the two per-image estimates makes of its shell sums, once they are gathered. Split out @@ -175,9 +175,19 @@ GlobalWilsonB CalcGlobalWilsonB(const std::vector &merged) { // Below ~4 A the Wilson plot is non-linear (bonding/solvent structure), so when the data extend to // lower resolution than that, restrict the fit to d <= 4 A - the standard Wilson-B convention. If // the whole dataset is coarser than 4 A, fall back to using all of it. + // But the window never gets narrower than the one a 3 A dataset fits over, 4-3 A in 1/d^2: with + // d_min just under 4 A the 4 A window is a sliver, and after the shells past the signal limit are + // dropped it held 0-2 shells. A narrower window is widened towards low resolution instead - into + // the curved part of the plot, which costs less than a slope over no range at all. constexpr double WILSON_LOW_RES_LIMIT_A = 4.0; - const float d_low = (d_max > WILSON_LOW_RES_LIMIT_A && d_min < WILSON_LOW_RES_LIMIT_A) - ? static_cast(WILSON_LOW_RES_LIMIT_A) : d_max; + constexpr double WILSON_MIN_WINDOW = 1.0 / (3.0 * 3.0) - 1.0 / (4.0 * 4.0); + float d_low = (d_max > WILSON_LOW_RES_LIMIT_A && d_min < WILSON_LOW_RES_LIMIT_A) + ? static_cast(WILSON_LOW_RES_LIMIT_A) : d_max; + const double widest_low_s2 = 1.0 / (static_cast(d_min) * d_min) - WILSON_MIN_WINDOW; + if (1.0 / (static_cast(d_low) * d_low) > widest_low_s2) + d_low = widest_low_s2 > 0.0 ? std::min(d_max, static_cast(1.0 / std::sqrt(widest_low_s2))) : d_max; + out.d_low_A = d_low; + out.d_high_A = d_min; const int nshells = 20; ResolutionShells shells(d_min, d_low, nshells); @@ -214,7 +224,7 @@ GlobalWilsonB CalcGlobalWilsonB(const std::vector &merged) { const double b = -2.0 * reg.slope; // ~ exp(-2 B s^2), s^2 = 1/(4 d^2), x = 1/d^2 out.n_shells = static_cast(x.size()); out.correlation = std::sqrt(std::clamp(static_cast(reg.r_square), 0.0, 1.0)); - if (std::isfinite(b) && b > 0.0 && b < WILSON_B_MAX) + if (std::isfinite(b) && b > 0.0) out.b = b; return out; } diff --git a/image_analysis/bragg_integration/CalcISigma.h b/image_analysis/bragg_integration/CalcISigma.h index 9e61767eb..6245661ad 100644 --- a/image_analysis/bragg_integration/CalcISigma.h +++ b/image_analysis/bragg_integration/CalcISigma.h @@ -26,6 +26,8 @@ struct GlobalWilsonB { double b = std::numeric_limits::quiet_NaN(); // Wilson B-factor estimate (A^2) double correlation = std::numeric_limits::quiet_NaN(); // |correlation| of the log-linear fit int n_shells = 0; // resolution shells actually used + double d_low_A = std::numeric_limits::quiet_NaN(); // the fit window, low-resolution end + double d_high_A = std::numeric_limits::quiet_NaN(); // and high-resolution end (A) }; GlobalWilsonB CalcGlobalWilsonB(const std::vector &merged); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index af996e7d6..54ec8d9eb 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -856,8 +856,9 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << in.overall; output << std::endl; if (std::isfinite(in.wilson_b) && in.wilson_b > 0.0) - output << fmt::format(" Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f})", - in.wilson_b, in.wilson_b_correlation) << std::endl; + output << fmt::format(" Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, fitted over {:.2f}-{:.2f} A)", + in.wilson_b, in.wilson_b_correlation, in.wilson_b_d_low_A, in.wilson_b_d_high_A) + << std::endl; if (std::isfinite(in.radiation_damage_delta_b)) output << fmt::format(" Radiation damage: relative B-factor change over run = {:+.2f} A^2 " "(first->last, {} batches)", diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 9ab944876..b8f2b797e 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -177,6 +177,8 @@ struct MergeStatistics { // Diagnostic only; not used in scaling. NaN when not determined. double wilson_b = NAN; double wilson_b_correlation = NAN; + double wilson_b_d_low_A = NAN; // the window the fit was made over + double wilson_b_d_high_A = NAN; // Radiation-damage monitor (rotation only): the relative Debye-Waller B change measured from the first // to the last frame of the run (A^2; positive = high-resolution intensity fades with dose = damage) and diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index db38ff87e..6c46f28eb 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1156,6 +1156,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("CC_ANOM", o.cc_anom, "{:.4f}")); if (std::isfinite(ms.wilson_b)) Add(s, KeyReal("WILSON_B", ms.wilson_b, "{:.2f}")); + if (std::isfinite(ms.wilson_b) && std::isfinite(ms.wilson_b_d_low_A)) + Add(s, KeyText("WILSON_B_RANGE_A", fmt::format("{:.2f} {:.2f}", ms.wilson_b_d_low_A, + ms.wilson_b_d_high_A), true)); // The error model in XDS's convention, so the numbers are directly comparable with a // CORRECT.LP. Kept in the default report because beamline monitoring scripts read them. Add(s, KeyReal("ERROR_MODEL_A", result.error_model_a, "{:.4f}")); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 6fd346d02..3e3a63028 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3358,9 +3358,12 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); sm.statistics.wilson_b = wilson.b; sm.statistics.wilson_b_correlation = wilson.correlation; + sm.statistics.wilson_b_d_low_A = wilson.d_low_A; + sm.statistics.wilson_b_d_high_A = wilson.d_high_A; if (std::isfinite(wilson.b) && wilson.b > 0.0) - logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", - wilson.b, wilson.correlation, wilson.n_shells); + logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells, fitted over " + "{:.2f}-{:.2f} A)", wilson.b, wilson.correlation, wilson.n_shells, wilson.d_low_A, + wilson.d_high_A); // Diffraction anisotropy. It describes the data and corrects no intensity; its tensor is // used only as the Wilson prior of the amplitudes, below. The tensor and its resolution signature come from diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9a8ea0e21..10f78955c 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2121,9 +2121,12 @@ static int RunRugnux(int argc, char **argv) { const GlobalWilsonB wilson = CalcGlobalWilsonB(merged_reflections); merged_statistics.wilson_b = wilson.b; merged_statistics.wilson_b_correlation = wilson.correlation; + merged_statistics.wilson_b_d_low_A = wilson.d_low_A; + merged_statistics.wilson_b_d_high_A = wilson.d_high_A; if (std::isfinite(wilson.b) && wilson.b > 0.0) - logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", - wilson.b, wilson.correlation, wilson.n_shells); + logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells, fitted over " + "{:.2f}-{:.2f} A)", wilson.b, wilson.correlation, wilson.n_shells, wilson.d_low_A, + wilson.d_high_A); } { diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 187c2cea7..aab3d5e9d 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -4,6 +4,7 @@ #include #include +#include "../image_analysis/bragg_integration/CalcISigma.h" #include "../image_analysis/scale_merge/ErrorModel.h" #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/scale_merge/Merge.h" @@ -555,3 +556,32 @@ TEST_CASE("ErrorModel_DeviceFitIsTheHostFit") { } } #endif + +// The dataset-wide Wilson B: a pure exp(-B / 2d^2) fall-off, read back over the fit window. A crystal +// diffracting to just under 4 A has almost nothing between its d_min and the 4 A convention, so the +// window is widened to the 1/d^2 width a 3 A dataset gets; and a high B is a measurement, not garbage. +TEST_CASE("Global Wilson B over a widened window, with no upper cap", "[wilson]") { + const auto merge = [](double b, double d_min) { + std::vector out; + for (int i = 0; i < 20000; ++i) { + // Uniform in 1/d^3 from 20 A to d_min, as reflections are. + const double s3 = (i + 0.5) / 20000.0 * (1.0 / std::pow(d_min, 3) - 1.0 / 8000.0) + 1.0 / 8000.0; + MergedReflection r; + r.h = i; r.k = 0; r.l = 0; + r.d = static_cast(1.0 / std::cbrt(s3)); + r.I = static_cast(1e6 * std::exp(-b / (2.0 * r.d * r.d))); + r.sigma = 0.01f * r.I; + out.push_back(r); + } + return out; + }; + const GlobalWilsonB high_res = CalcGlobalWilsonB(merge(40.0, 2.0)); + CHECK(high_res.b == Catch::Approx(40.0).epsilon(0.02)); + CHECK(high_res.d_low_A == Catch::Approx(4.0)); + + const GlobalWilsonB low_res = CalcGlobalWilsonB(merge(300.0, 3.9)); + CHECK(low_res.b == Catch::Approx(300.0).epsilon(0.02)); + CHECK(low_res.n_shells >= 3); + CHECK(1.0 / (3.9 * 3.9) - 1.0 / (low_res.d_low_A * low_res.d_low_A) + == Catch::Approx(1.0 / 9.0 - 1.0 / 16.0).epsilon(0.01)); +}