mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 72 Wiki Activity
Files
8874a788e68054c4445694f424eca00a2c6af459
Jungfraujoch/rugnux/ResultReport.cpp
T
leonarski_fandClaude Opus 5 7ce47bd4d1 Radiation damage: report a measurement, or nothing 
The monitor fitted each batch's relative-B on SINGLE observations -
ln(I_ref/I_obs) regressed on s^2, weighted by (I_obs/sigma)^2, with the
logarithm requiring I_obs > 0.  The observation therefore sits in the
response and in its own weight, and the positivity requirement keeps only
the upward half of the noise, so the estimate is biased downwards wherever
I/sigma approaches 1 and is unbounded in the limit.  Simulated: on a batch
with no relative-B at all and <I/sigma> = 0.3 it reads -28 A^2; on a batch
whose true relative-B is +30 A^2 it reads -29.  The bias grows with dose,
so it inverts the answer on exactly the data the number exists for.

That is not a corner case.  Re-measured on stored integrated intensities,
a 360 deg sweep obstructed over a 60 deg wedge - whose honest curve is flat
for 130 deg, dips over the wedge and comes back - printed a per-batch curve
saturated at -31.4 A^2 for twelve consecutive batches (the +-50 A^2 clamp
less the low-dose anchor, "no data here" reported as a measurement) under a
headline of -19 A^2 of radiation damage.  A deliberately dosed dataset
printed -39 A^2 where the honest measurement is about +34: the one crystal
with real damage got the sign wrong.  Eight of thirty-eight datasets
reported |dB| > 5 A^2 and their curves oscillate by tens of A^2.

So pool the observations into ten equal-occupancy resolution shells per
batch before taking the logarithm, and fit slope AND intercept over the
shell means, weighting each shell by its own pooled (I/sigma)^2.  A shell
mean is well determined where a single observation is not, it admits
negative intensities, and it carries the I/sigma that says whether the
batch can be measured at all.  The same simulations then reproduce the
truth to under 1 A^2 at every signal level.  The intercept keeps a batch
that is merely dimmer than the run - an attenuated beam, a mis-fitted frame
scale - out of the damage number: a batch mis-scaled by 2x read +12 A^2 of
"damage" without it and +0.05 with it.

A batch whose shells are too weak to fit is now absent from the curve,
printed as "-", instead of pinned to the clamp.  And the clamp itself is
now an argument of the solve rather than one shared constant: it guards
against divergence, and the correction keeps the bound it was tuned with,
but with the estimator fixed a heavily dosed crystal's honest relative-B
runs past it - the monitor pinned thirteen consecutive batches at +49 A^2,
which is the same defect in the other direction.  The monitor is given room
a real relative-B cannot reach and drops any batch that lands on it anyway.

The shells are laid inside the range the run actually diffracted to,
not across the whole merged range: a resolution limit taken from another
program or left generous spends most of an equal-occupancy grid on noise
and leaves a batch with too few shells to fit at all - on the battery that
silenced three crystals outright and cost two of them nineteen batches of
thirty-six.  Where the merged range already sits inside the signal the grid
is unchanged and so is every number.

The first->last headline
is reported only where a straight line explains at least half of the
curve's variance, or where the curve is flat to within a couple of A^2 and
the answer is simply "no damage"; otherwise there is no headline and the
report says the loss was not dose and points at the sweep-quality section.
The three shapes separate cleanly - progressive damage R^2 0.97, the
obstructed sweep 0.24, the clean control flat at +0.35 A^2.  And the label
now follows the sign: damage fades the high-resolution intensity, so only a
positive change is dose, where before any |dB| > 5 was called damage.

Report-only throughout - the monitor never touches corr, and the per-batch
curve's only consumer beyond the report is a sweep-quality field no reader
reads; classification runs on the per-frame scale and CC, and is unmoved.
The decay correction's global slope and the opt-in per-batch relative-B
share this estimator and are left alone here: they fold into the scale, so
Full 38-crystal rotation battery, twice (the second confirming the shell
placement), against a clean baseline at the same base:

  space groups   unchanged at 35/38
  merge metrics  move on three crystals only - the same three whose two-pass
                 lattice search takes a different branch on nearly every arm run
                 this session, one of which moves its own R_meas by 1.5 points on
                 thread count alone

A report-only change ought to be bit-identical and this is not quite, which is
worth saying plainly: the three crystals that move are the known unstable ones
and no space group moves, but "identical except where nothing is ever identical"
is a weaker statement than "identical", and the residue has not been chased to
ground.

The three validation cases behave as they must:

  60 deg beam-obstructed wedge, no decay   -23.02, labelled damage, twelve
                                           batches printing the clamp
                                        -> NOT_A_TREND, curve within 3 A^2, the
                                           two unmeasurable batches absent, and a
                                           pointer to the sweep-quality section
  genuine progressive damage               -39.46, sign inverted
                                        -> +94.50, monotone, corroborated by a
                                           per-image CC that falls 0.608 -> 0.159
                                           and never recovers
  clean control                            +0.18 -> +0.76, flat within 1 A^2

Across the battery the report now names four crystals as radiation-damaged
instead of ten; the other three are the two lowest-energy datasets and the
pink-beam one, each showing a monotone rise of about ten square Angstroms.

Sweep-quality classification is untouched, and the coupling that was assumed to
exist does not: rad_damage_b_batch reaches it through one field that is written
and never read. Ranges and reasons are identical on 35 of 38, the three that
differ by one to seven frames are the same unstable crystals, and the census of
stretches called radiation damage is one before and one after.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-12 08:13:06 +02:00

294 lines
17 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <fstream>
#include <sstream>
#include <spdlog/fmt/fmt.h>
#include "../common/GitInfo.h"
#include "../common/time_utc.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "ResultReport.h"
namespace {
// The version of this file format. Bumped when a key is renamed or removed, a table column moves,
// or a reason code changes meaning - a consumer can gate on it.
constexpr int REPORT_VERSION = 1;
const char *BANNER = " ******************************************************************************";
void Section(std::ostream &os, const std::string &title) {
os << "\n" << BANNER << "\n " << title << "\n" << BANNER << "\n\n";
}
// Every number a consumer might want is written as one of these, so it is one grep away.
template <class T> void Key(std::ostream &os, const char *key, const T &value) {
os << key << "= " << value << "\n";
}
std::string CellString(const UnitCell &c) {
return fmt::format("{:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}", c.a, c.b, c.c,
c.alpha, c.beta, c.gamma);
}
}
std::string RenderResultReport(const std::string &output_prefix,
const std::string &input_file,
const DiffractionExperiment &experiment,
const ProcessResult &result) {
std::ostringstream os;
const bool rotation = experiment.IsRotationIndexing();
const bool merged = result.has_merge_statistics;
std::vector<std::string> warnings = result.warnings;
os << BANNER << "\n"
<< " RUGNUX PROCESSING REPORT\n"
<< BANNER << "\n\n"
<< " What this run determined, written next to its other output. The `KEY= value` lines and\n"
<< " the tables below are a stable interface - a script greps them, and REPORT_VERSION says\n"
<< " when that interface last changed. Timing, rates and per-image progress are not here;\n"
<< " they are on stdout.\n\n";
Key(os, "REPORT_VERSION", REPORT_VERSION);
Key(os, "RUGNUX_VERSION", jfjoch_version());
if (!jfjoch_git_sha1().empty())
Key(os, "RUGNUX_GIT", jfjoch_git_sha1().substr(0, 6) + " " + jfjoch_git_date());
Key(os, "DATE", time_UTC(std::chrono::system_clock::now()));
Key(os, "INPUT_FILE", input_file);
Key(os, "OUTPUT_PREFIX", output_prefix);
// ---------------------------------------------------------------- 1. DATA SET
Section(os, "1. DATA SET");
Key(os, "EXPERIMENT_TYPE", rotation ? "ROTATION" : "STILLS");
Key(os, "IMAGES_PROCESSED", result.images_processed);
Key(os, "WAVELENGTH", fmt::format("{:.5f}", experiment.GetWavelength_A()));
if (const auto gonio = experiment.GetGoniometer()) {
Key(os, "OSCILLATION_RANGE", fmt::format("{:.4f}", gonio->GetIncrement_deg()));
Key(os, "STARTING_ANGLE", fmt::format("{:.3f}", gonio->GetStart_deg()));
const auto ax = gonio->GetAxis();
Key(os, "ROTATION_AXIS", fmt::format("{:.6f} {:.6f} {:.6f}", ax.x, ax.y, ax.z));
}
Key(os, "DETECTOR_DISTANCE", fmt::format("{:.3f}", result.used_distance_mm));
Key(os, "BEAM_CENTRE", fmt::format("{:.2f} {:.2f}", result.used_beam_x_pxl, result.used_beam_y_pxl));
os << "\n"
<< " The distance and beam centre above are the ones this result was integrated at, which on\n"
<< " a rotation run is the post-refined geometry rather than the values in the input file.\n";
if (result.pass_count > 1) {
os << "\n";
Key(os, "PASS", fmt::format("{} of {}", result.pass_number, result.pass_count));
Key(os, "PASS_DECISION", result.pass_decision);
os << "\n"
<< " A rotation run integrates twice: once at the geometry in the input file, then again at\n"
<< " the post-refined geometry. Every number in this report describes the pass named above,\n"
<< " whose files are " << output_prefix << ".*; the header-geometry pass was written to\n"
<< " " << output_prefix << "_01.* and is kept only for comparison.\n";
}
// ---------------------------------------------------------------- 2. INDEXING
Section(os, "2. INDEXING");
if (result.indexing_rate.has_value())
Key(os, "INDEXING_RATE", fmt::format("{:.4f}", result.indexing_rate.value()));
Key(os, "LATTICE_FOUND", (result.consensus_cell.has_value() ? "TRUE" : "FALSE"));
if (result.consensus_cell.has_value())
Key(os, "UNIT_CELL_CONSTANTS", CellString(*result.consensus_cell));
if (result.space_group_number.has_value())
Key(os, "SPACE_GROUP_NUMBER", result.space_group_number.value());
if (result.indexing_rate.value_or(0.0f) <= 0.0f)
warnings.emplace_back("No image indexed - no crystal lattice was determined from this dataset");
// ---------------------------------------------- 3. GEOMETRY POST-REFINEMENT
if (result.post_refine.has_value()) {
const auto &pr = *result.post_refine;
Section(os, "3. GEOMETRY POST-REFINEMENT");
os << " The rotation two-pass fits the detector distance and beam centre from the observed spot\n"
<< " positions, and the cell scale and rotation axis from the observed rocking angles. Each\n"
<< " step is committed only if it improves a held-out residual.\n\n";
Key(os, "POSTREFINE_EVENTS_USED", pr.events_used);
Key(os, "POSTREFINE_OBS_USED", pr.obs_used);
Key(os, "POSTREFINE_CELL_COMMITTED", pr.cell_refined ? "TRUE" : "FALSE");
Key(os, "POSTREFINE_DETECTOR_COMMITTED", pr.detector_refined ? "TRUE" : "FALSE");
Key(os, "POSTREFINE_DISTANCE", fmt::format("{:.3f} -> {:.3f}", pr.distance_before_mm,
pr.distance_after_mm));
Key(os, "POSTREFINE_BEAM_CENTRE", fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}",
pr.beam_x_before_px, pr.beam_y_before_px,
pr.beam_x_after_px, pr.beam_y_after_px));
Key(os, "GONIOMETER_ROTATION_SCALE", fmt::format("{:.5f}", pr.rotation_scale));
Key(os, "GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect ? "TRUE" : "FALSE");
os << "\n GONIOMETER_ROTATION_SCALE is the factor by which the stage actually turned relative to\n"
<< " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It drives\n"
<< " the second integration pass only when SUSPECT is TRUE - both cross-validated and outside\n"
<< " the tolerance - since a stage that is in fact well calibrated must be left alone. A\n"
<< " manual --rotation-scale replaces it and is applied to both passes.\n";
if (pr.rotation_scale_suspect)
warnings.emplace_back(fmt::format(
"The goniometer turned by a factor {:.5f} of the angles stored in the file - the "
"stage rotation looks mis-calibrated by {:+.2f}%. The correction was applied to this "
"run, but the fault is in the hardware and should be fixed there",
pr.rotation_scale, 100.0 * (pr.rotation_scale - 1.0)));
}
// ---------------------------------------------- 4. SPACE GROUP DETERMINATION
Section(os, "4. SPACE GROUP DETERMINATION");
if (result.space_group_search.has_value()) {
Key(os, "SPACE_GROUP_SEARCH", "DE_NOVO");
os << "\n" << SearchSpaceGroupResultToText(*result.space_group_search) << "\n";
} else if (result.space_group_number.has_value()) {
Key(os, "SPACE_GROUP_SEARCH", "FIXED");
os << "\n The space group was given, not determined here.\n";
} else {
Key(os, "SPACE_GROUP_SEARCH", "NONE");
os << "\n No space group was determined.\n";
}
// ---------------------------------------------------- 5. SCALING AND MERGING
Section(os, "5. SCALING AND MERGING");
if (!merged) {
Key(os, "MERGE", "NOT_PERFORMED");
os << "\n No scaling or merging was performed on this run, so there are no merging statistics, no\n"
<< " error model, and no sweep-quality diagnosis below. The integrated reflections are in\n"
<< " " << output_prefix << "_process.h5.\n";
} else {
const auto &o = result.merge_statistics.overall;
Key(os, "MERGE", "PERFORMED");
Key(os, "INCLUDE_RESOLUTION_RANGE", fmt::format("{:.3f} {:.3f}", o.d_max, o.d_min));
Key(os, "FRIEDELS_LAW", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE");
Key(os, "UNIQUE_REFLECTIONS", o.unique_reflections);
Key(os, "TOTAL_OBSERVATIONS", o.total_observations);
Key(os, "COMPLETENESS", o.possible_unique_reflections > 0
? fmt::format("{:.1f}", 100.0 * o.unique_reflections / o.possible_unique_reflections)
: std::string("nan"));
Key(os, "MULTIPLICITY", o.unique_reflections > 0
? fmt::format("{:.2f}", static_cast<double>(o.total_observations) / o.unique_reflections)
: std::string("nan"));
Key(os, "I_OVER_SIGMA", fmt::format("{:.2f}", o.mean_i_over_sigma));
Key(os, "R_MEAS", fmt::format("{:.4f}", o.r_meas));
Key(os, "CC_HALF", fmt::format("{:.4f}", o.cc_half));
Key(os, "SIGANO", fmt::format("{:.3f}", o.abs_diff_over_sigma_anomalous));
Key(os, "WILSON_B", fmt::format("{:.2f}", result.merge_statistics.wilson_b));
// The error model in XDS's convention, so the numbers are directly comparable with a CORRECT.LP.
Key(os, "ERROR_MODEL_A", fmt::format("{:.4f}", result.error_model_a));
Key(os, "ERROR_MODEL_B", fmt::format("{:.4e}", result.error_model_b));
Key(os, "ISA", fmt::format("{:.2f}", result.error_model_isa));
if (result.error_model_isa_asymptotic > 0.0)
Key(os, "ISA_ASYMPTOTIC", fmt::format("{:.2f}", result.error_model_isa_asymptotic));
Key(os, "REFERENCE_DATA_USED", result.has_reference ? "TRUE" : "FALSE");
// The shell table straight off the statistics rather than result.merge_statistics_text: that
// string also carries the twinning analysis and the advisories, which have sections of their own.
os << "\n ERROR_MODEL_A / ERROR_MODEL_B are in XDS's convention, sigma^2 = a*(sigma0^2 + b*I^2),\n"
<< " so ISA = 1/sqrt(a*b) means what CORRECT.LP's ISa means. ISA_ASYMPTOTIC, where present,\n"
<< " is the strong-reflection tier only.\n\n"
<< result.merge_statistics;
}
// --------------------------------------------------------------- 6. TWINNING
if (merged && result.twinning.l_test_pairs > 0) {
Section(os, "6. TWINNING");
Key(os, "TWINNING_SUSPECTED", result.twinning.twinning_suspected ? "TRUE" : "FALSE");
Key(os, "L_TEST_MEAN_ABS_L", fmt::format("{:.4f}", result.twinning.mean_abs_l));
Key(os, "L_TEST_MEAN_L_SQUARED", fmt::format("{:.4f}", result.twinning.mean_l_squared));
Key(os, "SECOND_MOMENT_I", fmt::format("{:.4f}", result.twinning.second_moment));
Key(os, "ESTIMATED_TWIN_FRACTION", fmt::format("{:.3f}", result.twinning.estimated_twin_fraction));
os << "\n" << TwinningAnalysisToText(result.twinning) << "\n";
if (result.twinning.twinning_suspected)
warnings.emplace_back(fmt::format(
"Twinning is indicated (<|L|> = {:.3f}, <I^2>/<I>^2 = {:.3f}, estimated twin "
"fraction {:.2f}) - refine against the merged data with care",
result.twinning.mean_abs_l, result.twinning.second_moment,
result.twinning.estimated_twin_fraction));
}
// ------------------------------------------------------- 7. RADIATION DAMAGE
if (!result.radiation_damage_text.empty()) {
Section(os, "7. RADIATION DAMAGE");
// 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";
}
// ------------------------------------------------------ 8. SWEEP QUALITY
const auto &sq = result.merge_statistics.sweep_quality;
Section(os, "8. SWEEP QUALITY");
os << " Stretches of the sweep over which the crystal delivered much less than the rest of the run.\n"
<< " REASON comes from a closed vocabulary, listed below so a consumer can tell an unknown code\n"
<< " from a missing one. SEVERITY is the fraction of the run's typical diffracting power missing\n"
<< " over the range (0 = as good as the run, 1 = nothing at all); SCALE and CC are the range's\n"
<< " mean per-image scale and CC-to-merge relative to the run median; INDEXED is the fraction of\n"
<< " the range's frames that were scaled at all. Nothing is excluded on the strength of this.\n\n";
Key(os, "SWEEP_QUALITY_STATUS", sq.measured ? "COMPUTED" : "NOT_COMPUTED");
Key(os, "SWEEP_QUALITY_COUNT", sq.ranges.size());
{
std::string codes;
for (int r = 0; r <= static_cast<int>(SweepQualityReason::RadiationDamage); ++r)
codes += (codes.empty() ? "" : " ")
+ std::string(SweepQualityReasonCode(static_cast<SweepQualityReason>(r)));
Key(os, "SWEEP_QUALITY_REASONS", codes);
}
if (sq.measured) {
Key(os, "SWEEP_ROTATION", fmt::format("{:.1f}", sq.sweep_deg));
Key(os, "FLUX_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.flux_peak_to_trough));
Key(os, "SCALE_MODULATION_PEAK_TO_TROUGH", fmt::format("{:.2f}", sq.modulation_peak_to_trough));
}
os << "\n"
<< " FIRST_IMAGE LAST_IMAGE N_IMAGES ROTATION REASON SEVERITY SCALE CC INDEXED\n"
<< " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
for (const auto &r : sq.ranges) {
os << fmt::format(" {:11d} {:11d} {:9d} {:8.1f} {:<20} {:8.2f} {:6.2f} {:6.2f} {:8.2f}\n",
r.first_image, r.last_image, r.last_image - r.first_image + 1, r.rotation_deg,
SweepQualityReasonCode(r.reason), r.severity, r.mean_relative_scale,
r.mean_relative_cc, r.indexed_fraction);
warnings.push_back(fmt::format(
"Frames {}-{} {} ({:.1f} deg, scale {:.2f} and CC {:.2f} of the run, {:.0f}% scaled)",
r.first_image, r.last_image, SweepQualityReasonText(r.reason), r.rotation_deg,
r.mean_relative_scale, r.mean_relative_cc, 100.0 * r.indexed_fraction));
}
os << " ----------- ----------- --------- -------- -------------------- -------- ------ ------ --------\n";
// --------------------------------------------------------------- 9. WARNINGS
if (result.cancelled)
warnings.emplace_back(fmt::format("Processing was cancelled after {} images - this report "
"describes an incomplete run", result.images_processed));
Section(os, "9. WARNINGS");
os << " Everything that needs a person's attention, one line each, marked so a script can find\n"
<< " them with a single grep for \"WARNING:\".\n\n";
Key(os, "WARNING_COUNT", warnings.size());
os << "\n";
for (const auto &w : warnings)
os << "WARNING: " << w << "\n";
if (warnings.empty())
os << " (none)\n";
os << "\n" << BANNER << "\n END OF REPORT\n" << BANNER << "\n";
return os.str();
}
void WriteResultReport(const std::string &output_prefix,
const std::string &input_file,
const DiffractionExperiment &experiment,
const ProcessResult &result,
Logger &logger) {
if (output_prefix.empty())
return; // "compute the statistics, persist nothing"
const std::string filename = output_prefix + "_report.txt";
// The report is unconditional, so it must never be the reason a run fails: a run that produced a
// good .mtz must survive an unwritable path or a full disk. Report the failure and carry on.
try {
std::ofstream file(filename);
file.exceptions(std::ios::failbit | std::ios::badbit);
file << RenderResultReport(output_prefix, input_file, experiment, result);
} catch (const std::exception &e) {
logger.Warning("Could not write the results report {}: {}", filename, e.what());
}
}
Reference in New Issue View Git Blame Copy Permalink