leonarski_fandClaude Opus 5
677ece7b59
The angles a rotation dataset stores are the COMMANDED ones, so a stage whose travel is miscalibrated leaves no trace in the header - every angle is self-consistently wrong. No existing parameter can absorb it either: the cell scale, the axis direction, the detector distance and the beam centre are all orthogonal to an error in rotation MAGNITUDE. So fit it as what it is - one scalar k, the ratio of the travel to the commanded angle - on the rocking events the geometry post-refinement already builds, after step A so the cell scale and the axis direction are fixed and k is the only free quantity. Two details decide whether the number means anything. The angle enters measured from the CENTRE of the sweep: the reference orientation was fitted against the commanded angles and has already absorbed their mean error, so measured from the goniometer's zero instead a constant missetting about the spindle leaks into k with a gain of <phi>/<phi^2>, which depends only on where the sweep happens to sit - on a short sweep starting near zero a 0.14 deg missetting fakes 1.4 % of k. Referred to the sweep centre that leak is identically zero at any width. And the robust loss is scaled to the scatter the events actually have, which varies by more than a decade between datasets, so any fixed constant is either inert or throws away real data. A stage fault is rare and a 1 % angle correction applied to a healthy dataset would damage it silently, so the correction is committed only when every test passes: at least 30 deg of sweep and 5000 events, |k-1| over 0.5 %, a misorientation of at least 0.5 deg at each end of the sweep, and the same k from every fifth of the sweep left out. The last test is not optional. A second lattice that dominates ONE END of a sweep - exactly what happens where the primary stops indexing - fakes a k that passes the other two, and the hkl-hash split used elsewhere in this file cannot see it, because both of its folds sit at the same angles and anything structured in phi survives in both. When it commits, the second pass re-integrates against the corrected angles. The pre-pass mosaicity is dropped with it: that is a width in degrees fitted against angles the second pass has just stopped using, and since the override can only ever raise the second pass's own estimate, carrying it over would hold the second pass at the rocking width the uncorrected angles produced - the correction half-applied. --rotation-scale asserts a known stage calibration by hand and overrides the fit. On the 38-crystal rotation battery the gate fires on exactly one dataset, at k = 1.01318 with 0.74 of that k surviving every fifth left out. The largest of the other 37 is 1.00211, which fails the end-error test; 34 of them sit below 1.0006. On the one that fires: R_meas 39.2 -> 23.9 % (XDS 37.1) CC1/2 86.5 -> 96.0 % (XDS 94.3) CC1/2 outer 1.4 -> 53.4 % (XDS 42.5) unique refl 40990 -> 41540 (XDS 41322) observations 74975 -> 103858 (XDS 129322) mosaicity 0.181 -> 0.159 deg which takes it from losing to XDS on R_meas, CC1/2 and outer-shell CC1/2 to beating it on all three, and the mosaicity drop is the inflation the uncorrected angles were producing. Its low-resolution R_meas is the one number that moves the wrong way, 12.0 -> 13.9 %, still well inside XDS's 18.3. No space group moves anywhere, and every other crystal's merge is unchanged beyond the two-pass loop's own jitter - measured here as the spread of the post-refined distance across arms that do not touch post-refinement at all, which is larger than anything this commit produces. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
288 lines
16 KiB
C++
288 lines
16 KiB
C++
// 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");
|
|
Key(os, "RADIATION_DAMAGE_RELATIVE_B", fmt::format("{:.2f}",
|
|
result.merge_statistics.radiation_damage_delta_b));
|
|
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());
|
|
}
|
|
}
|