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()); }