From f4620470383572c5c67f1dd8dc90d5437b34ea25 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:07:22 +0200 Subject: [PATCH 1/6] rugnux report: BEAM_CENTRE_OFF warning when the file's beam centre is wrong The pre-scan's beam-centre check and the post-refinement both measure where the direct beam lands, but the report only carried the centre the run used. Section 1 now sets the file's direct beam (FILE_DIRECT_BEAM) beside the background's (BACKGROUND_DIRECT_BEAM, with its sigma), the shift of the refined centre from the file's (BEAM_CENTRE_SHIFT) and the tolerance the geometry asks for (BEAM_CENTRE_TOLERANCE: BeamCenterNeed_pxl, the number the log's check prints, here on the run's own cell). A shift above the tolerance is a BEAM_CENTRE_OFF warning naming the file's, refined and background centres, the difference in px and mm and the value to correct the beamline to; a run that refined no centre is warned on the background alone where its fit resolves the difference (3 sigma). The summary gains a Beam centre line. The refined centre, not the background, is the trigger wherever one exists: on a battery set where XDS's refined direct beam and rugnux's agree to 1 px, the background fit is 12 px away from both. Report-only; no result changes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX_REPORT.md | 26 +++++++++++++-- rugnux/ReportDocument.h | 1 + rugnux/ResultReport.cpp | 67 ++++++++++++++++++++++++++++++++++++++ rugnux/Rugnux.cpp | 9 +++++ rugnux/Rugnux.h | 7 ++++ tests/ResultReportTest.cpp | 36 ++++++++++++++++++++ 6 files changed, 144 insertions(+), 2 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index ffaacd498..63e73672d 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -128,6 +128,28 @@ JFJOCH_DATASET_SETTINGS= {"beam_x_pxl": 2078.24, "beam_y_pxl": 2233.92, "detecto grep '^JFJOCH_DATASET_SETTINGS=' out_report.txt | cut -d' ' -f2- > geometry.json ``` +**The beam centre in the file.** Section 1 sets the file's centre beside the two the run measured: +**`FILE_DIRECT_BEAM=`** is where the input file puts the direct beam (its beam centre, under its own +detector tilt), **`BACKGROUND_DIRECT_BEAM=`** where the scattered background puts it, with the fit's +sigma (the pre-scan's beam-centre check), and **`BEAM_CENTRE_SHIFT=`** how far the centre the result was +integrated at (`DIRECT_BEAM=`, the refined one on a rotation run) lies from the file's, in pixels. +**`BEAM_CENTRE_TOLERANCE=`** is how far this crystal on this geometry lets the centre be wrong - the +displacement at which the accumulated reflections of the longest cell axis smear to 70% of their +height, the number the log's `Beam centre check` line states, here on the run's own cell. Where the +shift is larger, a **`BEAM_CENTRE_OFF`** warning gives the file's, the refined and the background +centre, the difference in pixels and millimetres, and the value to correct the beamline (or the +header) to: the run itself integrated at the refined centre, but every other program given the same +file starts from the wrong one. A run that refined no centre (stills) is warned on the background +alone, where its fit resolves the difference (three sigma) and the difference exceeds the tolerance, +worded as a check. The summary's `Beam centre` line states the outcome either way. + +``` +FILE_DIRECT_BEAM= 1564.31 1702.91 +BACKGROUND_DIRECT_BEAM= 1571.30 1699.85 0.62 +BEAM_CENTRE_SHIFT= 7.15 +BEAM_CENTRE_TOLERANCE= 3.49 +``` + **Which pass.** A rotation run integrates twice — once at the geometry in the input file (the first half of the sweep but at least 120°, `--prepass-fraction`), then the whole sweep again at the post-refined geometry — and can integrate a third time if a guard rejects the second pass. There is **one** report, for the pass that became the canonical output, and `PASS=` / @@ -156,13 +178,13 @@ and a run that produced a textbook data set read identically for their first thr `GONIO_SCALE`, `SPINDLE_CAP`, `ANISOTROPY`, `TWINNING`, `PSEUDO_TRANSLATION`, `LATTICE_TRANSLATION`, `MODEL_HAND`, `MODEL_NOT_VALIDATED`, `REFERENCE_MISMATCH`, `CANCELLED`, `RESOLUTION_FIT`, `FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`, `SCALING_NOT_CONVERGED`, `HARMONIC_CONTAMINATION`, - `SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `ICE_RINGS`, `POWDER_RINGS`. `NONE` when nothing fired. A code appears if and only if its + `SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `ICE_RINGS`, `POWDER_RINGS`, `BEAM_CENTRE_OFF`. `NONE` when nothing fired. A code appears if and only if its warning fired, so the flags and the `WARNING:` lines are two renderings of one list. Thresholds are deliberately low: a warning is a prompt to check something, not a verdict, and some of them fire on data that turn out fine — the closed type for machinery, the open sentence for a person. - Below them, one line each, the facts a reader needs before reading further: space group, cell, - mosaicity, the powder and ice rings, multiple lattices, resolution, completeness, signal, anomalous + mosaicity, the beam centre, the powder and ice rings, multiple lattices, resolution, completeness, signal, anomalous signal, anisotropy, pseudo-symmetry, twinning, supercell, radiation damage and the sweep. A line states a finding that is not a condition too - weak powder rings, a possible supercell no stronger than correct cells show, an arbitrary indexing choice (`Indexing choice`, a rotation sweep whose cell diff --git a/rugnux/ReportDocument.h b/rugnux/ReportDocument.h index eeee5e620..5f1361317 100644 --- a/rugnux/ReportDocument.h +++ b/rugnux/ReportDocument.h @@ -51,6 +51,7 @@ namespace PathologyCode { constexpr const char *SCALING_NOT_CONVERGED = "SCALING_NOT_CONVERGED"; // the per-frame scales hit their cap still moving constexpr const char *POWDER_RINGS = "POWDER_RINGS"; // non-ice powder rings hold many spots constexpr const char *ICE_RINGS = "ICE_RINGS"; // hexagonal-ice rings hold many spots + constexpr const char *BEAM_CENTRE_OFF = "BEAM_CENTRE_OFF"; // the file's beam centre is wrong for this geometry } // A closed type plus an open description. The type is what machinery switches on; the text is what a diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 7a752dae1..48b8847c7 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -42,6 +42,8 @@ namespace { // UNDECLARED_LATTICE_TRANSLATION_PCT and INDEXING_AMBIGUITY_OPERATORS; INDEXING_RATE is a // --developer key on rotation runs; POWDER_ICE_* / POWDER_NON_ICE_SPOT_FRACTION; // the ICE_RINGS and POWDER_RINGS flags. + // 16: FILE_DIRECT_BEAM / BACKGROUND_DIRECT_BEAM / BEAM_CENTRE_SHIFT / BEAM_CENTRE_TOLERANCE in + // section 1 and the BEAM_CENTRE_OFF flag - the file's beam centre against the refined one. // 8: every line that is not `KEY= value` and not blank starts with `#`, warnings included // (`# WARNING:`); the header gains RUGNUX_DOWNLOAD and BUILD_CXX_FLAGS, and RUGNUX_GIT is // the full build-time hash. @@ -408,6 +410,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, const bool rotation = experiment.IsRotationIndexing(); const bool merged = result.has_merge_statistics; doc.warnings = result.warnings; + // SUMMARY rows filled by the sections that measure them; empty leaves the row out. + std::string beam_centre_row; // ------------------------------------------------------------------ header block { @@ -543,6 +547,67 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, "integrated at is the best one reached, not a settled one (see PASS_DECISION)"); } + // The beam centre against the file's. The file's is read off the experiment the run was + // given - the header, with the header's tilt - and compared with where the result was + // integrated (the refined geometry) and with where the scattered background put the beam. The + // run corrects a wrong centre for itself; every other program reading the same file does not, + // which is why this is worth a warning however well the run went. The tolerance is the one the + // pre-scan's check prints (BeamCenterNeed_pxl), taken on the run's own cell. + { + const auto file_direct = experiment.GetDiffractionGeometry().GetDirectBeam_pxl(); + const double pixel_mm = experiment.GetPixelSize_mm(); + const double need = result.beam_centre_need_pxl; + const double moved = std::hypot(result.used_direct_beam_x_pxl - file_direct.first, + result.used_direct_beam_y_pxl - file_direct.second); + const auto &bg = result.background_beam_centre; + const double bg_moved = bg ? std::hypot(bg->beam_x_pxl - file_direct.first, + bg->beam_y_pxl - file_direct.second) : NAN; + Add(s, KeyText("FILE_DIRECT_BEAM", fmt::format("{:.2f} {:.2f}", file_direct.first, + file_direct.second))); + if (bg) + Add(s, KeyText("BACKGROUND_DIRECT_BEAM", fmt::format("{:.2f} {:.2f} {:.2f}", bg->beam_x_pxl, + bg->beam_y_pxl, bg->sigma_pxl))); + Add(s, KeyReal("BEAM_CENTRE_SHIFT", moved, "{:.2f}")); + if (std::isfinite(need)) + Add(s, KeyReal("BEAM_CENTRE_TOLERANCE", need, "{:.2f}")); + // The refined geometry is the evidence wherever the run refined one: the background fit can + // be pulled by an asymmetric background, and where both exist and disagree the spot + // positions are what the result rests on. Only a run that refined nothing (stills, no + // post-refinement) is warned on the background alone, and only where the fit resolves the + // difference at all (three of its own sigma, as the pre-scan's check reads it). + const bool bg_resolves = bg && bg_moved > 3.0 * bg->sigma_pxl; + const std::string bg_text = bg ? fmt::format("the scattered background at ({:.2f}, {:.2f}) +- {:.2f} px", + bg->beam_x_pxl, bg->beam_y_pxl, bg->sigma_pxl) + : std::string(); + if (std::isfinite(need) && moved > need) { + beam_centre_row = fmt::format("the file's is OFF by {:.1f} px ({:.2f} mm) - correct the " + "beamline value", moved, moved * pixel_mm); + Warn(doc, PathologyCode::BEAM_CENTRE_OFF, fmt::format( + "the beam centre in the input file is {:.1f} px ({:.2f} mm) from where this run found " + "the direct beam, more than the {:.1f} px this geometry needs it to be right to: the " + "file puts the direct beam at ({:.2f}, {:.2f}), the refined geometry at ({:.2f}, " + "{:.2f}){}. This run used the refined centre, but any other program given this file " + "starts from the wrong one - correct the beam centre at the beamline (or in the " + "file's header) to ({:.2f}, {:.2f})", + moved, moved * pixel_mm, need, file_direct.first, file_direct.second, + result.used_direct_beam_x_pxl, result.used_direct_beam_y_pxl, + bg ? ", " + bg_text : std::string(), + result.used_direct_beam_x_pxl, result.used_direct_beam_y_pxl)); + } else if (std::isfinite(need) && !result.post_refine && bg_resolves && bg_moved > need) { + beam_centre_row = fmt::format("the background puts it {:.1f} px ({:.2f} mm) from the file's - " + "check the beamline value", bg_moved, bg_moved * pixel_mm); + Warn(doc, PathologyCode::BEAM_CENTRE_OFF, fmt::format( + "the beam centre in the input file may be off: the file puts the direct beam at " + "({:.2f}, {:.2f}), {} - {:.1f} px ({:.2f} mm) away, more than the {:.1f} px this " + "geometry needs the centre to be right to. Nothing on this run refined the centre, so " + "this is the background's word alone: check the beam centre at the beamline and " + "correct it if this is confirmed", + file_direct.first, file_direct.second, bg_text, bg_moved, bg_moved * pixel_mm, need)); + } else if (std::isfinite(need)) + beam_centre_row = fmt::format("the file's is within the {:.1f} px this geometry needs ({:.1f} px " + "from the refined one)", need, moved); + } + Add(s, Blank()); // Where the geometry came from. This is the whole provenance story a user needs, because a // user acts on exactly one geometry: the one the result was produced with. @@ -2377,6 +2442,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, row("Indexing choice", fmt::format("alternative indexing {} - one sweep is indexed consistently, " "so it matters only against other data: index against those " "with -z or --model", result.indexing_ambiguity_operators)); + if (!beam_centre_row.empty()) + checked_row("Beam centre", {PathologyCode::BEAM_CENTRE_OFF}, beam_centre_row); if (!result.powder.rings_q_recipA.empty() || result.consensus_cell) { const auto &pw = result.powder; checked_row("Powder", {PathologyCode::POWDER_RINGS}, pw.rings_q_recipA.empty() ? std::string("no rings detected") diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 804f78bdd..bfa4dd921 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -434,8 +434,17 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // on the rotation two-pass that is under a third of what the run actually took, the pre-scan and the // first pass being invisible in it. Time the whole thing here, where every pass is inside. const auto run_start = std::chrono::steady_clock::now(); + DiffractionExperiment header_experiment = experiment_; ProcessResult result = RunAllPasses(observer); result.total_time_s = std::chrono::duration(std::chrono::steady_clock::now() - run_start).count(); + // For the report's beam-centre warning: what the background measured, and what this crystal asks + // of the centre - at the file's geometry, as the pre-scan's check asks it, but on the cell the run + // determined rather than the long axis that check has to assume. + JoinBeamCenterCapture(); + result.background_beam_centre = measured_beam_center_; + if (result.consensus_cell) + header_experiment.SetUnitCell(*result.consensus_cell); + result.beam_centre_need_pxl = BeamCenterNeed_pxl(header_experiment); return result; } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4df80b6b8..2c8068da2 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -530,6 +530,13 @@ struct ProcessResult { // did not post-refine (stills, single pass, or the refinement did not converge). std::optional post_refine; + // The beam centre the scattered background puts the direct beam at (the pre-scan's check), and + // how far the centre may be wrong on this crystal (BeamCenterNeed_pxl, with the run's cell where + // it has one). Report only: the file's centre is read off the experiment the report is given, the + // refined one off used_direct_beam_*. Empty / NaN where the check did not run. + std::optional background_beam_centre; + float beam_centre_need_pxl = NAN; + // Which pass this result is (1-based) and how many ran, plus one line saying which pass became the // canonical output and on what evidence. pass_count == 1 outside the rotation two-pass. int pass_number = 1; diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 1f964ed5c..010930e83 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -1200,3 +1200,39 @@ TEST_CASE("ResultReport_Checks", "[Diagnostics]") { CHECK(std::ranges::any_of(c.codes, [&](const std::string &code) { return std::ranges::find(doc.pathology_flags, code) != doc.pathology_flags.end(); })); } + +TEST_CASE("ResultReport_BeamCentreOff", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + x.BeamX_pxl(500.0f).BeamY_pxl(400.0f); + const auto file = x.GetDiffractionGeometry().GetDirectBeam_pxl(); + ProcessResult result; + result.post_refine = PostRefineResult{}; + result.beam_centre_need_pxl = 2.0f; + result.background_beam_centre = BeamCenterEstimate{.beam_x_pxl = static_cast(file.first) + 5.0f, + .beam_y_pxl = static_cast(file.second), + .sigma_pxl = 0.5f}; + + // Refined within the tolerance: the background's word is not enough where the centre was refined. + result.used_direct_beam_x_pxl = file.first + 1.0; + result.used_direct_beam_y_pxl = file.second; + auto text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("BEAM_CENTRE_OFF") == std::string::npos); + CHECK(text.find("\nBEAM_CENTRE_SHIFT= 1.00\n") != std::string::npos); + CHECK(text.find("\nBEAM_CENTRE_TOLERANCE= 2.00\n") != std::string::npos); + + // Refined 4 px off the file's centre: warned, with both centres, px and mm, and the advice. + result.used_direct_beam_x_pxl = file.first; + result.used_direct_beam_y_pxl = file.second + 4.0; + text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("\nPATHOLOGY_FLAGS= BEAM_CENTRE_OFF\n") != std::string::npos); + CHECK(text.find("4.0 px (0.30 mm)") != std::string::npos); + CHECK(text.find("correct the beam centre at the beamline") != std::string::npos); + + // Nothing refined (stills): the background alone warns, worded as a check. + result.post_refine.reset(); + result.used_direct_beam_x_pxl = file.first; + result.used_direct_beam_y_pxl = file.second; + text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("\nPATHOLOGY_FLAGS= BEAM_CENTRE_OFF\n") != std::string::npos); + CHECK(text.find("the background's word alone") != std::string::npos); +} From 2632b77eb5203c11bdeda1cf2cb53698b879f925 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:09:26 +0200 Subject: [PATCH 2/6] rugnux report: count overloads and overlapping spots, warn with advice Overloads: the integrated reflections whose signal disk held a saturated pixel are counted on the integration outcome before the merge, over the whole range and over the reflections coarser than 4 A, where they concentrate (log line "Overloads:", report OVERLOADED_REFLECTIONS*). An OVERLOADS warning fires from 0.05% of all, or 0.25% of the low-resolution reflections, and recommends attenuation / a shorter exposure. Measured: a strongly diffracting in-house sweep reads 0.017% / 0.42% (warned), an ordinary sweep with a few hundred overloaded partials 0.001% / 0.012%. Overlaps: both integration engines count the signal disks a nearer neighbour's prediction takes pixels from (the overlap ownership the default --overlap exclude already builds), and those left under --overlap-minpk of their own - BraggIntegrationCounts::overlapped / overlapped_major (log line "Overlapping spots:", report OVERLAPPED_REFLECTIONS*). An OVERLAPS warning fires from 1% in the second and recommends a longer distance, finer slicing or a smaller beam; on 10 battery sets it fires on the two long-axis ones (8.3% and 13.4%), the rest read 0.0-0.44%. Summary lines for both; documented in RUGNUX_REPORT.md. Counters only; nothing reads them. p.mtz md5-identical to the rc175 battery on 10 sets. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX_REPORT.md | 38 ++++++++++- image_analysis/MXAnalysisWithoutFPGA.cpp | 4 ++ .../BraggIntegrationEngine.h | 6 ++ .../BraggIntegrationEngineCPU.cpp | 5 ++ .../BraggIntegrationEngineGPU.cu | 12 +++- rugnux/ReportDocument.h | 2 + rugnux/ResultReport.cpp | 64 ++++++++++++++++++- rugnux/Rugnux.h | 15 +++++ rugnux/RugnuxImagePass.cpp | 14 +++- rugnux/RugnuxScaleMerge.cpp | 26 ++++++++ tests/ResultReportTest.cpp | 27 ++++++++ 11 files changed, 207 insertions(+), 6 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 63e73672d..d992bf483 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -178,13 +178,13 @@ and a run that produced a textbook data set read identically for their first thr `GONIO_SCALE`, `SPINDLE_CAP`, `ANISOTROPY`, `TWINNING`, `PSEUDO_TRANSLATION`, `LATTICE_TRANSLATION`, `MODEL_HAND`, `MODEL_NOT_VALIDATED`, `REFERENCE_MISMATCH`, `CANCELLED`, `RESOLUTION_FIT`, `FLIGHT_PATH`, `GEOMETRY_NOT_CONVERGED`, `SCALING_NOT_CONVERGED`, `HARMONIC_CONTAMINATION`, - `SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `ICE_RINGS`, `POWDER_RINGS`, `BEAM_CENTRE_OFF`. `NONE` when nothing fired. A code appears if and only if its + `SUPERCELL_POSSIBLE`, `MULTIPLE_LATTICES`, `ICE_RINGS`, `POWDER_RINGS`, `BEAM_CENTRE_OFF`, `OVERLOADS`, `OVERLAPS`. `NONE` when nothing fired. A code appears if and only if its warning fired, so the flags and the `WARNING:` lines are two renderings of one list. Thresholds are deliberately low: a warning is a prompt to check something, not a verdict, and some of them fire on data that turn out fine — the closed type for machinery, the open sentence for a person. - Below them, one line each, the facts a reader needs before reading further: space group, cell, - mosaicity, the beam centre, the powder and ice rings, multiple lattices, resolution, completeness, signal, anomalous + mosaicity, the beam centre, the powder and ice rings, overloads, overlapping spots, multiple lattices, resolution, completeness, signal, anomalous signal, anisotropy, pseudo-symmetry, twinning, supercell, radiation damage and the sweep. A line states a finding that is not a condition too - weak powder rings, a possible supercell no stronger than correct cells show, an arbitrary indexing choice (`Indexing choice`, a rotation sweep whose cell @@ -519,6 +519,40 @@ scaling. The summary has an `Ice` line, and where the merge found ice an `ICE_RI look at the ice-ring shells (one warning line, whichever of the two measurements raised it). The measurement is described in [CPU/GPU data analysis ▸ Resolution and ice-ring handling](CPU_DATA_ANALYSIS_IMAGE.md#33-resolution-and-ice-ring-handling). +## Overloads and overlapping spots + +Two ways the measurement, rather than the crystal, cost reflections - each with what to change at the +next collection. Both are in section 4, both have a summary line, and both warn at a low threshold. + +``` +OVERLOADED_REFLECTIONS= 1740 +OVERLOADED_REFLECTIONS_PCT= 0.017 +OVERLOADED_REFLECTIONS_LOW_RES_PCT= 0.418 +OVERLAPPED_REFLECTIONS_PCT= 11.10 +OVERLAPPED_REFLECTIONS_MAJOR_PCT= 0.44 +``` + +(a strongly diffracting lysozyme sweep: an `OVERLOADS` warning, no `OVERLAPS` one) + +**`OVERLOADED_REFLECTIONS=`** counts the integrated reflections - every frame's partial, before the +merge - with a saturated pixel in the signal disk (or one the detector could not count on that frame), +and **`OVERLOADED_REFLECTIONS_PCT=`** their share. They are not measurements: the merge leaves them out +(a rotation merge drops the whole rocking event, `OBSERVATIONS_REJECTED_OVERLOAD=` in section 3), and +they are the strongest reflections the crystal gave, so they concentrate at low resolution; +**`OVERLOADED_REFLECTIONS_LOW_RES_PCT=`** is their share among the reflections coarser than 4 Å. An +**`OVERLOADS`** warning fires from 0.05% of all integrated reflections or 0.25% of those coarser than 4 Å, +and recommends attenuating the beam or shortening the exposure (or a second, attenuated low-resolution +sweep). + +**`OVERLAPPED_REFLECTIONS_PCT=`** is the share of the predicted reflections whose signal disk a nearer +neighbour's prediction takes pixels from, as the integrator's overlap treatment (`--overlap`, default +`exclude`) assigns them, and **`OVERLAPPED_REFLECTIONS_MAJOR_PCT=`** the share left with less than +`--overlap-minpk` (default 75%) of the disk its own - the ones XDS's MINPK would reject. Under +`exclude` the shared pixels are left out of the profile fit, which keeps the reflection at a cost in +precision. An **`OVERLAPS`** warning fires from 1% in the second number and recommends a longer +detector distance, finer rotation slicing or a smaller beam. Neither key is written with +`--overlap off`, which assigns no pixel to anyone. + ## Index-2 superstructure On rotation data, section 4 reports whether the lattice the run adopted has intensity at half-integer diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 6933bd473..566782e62 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -201,6 +201,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, probe_counts.bkg_starved += after.bkg_starved - before.bkg_starved; probe_counts.bkg_starved_by_neighbour += after.bkg_starved_by_neighbour - before.bkg_starved_by_neighbour; probe_counts.profile_fallback += after.profile_fallback - before.profile_fallback; + probe_counts.overlapped += after.overlapped - before.overlapped; + probe_counts.overlapped_major += after.overlapped_major - before.overlapped_major; return ret; }; @@ -439,6 +441,8 @@ BraggIntegrationCounts MXAnalysisWithoutFPGA::BraggCounts() const { c.bkg_starved -= probe_counts.bkg_starved; c.bkg_starved_by_neighbour -= probe_counts.bkg_starved_by_neighbour; c.profile_fallback -= probe_counts.profile_fallback; + c.overlapped -= probe_counts.overlapped; + c.overlapped_major -= probe_counts.overlapped_major; return c; } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 981af9e6a..bc62af4b1 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -134,6 +134,12 @@ struct BraggIntegrationCounts { // costs data adds 4.1 percentage points on top of a floor of 0.02%. uint64_t bkg_starved_by_neighbour = 0; uint64_t profile_fallback = 0; // the profile fit disagreed with the box sum; the box sum was kept + // Signal disks a nearer neighbour's prediction takes pixels from (the overlap ownership, so counted + // only while an overlap treatment is on), and of those the ones left with less than + // overlap_min_peak of the disk their own - what --overlap reject drops, and XDS's MINPK. Reported + // only (the crowding of the pattern); nothing reads them. + uint64_t overlapped = 0; + uint64_t overlapped_major = 0; }; // One reflection's extracted intensity, produced by the derived engine and turned into a diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 6cc8b8287..3f57f5744 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -269,6 +269,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } } + if (use_owner && n_own < n_disk) { + ++counts.overlapped; + if (n_own < overlap_min_peak * n_disk) ++counts.overlapped_major; + } + int n_bkg_used = n_bkg; // pixels behind the FINAL background value (trim/clip shrink it) // A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the // reflection outright. A profile fit does not need it to: the fit is the amplitude of a diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 2acdbd862..0f458ce6e 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -49,7 +49,9 @@ constexpr int MOM_STRIDE = 3; constexpr int COUNT_BKG_STARVED = 0; constexpr int COUNT_BKG_STARVED_NEIGHBOUR = 1; constexpr int COUNT_PROFILE_FALLBACK = 2; -constexpr int COUNT_SLOTS = 3; +constexpr int COUNT_OVERLAPPED = 3; +constexpr int COUNT_OVERLAPPED_MAJOR = 4; +constexpr int COUNT_SLOTS = 5; // --- Mark the r2 signal region of every predicted reflection (race-free: one launch writes one // value). The region is the INNER stencil ellipse in the neighbour's own frame; see the CPU engine. @@ -277,6 +279,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, atomicAdd(counts + COUNT_BKG_STARVED, 1ull); if ((s_full || p.partial_ok) && s_nclean <= 5 && s_nring > 5) atomicAdd(counts + COUNT_BKG_STARVED_NEIGHBOUR, 1ull); + if (p.overlap && s_nown < s_ndisk) { + atomicAdd(counts + COUNT_OVERLAPPED, 1ull); + if (s_nown < p.minpk * s_ndisk) atomicAdd(counts + COUNT_OVERLAPPED_MAJOR, 1ull); + } s_bkg = s_accept ? ((double) (long long) s_bkgsum / (double) s_nbkg) : 0.0; s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0)); // The pixel is an integer, so comparing it against the floor of the threshold accepts @@ -1083,11 +1089,13 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu BraggIntegrationCounts BraggIntegrationEngineGPU::Counts() const { BraggIntegrationCounts out = counts; // `predicted` is counted on the host, as on the CPU engine - unsigned long long host[COUNT_SLOTS] = {0, 0, 0}; + unsigned long long host[COUNT_SLOTS] = {}; cuda_err(cudaMemcpyAsync(host, d_counts, sizeof(host), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaStreamSynchronize(*stream)); out.bkg_starved = host[COUNT_BKG_STARVED]; out.bkg_starved_by_neighbour = host[COUNT_BKG_STARVED_NEIGHBOUR]; out.profile_fallback = host[COUNT_PROFILE_FALLBACK]; + out.overlapped = host[COUNT_OVERLAPPED]; + out.overlapped_major = host[COUNT_OVERLAPPED_MAJOR]; return out; } diff --git a/rugnux/ReportDocument.h b/rugnux/ReportDocument.h index 5f1361317..6a171058c 100644 --- a/rugnux/ReportDocument.h +++ b/rugnux/ReportDocument.h @@ -52,6 +52,8 @@ namespace PathologyCode { constexpr const char *POWDER_RINGS = "POWDER_RINGS"; // non-ice powder rings hold many spots constexpr const char *ICE_RINGS = "ICE_RINGS"; // hexagonal-ice rings hold many spots constexpr const char *BEAM_CENTRE_OFF = "BEAM_CENTRE_OFF"; // the file's beam centre is wrong for this geometry + constexpr const char *OVERLOADS = "OVERLOADS"; // saturated pixels took reflections + constexpr const char *OVERLAPS = "OVERLAPS"; // neighbouring spots share signal pixels } // A closed type plus an open description. The type is what machinery switches on; the text is what a diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 48b8847c7..eca796bb3 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -100,6 +100,17 @@ namespace { // fired on were below it (0.4-0.6% of the sweep). constexpr double SWEEP_GAPS_WARNING_FRACTION = 0.01; + // Overloads are warned about from this share of the integrated reflections, or of those coarser than + // OVERLOAD_LOW_RES_A, where a handful already takes the strongest low-resolution terms. Low on purpose: + // the strong-diffraction in-house sweep reads 0.017% and 0.42%, and its merge drops 0.12% of its + // rocking events; an ordinary sweep with a few hundred overloaded partials reads 0.001% and 0.012%. + constexpr double OVERLOAD_WARNING_PCT = 0.05; + constexpr double OVERLOAD_LOW_RES_WARNING_PCT = 0.25; + + // Overlapping spots are warned about from this share of the predicted reflections left with less + // than --overlap-minpk of their profile their own - the ones XDS's MINPK would reject. + constexpr double OVERLAP_WARNING_PCT = 1.0; + // Whether <|L|> - on the merge reported, or on the one the space-group search was given - reads // outside its physical range. Then no twin verdict is read from it, in either direction. bool TwinStatisticsAbnormal(const ProcessResult &r) { @@ -411,7 +422,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, const bool merged = result.has_merge_statistics; doc.warnings = result.warnings; // SUMMARY rows filled by the sections that measure them; empty leaves the row out. - std::string beam_centre_row; + std::string beam_centre_row, overload_row, overlap_row; // ------------------------------------------------------------------ header block { @@ -1517,6 +1528,53 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, Blank()); } + // ---- overloads and overlapping spots: two ways the detector or the pattern took reflections + // the crystal gave, each with what to change at the next collection + if (result.integrated_reflections > 0) { + const double pct = 100.0 * result.overloaded_reflections / result.integrated_reflections; + const double pct_low = result.integrated_reflections_low_res > 0 + ? 100.0 * result.overloaded_reflections_low_res / result.integrated_reflections_low_res : 0.0; + Add(s, KeyInt("OVERLOADED_REFLECTIONS", result.overloaded_reflections)); + Add(s, KeyReal("OVERLOADED_REFLECTIONS_PCT", pct, "{:.3f}")); + Add(s, KeyReal("OVERLOADED_REFLECTIONS_LOW_RES_PCT", pct_low, "{:.3f}")); + overload_row = fmt::format("{:.3f}% of the integrated reflections, {:.3f}% of those coarser than {:.0f} A", + pct, pct_low, OVERLOAD_LOW_RES_A); + if (pct >= OVERLOAD_WARNING_PCT || pct_low >= OVERLOAD_LOW_RES_WARNING_PCT) + Warn(doc, PathologyCode::OVERLOADS, fmt::format( + "Overloads: {} integrated reflections ({:.3f}%; {:.2f}% of those coarser than {:.0f} A) held " + "a saturated pixel and are not measurements - the merge leaves them out, and they are the " + "strongest reflections the crystal gave. Attenuate the beam or shorten the exposure (or " + "collect a second, attenuated low-resolution sweep)", + result.overloaded_reflections, pct, pct_low, OVERLOAD_LOW_RES_A)); + } + if (result.predicted_reflections > 0 + && experiment.GetBraggIntegrationSettings().GetOverlap() != OverlapMode::Off) { + const double pct = 100.0 * static_cast(result.overlapped_reflections) + / static_cast(result.predicted_reflections); + const double pct_major = 100.0 * static_cast(result.overlapped_reflections_major) + / static_cast(result.predicted_reflections); + const double minpk = 100.0 * experiment.GetBraggIntegrationSettings().GetOverlapMinPeak(); + Add(s, KeyReal("OVERLAPPED_REFLECTIONS_PCT", pct, "{:.2f}")); + Add(s, KeyReal("OVERLAPPED_REFLECTIONS_MAJOR_PCT", pct_major, "{:.2f}")); + overlap_row = fmt::format("{:.2f}% of the predicted reflections share signal pixels, {:.2f}% keep under " + "{:.0f}% of their profile", pct, pct_major, minpk); + if (pct_major >= OVERLAP_WARNING_PCT) + Warn(doc, PathologyCode::OVERLAPS, fmt::format( + "Overlapping spots: {:.2f}% of the predicted reflections keep less than {:.0f}% of their " + "profile as their own ({:.2f}% share signal pixels with a neighbour) - the shared pixels " + "are left out of the fit, which costs precision. Move the detector back, slice the " + "rotation finer or use a smaller beam", pct_major, minpk, pct)); + } + if (result.integrated_reflections > 0 || result.predicted_reflections > 0) { + Add(s, Blank()); + Add(s, Prose(" OVERLOADED_REFLECTIONS counts the integrated reflections (every frame's partial) with a\n" + " saturated pixel in the signal disk, OVERLOADED_REFLECTIONS_LOW_RES_PCT their share among\n" + " those coarser than 4 A, where they concentrate. OVERLAPPED_REFLECTIONS_PCT is the share\n" + " of the predicted reflections whose signal disk a nearer neighbour takes pixels from,\n" + " OVERLAPPED_REFLECTIONS_MAJOR_PCT the share left with less than --overlap-minpk of it.")); + Add(s, Blank()); + } + // ---- translational pseudo-symmetry // Printed immediately before the twinning block, and in the same section, because the two // interact: a pseudo-translation that is not a half-integer vector biases <|L|> upwards and @@ -2444,6 +2502,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, "with -z or --model", result.indexing_ambiguity_operators)); if (!beam_centre_row.empty()) checked_row("Beam centre", {PathologyCode::BEAM_CENTRE_OFF}, beam_centre_row); + if (!overload_row.empty()) + checked_row("Overloads", {PathologyCode::OVERLOADS}, overload_row); + if (!overlap_row.empty()) + checked_row("Overlaps", {PathologyCode::OVERLAPS}, overlap_row); if (!result.powder.rings_q_recipA.empty() || result.consensus_cell) { const auto &pw = result.powder; checked_row("Powder", {PathologyCode::POWDER_RINGS}, pw.rings_q_recipA.empty() ? std::string("no rings detected") diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2c8068da2..fe821b0ce 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -537,6 +537,17 @@ struct ProcessResult { std::optional background_beam_centre; float beam_centre_need_pxl = NAN; + // Overloads: integrated reflections (partial observations, every frame) a saturated pixel of the + // signal disk made unmeasurable (Reflection::overloaded), over the whole range and over the + // reflections coarser than OVERLOAD_LOW_RES_A, where they concentrate. Counted on the integration + // outcome before the merge; zero where nothing was integrated. + int64_t integrated_reflections = 0, overloaded_reflections = 0; + int64_t integrated_reflections_low_res = 0, overloaded_reflections_low_res = 0; + // Overlapping spots, as the integrator's overlap treatment measured them (BraggIntegrationCounts): + // predicted reflections whose signal disk a nearer neighbour takes pixels from, and of those the + // ones left with less than --overlap-minpk of it. Zero with --overlap off. + uint64_t predicted_reflections = 0, overlapped_reflections = 0, overlapped_reflections_major = 0; + // Which pass this result is (1-based) and how many ran, plus one line saying which pass became the // canonical output and on what evidence. pass_count == 1 outside the rotation two-pass. int pass_number = 1; @@ -568,6 +579,10 @@ struct ProcessResult { std::optional spindle_lost_unique_fraction; }; +// Overloads are counted separately over the reflections coarser than this (A) - the Wilson-plot +// convention's 4 A (CalcGlobalWilsonB), where the strong low-resolution reflections that saturate sit. +constexpr float OVERLOAD_LOW_RES_A = 4.0f; + // The rotation two-pass quality guard: whether the refined-geometry pass (pass2) is worse than the // header-geometry pass (pass1), judged on their search merges. Returns why, in the words the run // reports it in, or an empty string when the refined pass stands. diff --git a/rugnux/RugnuxImagePass.cpp b/rugnux/RugnuxImagePass.cpp index cf698fec3..0f7e3e64a 100644 --- a/rugnux/RugnuxImagePass.cpp +++ b/rugnux/RugnuxImagePass.cpp @@ -127,7 +127,7 @@ void Rugnux::ImagePass(PipelineLocals &p) { // finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is // done once per worker at the end rather than per image. std::atomic bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0, - bragg_fallback = 0; + bragg_fallback = 0, bragg_overlapped = 0, bragg_overlapped_major = 0; // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of @@ -300,6 +300,8 @@ void Rugnux::ImagePass(PipelineLocals &p) { bragg_starved += bragg.bkg_starved; bragg_starved_nb += bragg.bkg_starved_by_neighbour; bragg_fallback += bragg.profile_fallback; + bragg_overlapped += bragg.overlapped; + bragg_overlapped_major += bragg.overlapped_major; }; if (observer) @@ -475,6 +477,16 @@ void Rugnux::ImagePass(PipelineLocals &p) { experiment_.GetBraggIntegrationSettings().GetR3(), bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb), frac(bragg_fallback)); + // Overlapping spots: signal disks a nearer neighbour's prediction takes pixels from, as the + // overlap treatment decided ownership. For the report (OVERLAPPED_REFLECTIONS). + result.predicted_reflections = bragg_predicted; + result.overlapped_reflections = bragg_overlapped; + result.overlapped_reflections_major = bragg_overlapped_major; + if (experiment_.GetBraggIntegrationSettings().GetOverlap() != OverlapMode::Off) + logger.Info("Overlapping spots: {:.2f}% of the predicted reflections share signal pixels with a " + "nearer neighbour, {:.2f}% keep less than {:.0f}% of their profile as their own", + frac(bragg_overlapped), frac(bragg_overlapped_major), + 100.0 * experiment_.GetBraggIntegrationSettings().GetOverlapMinPeak()); } // Every image failing is a total failure, not a run that produced nothing: it used to be reported diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 823d15df6..6c205a0d0 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -298,6 +298,32 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { phase("Scaling and merging"); Step(geometry_prepass ? "Scaling and post-refinement" : "Scaling and merging"); + // Overloads, for the report: the integrated reflections a saturated signal pixel made unmeasurable. + // Counted before anything below drops or flags one. The merge leaves them out (a rotation merge + // drops the whole rocking event) - they are the strongest reflections, mostly at low resolution, + // and what a model's low-resolution terms rest on. + { + for (const auto &io : indexer->GetIntegrationOutcome()) + for (const auto &r : io.reflections) { + const bool low_res = r.d >= OVERLOAD_LOW_RES_A; + ++result.integrated_reflections; + result.integrated_reflections_low_res += low_res; + if (r.overloaded) { + ++result.overloaded_reflections; + result.overloaded_reflections_low_res += low_res; + } + } + if (result.integrated_reflections > 0) + logger.Info("Overloads: {} of {} integrated reflections ({:.3f}%) held a saturated pixel, {:.3f}% of " + "the {} coarser than {:.0f} A", result.overloaded_reflections, + result.integrated_reflections, + 100.0 * result.overloaded_reflections / result.integrated_reflections, + result.integrated_reflections_low_res > 0 + ? 100.0 * result.overloaded_reflections_low_res / result.integrated_reflections_low_res + : 0.0, + result.integrated_reflections_low_res, OVERLOAD_LOW_RES_A); + } + // Ice-ring handling (--detect-ice-rings): flag reflections sitting on a hexagonal-ice powder // ring. Their integrated intensity is contaminated by the strong, variable ice background, so // they drag the per-image scale fit. Flag them here so scaling (the per-image G fit and the diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 010930e83..cba64eef1 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -1236,3 +1236,30 @@ TEST_CASE("ResultReport_BeamCentreOff", "[Diagnostics]") { CHECK(text.find("\nPATHOLOGY_FLAGS= BEAM_CENTRE_OFF\n") != std::string::npos); CHECK(text.find("the background's word alone") != std::string::npos); } + +TEST_CASE("ResultReport_OverloadsAndOverlaps", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + ProcessResult result; + result.integrated_reflections = 100000; + result.integrated_reflections_low_res = 5000; + result.overloaded_reflections = 5; + result.overloaded_reflections_low_res = 5; + result.predicted_reflections = 100000; + result.overlapped_reflections = 5000; + result.overlapped_reflections_major = 500; + + auto text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("\nOVERLOADED_REFLECTIONS= 5\n") != std::string::npos); + CHECK(text.find("\nOVERLOADED_REFLECTIONS_LOW_RES_PCT= 0.100\n") != std::string::npos); + CHECK(text.find("\nOVERLAPPED_REFLECTIONS_MAJOR_PCT= 0.50\n") != std::string::npos); + CHECK(text.find("\nPATHOLOGY_FLAGS= NONE\n") != std::string::npos); + + // 1.2% of the low-resolution reflections overloaded, 1.5% left under MINPK: both warn. + result.overloaded_reflections = 60; + result.overloaded_reflections_low_res = 60; + result.overlapped_reflections_major = 1500; + text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("\nPATHOLOGY_FLAGS= OVERLOADS OVERLAPS\n") != std::string::npos); + CHECK(text.find("Attenuate the beam") != std::string::npos); + CHECK(text.find("Move the detector back") != std::string::npos); +} From 244f3c0f9d594dae5c6f83f104481bbf31bb1304 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:15:05 +0200 Subject: [PATCH 3/6] rugnux: write the refined geometry as _XDS.INP A rotation run now writes the geometry it integrated at as an XDS input beside the report: ORGX/ORGY (the PONI, 1-based), DETECTOR_DISTANCE, the detector axes carrying the tilt, wavelength, the rotation axis the run integrated at (ProcessResult::used_rotation_axis, after any sign rescue), oscillation (with an adopted rotation scale), starting angle/frame, data range, NX/NY/QX/QY, OVERLOAD, sensor thickness, space group and cell (left commented where the group is not in its number's reference setting), FRIEDEL'S_LAW per -A and the merged resolution range. HDF5 inputs get the master file as an x_??????.h5 template and a LIB= line to point at the reader plugin; a CBF series gets its ???? template and its own frame numbers. The report names the file (XDS_INPUT=). Conventions per Kabsch (2010) Acta Cryst. D66, 125-132. The rotation axis goes out negated: rugnux turns the crystal by -phi about its axis (IndexAndRefine), XDS by +phi about ROTATION_AXIS. Checked with XDS IDXREF on two exports (one CBF, one HDF5 via durin): the negated axis indexes 99.5% / 14% of the spots, rugnux's own sign 37% / 0.04%. On the CBF set XDS then ran DEFPIX/INTEGRATE/CORRECT from the file unchanged (apart from the job list): 99.8% complete, ISa 23.3, direct beam within 0.15 px and cell within 0.07 A of rugnux's. Output only; no result changes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX.md | 3 +- docs/RUGNUX_REPORT.md | 3 + docs/RUGNUX_TUTORIAL.md | 14 ++++ rugnux/CMakeLists.txt | 2 + rugnux/ResultReport.cpp | 3 + rugnux/Rugnux.cpp | 2 + rugnux/Rugnux.h | 5 ++ rugnux/XDSInput.cpp | 167 +++++++++++++++++++++++++++++++++++++ rugnux/XDSInput.h | 29 +++++++ rugnux/rugnux_cli.cpp | 5 ++ tests/ResultReportTest.cpp | 34 ++++++++ 11 files changed, 266 insertions(+), 1 deletion(-) create mode 100644 rugnux/XDSInput.cpp create mode 100644 rugnux/XDSInput.h diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 8b3198a35..2b8c7f7e0 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -52,7 +52,7 @@ that, or to lift the per-image loop's default ceiling of 16 workers per GPU. Nothing more is needed to pick the workflow: a dataset carrying a **goniometer axis** is processed as a rotation sweep, one without as **independent stills**, and scaling and merging run by default -in both. A rotation run that merges — the default — leaves eight files next to each other: +in both. A rotation run that merges — the default — leaves nine files next to each other: ``` myrun.mtz merged intensities + French-Wilson amplitudes, for CCP4 / phenix @@ -63,6 +63,7 @@ myrun_unmerged.mtz every observation before scaling, for pointless / aimless / myrun_report.txt what the run determined: cell, space group, statistics, warnings myrun_plot.txt one row per image, for plotting how the crystal behaved over the sweep myrun_detector.jpg the detector with the pixel mask and the beam-stop shadow drawn on it +myrun_XDS.INP the geometry the run refined, as an XDS input - a starting point for XDS myrun_log.txt the full timestamped log of the run (debug level) ``` diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index d992bf483..13c2dad66 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -150,6 +150,9 @@ BEAM_CENTRE_SHIFT= 7.15 BEAM_CENTRE_TOLERANCE= 3.49 ``` +**`XDS_INPUT=`** in section 1 names `_XDS.INP`, the same geometry written as an XDS input +(rotation runs; see the tutorial's list of output files). + **Which pass.** A rotation run integrates twice — once at the geometry in the input file (the first half of the sweep but at least 120°, `--prepass-fraction`), then the whole sweep again at the post-refined geometry — and can integrate a third time if a guard rejects the second pass. There is **one** report, for the pass that became the canonical output, and `PASS=` / diff --git a/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index 213fe8db5..8f68b2acd 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -331,6 +331,20 @@ Absorption is corrected only by the empirical surface a rotation run fits and ke is missing. A shadow drawn where the background is unbroken, or a dark patch the run left uncoloured, is visible at a glance; read it against the `Beam stop shadow: N pixels ...` line in the log. +- `_XDS.INP` — **the geometry the run ended on, as an XDS input**: on a rotation run the + post-refined beam centre (`ORGX`/`ORGY`, the foot of the perpendicular, numbered from 1 as XDS + numbers pixels), distance, detector axes (`DIRECTION_OF_DETECTOR_X-AXIS`/`Y-AXIS`, which carry the + tilt), wavelength, rotation axis (the one the run integrated at, after any sign rescue, negated: XDS turns + the crystal the other way about the axis it is given), oscillation + (scaled where the run adopted a goniometer rotation scale), starting angle and data range, plus the + space group and cell (left commented out where the group is not in the reference setting of its + number, which is the only setting a number names), `FRIEDEL'S_LAW` per `-A` and the merged + resolution range. For HDF5 the template is the master file and the `LIB=` line for the reader + plugin is left for the user to point at durin or neggia; the file's own pixel mask reaches XDS + through the plugin, the beam-stop shadow rugnux detected does not (XDS finds its own in `DEFPIX`). + The conventions are XDS's own (W. Kabsch, "XDS" (2010), Acta Cryst. D66, 125-132 + [doi:10.1107/S0907444909047337](https://doi.org/10.1107/S0907444909047337)). It is a starting + point - XDS refines the geometry again - and the report names it as **`XDS_INPUT=`** in section 1. Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the console. By default the written resolution is trimmed automatically where CC1/2 falls off diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 6b437d83b..5db9f5764 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -23,6 +23,8 @@ ADD_LIBRARY(JFJochRugnux STATIC DiagnosticOutput.h BrokerFeed.cpp BrokerFeed.h + XDSInput.cpp + XDSInput.h ) TARGET_LINK_LIBRARIES(JFJochRugnux JFJochReader JFJochImageAnalysis JFJochWriter JFJochModelWriter JFJochPreview gemmi fftw3f diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index eca796bb3..8cfc7353d 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -548,6 +548,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (result.spot_resolution_estimate_A.has_value()) Add(s, KeyReal("SPOT_RESOLUTION_ESTIMATE", *result.spot_resolution_estimate_A, "{:.2f}")); + // The same geometry once more, as an XDS input (rugnux/XDSInput.h). + if (!result.xds_input_path.empty()) + Add(s, KeyText("XDS_INPUT", result.xds_input_path)); if (result.pass_count > 1) { Add(s, KeyText("PASS", fmt::format("{} of {}", result.pass_number, result.pass_count))); Add(s, KeyText("PASS_DECISION", result.pass_decision, true)); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index bfa4dd921..33dd36c9f 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1075,6 +1075,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.used_direct_beam_x_pxl = direct.first; result.used_direct_beam_y_pxl = direct.second; } + if (const auto g = experiment_.GetGoniometer()) + result.used_rotation_axis = g->GetAxis(); if (census_future.valid()) { TimingMark mark(logger, "leftover-lattice census wait"); result.leftover_lattices = census_future.get(); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index fe821b0ce..cbe0bafbd 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -301,6 +301,9 @@ struct ProcessResult { // rotations, which is a geometry no pass ever ran at. std::array used_detector_tilt_deg{}; double used_direct_beam_x_pxl = 0.0, used_direct_beam_y_pxl = 0.0; + // The rotation axis the same pass integrated at - the sign rescue can adopt the opposite of the + // file's (see LogXDSGeometry). Empty on stills. + std::optional used_rotation_axis; std::optional consensus_cell; // The centring `consensus_cell` is written in; empty when nothing named one (the stills path // averages per-image cells and has no lattice search behind them). A conventional cell says @@ -326,6 +329,8 @@ struct ProcessResult { std::optional lattice_angle_departure_deg; MeanProcessingTime mean_processing_time{}; std::optional written_master_path; + // _XDS.INP, where the CLI wrote the run's geometry as an XDS input (rugnux/XDSInput.h). + std::string xds_input_path; std::string merge_statistics_text; // populated when scaling/merging ran // Structured merge statistics (per-shell + overall), populated when scaling/merging ran. ISa is diff --git a/rugnux/XDSInput.cpp b/rugnux/XDSInput.cpp new file mode 100644 index 000000000..e8361c8f8 --- /dev/null +++ b/rugnux/XDSInput.cpp @@ -0,0 +1,167 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include +#include +#include + +#include + +#include "../common/GitInfo.h" +#include "../common/JFJochMath.h" +#include "XDSInput.h" + +std::string RenderXDSInput(const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + int64_t first_image, int64_t end_image) { + // The geometry the result was integrated at, on a copy of the file's experiment. + DiffractionExperiment x = experiment; + if (std::isfinite(result.used_beam_x_pxl) && std::isfinite(result.used_distance_mm)) { + constexpr double RAD = PI / 180.0; + const auto &t = result.used_detector_tilt_deg; + x.BeamX_pxl(result.used_beam_x_pxl).BeamY_pxl(result.used_beam_y_pxl) + .DetectorDistance_mm(result.used_distance_mm) + .PoniRot1_rad(static_cast(t[0] * RAD)).PoniRot2_rad(static_cast(t[1] * RAD)) + .PoniRot3_rad(static_cast(t[2] * RAD)); + } + const auto geom = x.GetDiffractionGeometry(); + const float beam_x = geom.GetBeamX_pxl(); + const float beam_y = geom.GetBeamY_pxl(); + const float pixel_mm = geom.GetPixelSize_mm(); + + // Following the XDS conventions, Kabsch (2010) Acta Cryst. D66, 125-132: XDS places a pixel at + // x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis) + // in a laboratory frame that is ours (x along a column step, y along a row step, z along the + // beam), with pixels numbered from 1 and ORGX/ORGY the foot of the perpendicular from the crystal - + // our PONI (see LogXDSGeometry in Rugnux.cpp and docs/DETECTOR_GEOMETRY.md). So the tilt goes out + // as the two detector axes, taken as differences of LabCoord. + const Coord origin = geom.LabCoord(beam_x, beam_y); + const Coord x_axis = ((geom.LabCoord(beam_x + 1.0f, beam_y) - origin) / pixel_mm).Normalize(); + const Coord y_axis = ((geom.LabCoord(beam_x, beam_y + 1.0f) - origin) / pixel_mm).Normalize(); + + const bool hdf5 = input_file.ends_with(".h5") || input_file.ends_with(".nxs"); + const bool cbf = input_file.ends_with(".cbf") || input_file.ends_with(".cbf.gz"); + std::string input_path = input_file; + try { + input_path = std::filesystem::absolute(input_file).string(); + } catch (const std::exception &) { + } + + // A CBF sweep is one file per frame: XDS takes a template with ? for the frame number - the last run + // of digits in the name - and numbers the frames by it, so the data range starts at the given file's. + std::string template_path = input_path; + int64_t first_frame = 1; + // The HDF5 plugins take the master file's name with "master" as ??????, as XDS has it for + // an EIGER data set. + if (hdf5 && input_path.ends_with("_master.h5")) + template_path.replace(template_path.size() - 9, 6, "??????"); + if (cbf) { + const size_t name = template_path.find_last_of('/') + 1; + const size_t end = template_path.find_last_of("0123456789"); + if (end != std::string::npos && end >= name) { + size_t begin = end; + while (begin > name && std::isdigit(static_cast(template_path[begin - 1]))) + --begin; + first_frame = std::stoll(template_path.substr(begin, end - begin + 1)); + template_path.replace(begin, end - begin + 1, std::string(end - begin + 1, '?')); + } + } + + std::string s; + s += fmt::format("! XDS.INP written by rugnux {} from {}\n", jfjoch_version(), input_path); + s += "! The detector geometry is the one rugnux integrated at - on a rotation run the post-refined\n" + "! one - not the file's; XDS refines it again from here. Check before use: this file is a\n" + "! starting point, not a validated XDS setup.\n"; + s += "JOB= XYCORR INIT COLSPOT IDXREF DEFPIX INTEGRATE CORRECT\n\n"; + + if (hdf5) + s += "! HDF5 is read through a plugin: uncomment and point LIB= at durin-plugin.so (or the\n" + "! Dectris neggia plugin), which opens the _master.h5 the template names with ??????.\n" + "! The plugin applies the file's own pixel mask.\n" + "!LIB= /path/to/durin-plugin.so\n"; + else if (!cbf) + s += "! rugnux read this input itself; XDS may need the frames converted, and the template\n" + "! below written with ?????? in place of the frame number.\n"; + s += fmt::format("NAME_TEMPLATE_OF_DATA_FRAMES= {}\n", template_path); + s += fmt::format("DATA_RANGE= {} {}\n", first_frame + first_image, first_frame + end_image - 1); + s += fmt::format("SPOT_RANGE= {} {}\n\n", first_frame + first_image, first_frame + end_image - 1); + + s += fmt::format("DETECTOR= {} MINIMUM_VALID_PIXEL_VALUE= 0 OVERLOAD= {}\n", + cbf ? "PILATUS" : "EIGER", SaturationValueFromLimit(experiment.GetSaturationLimit())); + s += fmt::format("NX= {} NY= {} QX= {:.6f} QY= {:.6f}\n", + experiment.GetXPixelsNumConv(), experiment.GetYPixelsNumConv(), pixel_mm, pixel_mm); + const float thickness_um = experiment.GetDetectorSetup().GetSensorThickness_um(); + if (thickness_um > 0.0f) + s += fmt::format("SENSOR_THICKNESS= {:.3f}\n", thickness_um / 1000.0f); + s += fmt::format("ORGX= {:.2f} ORGY= {:.2f}\n", beam_x + 1.0f, beam_y + 1.0f); + s += fmt::format("DETECTOR_DISTANCE= {:.3f}\n", geom.GetDetectorDistance_mm()); + s += fmt::format("DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}\n", x_axis.x, x_axis.y, x_axis.z); + s += fmt::format("DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}\n", y_axis.x, y_axis.y, y_axis.z); + s += "INCIDENT_BEAM_DIRECTION= 0 0 1\n"; + s += fmt::format("X-RAY_WAVELENGTH= {:.6f}\n", geom.GetWavelength_A()); + s += "! Polarization is not in the file; a typical synchrotron:\n" + "!FRACTION_OF_POLARIZATION= 0.99 POLARIZATION_PLANE_NORMAL= 0 1 0\n\n"; + + if (const auto gonio = experiment.GetGoniometer()) { + // Negated: rugnux turns the crystal by -phi about its axis (IndexAndRefine takes the lattice + // to an image with GetTransformationAngle(-angle)), XDS by a right-handed +phi about + // ROTATION_AXIS. Checked with IDXREF on the export: 99.5% of the spots index with the + // negated axis, 37% (the chance floor of a mis-signed spindle) with rugnux's own. + const Coord axis = -result.used_rotation_axis.value_or(gonio->GetAxis()); + // The oscillation the run integrated at: the file's, scaled where the run adopted a goniometer + // rotation scale (as OSCILLATION_RANGE in the report). + const double scale = result.post_refine && result.post_refine->rotation_scale_suspect + ? result.post_refine->rotation_scale : 1.0; + // + 0.0f: no "-0.000000" from a negated zero. + s += fmt::format("ROTATION_AXIS= {:.6f} {:.6f} {:.6f}\n", axis.x + 0.0f, axis.y + 0.0f, axis.z + 0.0f); + s += fmt::format("OSCILLATION_RANGE= {:.6f}\n", gonio->GetIncrement_deg() * scale); + s += fmt::format("STARTING_FRAME= {}\n", first_frame); + s += fmt::format("STARTING_ANGLE= {:.4f}\n\n", gonio->GetStart_deg()); + } + + if (result.space_group && result.consensus_cell) { + const auto &c = *result.consensus_cell; + // A number names only the reference setting; a cell in another one would be read wrongly. + const char *prefix = result.space_group->is_reference_setting() ? "" : "!"; + if (*prefix) + s += fmt::format("! The run's space group {} is not in the reference setting of its number;\n" + "! XDS would read the cell below in the wrong setting, so both are left out.\n", + result.space_group->xhm()); + s += fmt::format("{}SPACE_GROUP_NUMBER= {}\n", prefix, result.space_group->number); + s += fmt::format("{}UNIT_CELL_CONSTANTS= {:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}\n", prefix, + c.a, c.b, c.c, c.alpha, c.beta, c.gamma); + } + s += fmt::format("FRIEDEL'S_LAW= {}\n", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE"); + if (result.has_merge_statistics) + s += fmt::format("INCLUDE_RESOLUTION_RANGE= {:.3f} {:.3f}\n", + result.merge_statistics.overall.d_max, result.merge_statistics.overall.d_min); + s += "\n! rugnux also masked the beam-stop shadow it detected (_detector.jpg); XDS finds its\n" + "! own in DEFPIX, or add UNTRUSTED_RECTANGLE= / UNTRUSTED_QUADRILATERAL= lines here.\n"; + return s; +} + +std::string WriteXDSInput(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + int64_t first_image, int64_t end_image, + Logger &logger) { + if (output_prefix.empty() || !experiment.GetGoniometer()) + return {}; + const std::string filename = output_prefix + "_XDS.INP"; + try { + std::ofstream f(filename); + if (!f) + throw std::runtime_error("cannot open for writing"); + f << RenderXDSInput(input_file, experiment, result, first_image, end_image); + if (!f) + throw std::runtime_error("write failed"); + } catch (const std::exception &e) { + logger.Warning("Could not write the XDS input {}: {}", filename, e.what()); + return {}; + } + return filename; +} diff --git a/rugnux/XDSInput.h b/rugnux/XDSInput.h new file mode 100644 index 000000000..393cfd7ef --- /dev/null +++ b/rugnux/XDSInput.h @@ -0,0 +1,29 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../common/DiffractionExperiment.h" +#include "../common/Logger.h" +#include "Rugnux.h" + +// _XDS.INP: the geometry the run integrated at - on a rotation run the post-refined one - as an +// XDS input, so the same data can be processed by XDS from where rugnux ended rather than from the +// header. `experiment` is the one the run was given (the file's), `result` carries what the run +// refined; the images are the file's 0-based [first_image, end_image). +std::string RenderXDSInput(const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + int64_t first_image, int64_t end_image); + +// Writes _XDS.INP and returns its name; returns an empty string, writing nothing, on a +// run with no output prefix or no rotation, and logs and swallows a failure to write. +std::string WriteXDSInput(const std::string &output_prefix, + const std::string &input_file, + const DiffractionExperiment &experiment, + const ProcessResult &result, + int64_t first_image, int64_t end_image, + Logger &logger); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9e74ef787..4e9d30bc5 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -29,6 +29,7 @@ #include "../reader/JFJochSMVReader.h" #include "../reader/JFJochStartMessageReader.h" #include "BrokerFeed.h" +#include "XDSInput.h" #include "../writer/HDF5NXmx.h" #include #include @@ -3063,6 +3064,10 @@ static int RunRugnux(int argc, char **argv) { provenance.wall_time_s = std::chrono::duration( std::chrono::steady_clock::now() - invocation_start).count(); TimingMark report_mark(logger, "results report"); + // The geometry the run ended on, as an XDS input beside the report, which names it. + if (rotation_indexing) + result.xds_input_path = WriteXDSInput(output_prefix, input_file, experiment, result, start_image, + end_image, logger); const ReportDocument doc = WriteResultReport(output_prefix, input_file, experiment, result, logger, provenance); report_mark.Stop(); diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index cba64eef1..64d68342b 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -7,6 +7,7 @@ #include "../image_analysis/scale_merge/Merge.h" #include "../reader/JFJochHDF5Reader.h" #include "../rugnux/ResultReport.h" +#include "../rugnux/XDSInput.h" #include "../common/CUDAWrapper.h" #include "../writer/FileWriter.h" @@ -1263,3 +1264,36 @@ TEST_CASE("ResultReport_OverloadsAndOverlaps", "[Diagnostics]") { CHECK(text.find("Attenuate the beam") != std::string::npos); CHECK(text.find("Move the detector back") != std::string::npos); } + +TEST_CASE("XDSInput_Render", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + x.BeamX_pxl(500.0f).BeamY_pxl(400.0f).DetectorDistance_mm(100.0f); + x.Goniometer(GoniometerAxis("omega", 10.0f, 0.1f, Coord(1, 0, 0), {})); + ProcessResult result; + // The refined geometry, not the file's, goes out - 1-based, untilted here. + result.used_beam_x_pxl = 510.0f; + result.used_beam_y_pxl = 395.5f; + result.used_distance_mm = 120.0f; + result.used_rotation_axis = Coord(-1, 0, 0); + + const auto text = RenderXDSInput("/data/x_master.h5", x, result, 0, 900); + CHECK(text.find("\nORGX= 511.00 ORGY= 396.50\n") != std::string::npos); + CHECK(text.find("\nDETECTOR_DISTANCE= 120.000\n") != std::string::npos); + CHECK(text.find("\nDIRECTION_OF_DETECTOR_X-AXIS= 1.000000 0.000000 0.000000\n") != std::string::npos); + // XDS turns the crystal the other way about the axis it is given. + CHECK(text.find("\nROTATION_AXIS= 1.000000 0.000000 0.000000\n") != std::string::npos); + CHECK(text.find("\nNAME_TEMPLATE_OF_DATA_FRAMES= /data/x_??????.h5\n") != std::string::npos); + CHECK(text.find("\nDATA_RANGE= 1 900\n") != std::string::npos); + CHECK(text.find("\nSTARTING_ANGLE= 10.0000\n") != std::string::npos); + CHECK(text.find("\nOSCILLATION_RANGE= 0.100000\n") != std::string::npos); +} + +TEST_CASE("XDSInput_CBFTemplate", "[Diagnostics]") { + DiffractionExperiment x(DetJF(1)); + x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(1, 0, 0), {})); + const auto text = RenderXDSInput("/data/run_2_0101.cbf", x, ProcessResult{}, 0, 50); + CHECK(text.find("\nNAME_TEMPLATE_OF_DATA_FRAMES= /data/run_2_????.cbf\n") != std::string::npos); + CHECK(text.find("\nDATA_RANGE= 101 150\n") != std::string::npos); + CHECK(text.find("\nSTARTING_FRAME= 101\n") != std::string::npos); + CHECK(text.find("\nDETECTOR= PILATUS") != std::string::npos); +} From 7624bbe12e2d3d37a271b534b0081ec5ec6a6da7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:16:37 +0200 Subject: [PATCH 4/6] battery: paper_plots.py draws the other figures when no run has REFMAC R-free The REFMAC check is opt-in, so a battery run without --model-check has no R-free pairs, and the axis range min()/max() over them crashed the script before any figure was drawn. The R-free figure is now its own gnuplot block, left out (with a line in summary.txt saying why) when there is nothing to draw; d_min and the time figures are drawn as before. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- tools/battery/paper_plots.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/tools/battery/paper_plots.py b/tools/battery/paper_plots.py index fbc2a3d6c..5229d5c3e 100644 --- a/tools/battery/paper_plots.py +++ b/tools/battery/paper_plots.py @@ -21,6 +21,7 @@ Left figure: R-free of the deposited model after the battery's REFMAC check on t (x) and on Rugnux's data (y). By default the refined protocol (model_check.py: the same model refined against each data set at its own resolution, R-free over the free reflections both data sets share); --rfree first-cycle draws the first-cycle rigid-body numbers at the common resolution instead. +A run without REFMAC values (the check is opt-in, --model-check) gets no R-free figure; the rest are drawn. Right figure: the deposited high-resolution limit (x) against Rugnux's own cut (y), log axes so the low-resolution entries show. Small molecules (no deposited model) and sets whose verdict is `fail` (a wrong lattice: nothing to compare) are left out. @@ -75,7 +76,10 @@ set style line 2 pt 7 ps 0.35 lc rgb "#0072b2" # counts: names left-aligned, numbers right-aligned in a column beside them set label 11 "Structures\nFacilities\nBeamlines" at graph 0.05, graph 0.93 left front set label 12 "%(n)d\n%(nf)d\n%(nb)d" at graph 0.31, graph 0.93 right front +""" +# Drawn only where the runs carry REFMAC R-free values; the battery's REFMAC check is opt-in. +GNUPLOT_RFREE = r""" # --- R-free set output "%(out)s/rfree.pdf" set xlabel "{/:Italic R}_{free}, depositor data" @@ -86,7 +90,9 @@ set format x "%%.1f"; set format y "%%.1f" set label 2 "Within ±%(band).2f: %(within)d/%(n)d" at graph 0.95, graph 0.07 right front plot x ls 1, "%(out)s/rfree.dat" using 1:2 ls 2 unset label 2 +""" +GNUPLOT_REST = r""" # --- high-resolution limit set output "%(out)s/dmin.pdf" set logscale xy @@ -255,7 +261,8 @@ def main(): else "set key off") params = dict(out=out, tkey=tkey, n=n, nf=nf, nb=nb, dn=dn, dnf=dnf, dnb=dnb, within=within, band=RFREE_BAND, deeper=deeper, shallower=sum(y > x for x, y, _ in dmin), - rlo=0.05 * math.floor((min(values) - 0.01) / 0.05), rhi=0.05 * math.ceil((max(values) + 0.01) / 0.05), + rlo=0.05 * math.floor((min(values, default=0.2) - 0.01) / 0.05), + rhi=0.05 * math.ceil((max(values, default=0.3) + 0.01) / 0.05), dlo=0.5, dhi=8.0, bin=TIME_BIN, thi=TIME_BIN * math.ceil(max(secs) / TIME_BIN), tyhi=10 * math.ceil((tallest + 1) / 10), tmed=tmed, tn=len(times), @@ -263,9 +270,9 @@ def main(): zyhi=5 * math.ceil((max(collections.Counter(math.floor(t / ZOOM_BIN) for t in secs if t < ZOOM_MAX).values()) + 1) / 5)) gp = os.path.join(out, "figures.gp") - open(gp, "w").write(GNUPLOT % params) + open(gp, "w").write((GNUPLOT + (GNUPLOT_RFREE if rfree else "") + GNUPLOT_REST) % params) subprocess.run(["gnuplot", gp], check=True) - for fig in ("rfree", "dmin", "time", "time_zoom"): # the PNG is the PDF rasterised at 600 d.p.i. + for fig in (["rfree"] if rfree else []) + ["dmin", "time", "time_zoom"]: # the PNG is the PDF rasterised at 600 d.p.i. subprocess.run(["pdftoppm", "-r", "600", "-png", "-singlefile", os.path.join(out, fig + ".pdf"), os.path.join(out, fig)], check=True) @@ -276,8 +283,10 @@ def main(): fh.write(f"run {path}\n label {m.get('label')} binary git {b.get('source_head')} " f"dirty {b.get('source_dirty')} sha256 {b.get('sha256')}\n flags {b.get('flags')} " f"runner git {m.get('runner_git')} dirty {m.get('runner_dirty')}\n") - fh.write(f"R-free: {n} structures, {nf} facilities, {nb} beamlines; within +-{RFREE_BAND}: {within}\n" - f"d_min: {dn} structures, {dnf} facilities, {dnb} beamlines; Rugnux reaches further on " + fh.write((f"R-free: {n} structures, {nf} facilities, {nb} beamlines; within +-{RFREE_BAND}: {within}\n" + if rfree else "R-free: no REFMAC values in these runs (run the battery with --model-check); " + "no rfree figure\n") + + f"d_min: {dn} structures, {dnf} facilities, {dnb} beamlines; Rugnux reaches further on " f"{deeper}\n") queue_wait = max(time_rows[s]["elapsed_s"] - t for s, _, t, _, _ in times) fh.write(f"time: {len(times)} sets of arm(s) {', '.join(a.time_arms)}; field rugnux_wall_s " From 7ed15441363246cf81958b49adba84a5f554da1f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 18:57:06 +0200 Subject: [PATCH 5/6] rugnux report: BEAM_CENTRE_OFF only beyond three times the tolerance At one times BeamCenterNeed_pxl the warning fired on 107 of the 250 battery sets - most headers are a few pixels off. It now fires beyond three times the tolerance; between one and three times the summary's Beam centre line states the offset and recommends correcting the beamline value, without a warning. The background-only branch (runs that refined no centre) uses the same bound. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX_REPORT.md | 8 +++++--- rugnux/ResultReport.cpp | 24 +++++++++++++++++------- tests/ResultReportTest.cpp | 12 +++++++++--- 3 files changed, 31 insertions(+), 13 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 13c2dad66..162a97c74 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -136,12 +136,14 @@ integrated at (`DIRECT_BEAM=`, the refined one on a rotation run) lies from the **`BEAM_CENTRE_TOLERANCE=`** is how far this crystal on this geometry lets the centre be wrong - the displacement at which the accumulated reflections of the longest cell axis smear to 70% of their height, the number the log's `Beam centre check` line states, here on the run's own cell. Where the -shift is larger, a **`BEAM_CENTRE_OFF`** warning gives the file's, the refined and the background +shift is more than three times that, a **`BEAM_CENTRE_OFF`** warning gives the file's, the refined and the background centre, the difference in pixels and millimetres, and the value to correct the beamline (or the header) to: the run itself integrated at the refined centre, but every other program given the same file starts from the wrong one. A run that refined no centre (stills) is warned on the background -alone, where its fit resolves the difference (three sigma) and the difference exceeds the tolerance, -worded as a check. The summary's `Beam centre` line states the outcome either way. +alone, where its fit resolves the difference (three sigma) and the difference exceeds three times the +tolerance, worded as a check. Between one and three times the tolerance the summary's `Beam centre` +line says so and recommends correcting the beamline value, without a warning - most headers are a few +pixels off, and at one times the tolerance the warning fired on 107 of 250 battery sets. ``` FILE_DIRECT_BEAM= 1564.31 1702.91 diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 8cfc7353d..4be52ae86 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -100,6 +100,11 @@ namespace { // fired on were below it (0.4-0.6% of the sweep). constexpr double SWEEP_GAPS_WARNING_FRACTION = 0.01; + // The beam centre is warned about beyond this many times what the geometry needs of it + // (BEAM_CENTRE_TOLERANCE); between one and this many times the summary states it without a warning. + // At one times the tolerance 107 of the 250 battery sets fired - most headers are a few pixels off. + constexpr double BEAM_CENTRE_WARNING_NEEDS = 3.0; + // Overloads are warned about from this share of the integrated reflections, or of those coarser than // OVERLOAD_LOW_RES_A, where a handful already takes the strongest low-resolution terms. Low on purpose: // the strong-diffraction in-house sweep reads 0.017% and 0.42%, and its merge drops 0.12% of its @@ -593,13 +598,14 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, const std::string bg_text = bg ? fmt::format("the scattered background at ({:.2f}, {:.2f}) +- {:.2f} px", bg->beam_x_pxl, bg->beam_y_pxl, bg->sigma_pxl) : std::string(); - if (std::isfinite(need) && moved > need) { + const double warn_at = BEAM_CENTRE_WARNING_NEEDS * need; + if (std::isfinite(need) && moved > warn_at) { beam_centre_row = fmt::format("the file's is OFF by {:.1f} px ({:.2f} mm) - correct the " "beamline value", moved, moved * pixel_mm); Warn(doc, PathologyCode::BEAM_CENTRE_OFF, fmt::format( "the beam centre in the input file is {:.1f} px ({:.2f} mm) from where this run found " - "the direct beam, more than the {:.1f} px this geometry needs it to be right to: the " - "file puts the direct beam at ({:.2f}, {:.2f}), the refined geometry at ({:.2f}, " + "the direct beam, more than three times the {:.1f} px this geometry needs it to be right " + "to: the file puts the direct beam at ({:.2f}, {:.2f}), the refined geometry at ({:.2f}, " "{:.2f}){}. This run used the refined centre, but any other program given this file " "starts from the wrong one - correct the beam centre at the beamline (or in the " "file's header) to ({:.2f}, {:.2f})", @@ -607,17 +613,21 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, result.used_direct_beam_x_pxl, result.used_direct_beam_y_pxl, bg ? ", " + bg_text : std::string(), result.used_direct_beam_x_pxl, result.used_direct_beam_y_pxl)); - } else if (std::isfinite(need) && !result.post_refine && bg_resolves && bg_moved > need) { + } else if (std::isfinite(need) && !result.post_refine && bg_resolves && bg_moved > warn_at) { beam_centre_row = fmt::format("the background puts it {:.1f} px ({:.2f} mm) from the file's - " "check the beamline value", bg_moved, bg_moved * pixel_mm); Warn(doc, PathologyCode::BEAM_CENTRE_OFF, fmt::format( "the beam centre in the input file may be off: the file puts the direct beam at " - "({:.2f}, {:.2f}), {} - {:.1f} px ({:.2f} mm) away, more than the {:.1f} px this " - "geometry needs the centre to be right to. Nothing on this run refined the centre, so " + "({:.2f}, {:.2f}), {} - {:.1f} px ({:.2f} mm) away, more than three times the {:.1f} px " + "this geometry needs the centre to be right to. Nothing on this run refined the centre, so " "this is the background's word alone: check the beam centre at the beamline and " "correct it if this is confirmed", file_direct.first, file_direct.second, bg_text, bg_moved, bg_moved * pixel_mm, need)); - } else if (std::isfinite(need)) + } else if (std::isfinite(need) && moved > need) + beam_centre_row = fmt::format("the file's is {:.1f} px ({:.2f} mm) from the refined one, more than " + "the {:.1f} px this geometry needs but under three times it - worth " + "correcting at the beamline", moved, moved * pixel_mm, need); + else if (std::isfinite(need)) beam_centre_row = fmt::format("the file's is within the {:.1f} px this geometry needs ({:.1f} px " "from the refined one)", need, moved); } diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index 64d68342b..ff0e11731 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -1209,7 +1209,7 @@ TEST_CASE("ResultReport_BeamCentreOff", "[Diagnostics]") { ProcessResult result; result.post_refine = PostRefineResult{}; result.beam_centre_need_pxl = 2.0f; - result.background_beam_centre = BeamCenterEstimate{.beam_x_pxl = static_cast(file.first) + 5.0f, + result.background_beam_centre = BeamCenterEstimate{.beam_x_pxl = static_cast(file.first) + 7.0f, .beam_y_pxl = static_cast(file.second), .sigma_pxl = 0.5f}; @@ -1221,12 +1221,18 @@ TEST_CASE("ResultReport_BeamCentreOff", "[Diagnostics]") { CHECK(text.find("\nBEAM_CENTRE_SHIFT= 1.00\n") != std::string::npos); CHECK(text.find("\nBEAM_CENTRE_TOLERANCE= 2.00\n") != std::string::npos); - // Refined 4 px off the file's centre: warned, with both centres, px and mm, and the advice. + // 4 px off: over the tolerance but under three times it - stated, not warned. result.used_direct_beam_x_pxl = file.first; result.used_direct_beam_y_pxl = file.second + 4.0; text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("BEAM_CENTRE_OFF") == std::string::npos); + CHECK(text.find("under three times it") != std::string::npos); + + // Refined 8 px off the file's centre: warned, with both centres, px and mm, and the advice. + result.used_direct_beam_y_pxl = file.second + 8.0; + text = RenderResultReport("p", "in.h5", x, result); CHECK(text.find("\nPATHOLOGY_FLAGS= BEAM_CENTRE_OFF\n") != std::string::npos); - CHECK(text.find("4.0 px (0.30 mm)") != std::string::npos); + CHECK(text.find("8.0 px (0.60 mm)") != std::string::npos); CHECK(text.find("correct the beam centre at the beamline") != std::string::npos); // Nothing refined (stills): the background alone warns, worded as a check. From c0bcf1a7a9217be4e911b84baa93a71b628fdd18 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 19:09:27 +0200 Subject: [PATCH 6/6] rugnux report: XDS geometry keywords as one report section, not p_XDS.INP The refined geometry is no longer written as a whole p_XDS.INP (and the report's XDS_INPUT= key is gone). Instead the report carries one section, XDS GEOMETRY KEYWORDS, with the XDS.INP keywords that carry the geometry the result was integrated at: ORGX/ORGY (PONI, 1-based), DETECTOR_DISTANCE, the detector axes where the detector is tilted, INCIDENT_BEAM_DIRECTION with X-RAY_WAVELENGTH, and ROTATION_AXIS in XDS's sense (rugnux's axis negated, as verified with XDS IDXREF). No cell or space group: XDS determines those. The lines carry a "# XDS.INP> " prefix, so sed lifts them ready to paste. The same keywords were logged by LogXDSGeometry at the start of every pass, with the axis in rugnux's sense; that printing is removed, so the report section is the one place they are stated. Docs updated (RUGNUX_REPORT, RUGNUX_INTEGRATION, RUGNUX, RUGNUX_TUTORIAL). Output only; p.mtz md5-identical on 3 battery sets. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/RUGNUX.md | 3 +- docs/RUGNUX_INTEGRATION.md | 24 +++--- docs/RUGNUX_REPORT.md | 28 ++++++- docs/RUGNUX_TUTORIAL.md | 14 ---- rugnux/CMakeLists.txt | 2 - rugnux/ResultReport.cpp | 64 +++++++++++++- rugnux/Rugnux.cpp | 68 --------------- rugnux/Rugnux.h | 4 +- rugnux/XDSInput.cpp | 167 ------------------------------------- rugnux/XDSInput.h | 29 ------- rugnux/rugnux_cli.cpp | 5 -- tests/ResultReportTest.cpp | 32 +++---- 12 files changed, 113 insertions(+), 327 deletions(-) delete mode 100644 rugnux/XDSInput.cpp delete mode 100644 rugnux/XDSInput.h diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 2b8c7f7e0..8b3198a35 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -52,7 +52,7 @@ that, or to lift the per-image loop's default ceiling of 16 workers per GPU. Nothing more is needed to pick the workflow: a dataset carrying a **goniometer axis** is processed as a rotation sweep, one without as **independent stills**, and scaling and merging run by default -in both. A rotation run that merges — the default — leaves nine files next to each other: +in both. A rotation run that merges — the default — leaves eight files next to each other: ``` myrun.mtz merged intensities + French-Wilson amplitudes, for CCP4 / phenix @@ -63,7 +63,6 @@ myrun_unmerged.mtz every observation before scaling, for pointless / aimless / myrun_report.txt what the run determined: cell, space group, statistics, warnings myrun_plot.txt one row per image, for plotting how the crystal behaved over the sweep myrun_detector.jpg the detector with the pixel mask and the beam-stop shadow drawn on it -myrun_XDS.INP the geometry the run refined, as an XDS input - a starting point for XDS myrun_log.txt the full timestamped log of the run (debug level) ``` diff --git a/docs/RUGNUX_INTEGRATION.md b/docs/RUGNUX_INTEGRATION.md index 542f5d8c5..f97e653f9 100644 --- a/docs/RUGNUX_INTEGRATION.md +++ b/docs/RUGNUX_INTEGRATION.md @@ -504,23 +504,27 @@ sweep that does not have them will get this far and no further. ## Comparing the geometry with XDS -Every run logs the detector geometry a second time in **XDS's convention**, so it can be read -straight across against the `IDXREF.LP` / `CORRECT.LP` of an XDS run on the same data: +The report states the geometry the result was integrated at in **XDS's convention**, in its +`XDS GEOMETRY KEYWORDS` section (see [the results report](RUGNUX_REPORT.md)), so it can be read +straight across against the `IDXREF.LP` / `CORRECT.LP` of an XDS run on the same data, or pasted into +`XDS.INP`: ``` -XDS convention: ORGX= 1091.00 ORGY= 1137.00 DETECTOR_DISTANCE= 75.0000 -XDS convention: DIRECTION_OF_DETECTOR_X-AXIS= 1.000000 0.000000 0.000000 -XDS convention: DIRECTION_OF_DETECTOR_Y-AXIS= 0.000000 1.000000 0.000000 -XDS convention: INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= 1.000000 QX= QY= 0.075000 -XDS convention: ROTATION_AXIS= -1.000000 0.000000 0.000000 +# XDS.INP> ORGX= 1091.00 ORGY= 1137.00 +# XDS.INP> DETECTOR_DISTANCE= 75.000 +# XDS.INP> INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= 1.000000 +# XDS.INP> ROTATION_AXIS= 1.000000 0.000000 0.000000 ``` -XDS is never given this geometry — the [XDS plugin](SOFTWARE_INTEGRATION.md) supplies image data -only, and XDS refines its own from `XDS.INP` — which is what makes the comparison worth having. The +The [XDS plugin](SOFTWARE_INTEGRATION.md) supplies image data only, so XDS refines its own geometry +from `XDS.INP` - which is what makes the comparison worth having. The two laboratory frames coincide (x along increasing detector column, y along increasing row, z along the beam), so the numbers are directly comparable, and a tilt appears as the two detector axis vectors rather than as angles, which is how XDS reports it after refinement. Two things to keep in mind: **`ORGX`/`ORGY` are 1-based**, because XDS counts pixels from 1 and Jungfraujoch from 0; and they are the **PONI**, the same quantity Jungfraujoch's beam centre is — so no correction is needed — -but not the direct beam once the detector is tilted (see above). +but not the direct beam once the detector is tilted (see above). **`ROTATION_AXIS` is rugnux's axis +negated**: rugnux turns the crystal by -phi about its axis, XDS by +phi, and with rugnux's own sign +`IDXREF` indexes only the chance floor of a mis-signed spindle. The detector axes appear only where the +detector is tilted. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 162a97c74..cf20d76cf 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -152,8 +152,32 @@ BEAM_CENTRE_SHIFT= 7.15 BEAM_CENTRE_TOLERANCE= 3.49 ``` -**`XDS_INPUT=`** in section 1 names `_XDS.INP`, the same geometry written as an XDS input -(rotation runs; see the tutorial's list of output files). +**`XDS GEOMETRY KEYWORDS`**, the section after section 1, is the same geometry as the `XDS.INP` +keywords that carry it, ready to paste into an XDS run: `ORGX`/`ORGY` (the PONI, numbered from 1 as XDS +numbers pixels), `DETECTOR_DISTANCE`, the detector axes `DIRECTION_OF_DETECTOR_X-AXIS`/`Y-AXIS` where the +detector is tilted (they carry the tilt), `INCIDENT_BEAM_DIRECTION` with `X-RAY_WAVELENGTH`, and +`ROTATION_AXIS` - the axis the run integrated at, **negated**, because XDS turns the crystal the other way +about the axis it is given (checked with XDS `IDXREF`). Geometry only: XDS determines the cell and the +space group itself. The lines carry a fixed prefix, so + +```bash +sed -n 's/^# XDS.INP> //p' myrun_report.txt +``` + +prints them as XDS reads them: + +``` +ORGX= 749.15 ORGY= 817.32 +DETECTOR_DISTANCE= 175.668 +DIRECTION_OF_DETECTOR_X-AXIS= 1.000000 0.000000 0.000444 +DIRECTION_OF_DETECTOR_Y-AXIS= 0.000000 1.000000 -0.000976 +INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= 0.918400 +ROTATION_AXIS= 1.000000 0.000000 0.000000 +``` + +The conventions are XDS's own (W. Kabsch, "XDS" (2010), Acta Cryst. D66, 125-132 +[doi:10.1107/S0907444909047337](https://doi.org/10.1107/S0907444909047337)); this section is the one +place the run states them. **Which pass.** A rotation run integrates twice — once at the geometry in the input file (the first half of the sweep but at least 120°, `--prepass-fraction`), then the whole sweep again at the post-refined geometry — and can integrate a third time if a guard rejects the second pass. diff --git a/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index 8f68b2acd..213fe8db5 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -331,20 +331,6 @@ Absorption is corrected only by the empirical surface a rotation run fits and ke is missing. A shadow drawn where the background is unbroken, or a dark patch the run left uncoloured, is visible at a glance; read it against the `Beam stop shadow: N pixels ...` line in the log. -- `_XDS.INP` — **the geometry the run ended on, as an XDS input**: on a rotation run the - post-refined beam centre (`ORGX`/`ORGY`, the foot of the perpendicular, numbered from 1 as XDS - numbers pixels), distance, detector axes (`DIRECTION_OF_DETECTOR_X-AXIS`/`Y-AXIS`, which carry the - tilt), wavelength, rotation axis (the one the run integrated at, after any sign rescue, negated: XDS turns - the crystal the other way about the axis it is given), oscillation - (scaled where the run adopted a goniometer rotation scale), starting angle and data range, plus the - space group and cell (left commented out where the group is not in the reference setting of its - number, which is the only setting a number names), `FRIEDEL'S_LAW` per `-A` and the merged - resolution range. For HDF5 the template is the master file and the `LIB=` line for the reader - plugin is left for the user to point at durin or neggia; the file's own pixel mask reaches XDS - through the plugin, the beam-stop shadow rugnux detected does not (XDS finds its own in `DEFPIX`). - The conventions are XDS's own (W. Kabsch, "XDS" (2010), Acta Cryst. D66, 125-132 - [doi:10.1107/S0907444909047337](https://doi.org/10.1107/S0907444909047337)). It is a starting - point - XDS refines the geometry again - and the report names it as **`XDS_INPUT=`** in section 1. Merged statistics (⟨I/σ⟩, CC1/2, completeness, …), the error model and timing are printed to the console. By default the written resolution is trimmed automatically where CC1/2 falls off diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 5db9f5764..6b437d83b 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -23,8 +23,6 @@ ADD_LIBRARY(JFJochRugnux STATIC DiagnosticOutput.h BrokerFeed.cpp BrokerFeed.h - XDSInput.cpp - XDSInput.h ) TARGET_LINK_LIBRARIES(JFJochRugnux JFJochReader JFJochImageAnalysis JFJochWriter JFJochModelWriter JFJochPreview gemmi fftw3f diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 4be52ae86..3b563f951 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -43,7 +43,8 @@ namespace { // --developer key on rotation runs; POWDER_ICE_* / POWDER_NON_ICE_SPOT_FRACTION; // the ICE_RINGS and POWDER_RINGS flags. // 16: FILE_DIRECT_BEAM / BACKGROUND_DIRECT_BEAM / BEAM_CENTRE_SHIFT / BEAM_CENTRE_TOLERANCE in - // section 1 and the BEAM_CENTRE_OFF flag - the file's beam centre against the refined one. + // section 1 and the BEAM_CENTRE_OFF flag - the file's beam centre against the refined one; + // the XDS GEOMETRY KEYWORDS section; OVERLOADED_* / OVERLAPPED_* and the OVERLOADS / OVERLAPS flags. // 8: every line that is not `KEY= value` and not blank starts with `#`, warnings included // (`# WARNING:`); the header gains RUGNUX_DOWNLOAD and BUILD_CXX_FLAGS, and RUGNUX_GIT is // the full build-time hash. @@ -553,9 +554,6 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (result.spot_resolution_estimate_A.has_value()) Add(s, KeyReal("SPOT_RESOLUTION_ESTIMATE", *result.spot_resolution_estimate_A, "{:.2f}")); - // The same geometry once more, as an XDS input (rugnux/XDSInput.h). - if (!result.xds_input_path.empty()) - Add(s, KeyText("XDS_INPUT", result.xds_input_path)); if (result.pass_count > 1) { Add(s, KeyText("PASS", fmt::format("{} of {}", result.pass_number, result.pass_count))); Add(s, KeyText("PASS_DECISION", result.pass_decision, true)); @@ -682,6 +680,64 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, doc.sections.push_back(std::move(s)); } + // ------------------------------------------------- XDS GEOMETRY KEYWORDS + // The geometry the result was integrated at - on a rotation run the post-refined one - as the XDS.INP + // keywords that carry it, so it can be pasted into an XDS run rather than retyped. Geometry only: + // XDS determines the cell and the space group itself. Comment lines, like all prose here, with a + // fixed prefix so `sed -n 's/^# XDS.INP> //p'` lifts them. + // + // Following the XDS conventions, Kabsch (2010) Acta Cryst. D66, 125-132: XDS places a pixel at + // x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis) + // in a laboratory frame that is ours (x along a column step, y along a row step, z along the beam), + // with pixels numbered from 1 and ORGX/ORGY the foot of the perpendicular from the crystal - our PONI + // (docs/DETECTOR_GEOMETRY.md). A tilt therefore goes out as the two detector axes, taken as + // differences of LabCoord. + // + // ROTATION_AXIS is rugnux's NEGATED: rugnux turns the crystal by -phi about its axis (IndexAndRefine + // takes the lattice to an image with GetTransformationAngle(-angle)), XDS by a right-handed +phi about + // ROTATION_AXIS. Checked with IDXREF on two sets given this geometry: the negated axis indexes 99.5% + // and 14% of the spots, rugnux's own sign 37% and 0.04% (the chance floor of a mis-signed spindle). + // The axis is the one the run integrated at, after any sign rescue. + if (std::isfinite(result.used_beam_x_pxl) && std::isfinite(result.used_distance_mm)) { + constexpr double RAD = PI / 180.0; + const auto &t = result.used_detector_tilt_deg; + DiffractionExperiment x = experiment; + x.BeamX_pxl(result.used_beam_x_pxl).BeamY_pxl(result.used_beam_y_pxl) + .DetectorDistance_mm(result.used_distance_mm) + .PoniRot1_rad(static_cast(t[0] * RAD)).PoniRot2_rad(static_cast(t[1] * RAD)) + .PoniRot3_rad(static_cast(t[2] * RAD)); + const auto geom = x.GetDiffractionGeometry(); + const float bx = geom.GetBeamX_pxl(), by = geom.GetBeamY_pxl(), pixel_mm = geom.GetPixelSize_mm(); + const Coord origin = geom.LabCoord(bx, by); + const Coord x_axis = ((geom.LabCoord(bx + 1.0f, by) - origin) / pixel_mm).Normalize(); + const Coord y_axis = ((geom.LabCoord(bx, by + 1.0f) - origin) / pixel_mm).Normalize(); + + std::string k; + const auto line = [&k](const std::string &l) { k += " XDS.INP> " + l + "\n"; }; + line(fmt::format("ORGX= {:.2f} ORGY= {:.2f}", bx + 1.0f, by + 1.0f)); + line(fmt::format("DETECTOR_DISTANCE= {:.3f}", geom.GetDetectorDistance_mm())); + if (t[0] != 0.0 || t[1] != 0.0 || t[2] != 0.0) { + line(fmt::format("DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}", x_axis.x, x_axis.y, x_axis.z)); + line(fmt::format("DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}", y_axis.x, y_axis.y, y_axis.z)); + } + line(fmt::format("INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= {:.6f}", geom.GetWavelength_A())); + if (const auto gonio = experiment.GetGoniometer()) { + const Coord axis = -result.used_rotation_axis.value_or(gonio->GetAxis()); + // + 0.0f: no "-0.000000" from a negated zero. + line(fmt::format("ROTATION_AXIS= {:.6f} {:.6f} {:.6f}", axis.x + 0.0f, axis.y + 0.0f, axis.z + 0.0f)); + } + k.pop_back(); + + ReportSection s; + s.title = "XDS GEOMETRY KEYWORDS"; + Add(s, Prose(" The geometry this result was integrated at, as XDS.INP keywords - ORGX/ORGY numbered from\n" + " 1, ROTATION_AXIS in XDS's sense of rotation. Geometry only; XDS determines the cell and\n" + " space group itself. Lift them with: sed -n 's/^# XDS.INP> //p' ")); + Add(s, Blank()); + Add(s, Prose(k)); + doc.sections.push_back(std::move(s)); + } + // ------------------------------------------------------------- 2. CRYSTAL { ReportSection s; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 33dd36c9f..41f94a17f 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -309,73 +309,6 @@ Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment, } } -// Report the detector geometry the way XDS states it, so a run can be compared line by line against -// the IDXREF.LP / CORRECT.LP of XDS on the same data. XDS is never given this geometry - the durin -// plugin hands it image data only (plugin_get_header returns dimensions, bytes per pixel, pixel size -// and frame count, nothing else) and XDS refines its own from XDS.INP - which is precisely why -// having ours in the same convention is worth printing. -// -// The two laboratory frames coincide: x along increasing detector column, y along increasing row, -// z along the beam. XDS places a pixel at -// x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis) -// and DiffractionGeometry::LabCoord is that same expression with X_axis = poni_rot*(1,0,0) and -// Y_axis = poni_rot*(0,1,0). So the tilt goes out as the two axis vectors rather than as angles, -// which is also how XDS reports it after refinement. Taking the axes as differences of LabCoord -// keeps them in step with whatever the geometry currently is, tilt included. -// -// ORGX/ORGY are 1-based - XDS counts pixels from 1 - and are the foot of the perpendicular from the -// crystal rather than the direct beam, which is what our beam centre is too -// (docs/DETECTOR_GEOMETRY.md), so they map across without correction. -// -// ROTATION_AXIS is printed once per pass, and it is not always the vector stored in the file: the -// sign rescue further down adopts the opposite sign when the file's own indexes almost nothing, and -// the following pass then prints what it adopted. On the corpus the first and last line of a run -// disagree in sign on 12 of 115 rotation runs, so the LAST one - the axis the run actually -// integrated at - is the one to read. -// -// That sign is not reliably the one XDS wants - but it is not a constant either, so do NOT negate -// this line on principle. It is paired with the detector axis directions printed just above it and -// has to be read together with them: two rotation datasets whose detector axes are both +x/+y wanted -// the negation, while a third whose axes are negative wanted the sign as printed. The measure is the -// fraction of COLSPOT spots IDXREF can index - 8.9% as printed against 92.6% negated, and 1060/10080 -// against 9624/10080; the other way round, 43.9% as printed against 13.5%. -// -// Read the HIGH member of each pair. That statistic has a chance floor which reaches about 15% on a -// dense reciprocal lattice - a null with the frame numbers permuted scored 14.9% against a real -// 15.4% - so a low value means "IDXREF gave up", not a measured near-zero: both wrong-sign runs -// above stopped on INSUFFICIENT PERCENTAGE (< 50%) OF INDEXED REFLECTIONS, one of them after -// refining the right cell and rejecting 9020 spots as too far from their ideal positions, which is -// the signature of a mis-signed spindle. So try both signs and keep the one that indexes; a sign -// carried over from another file is worth nothing. -// -// One of those two runs read a file our own converter had written, whose axis default was the -// negative sign. It therefore says what XDS does with a given sign, which is what this comment is -// about, and nothing about what a foreign file stores. -static void LogXDSGeometry(const DiffractionExperiment &experiment) { - Logger logger("Rugnux"); - - const auto geom = experiment.GetDiffractionGeometry(); - const float beam_x = geom.GetBeamX_pxl(); - const float beam_y = geom.GetBeamY_pxl(); - const float pixel_mm = geom.GetPixelSize_mm(); - - const Coord origin = geom.LabCoord(beam_x, beam_y); - const Coord x_axis = ((geom.LabCoord(beam_x + 1.0f, beam_y) - origin) / pixel_mm).Normalize(); - const Coord y_axis = ((geom.LabCoord(beam_x, beam_y + 1.0f) - origin) / pixel_mm).Normalize(); - - logger.Info("XDS convention: ORGX= {:.2f} ORGY= {:.2f} DETECTOR_DISTANCE= {:.4f}", - beam_x + 1.0f, beam_y + 1.0f, geom.GetDetectorDistance_mm()); - logger.Info("XDS convention: DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}", - x_axis.x, x_axis.y, x_axis.z); - logger.Info("XDS convention: DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}", - y_axis.x, y_axis.y, y_axis.z); - logger.Info("XDS convention: INCIDENT_BEAM_DIRECTION= 0 0 1 X-RAY_WAVELENGTH= {:.6f} QX= QY= {:.6f}", - geom.GetWavelength_A(), pixel_mm); - if (const auto goniometer = experiment.GetGoniometer()) - logger.Info("XDS convention: ROTATION_AXIS= {:.6f} {:.6f} {:.6f}", - goniometer->GetAxis().x, goniometer->GetAxis().y, goniometer->GetAxis().z); -} - namespace rugnux_internal { // The FFT recovers a basis vector only up to fft_max_unit_cell_A: FFTIndexer sizes its projected // histogram from that value and the transform's last usable bin IS that length. So a reference cell @@ -787,7 +720,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b : calibration ? "powder calibration" : "azimuthal integration"; logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]", images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name); - LogXDSGeometry(experiment_); if (observer) observer->OnPhase(full ? "Full analysis" : calibration ? "Powder calibration" : "Azimuthal integration"); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index cbe0bafbd..55a34ad6e 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -302,7 +302,7 @@ struct ProcessResult { std::array used_detector_tilt_deg{}; double used_direct_beam_x_pxl = 0.0, used_direct_beam_y_pxl = 0.0; // The rotation axis the same pass integrated at - the sign rescue can adopt the opposite of the - // file's (see LogXDSGeometry). Empty on stills. + // file's. The report's XDS keyword block reads it. Empty on stills. std::optional used_rotation_axis; std::optional consensus_cell; // The centring `consensus_cell` is written in; empty when nothing named one (the stills path @@ -329,8 +329,6 @@ struct ProcessResult { std::optional lattice_angle_departure_deg; MeanProcessingTime mean_processing_time{}; std::optional written_master_path; - // _XDS.INP, where the CLI wrote the run's geometry as an XDS input (rugnux/XDSInput.h). - std::string xds_input_path; std::string merge_statistics_text; // populated when scaling/merging ran // Structured merge statistics (per-shell + overall), populated when scaling/merging ran. ISa is diff --git a/rugnux/XDSInput.cpp b/rugnux/XDSInput.cpp deleted file mode 100644 index e8361c8f8..000000000 --- a/rugnux/XDSInput.cpp +++ /dev/null @@ -1,167 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include -#include -#include -#include -#include - -#include - -#include "../common/GitInfo.h" -#include "../common/JFJochMath.h" -#include "XDSInput.h" - -std::string RenderXDSInput(const std::string &input_file, - const DiffractionExperiment &experiment, - const ProcessResult &result, - int64_t first_image, int64_t end_image) { - // The geometry the result was integrated at, on a copy of the file's experiment. - DiffractionExperiment x = experiment; - if (std::isfinite(result.used_beam_x_pxl) && std::isfinite(result.used_distance_mm)) { - constexpr double RAD = PI / 180.0; - const auto &t = result.used_detector_tilt_deg; - x.BeamX_pxl(result.used_beam_x_pxl).BeamY_pxl(result.used_beam_y_pxl) - .DetectorDistance_mm(result.used_distance_mm) - .PoniRot1_rad(static_cast(t[0] * RAD)).PoniRot2_rad(static_cast(t[1] * RAD)) - .PoniRot3_rad(static_cast(t[2] * RAD)); - } - const auto geom = x.GetDiffractionGeometry(); - const float beam_x = geom.GetBeamX_pxl(); - const float beam_y = geom.GetBeamY_pxl(); - const float pixel_mm = geom.GetPixelSize_mm(); - - // Following the XDS conventions, Kabsch (2010) Acta Cryst. D66, 125-132: XDS places a pixel at - // x_lab(i,j) = (i-ORGX)*QX*X_axis + (j-ORGY)*QY*Y_axis + DETECTOR_DISTANCE*(X_axis x Y_axis) - // in a laboratory frame that is ours (x along a column step, y along a row step, z along the - // beam), with pixels numbered from 1 and ORGX/ORGY the foot of the perpendicular from the crystal - - // our PONI (see LogXDSGeometry in Rugnux.cpp and docs/DETECTOR_GEOMETRY.md). So the tilt goes out - // as the two detector axes, taken as differences of LabCoord. - const Coord origin = geom.LabCoord(beam_x, beam_y); - const Coord x_axis = ((geom.LabCoord(beam_x + 1.0f, beam_y) - origin) / pixel_mm).Normalize(); - const Coord y_axis = ((geom.LabCoord(beam_x, beam_y + 1.0f) - origin) / pixel_mm).Normalize(); - - const bool hdf5 = input_file.ends_with(".h5") || input_file.ends_with(".nxs"); - const bool cbf = input_file.ends_with(".cbf") || input_file.ends_with(".cbf.gz"); - std::string input_path = input_file; - try { - input_path = std::filesystem::absolute(input_file).string(); - } catch (const std::exception &) { - } - - // A CBF sweep is one file per frame: XDS takes a template with ? for the frame number - the last run - // of digits in the name - and numbers the frames by it, so the data range starts at the given file's. - std::string template_path = input_path; - int64_t first_frame = 1; - // The HDF5 plugins take the master file's name with "master" as ??????, as XDS has it for - // an EIGER data set. - if (hdf5 && input_path.ends_with("_master.h5")) - template_path.replace(template_path.size() - 9, 6, "??????"); - if (cbf) { - const size_t name = template_path.find_last_of('/') + 1; - const size_t end = template_path.find_last_of("0123456789"); - if (end != std::string::npos && end >= name) { - size_t begin = end; - while (begin > name && std::isdigit(static_cast(template_path[begin - 1]))) - --begin; - first_frame = std::stoll(template_path.substr(begin, end - begin + 1)); - template_path.replace(begin, end - begin + 1, std::string(end - begin + 1, '?')); - } - } - - std::string s; - s += fmt::format("! XDS.INP written by rugnux {} from {}\n", jfjoch_version(), input_path); - s += "! The detector geometry is the one rugnux integrated at - on a rotation run the post-refined\n" - "! one - not the file's; XDS refines it again from here. Check before use: this file is a\n" - "! starting point, not a validated XDS setup.\n"; - s += "JOB= XYCORR INIT COLSPOT IDXREF DEFPIX INTEGRATE CORRECT\n\n"; - - if (hdf5) - s += "! HDF5 is read through a plugin: uncomment and point LIB= at durin-plugin.so (or the\n" - "! Dectris neggia plugin), which opens the _master.h5 the template names with ??????.\n" - "! The plugin applies the file's own pixel mask.\n" - "!LIB= /path/to/durin-plugin.so\n"; - else if (!cbf) - s += "! rugnux read this input itself; XDS may need the frames converted, and the template\n" - "! below written with ?????? in place of the frame number.\n"; - s += fmt::format("NAME_TEMPLATE_OF_DATA_FRAMES= {}\n", template_path); - s += fmt::format("DATA_RANGE= {} {}\n", first_frame + first_image, first_frame + end_image - 1); - s += fmt::format("SPOT_RANGE= {} {}\n\n", first_frame + first_image, first_frame + end_image - 1); - - s += fmt::format("DETECTOR= {} MINIMUM_VALID_PIXEL_VALUE= 0 OVERLOAD= {}\n", - cbf ? "PILATUS" : "EIGER", SaturationValueFromLimit(experiment.GetSaturationLimit())); - s += fmt::format("NX= {} NY= {} QX= {:.6f} QY= {:.6f}\n", - experiment.GetXPixelsNumConv(), experiment.GetYPixelsNumConv(), pixel_mm, pixel_mm); - const float thickness_um = experiment.GetDetectorSetup().GetSensorThickness_um(); - if (thickness_um > 0.0f) - s += fmt::format("SENSOR_THICKNESS= {:.3f}\n", thickness_um / 1000.0f); - s += fmt::format("ORGX= {:.2f} ORGY= {:.2f}\n", beam_x + 1.0f, beam_y + 1.0f); - s += fmt::format("DETECTOR_DISTANCE= {:.3f}\n", geom.GetDetectorDistance_mm()); - s += fmt::format("DIRECTION_OF_DETECTOR_X-AXIS= {:.6f} {:.6f} {:.6f}\n", x_axis.x, x_axis.y, x_axis.z); - s += fmt::format("DIRECTION_OF_DETECTOR_Y-AXIS= {:.6f} {:.6f} {:.6f}\n", y_axis.x, y_axis.y, y_axis.z); - s += "INCIDENT_BEAM_DIRECTION= 0 0 1\n"; - s += fmt::format("X-RAY_WAVELENGTH= {:.6f}\n", geom.GetWavelength_A()); - s += "! Polarization is not in the file; a typical synchrotron:\n" - "!FRACTION_OF_POLARIZATION= 0.99 POLARIZATION_PLANE_NORMAL= 0 1 0\n\n"; - - if (const auto gonio = experiment.GetGoniometer()) { - // Negated: rugnux turns the crystal by -phi about its axis (IndexAndRefine takes the lattice - // to an image with GetTransformationAngle(-angle)), XDS by a right-handed +phi about - // ROTATION_AXIS. Checked with IDXREF on the export: 99.5% of the spots index with the - // negated axis, 37% (the chance floor of a mis-signed spindle) with rugnux's own. - const Coord axis = -result.used_rotation_axis.value_or(gonio->GetAxis()); - // The oscillation the run integrated at: the file's, scaled where the run adopted a goniometer - // rotation scale (as OSCILLATION_RANGE in the report). - const double scale = result.post_refine && result.post_refine->rotation_scale_suspect - ? result.post_refine->rotation_scale : 1.0; - // + 0.0f: no "-0.000000" from a negated zero. - s += fmt::format("ROTATION_AXIS= {:.6f} {:.6f} {:.6f}\n", axis.x + 0.0f, axis.y + 0.0f, axis.z + 0.0f); - s += fmt::format("OSCILLATION_RANGE= {:.6f}\n", gonio->GetIncrement_deg() * scale); - s += fmt::format("STARTING_FRAME= {}\n", first_frame); - s += fmt::format("STARTING_ANGLE= {:.4f}\n\n", gonio->GetStart_deg()); - } - - if (result.space_group && result.consensus_cell) { - const auto &c = *result.consensus_cell; - // A number names only the reference setting; a cell in another one would be read wrongly. - const char *prefix = result.space_group->is_reference_setting() ? "" : "!"; - if (*prefix) - s += fmt::format("! The run's space group {} is not in the reference setting of its number;\n" - "! XDS would read the cell below in the wrong setting, so both are left out.\n", - result.space_group->xhm()); - s += fmt::format("{}SPACE_GROUP_NUMBER= {}\n", prefix, result.space_group->number); - s += fmt::format("{}UNIT_CELL_CONSTANTS= {:.3f} {:.3f} {:.3f} {:.3f} {:.3f} {:.3f}\n", prefix, - c.a, c.b, c.c, c.alpha, c.beta, c.gamma); - } - s += fmt::format("FRIEDEL'S_LAW= {}\n", experiment.GetScalingSettings().GetMergeFriedel() ? "TRUE" : "FALSE"); - if (result.has_merge_statistics) - s += fmt::format("INCLUDE_RESOLUTION_RANGE= {:.3f} {:.3f}\n", - result.merge_statistics.overall.d_max, result.merge_statistics.overall.d_min); - s += "\n! rugnux also masked the beam-stop shadow it detected (_detector.jpg); XDS finds its\n" - "! own in DEFPIX, or add UNTRUSTED_RECTANGLE= / UNTRUSTED_QUADRILATERAL= lines here.\n"; - return s; -} - -std::string WriteXDSInput(const std::string &output_prefix, - const std::string &input_file, - const DiffractionExperiment &experiment, - const ProcessResult &result, - int64_t first_image, int64_t end_image, - Logger &logger) { - if (output_prefix.empty() || !experiment.GetGoniometer()) - return {}; - const std::string filename = output_prefix + "_XDS.INP"; - try { - std::ofstream f(filename); - if (!f) - throw std::runtime_error("cannot open for writing"); - f << RenderXDSInput(input_file, experiment, result, first_image, end_image); - if (!f) - throw std::runtime_error("write failed"); - } catch (const std::exception &e) { - logger.Warning("Could not write the XDS input {}: {}", filename, e.what()); - return {}; - } - return filename; -} diff --git a/rugnux/XDSInput.h b/rugnux/XDSInput.h deleted file mode 100644 index 393cfd7ef..000000000 --- a/rugnux/XDSInput.h +++ /dev/null @@ -1,29 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include - -#include "../common/DiffractionExperiment.h" -#include "../common/Logger.h" -#include "Rugnux.h" - -// _XDS.INP: the geometry the run integrated at - on a rotation run the post-refined one - as an -// XDS input, so the same data can be processed by XDS from where rugnux ended rather than from the -// header. `experiment` is the one the run was given (the file's), `result` carries what the run -// refined; the images are the file's 0-based [first_image, end_image). -std::string RenderXDSInput(const std::string &input_file, - const DiffractionExperiment &experiment, - const ProcessResult &result, - int64_t first_image, int64_t end_image); - -// Writes _XDS.INP and returns its name; returns an empty string, writing nothing, on a -// run with no output prefix or no rotation, and logs and swallows a failure to write. -std::string WriteXDSInput(const std::string &output_prefix, - const std::string &input_file, - const DiffractionExperiment &experiment, - const ProcessResult &result, - int64_t first_image, int64_t end_image, - Logger &logger); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 4e9d30bc5..9e74ef787 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -29,7 +29,6 @@ #include "../reader/JFJochSMVReader.h" #include "../reader/JFJochStartMessageReader.h" #include "BrokerFeed.h" -#include "XDSInput.h" #include "../writer/HDF5NXmx.h" #include #include @@ -3064,10 +3063,6 @@ static int RunRugnux(int argc, char **argv) { provenance.wall_time_s = std::chrono::duration( std::chrono::steady_clock::now() - invocation_start).count(); TimingMark report_mark(logger, "results report"); - // The geometry the run ended on, as an XDS input beside the report, which names it. - if (rotation_indexing) - result.xds_input_path = WriteXDSInput(output_prefix, input_file, experiment, result, start_image, - end_image, logger); const ReportDocument doc = WriteResultReport(output_prefix, input_file, experiment, result, logger, provenance); report_mark.Stop(); diff --git a/tests/ResultReportTest.cpp b/tests/ResultReportTest.cpp index ff0e11731..5d60a7a0d 100644 --- a/tests/ResultReportTest.cpp +++ b/tests/ResultReportTest.cpp @@ -7,7 +7,6 @@ #include "../image_analysis/scale_merge/Merge.h" #include "../reader/JFJochHDF5Reader.h" #include "../rugnux/ResultReport.h" -#include "../rugnux/XDSInput.h" #include "../common/CUDAWrapper.h" #include "../writer/FileWriter.h" @@ -1271,35 +1270,26 @@ TEST_CASE("ResultReport_OverloadsAndOverlaps", "[Diagnostics]") { CHECK(text.find("Move the detector back") != std::string::npos); } -TEST_CASE("XDSInput_Render", "[Diagnostics]") { +TEST_CASE("ResultReport_XDSGeometryKeywords", "[Diagnostics]") { DiffractionExperiment x(DetJF(1)); x.BeamX_pxl(500.0f).BeamY_pxl(400.0f).DetectorDistance_mm(100.0f); x.Goniometer(GoniometerAxis("omega", 10.0f, 0.1f, Coord(1, 0, 0), {})); ProcessResult result; + result.consensus_cell = UnitCell{.a = 79.0f, .b = 79.0f, .c = 38.0f, + .alpha = 90.0f, .beta = 90.0f, .gamma = 90.0f}; // The refined geometry, not the file's, goes out - 1-based, untilted here. result.used_beam_x_pxl = 510.0f; result.used_beam_y_pxl = 395.5f; result.used_distance_mm = 120.0f; result.used_rotation_axis = Coord(-1, 0, 0); - const auto text = RenderXDSInput("/data/x_master.h5", x, result, 0, 900); - CHECK(text.find("\nORGX= 511.00 ORGY= 396.50\n") != std::string::npos); - CHECK(text.find("\nDETECTOR_DISTANCE= 120.000\n") != std::string::npos); - CHECK(text.find("\nDIRECTION_OF_DETECTOR_X-AXIS= 1.000000 0.000000 0.000000\n") != std::string::npos); + const auto text = RenderResultReport("p", "in.h5", x, result); + CHECK(text.find("XDS GEOMETRY KEYWORDS") != std::string::npos); + CHECK(text.find("\n# XDS.INP> ORGX= 511.00 ORGY= 396.50\n") != std::string::npos); + CHECK(text.find("\n# XDS.INP> DETECTOR_DISTANCE= 120.000\n") != std::string::npos); // XDS turns the crystal the other way about the axis it is given. - CHECK(text.find("\nROTATION_AXIS= 1.000000 0.000000 0.000000\n") != std::string::npos); - CHECK(text.find("\nNAME_TEMPLATE_OF_DATA_FRAMES= /data/x_??????.h5\n") != std::string::npos); - CHECK(text.find("\nDATA_RANGE= 1 900\n") != std::string::npos); - CHECK(text.find("\nSTARTING_ANGLE= 10.0000\n") != std::string::npos); - CHECK(text.find("\nOSCILLATION_RANGE= 0.100000\n") != std::string::npos); -} - -TEST_CASE("XDSInput_CBFTemplate", "[Diagnostics]") { - DiffractionExperiment x(DetJF(1)); - x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(1, 0, 0), {})); - const auto text = RenderXDSInput("/data/run_2_0101.cbf", x, ProcessResult{}, 0, 50); - CHECK(text.find("\nNAME_TEMPLATE_OF_DATA_FRAMES= /data/run_2_????.cbf\n") != std::string::npos); - CHECK(text.find("\nDATA_RANGE= 101 150\n") != std::string::npos); - CHECK(text.find("\nSTARTING_FRAME= 101\n") != std::string::npos); - CHECK(text.find("\nDETECTOR= PILATUS") != std::string::npos); + CHECK(text.find("\n# XDS.INP> ROTATION_AXIS= 1.000000 0.000000 0.000000\n") != std::string::npos); + // Untilted: no detector axes; and never the cell or the space group. + CHECK(text.find("XDS.INP> DIRECTION_OF_DETECTOR") == std::string::npos); + CHECK(text.find("XDS.INP> UNIT_CELL") == std::string::npos); }