From 307987c86510a9eb97b794a450b1bc2fe5249e1b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 4 Sep 2026 10:41:50 +0200 Subject: [PATCH] rugnux: the mounting's cost is computed exactly from the measured group, not guessed from an angle Offline, a merged rotation run has what a still lacks - a determined point group and an exact indexed orientation - so the run-level number no longer needs the pessimistic presumed-diad bound, and it no longer uses the nearest-axis angle either. That 15-deg warning heuristic was wrong in both directions: an aligned in-plane 2-fold of a dihedral group is repaired by the principal axis, a cubic group is never severe in any orientation, and a lone diad perpendicular to the spindle is severe with no axis anywhere near the spindle at all. The group's proper rotations are applied to the sweep's blind double cone in the crystal's actual orientation, and what no operator maps onto measured territory is counted, weighted by each shell's own cone width so the result is a fraction of unique reflections to this run's resolution limit. Friedel and the improper operators need no separate handling - the cone and the measured region are both inversion-symmetric. The number is machine-readable on purpose: SPINDLE_LOST_UNIQUE_FRACTION in the report (0-1, a bare number a pipeline can act on) and /entry/MX/spindleLostUniqueFraction in the master, with the warning prose only on top of it, fired when the group recovers less than half the cone's content. REPORT_VERSION stays 7: the format's own rule is that adding a key does not move it. This also settles what the nearest-axis keys hedged: with the measured group the mounting is cleared or convicted exactly, so their documentation now calls them descriptive and points at the new key for the verdict. Verified against Monte Carlo: P1 loses 2.0% of unique reflections at theta_max = 15 deg with nothing repaired; a lone diad on or perpendicular to the spindle repairs nothing; an axis of order >= 3 perpendicular to the spindle repairs everything; 622 with an in-plane diad on the spindle loses nothing; cubic loses nothing in any orientation. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- common/JFJochMessages.h | 5 + docs/CBOR.md | 1 + docs/HDF5.md | 1 + docs/RUGNUX_REPORT.md | 20 ++-- frame_serialize/CBORStream2Deserializer.cpp | 2 + frame_serialize/CBORStream2Serializer.cpp | 1 + rugnux/CMakeLists.txt | 2 + rugnux/ResultReport.cpp | 8 ++ rugnux/Rugnux.cpp | 64 +++++++---- rugnux/Rugnux.h | 6 ++ rugnux/SpindleCuspLoss.cpp | 111 ++++++++++++++++++++ rugnux/SpindleCuspLoss.h | 36 +++++++ tests/CBORTest.cpp | 2 + tests/CMakeLists.txt | 1 + tests/SpindleCuspLossTest.cpp | 93 ++++++++++++++++ writer/HDF5NXmx.cpp | 4 + 16 files changed, 331 insertions(+), 26 deletions(-) create mode 100644 rugnux/SpindleCuspLoss.cpp create mode 100644 rugnux/SpindleCuspLoss.h create mode 100644 tests/SpindleCuspLossTest.cpp diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index c75e9dc45..c4dba77a5 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -380,6 +380,11 @@ struct EndMessage { // Run mean of the per-image spindle_blind_fraction, over the frames that had a value; absent // when none did. Written to /entry/MX/spindleBlindFractionMean. std::optional spindle_blind_fraction; + // The exact run-level counterpart, offline only (rugnux): the fraction (0-1) of unique + // reflections to the run's resolution limit that the MEASURED point group could not recover + // from the sweep's blind cone, in the crystal's indexed orientation. Not the per-image + // worst-case bound. Written to /entry/MX/spindleLostUniqueFraction. + std::optional spindle_lost_unique_fraction; std::optional end_date; diff --git a/docs/CBOR.md b/docs/CBOR.md index 6c4d3e36a..536a11803 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -303,6 +303,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | max_receiver_delay | uint64 | Internal performance of Jungfraujoch | | | bkg_estimate | float | Mean background estimate for the whole run | | | spindle_blind_fraction | float | Run mean of the per-image spindle_blind_fraction, over the frames that had one | | +| spindle_lost_unique_fraction | float | Fraction (0-1) of unique reflections the mounting made unmeasurable, exact under the measured point group; offline (rugnux) only | | | indexing_rate | float | Mean indexing rate for the whole run | | | unit_cell | object (optional) | Unit cell of the system, based on the actual experiment: a, b, c \[angstrom\] and alpha, beta, gamma \[degree\] | | | rotation_lattice_type | object | Bravais lattice classification of the total rotation solution over the run, if available; same schema as `lattice_type` | | diff --git a/docs/HDF5.md b/docs/HDF5.md index 3442b2b22..52a46f29c 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -369,6 +369,7 @@ variants. | `imageIndexedMean` | | mean indexing rate over the run | | `bkgEstimateMean` | photons | mean background over the run | | `spindleBlindFractionMean` | fraction (0-1) | mean `spindleBlindFraction` over the frames that had one | +| `spindleLostUniqueFraction` | fraction (0-1) | unique reflections (to the run's resolution limit) the mounting made unmeasurable, exact under the measured point group and indexed orientation; offline (rugnux) only | | `iceRingScoreMean` | ratio | mean `iceRingScore` over the run — the single "how icy was this dataset" number (1 = no ice) | | `indexedLatticeCount` | | per-image lattice count summary (master). *Note: data files use `indexingLatticeCount`; readers accept either.* | | `reindexMatrix` | | change of basis from the setting the per-image data are in to the setting of `/entry/sample/unit_cell` (`[9]`, `int32`, flattened 3×3, row major) — see below | diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 9f047eccc..52e5ac65d 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -243,13 +243,19 @@ nothing else. `SPINDLE_SYMMETRY_AXIS_ANGLE_DEG=` and `SPINDLE_SYMMETRY_AXIS_ORDER=` in section 4 say how the crystal sat on the goniometer: the angle between the spindle and the nearest symmetry axis, and that axis's -order. Below about 15 deg the axis maps the sweep's blind cone onto itself and the reflections in -that cone stay missing however long the sweep runs, which is why the same number also raises a -warning there. It is written on every rotation run that determined a space group, including the -ordinary well-mounted case, so an incomplete cusp can be attributed to the mounting. A large angle -does not on its own clear it: an axis within the same margin of perpendicular to the spindle carries -the blind cone onto its opposite half, which is equally unmeasured, and only the nearest axis is -reported. (New in `REPORT_VERSION= 7`.) +order. They are descriptive: neither convicts nor clears the mounting on its own, because an aligned +axis of any order maps the sweep's blind cone onto itself while an axis near perpendicular does the +same only when it is a lone 2-fold, and only the nearest axis is reported. (New in +`REPORT_VERSION= 7`.) + +`SPINDLE_LOST_UNIQUE_FRACTION=` in section 4 is the exact verdict the angle cannot give: the fraction +(0-1, so 0.0300 means 3%) of unique reflections, to this run's resolution limit, that the mounting +made unmeasurable - the part of the sweep's blind double cone that no operator of the measured point +group maps onto measured territory, computed in the crystal's actual indexed orientation. 0.0000 +means the mounting cost nothing; the run warns when the group recovers less than half of the cone's +content. The same number is written to the master file as `/entry/MX/spindleLostUniqueFraction`, so a +pipeline can read it from either output without parsing prose. Written on every rotation run that +determined a space group and merged reflections. `SPACE_GROUP_ENANTIOMORPH=` in section 4 reads **`ASSUMED_FROM_MODEL`** when the hand written in the files is the model's. Assumed, not determined: merged intensities cannot see the hand at all — |F| is diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index bb2daba1e..367e85615 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -1479,6 +1479,8 @@ namespace { message.bkg_estimate = GetCBORFloat(value); else if (key == "spindle_blind_fraction") message.spindle_blind_fraction = GetCBORFloat(value); + else if (key == "spindle_lost_unique_fraction") + message.spindle_lost_unique_fraction = GetCBORFloat(value); else if (key == "indexing_rate") message.indexing_rate = GetCBORFloat(value); else if (key == "indexed_lattice_count") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 1ee9bed6b..9a86d4091 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -802,6 +802,7 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "indexing_rate", message.indexing_rate); CBOR_ENC(mapEncoder, "bkg_estimate", message.bkg_estimate); CBOR_ENC(mapEncoder, "spindle_blind_fraction", message.spindle_blind_fraction); + CBOR_ENC(mapEncoder, "spindle_lost_unique_fraction", message.spindle_lost_unique_fraction); CBOR_ENC(mapEncoder, "rotation_lattice_type", message.rotation_lattice_type); if (message.rotation_lattice.has_value()) diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index a8d0e7c24..8883c0df2 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -18,6 +18,8 @@ ADD_LIBRARY(Rugnux STATIC ResultReport.h SpotWidth.cpp SpotWidth.h + SpindleCuspLoss.cpp + SpindleCuspLoss.h WriteModel.cpp WriteModel.h ) diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index c77dffd06..15b9c4875 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -314,6 +314,14 @@ std::string RenderResultReport(const std::string &output_prefix, fmt::format("{:.1f}", *result.spindle_symmetry_axis_deg)); Key(os, "SPINDLE_SYMMETRY_AXIS_ORDER", std::to_string(result.spindle_symmetry_axis_order)); } + // The exact verdict on the mounting, which the angle above cannot give on its own: the + // fraction (0-1) of unique reflections, to this run's resolution limit, that the measured + // point group could not recover from the sweep's blind cone. 0.0000 means the mounting + // cost nothing. Machine-readable on purpose - a downstream pipeline decides on this + // number, not on the warning prose. + if (result.spindle_lost_unique_fraction.has_value()) + Key(os, "SPINDLE_LOST_UNIQUE_FRACTION", + fmt::format("{:.4f}", *result.spindle_lost_unique_fraction)); os << "\n SPACE_GROUP_ALTERNATIVES names every group these data cannot separate from the one\n" << " adopted (enantiomorphic partners, origin-ambiguous pairs, groups a gap in the data\n" diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 469233cb6..4b1a6fd9c 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -5,6 +5,7 @@ #include "Rugnux.h" #include "ModelValidation.h" #include "WriteModel.h" +#include "SpindleCuspLoss.h" #include "SpotWidth.h" #include @@ -4702,38 +4703,63 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // search itself made, so the text cannot claim "no twin law exists" on its own say-so. result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; - // Symmetry axis vs spindle. Reported always on rotation data (it is a property of how the - // crystal was mounted, which the user can change), warned about when the two nearly coincide. + // How the crystal sat on the spindle. Reported always on rotation data (it is a property of + // the mounting, which the user can change). The angle to the nearest symmetry axis is kept + // as a descriptive key, but the WARNING is decided by the exact orbit computation: the old + // any-proper-axis-within-15-deg rule was wrong in both directions - an aligned in-plane + // 2-fold of a dihedral group is repaired by the principal axis, a cubic group is never + // severe in any orientation, and a lone diad PERPENDICULAR to the spindle is severe while + // no axis is anywhere near it. if (experiment_.IsRotationIndexing() && twin_sg && end_msg.rotation_lattice.has_value()) { if (const auto gonio = experiment_.GetGoniometer()) { const auto closest = ClosestSymmetryAxisToSpindle(*twin_sg, *end_msg.rotation_lattice, gonio->GetAxis()); if (closest.has_value()) { - // Practical bound, not a derived one: the blind cone's half-angle is the maximum - // Bragg angle (~15 deg for 2 A data at 1 A), and measured on this battery a 13.6 deg - // case shows the loss while a 16.2 deg one is 99.7% complete. - constexpr double WARN_DEG = 15.0; const auto [angle, order] = *closest; result.spindle_symmetry_axis_deg = angle; result.spindle_symmetry_axis_order = order; - if (angle < WARN_DEG) { + stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at " + << std::fixed << std::setprecision(1) << angle << " deg\n"; + } + + // The measured resolution of THIS run sets the cone: past merging the question is + // what these data lost, not what the detector could have reached. + double d_min = 0; + for (const auto &r : sm.merged) + if (std::isfinite(r.d) && r.d > 0 && (d_min == 0 || r.d < d_min)) + d_min = r.d; + if (const auto loss = SpindleUnrepairedFraction( + *twin_sg, *end_msg.rotation_lattice, gonio->GetAxis(), + experiment_.GetWavelength_A(), d_min)) { + result.spindle_lost_unique_fraction = loss->lost_unique_fraction; + end_msg.spindle_lost_unique_fraction = + static_cast(loss->lost_unique_fraction); + stats_text << "Unique reflections the mounting made unmeasurable: " + << std::fixed << std::setprecision(2) + << 100.0 * loss->lost_unique_fraction << "% (of the " + << std::setprecision(2) << 100.0 * (1.0 - std::cos(loss->theta_max_deg * PI / 180.0)) + << "% a sweep's blind cone holds at " << std::setprecision(2) + << d_min << " A)\n\n"; + // Warn when the point group recovers less than half the cone - the same half + // that fixes the per-image trigger threshold. The loss is a coherent cap about + // the spindle direction, which costs a map more than the same percentage lost + // at random. + if (loss->cone_fraction >= 0.5) { const std::string msg = fmt::format( - "The crystal's {}-fold axis is only {:.1f} deg from the spindle. A rotation sweep " - "never records the reflections whose reciprocal vector lies within the Bragg angle " - "of the spindle, and symmetry normally supplies them from an equivalent elsewhere; " - "it cannot here, because that axis maps the blind region onto itself. Those " - "reflections stay missing however long the sweep runs. The loss is confined to " - "that cone rather than spread over the data, so overall completeness may still " - "look reasonable - check it near the spindle direction. A second sweep on a " - "different axis, or re-mounting, recovers them.", - order, angle); + "The mounting makes {:.1f}% of unique reflections unmeasurable: the " + "measured point group cannot map that part of the sweep's blind cone " + "onto measured territory, however long the sweep runs. The loss is a " + "coherent cap about the spindle direction rather than a scatter, so " + "overall completeness may still look reasonable - check it near the " + "spindle. A second sweep about a different axis, or re-mounting, " + "recovers it.", + 100.0 * loss->lost_unique_fraction); logger.Warning("{}", msg); stats_text << " !! " << msg << "\n\n"; result.warnings.push_back(msg); - } else { - stats_text << "Closest symmetry axis to the spindle: " << order << "-fold at " - << std::fixed << std::setprecision(1) << angle << " deg\n\n"; } + } else { + stats_text << "\n"; } } } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index bc56961a9..f1c43e3c5 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -324,6 +324,12 @@ struct ProcessResult { // warning. Absent on stills, or where no space group was determined. std::optional spindle_symmetry_axis_deg; int spindle_symmetry_axis_order = 0; + // Fraction of unique reflections (to the processing resolution limit) the mounting made + // unmeasurable: the part of the sweep's blind double cone that no operator of the measured + // point group maps onto measured territory, in the crystal's actual indexed orientation. + // Exact, not the per-image worst-case bound. 0 = symmetry (or the mounting) recovered + // everything a sweep can lose. + std::optional spindle_lost_unique_fraction; }; // How far a cell stands from the metric its space group requires. A group's own rotations must leave diff --git a/rugnux/SpindleCuspLoss.cpp b/rugnux/SpindleCuspLoss.cpp new file mode 100644 index 000000000..f3adcc027 --- /dev/null +++ b/rugnux/SpindleCuspLoss.cpp @@ -0,0 +1,111 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include "SpindleCuspLoss.h" +#include "../common/JFJochMath.h" + +std::optional SpindleUnrepairedFraction(const gemmi::SpaceGroup &sg, + const CrystalLattice &lattice, + const Coord &spindle, + double wavelength_A, + double d_min_A) { + const double spindle_len = spindle.Length(); + if (spindle_len < 1e-9 || wavelength_A <= 0 || d_min_A <= 0) + return {}; + const Coord axis = spindle * static_cast(1.0 / spindle_len); + + const double sin_tm = std::min(1.0, wavelength_A / (2.0 * d_min_A)); + const double cos_tm = std::sqrt(std::max(0.0, 1.0 - sin_tm * sin_tm)); + if (sin_tm <= 0) + return {}; + + // Proper rotations only, identity included. Friedel and the improper operators need no separate + // treatment: the blind double cone and the measured region are both inversion-symmetric, so -R + // lands a point inside or outside exactly as R does. + std::vector, 3>> rot; + for (const auto &op : sg.operations().derive_symmorphic().sym_ops) { + const auto &w = op.rot; + const double det = + static_cast(w[0][0]) * (w[1][1] * w[2][2] - w[1][2] * w[2][1]) + - static_cast(w[0][1]) * (w[1][0] * w[2][2] - w[1][2] * w[2][0]) + + static_cast(w[0][2]) * (w[1][0] * w[2][1] - w[1][1] * w[2][0]); + if (det <= 0) + continue; + if (op.rot == gemmi::Op::identity().rot) + continue; // handled exactly below, without the float round trip through the basis + std::array, 3> m{}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + m[i][j] = static_cast(w[i][j]) / gemmi::Op::DEN; + rot.push_back(m); + } + + const Coord a = lattice.Vec0(), b = lattice.Vec1(), c = lattice.Vec2(); + const Coord as = lattice.Astar(), bs = lattice.Bstar(), cs = lattice.Cstar(); + + // Deterministic spiral over ONE lobe of the blind cone; the other lobe contributes identically + // (the orbit of -p is the negated orbit of p, and only |cos| to the spindle enters). Only cone + // directions can be lost - the identity is in every orbit - so nothing outside is sampled. + constexpr int N_DIRECTIONS = 8192; + const double golden_angle = PI * (3.0 - std::sqrt(5.0)); + + // Any unit vector perpendicular to the spindle, to open the cap around it. + const Coord seed = std::fabs(axis.x) < 0.9f ? Coord(1, 0, 0) : Coord(0, 1, 0); + const Coord e1 = (axis % seed).Normalize(); + const Coord e2 = axis % e1; + + double lost_sum = 0, p1_sum = 0; + for (int i = 0; i < N_DIRECTIONS; i++) { + const double z = cos_tm + (1.0 - cos_tm) * (static_cast(i) + 0.5) / N_DIRECTIONS; + const double r = std::sqrt(std::max(0.0, 1.0 - z * z)); + const double phi = golden_angle * i; + const Coord p = axis * static_cast(z) + + e1 * static_cast(r * std::cos(phi)) + + e2 * static_cast(r * std::sin(phi)); + + // h = (a.p, b.p, c.p) are continuous Miller coordinates; an operator sends the reflection + // h to h' = h W (gemmi's row convention) and h' returns to Cartesian through the reciprocal + // basis. No matrix inversion: direct and reciprocal bases are mutually inverse. + const double h0 = a * p, h1 = b * p, h2 = c * p; + + // The orbit's largest folded angle to the spindle, as its sine (the fold to [0, 90] makes + // the sine monotone). A point at direction p and radius q is unmeasured iff every image of + // its orbit is inside the cone of ITS OWN shell, i.e. iff q > 2 sin(a_max) / lambda - so + // the direction loses the radial fraction 1 - (sin a_max / sin theta_max)^3 of its unique + // reflections, reciprocal volume being what unique reflections are proportional to. + // The identity's own angle comes straight from z, exactly; only the non-trivial + // operators go through the basis. + const double sin_id = std::sqrt(std::max(0.0, 1.0 - z * z)); + double max_sin = sin_id; + for (const auto &m : rot) { + const double g0 = h0 * m[0][0] + h1 * m[1][0] + h2 * m[2][0]; + const double g1 = h0 * m[0][1] + h1 * m[1][1] + h2 * m[2][1]; + const double g2 = h0 * m[0][2] + h1 * m[1][2] + h2 * m[2][2]; + const Coord s = as * static_cast(g0) + bs * static_cast(g1) + + cs * static_cast(g2); + const double len = s.Length(); + if (len < 1e-12) + continue; + const double cos_a = std::min(1.0, std::fabs(s * axis) / len); + max_sin = std::max(max_sin, std::sqrt(std::max(0.0, 1.0 - cos_a * cos_a))); + } + if (max_sin < sin_tm) { + const double x = max_sin / sin_tm; + lost_sum += 1.0 - x * x * x; + } + const double y = sin_id / sin_tm; + p1_sum += 1.0 - y * y * y; + } + + SpindleCuspLoss ret; + ret.theta_max_deg = std::asin(sin_tm) * 180.0 / PI; + // The double cone is (1 - cos theta_max) of the sphere, and the sample covers exactly the cone, + // so the cone average rescales by that solid angle to a fraction of ALL unique reflections. + ret.lost_unique_fraction = (1.0 - cos_tm) * lost_sum / N_DIRECTIONS; + ret.cone_fraction = p1_sum > 0 ? lost_sum / p1_sum : 0.0; + return ret; +} diff --git a/rugnux/SpindleCuspLoss.h b/rugnux/SpindleCuspLoss.h new file mode 100644 index 000000000..d5fd0452f --- /dev/null +++ b/rugnux/SpindleCuspLoss.h @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include + +#include "../common/Coord.h" +#include "../common/CrystalLattice.h" + +// The exact, orbit-based counterpart of the per-image spindle severity. A still has to assume the +// worst - that the nearest plausible lattice row is a lone 2-fold - but by the time a rotation run +// is merged the point group and the indexed orientation are KNOWN, so nothing needs guessing: apply +// the group's proper rotations to the sweep's blind double cone in the crystal's actual orientation +// and count what no operator maps onto measured territory. A unique reflection is lost exactly when +// its whole orbit stays inside the cone (Friedel never helps: the cone is double-sided, so the +// inversion maps it onto itself). +struct SpindleCuspLoss { + // Fraction of all unique reflections to d_min that the mounting made unmeasurable. Weighted by + // reciprocal volume: each shell has its own, narrower cone theta(d) = asin(lambda/2d), so a + // direction at angle a from the spindle is only blind past q = 2 sin(a) / lambda. + double lost_unique_fraction = 0; + // The same loss as a fraction of the blind cone's own unique content: 1 = the group repaired + // nothing (P1, or a lone diad on or perpendicular to the spindle), 0 = symmetry recovers the + // whole cone. This is what severity means independent of how wide the cone happens to be. + double cone_fraction = 0; + double theta_max_deg = 0; +}; + +std::optional SpindleUnrepairedFraction(const gemmi::SpaceGroup &sg, + const CrystalLattice &lattice, + const Coord &spindle, + double wavelength_A, + double d_min_A); diff --git a/tests/CBORTest.cpp b/tests/CBORTest.cpp index 638506425..0b76e1726 100644 --- a/tests/CBORTest.cpp +++ b/tests/CBORTest.cpp @@ -498,6 +498,7 @@ TEST_CASE("CBORSerialize_End", "[CBOR]") { .max_receiver_delay = 3456, .efficiency = 0.99, .spindle_blind_fraction = 0.31f, + .spindle_lost_unique_fraction = 0.021f, .end_date = "ccc", .run_name = "bla5", .run_number = 45676782, @@ -537,6 +538,7 @@ TEST_CASE("CBORSerialize_End", "[CBOR]") { CHECK(output_message.rotation_lattice->GetUnitCell().c == Catch::Approx(60.0)); REQUIRE(output_message.spindle_blind_fraction == message.spindle_blind_fraction); + REQUIRE(output_message.spindle_lost_unique_fraction == message.spindle_lost_unique_fraction); // The per-image vector holds NaN where a frame had no value (CANNOT SAY); the hole must // survive the round trip as a hole, not as a number. REQUIRE(output_message.v_spindle_blind_fraction.size() == 3); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 299a6a503..bb6e33adf 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -83,6 +83,7 @@ ADD_EXECUTABLE(jfjoch_test TimeTest.cpp RotationIndexerTest.cpp SpindleBlindFractionTest.cpp + SpindleCuspLossTest.cpp TopPixelsTest.cpp HKLKeyTest.cpp TCPImagePusherTest.cpp diff --git a/tests/SpindleCuspLossTest.cpp b/tests/SpindleCuspLossTest.cpp new file mode 100644 index 000000000..162476764 --- /dev/null +++ b/tests/SpindleCuspLossTest.cpp @@ -0,0 +1,93 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "../rugnux/SpindleCuspLoss.h" + +using Catch::Matchers::WithinAbs; + +namespace { + // lambda = 1 A against d_min = 1 / (2 sin 15 deg) puts theta_max at exactly 15 deg. + constexpr double WVL = 1.0; + const double D_MIN = 1.0 / (2.0 * std::sin(15.0 * M_PI / 180.0)); + + Coord Perpendicular(const Coord &v) { + const Coord seed = std::fabs(v.Normalize().x) < 0.9f ? Coord(1, 0, 0) : Coord(0, 1, 0); + return (v % seed).Normalize(); + } +} + +TEST_CASE("SpindleCuspLoss_Orbit", "[Rugnux][Spindle]") { + SECTION("P1 repairs nothing, anywhere") { + const CrystalLattice latt(UnitCell{40, 50, 60, 83, 95, 102}); + const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 1"), latt, + Coord(0.3f, -0.5f, 0.8f), WVL, D_MIN); + REQUIRE(r.has_value()); + CHECK_THAT(r->cone_fraction, WithinAbs(1.0, 1e-9)); + // Radially weighted double-cone content at theta_max = 15 deg (verified by Monte Carlo). + CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0203, 0.0005)); + CHECK_THAT(r->theta_max_deg, WithinAbs(15.0, 1e-6)); + } + + SECTION("a lone diad on the spindle, and one perpendicular to it, both lose the whole cone") { + const CrystalLattice latt(UnitCell{40, 50, 60, 90, 100, 90}); // P2, unique axis b + const auto &sg = *gemmi::find_spacegroup_by_name("P 2"); + const Coord diad = latt.Vec1(); // b is perpendicular to a and c, so it IS the axis + const auto on = SpindleUnrepairedFraction(sg, latt, diad, WVL, D_MIN); + REQUIRE(on.has_value()); + CHECK_THAT(on->cone_fraction, WithinAbs(1.0, 1e-4)); + const auto perp = SpindleUnrepairedFraction(sg, latt, Perpendicular(diad), WVL, D_MIN); + REQUIRE(perp.has_value()); + CHECK_THAT(perp->cone_fraction, WithinAbs(1.0, 1e-4)); + } + + SECTION("the same diad at 60 deg from the spindle loses nothing") { + const CrystalLattice latt(UnitCell{40, 50, 60, 90, 100, 90}); + const Coord diad = latt.Vec1().Normalize(); + const Coord spindle = diad * 0.5f + Perpendicular(diad) * static_cast(std::sqrt(0.75)); + const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 2"), latt, + spindle, WVL, D_MIN); + REQUIRE(r.has_value()); + CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9)); + } + + SECTION("an axis of order >= 3 perpendicular to the spindle repairs the cone completely") { + // The situation the per-image worst-case bound must flag but the known group clears. + const CrystalLattice latt(UnitCell{50, 50, 60, 90, 90, 120}); + const Coord three_fold = latt.Vec2(); // c, the 3-fold of P3 + const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 3"), latt, + Perpendicular(three_fold), WVL, D_MIN); + REQUIRE(r.has_value()); + CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9)); + } + + SECTION("a dihedral group repairs an in-plane diad the warning heuristic used to flag") { + // 622 with an in-plane 2-fold exactly on the spindle: the old any-axis-within-15-deg rule + // warned here, but the principal 6-fold maps the cone off itself and nothing is lost. + const CrystalLattice latt(UnitCell{50, 50, 60, 90, 90, 120}); + const auto r = SpindleUnrepairedFraction(*gemmi::find_spacegroup_by_name("P 6 2 2"), latt, + latt.Vec0(), WVL, D_MIN); + REQUIRE(r.has_value()); + CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9)); + } + + SECTION("a cubic group is never severe, whatever the mounting") { + const CrystalLattice latt(UnitCell{60, 60, 60, 90, 90, 90}); + const auto &sg = *gemmi::find_spacegroup_by_name("P 4 3 2"); + for (const auto &spindle : {Coord(1, 0, 0), Coord(0.3f, -0.5f, 0.8f), Coord(1, 1, 1)}) { + const auto r = SpindleUnrepairedFraction(sg, latt, spindle, WVL, D_MIN); + REQUIRE(r.has_value()); + CHECK_THAT(r->lost_unique_fraction, WithinAbs(0.0, 1e-9)); + } + } + + SECTION("no spindle, no wavelength or no resolution gives no answer") { + const CrystalLattice latt(UnitCell{40, 50, 60, 90, 90, 90}); + const auto &sg = *gemmi::find_spacegroup_by_name("P 1"); + CHECK_FALSE(SpindleUnrepairedFraction(sg, latt, Coord(0, 0, 0), WVL, D_MIN).has_value()); + CHECK_FALSE(SpindleUnrepairedFraction(sg, latt, Coord(0, 0, 1), 0.0, D_MIN).has_value()); + CHECK_FALSE(SpindleUnrepairedFraction(sg, latt, Coord(0, 0, 1), WVL, 0.0).has_value()); + } +} diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index a0b1be62b..c37e745b0 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1182,6 +1182,10 @@ void NXmx::Finalize(const EndMessage &end) { if (end.spindle_blind_fraction) { SaveScalar(*hdf5_file, "/entry/MX/spindleBlindFractionMean", end.spindle_blind_fraction.value()); } + if (end.spindle_lost_unique_fraction) { + SaveScalar(*hdf5_file, "/entry/MX/spindleLostUniqueFraction", + end.spindle_lost_unique_fraction.value()); + } if (end.ice_ring_score_mean) { SaveScalar(*hdf5_file, "/entry/MX/iceRingScoreMean", end.ice_ring_score_mean.value()); }