diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index bf2805f8..14ecdecd 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 * Scaling: rotation merging now also fits a correction surface over **rotation angle × detector position**, which reaches the absorption left behind when the diffracting volume drifts through the beam; like the other surfaces it is cross-validated and applied only where it improves reflections it was not fitted on. * Scaling: the rotation correction surfaces (decay, absorption, detector-plane modulation) no longer pull each cell's factor towards zero in proportion to that cell's own noise, which had put a spurious centre-to-edge ramp on runs whose outer resolution shells carry no signal. +* rugnux: The radiation-damage relative *B*-factor is now measured on resolution-shell means instead of single observations, a batch whose data cannot support a measurement prints `-` rather than a value, and no first→last number is reported for a curve that no straight line describes — a disturbance confined to part of a sweep is no longer reported as dose. * rugnux: A goniometer that turned further than the angles stored in the file — which are the commanded ones — is now measured and **corrected in the second rotation pass**, or set by hand with `--rotation-scale `; the fitted correction is applied only when it exceeds 0.5 %, moves each end of the sweep by at least 0.5°, and comes out the same with any fifth of the sweep left out. * rugnux: Every `mx` and `scale` run now writes `_report.txt`, a results report modelled on XDS's `CORRECT.LP` — `KEY= value` lines, fixed-width tables and `WARNING:` sentences covering indexing, geometry post-refinement, the space-group decision, merging, twinning, radiation damage and the stretches of the sweep over which the crystal delivered much less than the rest of the run. * Space-group search: a screw axis is now claimed from how decisively its predicted-absent reflections are weaker than the rest of their own axial row, instead of from a minimum count of them, so a screw survives a sweep that recorded few axial reflections and is refused on a row too weak to decide either way. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 2a12663e..ee40e70b 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -758,6 +758,8 @@ Each surface is **cross-validated**: fitted on even-numbered frames and kept onl **Radiation-damage report (rotation, report-only).** Independently of whether any decay correction is applied, rugnux measures and reports the relative Debye–Waller $B$ across the sweep: the per-image scale's correlation to the merge and the per-image mosaicity versus frame (dose), together with a per-batch relative-$B$ curve whose first→last change is a single headline number (measured before any decay correction, against the least-damaged early wedge). It is written to the log and to the merged mmCIF as a data-quality-vs-dose diagnostic and **never** alters the merged intensities — distinct from the decay correction above, which does fold into the scale. +Each batch's $B$ is fitted on **resolution-shell means**, not on single observations: $\ln(I_\mathrm{ref}/I_\mathrm{obs})$ of one weak observation is unbounded and biased — the observation appears in the response and in its own weight, and the logarithm needs $I_\mathrm{obs} > 0$, which keeps only the upward half of the noise — so where a stretch of the sweep holds little signal that bias alone reaches tens of Ų, and on genuinely decayed data it grows with dose until it reverses the sign of the answer. The shell fit carries an intercept as well as a slope, so a batch that is merely *dimmer* than the run (an attenuated beam, a mis-fitted frame scale) is not reported as damage. A batch whose shells are too weak to fit, or whose solved value reaches the numerical bound the smoothing solve clamps to, is reported as **absent** rather than as a number, and the first→last headline is reported only where a straight line describes the curve: radiation damage is progressive, so a curve that dips and recovers is a disturbance rather than dose, and is left to the sweep-quality report to name. + ### 10.7 R-free test-set flags A fraction of the unique reflections (`rfree_fraction`, default 0.05) is flagged as a **free (test) set**, written to the output (MTZ `FreeR_flag`, mmCIF `_refln.status_free`, a text-HKL column) for model validation (§14) and for downstream refinement. The flag is a pure function of the reflection's **Friedel-merged (Laue) ASU index**, which gives three properties: diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 71e661f1..5e9fc34b 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -330,7 +330,11 @@ rather than a source of added noise, which is why they are safe to leave on. Independently of any correction, a rotation run prints a **radiation-damage report** — the per-image scale correlation-to-merge and mosaicity versus dose, and the relative *B*-factor change over the run (first→last) together with a per-batch relative-*B* curve, also written to the merged mmCIF. It is a -data-quality-vs-dose diagnostic and never alters the merged intensities. +data-quality-vs-dose diagnostic and never alters the merged intensities. A batch whose data cannot +support a measurement prints `-` instead of a value, and the first→last number is printed only where a +straight line describes the curve — damage is progressive, so a curve that dips and recovers is a +disturbance of the sweep, not dose, and the report says so and points at the sweep-quality section +(`RADIATION_DAMAGE_RELATIVE_B= NOT_A_TREND`). ### Still / serial data diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 0074937d..4635a92a 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -94,7 +94,9 @@ struct MergeStatistics { // to the last frame of the run (A^2; positive = high-resolution intensity fades with dose = damage) and // the per-batch relative-B curve it was derived from. Measured before any decay/relative-B correction is // applied, so it reports how much radiation damage was present. Diagnostic; NaN / empty for stills or - // when not determined. batch_deg is the rotation width per batch of the curve. + // when not determined. batch_deg is the rotation width per batch of the curve. A batch the data cannot + // measure is NaN in the curve, and delta_b is NaN when the curve is not a trend a single number + // summarises - damage is progressive, so a curve that is not is telling the user about something else. double radiation_damage_delta_b = NAN; std::vector radiation_damage_b_batch; double radiation_damage_batch_deg = 0.0; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index cd303138..e7d7c2ff 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -1011,17 +1011,25 @@ void RotationScaleMerge::RefineDecay(int n_groups) { 0.5 * slope * n_frames, 0.5 * slope); } +// A^2, the bound the Gauss-Seidel solve clamps each batch's relative-B to. It is a numerical guard against +// divergence, never an answer: a run of identical values in the solved curve means the solve reached the +// bound, not that the crystal did anything. The CORRECTION keeps the tight bound it was tuned with; the +// MONITOR is given room a real relative-B cannot reach, and drops any batch that lands on it anyway - a +// number the user reads as radiation damage has to be a measurement. +constexpr double RELATIVE_B_MAX = 50.0; +constexpr double RELATIVE_B_MAX_MONITOR = 200.0; + std::vector RotationScaleMerge::SolveCurvatureSmoothedB(const std::vector &num, - const std::vector &den) const { - // Solve for a SMOOTH per-batch relative-B from the per-batch normal-equation data (num_c = sum w s^2 y, - // den_c = sum w s^4). A per-batch regression alone (b_c = num_c/den_c) over-fits: a genuine relative-B is + const std::vector &den, + double b_max) const { + // Solve for a SMOOTH per-batch relative-B from the per-batch normal equations for b (num_c, den_c, such + // that the batch on its own reads b_c = num_c/den_c). A per-batch regression alone over-fits: a relative-B is // a low-curvature function of dose / orientation, but an unconstrained fit chases resolution-correlated // noise into a rough, oscillating curve. So minimise data-fidelity PLUS a second-difference (curvature) // penalty sum_c den_c (b_c - num_c/den_c)^2 + mu * sum_c (b_{c-1} - 2 b_c + b_{c+1})^2, which leaves a // straight RAMP untouched but suppresses oscillation, solving the small SPD normal system A b = num, // A = diag(den + lambda) + mu * D^T D, by Gauss-Seidel (n_batch <= a few tens). Not anchored - callers set // the gauge (zero-mean for the correction, low-dose reference for the monitor). - constexpr double RELATIVE_B_MAX = 50.0; // A^2, clamp a single batch's relative-B (generous guard) constexpr double RELATIVE_B_CURVATURE_K = 8.0; // curvature penalty weight (x the median batch leverage) const int n_batch = static_cast(num.size()); std::vector dsorted(den); @@ -1045,7 +1053,7 @@ std::vector RotationScaleMerge::SolveCurvatureSmoothedB(const std::vecto double s = num[c]; for (int j = 0; j < n_batch; ++j) if (j != c) s -= Aat(c, j) * b[j]; - b[c] = std::clamp(s / Aat(c, c), -RELATIVE_B_MAX, RELATIVE_B_MAX); + b[c] = std::clamp(s / Aat(c, c), -b_max, b_max); } return b; } @@ -1082,7 +1090,7 @@ std::vector RotationScaleMerge::FitRelativeBCurve(int n_groups, int n_ba const double s2 = s2_of(o.d), y = std::log(Iref / Is); num[batch_of(o)] += w * s2 * y; den[batch_of(o)] += w * s2 * s2; } - std::vector b = SolveCurvatureSmoothedB(num, den); + std::vector b = SolveCurvatureSmoothedB(num, den, RELATIVE_B_MAX); double dw = 0.0, dbw = 0.0; for (int c = 0; c < n_batch; ++c) { dw += den[c]; dbw += den[c] * b[c]; } const double bmean = dw > 0.0 ? dbw / dw : 0.0; @@ -1099,6 +1107,8 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { // A whole-run-mean reference inverts on heavily damaged data (the surviving late high-res reflections // look strong against a mean dominated by the frames where they had already faded). Stores the first->last // relative-B change (the headline "B-factor difference") and the per-batch curve; never touches corr. + // A batch the data cannot measure reads NaN, and the headline is left NaN unless a straight line + // actually describes the curve - what the user reads has to be a measurement or nothing. rad_damage_delta_b = std::numeric_limits::quiet_NaN(); rad_damage_b_batch.clear(); rad_damage_batch_deg = 0.0; @@ -1128,34 +1138,143 @@ void RotationScaleMerge::MeasureRadiationDamageB(int n_groups) { const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); sw0[o.group] += w; swI0[o.group] += w * static_cast(o.I) * o.corr; } - // Per batch: slope of ln(I_ref0 / I_obs) vs s^2, over groups also seen in the low-dose reference. - std::vector num(n_batch, 0.0), den(n_batch, 0.0); + // Resolution shells of equal occupancy in s^2, from a strided sample of the fulls - the per-batch fit + // regresses on the shell means, and an equal-WIDTH grid leaves the low-resolution shells (the ones that + // still carry signal in a weak batch) with almost no reflections. + constexpr int MONITOR_SHELLS = 10; + constexpr double MONITOR_MIN_SHELL_ISIGMA = 3.0; // pooled I/sigma a shell needs to be worth fitting + std::vector sample; + for (size_t i = 0; i < fulls.size(); i += 1 + fulls.size() / 200000) + if (usable(fulls[i])) sample.push_back(s2_of(fulls[i].d)); + if (static_cast(sample.size()) < MONITOR_SHELLS) + return; + std::sort(sample.begin(), sample.end()); + auto edges_over = [&](size_t n) { + std::vector e(MONITOR_SHELLS - 1); + for (int k = 1; k < MONITOR_SHELLS; ++k) e[k - 1] = sample[n * k / MONITOR_SHELLS]; + return e; + }; + std::vector edge = edges_over(sample.size()); + auto shell_of = [&](double s2) { + return static_cast(std::upper_bound(edge.begin(), edge.end(), s2) - edge.begin()); + }; + // The grid has to sit inside the range the crystal actually diffracted to. A merged resolution range + // can reach well past the signal - a limit taken from another program, a generous default - and shells + // spread over all of it would spend most of the grid on noise, leaving a batch with too few to fit at + // all. Pooled over the whole run a shell either has signal or it has none, so walk out to the first + // that has none and lay the shells inside that. + std::vector rw(MONITOR_SHELLS, 0.0), rwI(MONITOR_SHELLS, 0.0); + for (const auto &o : fulls) { + if (!usable(o)) continue; + const double sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); + const int k = shell_of(s2_of(o.d)); + rw[k] += w; rwI[k] += w * static_cast(o.I) * o.corr; + } + int n_live = 0; + while (n_live < MONITOR_SHELLS && rw[n_live] > 0.0 + && rwI[n_live] / std::sqrt(rw[n_live]) >= MONITOR_MIN_SHELL_ISIGMA) + ++n_live; + if (n_live == 0) // nothing in the run is strong enough to measure a resolution-dependent fade on + return; + if (n_live < MONITOR_SHELLS) + edge = edges_over(std::upper_bound(sample.begin(), sample.end(), edge[n_live - 1]) - sample.begin()); + + // Per (batch, shell): the inverse-variance pooled observation and its low-dose reference. The pooling + // has to happen BEFORE the logarithm. ln(I_ref / I_obs) of a SINGLE observation is unbounded and biased + // downwards as soon as I/sigma approaches 1 - the observation sits in the response and in its own weight, + // and the I_obs > 0 the logarithm needs keeps only the upward half of the noise. Simulated on a batch + // with no relative-B at all, that bias alone reads -28 A^2 at = 0.3 and pins the fit to the + // clamp; on a genuinely decayed batch it grows with dose until it reverses the sign of the answer + // (-29 A^2 where the truth is +30). Pooling first holds the same simulations to under 1 A^2. A shell + // mean is well determined where a single observation is not, admits negative intensities, and carries + // its own I/sigma - which is what says whether the batch can be measured at all. + const int n_cell = n_batch * MONITOR_SHELLS; + std::vector sw(n_cell, 0.0), swI(n_cell, 0.0), swR(n_cell, 0.0), sws2(n_cell, 0.0); for (const auto &o : fulls) { if (!usable(o) || sw0[o.group] <= 0.0) continue; const double Iref = swI0[o.group] / sw0[o.group]; - const double Is = static_cast(o.I) * o.corr, sc = static_cast(o.sigma) * o.corr; - if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue; - const double w = (Is / sc) * (Is / sc); - const double s2 = s2_of(o.d), y = std::log(Iref / Is); - num[batch_of(o)] += w * s2 * y; den[batch_of(o)] += w * s2 * s2; + if (!std::isfinite(Iref)) continue; + const double s2 = s2_of(o.d), sc = static_cast(o.sigma) * o.corr, w = 1.0 / (sc * sc); + const int i = batch_of(o) * MONITOR_SHELLS + shell_of(s2); + sw[i] += w; swI[i] += w * static_cast(o.I) * o.corr; swR[i] += w * Iref; sws2[i] += w * s2; } - std::vector b = SolveCurvatureSmoothedB(num, den); - // Anchor so the reference batches read ~0 (b = added B relative to the low-dose start). + // Per batch: slope of ln( / ) against s^2 over the shells that hold signal, weighted by + // the shell's own (I/sigma)^2 since that is the precision of its log. Fitted WITH an intercept and the + // normal equations centred, so that only the resolution-DEPENDENT part reaches the curve: a batch that + // is simply dimmer than the run - an attenuated beam, a mis-fitted frame scale - moves the intercept, + // and calling that radiation damage is the error this monitor exists to avoid. A shell too weak to + // measure, and a batch with too few shells left to carry a slope, are absent rather than estimated. + constexpr int MONITOR_MIN_SHELLS = 4; // shells a batch needs for a slope and an intercept + std::vector num(n_batch, 0.0), den(n_batch, 0.0); + std::vector measured(n_batch, 0); + for (int c = 0; c < n_batch; ++c) { + double sw2 = 0, sx = 0, sy = 0, sxx = 0, sxy = 0; + int n_shell = 0; + for (int k = 0; k < MONITOR_SHELLS; ++k) { + const int i = c * MONITOR_SHELLS + k; + if (!(sw[i] > 0.0)) continue; + const double pI = swI[i] / sw[i], pR = swR[i] / sw[i], isig = pI * std::sqrt(sw[i]); + if (!(pI > 0.0) || !(pR > 0.0) || isig < MONITOR_MIN_SHELL_ISIGMA) continue; + const double s2 = sws2[i] / sw[i], y = std::log(pR / pI), w = isig * isig; + sw2 += w; sx += w * s2; sy += w * y; sxx += w * s2 * s2; sxy += w * s2 * y; + ++n_shell; + } + if (n_shell < MONITOR_MIN_SHELLS) + continue; + const double lever = sxx - sx * sx / sw2; + if (!(lever > 0.0)) + continue; + num[c] = sxy - sx * sy / sw2; den[c] = lever; measured[c] = 1; + } + std::vector b = SolveCurvatureSmoothedB(num, den, RELATIVE_B_MAX_MONITOR); + for (int c = 0; c < n_batch; ++c) + if (std::fabs(b[c]) >= RELATIVE_B_MAX_MONITOR) measured[c] = 0; + // Anchor so the reference batches read ~0 (b = added B relative to the low-dose start). With no + // measured batch at the start of the run there is nothing to measure the rest of it against. const int n_ref_batch = std::max(1, ref_frames / frames_per_batch); double ref_mean = 0.0; - for (int c = 0; c < n_ref_batch && c < n_batch; ++c) ref_mean += b[c]; - ref_mean /= std::min(n_ref_batch, n_batch); - for (auto &v : b) v -= ref_mean; - // Headline: the first->last relative-B change from a least-squares linear trend of the (smoothed) curve - - // robust to a single noisy end batch, unlike raw endpoints. Positive = high-res fades with dose = damage. - double sx = 0, sy = 0, sxx = 0, sxy = 0; - for (int c = 0; c < n_batch; ++c) { sx += c; sy += b[c]; sxx += double(c) * c; sxy += double(c) * b[c]; } - const double det = n_batch * sxx - sx * sx; - const double slope = det > 0.0 ? (n_batch * sxy - sx * sy) / det : 0.0; - rad_damage_delta_b = slope * (n_batch - 1); + int n_ref = 0; + for (int c = 0; c < n_ref_batch && c < n_batch; ++c) + if (measured[c]) { ref_mean += b[c]; ++n_ref; } + if (n_ref == 0) + return; + ref_mean /= n_ref; rad_damage_batch_deg = MONITOR_BATCH_DEG; rad_damage_b_batch.resize(n_batch); - for (int c = 0; c < n_batch; ++c) rad_damage_b_batch[c] = static_cast(b[c]); + for (int c = 0; c < n_batch; ++c) + rad_damage_b_batch[c] = measured[c] ? static_cast(b[c] - ref_mean) + : std::numeric_limits::quiet_NaN(); + + // Headline: the relative-B change across the measured part of the run, from a least-squares linear trend + // of the curve - but reported only where that line DESCRIBES the curve, explaining at least half of its + // variance. Radiation damage is progressive, so a real one is a trend; a disturbance confined to part of + // the sweep is not, and summarising it first->last invents a decay that never happened. A curve flat to + // within a couple of A^2 is reported whatever its R^2 - there the answer is "no damage" and the shape of + // the remaining wiggle is not a question. Positive = high-res fades with dose = damage. + constexpr double MONITOR_TREND_FLAT_B = 2.0; // A^2, scatter at which the curve is simply flat + constexpr double MONITOR_TREND_UNEXPLAINED = 0.5; // ... else the line must explain half the variance + double sx = 0, sy = 0, sxx = 0, sxy = 0, n = 0; + int first = -1, last = -1; + for (int c = 0; c < n_batch; ++c) { + if (!measured[c]) continue; + const double y = b[c] - ref_mean; + sx += c; sy += y; sxx += double(c) * c; sxy += double(c) * y; n += 1; + if (first < 0) first = c; + last = c; + } + const double det = n * sxx - sx * sx; + if (!(det > 0.0)) + return; + const double slope = (n * sxy - sx * sy) / det, icept = (sy - slope * sx) / n, ymean = sy / n; + double ssr = 0.0, sst = 0.0; + for (int c = 0; c < n_batch; ++c) + if (measured[c]) { + const double y = b[c] - ref_mean; + ssr += (y - icept - slope * c) * (y - icept - slope * c); + sst += (y - ymean) * (y - ymean); + } + if (ssr <= MONITOR_TREND_UNEXPLAINED * sst || ssr <= n * MONITOR_TREND_FLAT_B * MONITOR_TREND_FLAT_B) + rad_damage_delta_b = slope * (last - first); } namespace { @@ -1363,7 +1482,7 @@ void RotationScaleMerge::MeasureSweepQuality(const std::vector &partial if (!rad_damage_b_batch.empty() && rad_damage_batch_deg > 0.0) { const int c = std::min(rad_damage_b_batch.size() - 1, static_cast(f * osc_deg / rad_damage_batch_deg)); - sum_b += rad_damage_b_batch[c]; n_b += 1; + if (std::isfinite(rad_damage_b_batch[c])) { sum_b += rad_damage_b_batch[c]; n_b += 1; } } } r.rotation_deg = static_cast(n_in * osc_deg); diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index e0b1edcb..5801f254 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -260,11 +260,12 @@ private: // physical total-dB floor). RefineAbsorption fits a smooth factor over the diffracted-beam direction in // the goniometer frame (path-length / absorption; negligible at hard X-rays, matters at low energy). void RefineDecay(int n_groups); - // Solve a smooth per-batch relative-B from the per-batch normal equations (num_c=sum w s^2 y, - // den_c=sum w s^4): data-fidelity + a second-difference (curvature) penalty, by Gauss-Seidel. Returns the - // un-anchored curve; the caller sets the gauge. Shared by the correction and the radiation-damage monitor. + // Solve a smooth per-batch relative-B from the per-batch normal equations for b (num_c, den_c): + // data-fidelity + a second-difference (curvature) penalty, by Gauss-Seidel, each batch clamped to + // +-b_max. Returns the un-anchored curve; the caller sets the gauge and the clamp it can live with. + // Shared by the correction and the radiation-damage monitor. std::vector SolveCurvatureSmoothedB(const std::vector &num, - const std::vector &den) const; + const std::vector &den, double b_max) const; // Fit a smoothed per-batch relative-B curve (A^2 per batch) on the fulls over the ASU-group subset // {group&1==gparity} (gparity<0 = all): the weighted s^2 slope of ln(Iref/Iobs) per batch against a // subset-global reference, smoothed and zero-mean-anchored. Drives the per-batch correction. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 99bec2b6..a33d1736 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -205,8 +205,14 @@ std::string RenderResultReport(const std::string &output_prefix, // ------------------------------------------------------- 7. RADIATION DAMAGE if (!result.radiation_damage_text.empty()) { Section(os, "7. RADIATION DAMAGE"); - Key(os, "RADIATION_DAMAGE_RELATIVE_B", fmt::format("{:.2f}", - result.merge_statistics.radiation_damage_delta_b)); + // A number, or a word saying why there is none: NOT_A_TREND where the per-batch curve was measured + // but no straight line describes it (damage is progressive, so that curve is not dose), NOT_MEASURED + // where the monitor could not run at all. + const double db = result.merge_statistics.radiation_damage_delta_b; + Key(os, "RADIATION_DAMAGE_RELATIVE_B", + std::isfinite(db) ? fmt::format("{:.2f}", db) + : result.merge_statistics.radiation_damage_b_batch.empty() ? std::string("NOT_MEASURED") + : std::string("NOT_A_TREND")); os << "\n" << result.radiation_damage_text << "\n"; } diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index efae3708..e6d6c9ae 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2114,17 +2114,24 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b cc9 < cc0 - 0.05 ? " (falling: radiation damage / crystal decay)" : ""); } // Relative B-factor across the sweep (measured before any decay correction, against the - // low-dose start): the resolution-dependent complement to the CC/mosaicity read above. A large - // magnitude flags radiation damage; positive is the typical direction (high-res fades with dose). - if (std::isfinite(sm.statistics.radiation_damage_delta_b)) { - os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}", - sm.statistics.radiation_damage_delta_b, - std::fabs(sm.statistics.radiation_damage_delta_b) > 5.0 - ? " (significant dose-dependent scaling: radiation damage)" : ""); + // low-dose start): the resolution-dependent complement to the CC/mosaicity read above. + // Radiation damage fades the high-resolution intensity, so only a POSITIVE change is dose. + // A curve no straight line describes gets no first->last number at all, and a batch whose + // data could not measure one prints "-" rather than a value. + if (!sm.statistics.radiation_damage_b_batch.empty()) { + if (std::isfinite(sm.statistics.radiation_damage_delta_b)) + os << fmt::format("\n => relative B-factor change over run = {:+.2f} A^2 (first->last){}", + sm.statistics.radiation_damage_delta_b, + sm.statistics.radiation_damage_delta_b > 5.0 + ? " (significant dose-dependent scaling: radiation damage)" : ""); + else + os << "\n => relative B-factor: the per-batch curve below is not a trend, so no " + "first->last number describes it - dose does not come back, so whatever moved " + "here was not dose; see the sweep-quality report"; os << fmt::format("\n per-batch relative-B (A^2, {:.0f} deg/batch):", sm.statistics.radiation_damage_batch_deg); for (float bb : sm.statistics.radiation_damage_b_batch) - os << fmt::format(" {:.1f}", bb); + os << (std::isfinite(bb) ? fmt::format(" {:.1f}", bb) : std::string(" -")); } logger.Info("{}", os.str()); result.radiation_damage_text = os.str(); diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 9462eb5a..2fe6c2fe 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -109,6 +109,30 @@ TEST_CASE("ResultReport_RenderEmpty", "[Diagnostics]") { CHECK(not_merged_text.find("FIRST_IMAGE LAST_IMAGE") != std::string::npos); } +TEST_CASE("ResultReport_RadiationDamage", "[Diagnostics]") { + // RADIATION_DAMAGE_RELATIVE_B is a number only when there is one. A curve that was measured but + // that no straight line describes, and a monitor that could not run at all, are different answers, + // and a consumer has to be able to tell them apart - and both from a measured zero. + DiffractionExperiment x(DetJF(1)); + ProcessResult result; + result.has_merge_statistics = true; + result.radiation_damage_text = "per-batch relative-B"; + result.merge_statistics.radiation_damage_batch_deg = 10.0; + result.merge_statistics.radiation_damage_b_batch = {0.0f, 4.0f, NAN}; + + result.merge_statistics.radiation_damage_delta_b = 8.25; + CHECK(RenderResultReport("p", "in.h5", x, result).find("\nRADIATION_DAMAGE_RELATIVE_B= 8.25\n") + != std::string::npos); + + result.merge_statistics.radiation_damage_delta_b = NAN; + CHECK(RenderResultReport("p", "in.h5", x, result).find("\nRADIATION_DAMAGE_RELATIVE_B= NOT_A_TREND\n") + != std::string::npos); + + result.merge_statistics.radiation_damage_b_batch.clear(); + CHECK(RenderResultReport("p", "in.h5", x, result).find("\nRADIATION_DAMAGE_RELATIVE_B= NOT_MEASURED\n") + != std::string::npos); +} + TEST_CASE("SweepQuality_HDF5RoundTrip", "[HDF5][Full][Diagnostics]") { // The per-image codes have to survive the writer and come back out of the reader. Without an // assertion here the field can ship as all-zeros without anyone noticing.