diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.cpp b/image_analysis/geom_refinement/AssignSpotsToRings.cpp index 2bc6a0c87..ef501323d 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.cpp +++ b/image_analysis/geom_refinement/AssignSpotsToRings.cpp @@ -186,7 +186,12 @@ namespace { } } -std::vector CalculateXtalRings(const UnitCell &cell, ReflectionCondition condition, int hkl_max) { +namespace { +// The ring enumeration itself. `present` decides which hkl the lattice actually diffracts into - a +// fixed centring condition for the built-in standards, the space group's own absences for a cell the +// user supplied. +template +std::vector XtalRings(const UnitCell &cell, int hkl_max, Present present) { CrystalLattice latt(cell); Coord Astar = latt.Astar(); @@ -201,7 +206,7 @@ std::vector CalculateXtalRings(const UnitCell &cell, ReflectionCondition for (int k = -hkl_max; k <= hkl_max; k++) { for (int l = 0; l <= hkl_max; l++) { if (h == 0 && k == 0 && l == 0) continue; - if (!reflection_present(condition, h, k, l)) continue; + if (!present(h, k, l)) continue; auto p = Astar * h + Bstar * k + Cstar * l; float Q = 2.0f * PI * p.Length(); u.push_back(Q); @@ -217,6 +222,20 @@ std::vector CalculateXtalRings(const UnitCell &cell, ReflectionCondition return u; } +} // namespace + +std::vector CalculateXtalRings(const UnitCell &cell, ReflectionCondition condition, int hkl_max) { + return XtalRings(cell, hkl_max, [condition](int h, int k, int l) { + return reflection_present(condition, h, k, l); + }); +} + +std::vector CalculateXtalRings(const UnitCell &cell, const gemmi::SpaceGroup &sg, int hkl_max) { + const gemmi::GroupOps ops = sg.operations(); + return XtalRings(cell, hkl_max, [&ops](int h, int k, int l) { + return !ops.is_systematically_absent(gemmi::Op::Miller{h, k, l}); + }); +} std::vector CalculateCubicXtalRings(float a, int hkl_max) { diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.h b/image_analysis/geom_refinement/AssignSpotsToRings.h index 69cd79ef0..360c5871e 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.h +++ b/image_analysis/geom_refinement/AssignSpotsToRings.h @@ -4,6 +4,8 @@ #pragma once #include +#include + #include "../../common/UnitCell.h" #include "../../common/SpotToSave.h" #include "RingOptimizer.h" @@ -45,6 +47,14 @@ enum class ReflectionCondition { std::vector CalculateXtalRings(const UnitCell &cell, ReflectionCondition condition = ReflectionCondition::All, int hkl_max = 6); + +// The same, for a cell whose symmetry is known: the absences come from the space group itself rather +// than from the three conditions above, so a centring, a glide and a screw are all handled. This is what +// a user-supplied calibrant cell needs - the fit pairs the innermost OBSERVED ring with the innermost +// LISTED one, so a list that opens with a forbidden reflection scales the whole calibration by the ratio +// between them. +std::vector CalculateXtalRings(const UnitCell &cell, const gemmi::SpaceGroup &sg, + int hkl_max = 6); std::vector CalculateCubicXtalRings( float a, int hkl_max = 4); float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_pxl, float d_A); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index dd3465172..7eb344e93 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -98,7 +98,7 @@ void print_usage() { std::cout << std::endl; std::cout << " Calibration (--mode calibration)" << std::endl; - std::cout << " --calibrant Powder standard: lab6|agbh|ceo2|si|ice (default: lab6, case-insensitive)" << std::endl; + std::cout << " --calibrant Powder standard: lab6|agbh|ceo2|si|ice (default: lab6, case-insensitive). Overridden by -C: a unit cell given there IS the standard, with its rings enumerated from the cell and its absences from -S where that is given, or taken as primitive where it is not. -S supplies SYMMETRY absences only: a standard whose extinctions come from its atomic basis rather than from a symmetry element - silicon, whose 222 is absent in Fd-3m for that reason - is better named with --calibrant" << std::endl; std::cout << " --calibration How the rings are measured: rings|spots (default: rings). rings sums the (q x azimuth) azimuthal profile over every processed image and fits the ring arcs in it; spots pools the found spots and fits those. -s/-e/--stride select the images; rings defaults --azim-phi-bins to 32" << std::endl; std::cout << " --no-refine-tilt Do not refine the detector tilt: hold rot1/rot2 at the header value and fit only the beam centre and the distance. For a calibration handed to a program that cannot express a tilted detector, e.g. XDS" << std::endl; std::cout << std::endl; @@ -1975,7 +1975,27 @@ static int RunRugnux(int argc, char **argv) { config.mode = ProcessMode::Calibration; config.calibration_method = calibration_method; config.calibration_refine_tilt = calibration_refine_tilt; - config.calibrant_ring_q = CalibrantRings(calibrant); + // -C wins over --calibrant: a cell given on the command line IS the standard, and the built-in + // table is only the convenience for the five that have names. -S, where it is given, supplies + // the absences - without it the cell is taken as primitive, which for a centred standard would + // open the ring list with a reflection that is not there and scale the whole calibration. + std::string calibrant_label = calibrant; + if (fixed_reference_unit_cell.has_value()) { + config.calibrant_ring_q = space_group != nullptr + ? CalculateXtalRings(*fixed_reference_unit_cell, *space_group) + : CalculateXtalRings(*fixed_reference_unit_cell); + calibrant_label = fmt::format("cell {:.4f},{:.4f},{:.4f},{:.3f},{:.3f},{:.3f} ({})", + fixed_reference_unit_cell->a, fixed_reference_unit_cell->b, + fixed_reference_unit_cell->c, fixed_reference_unit_cell->alpha, + fixed_reference_unit_cell->beta, fixed_reference_unit_cell->gamma, + space_group != nullptr ? space_group->hm : "assumed primitive"); + if (config.calibrant_ring_q.empty()) { + logger.Error("The unit cell given with -C produces no rings"); + return 1; + } + } else { + config.calibrant_ring_q = CalibrantRings(calibrant); + } config.start_image = start_image; config.end_image = end_image; config.stride = image_stride; @@ -2008,7 +2028,7 @@ static int RunRugnux(int argc, char **argv) { experiment.MaxSpotCount(max_spot_count_override.value_or(RUGNUX_MAX_SPOT_COUNT)); logger.Info("Powder calibration: calibrant {} ({} rings), method {}, {} azimuthal bins, tilt {}", - calibrant, config.calibrant_ring_q.size(), + calibrant_label, config.calibrant_ring_q.size(), calibration_method == CalibrationMethod::Rings ? "rings" : "spots", experiment.GetAzimuthalIntegrationSettings().GetAzimuthalBinCount(), calibration_refine_tilt ? "refined" : "fixed"); diff --git a/tests/CalibrationTest.cpp b/tests/CalibrationTest.cpp index 1f6bda6c4..b525b062c 100644 --- a/tests/CalibrationTest.cpp +++ b/tests/CalibrationTest.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include "../common/Definitions.h" #include "../common/JFJochMath.h" @@ -168,3 +169,62 @@ TEST_CASE("Calibration_CrowdedStandardsNarrowTheWindow", "[DetGeomCalib]") { const auto lab6 = CalibrantRings("lab6"); CHECK(RingMatchWindow(lab6, 0, RING_MATCH_Q_RECIPA) == Catch::Approx(RING_MATCH_Q_RECIPA)); } + +// A cell given with -C takes its absences from -S. That path is independent of the hand-written +// ReflectionConditions the built-in table uses, so the two must agree where the standard's absences are +// a property of its SYMMETRY - which is what says the gemmi route is safe to hand a user's cell. +TEST_CASE("Calibration_SpaceGroupAbsencesMatchTheBuiltInConditions", "[DetGeomCalib]") { + struct Standard { std::string name; UnitCell cell; std::string hm; }; + const std::vector standards = { + {"lab6", UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90, 90, 90), "P m -3 m"}, + {"ceo2", UnitCell(5.4115, 5.4115, 5.4115, 90, 90, 90), "F m -3 m"} + }; + for (const auto &s : standards) { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(s.hm); + REQUIRE(sg != nullptr); + const auto from_sg = CalculateXtalRings(s.cell, *sg); + const auto from_table = CalibrantRings(s.name); + REQUIRE(from_sg.size() == from_table.size()); + for (size_t i = 0; i < from_sg.size(); ++i) + CHECK(from_sg[i] == Catch::Approx(from_table[i]).epsilon(1e-6)); + } +} + +// Silicon is the case where they must NOT agree, and it is worth pinning because it bounds what -C -S +// can do. Fd-3m's symmetry absences are only the F centring; silicon's further extinctions - 222 is the +// first - come from its two-atom basis, i.e. from the structure factor and not from any symmetry +// element, so gemmi cannot know them and reports 24 rings where the table's diamond condition gives 18. +// The extra ones are exactly the all-even reflections with h+k+l not a multiple of 4. They do not move +// the FIRST ring, so the calibration is not scaled wholesale - but they are rings carrying no intensity +// offered to the matcher in the middle of the list, which is why --calibrant si still exists. +TEST_CASE("Calibration_SpaceGroupCannotKnowStructureFactorAbsences", "[DetGeomCalib]") { + const UnitCell si(5.43102, 5.43102, 5.43102, 90, 90, 90); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("F d -3 m"); + REQUIRE(sg != nullptr); + const auto from_sg = CalculateXtalRings(si, *sg); + const auto from_table = CalibrantRings("si"); + CHECK(from_sg.size() > from_table.size()); + // The table's rings are a subset of the space group's - nothing is LOST by asking gemmi, only added. + for (const float q : from_table) { + const bool present = std::any_of(from_sg.begin(), from_sg.end(), + [q](float r) { return std::fabs(r - q) < 1e-4f; }); + CHECK(present); + } + // ...and the first ring, the one the distance is seeded from, is the same either way. + CHECK(from_sg.front() == Catch::Approx(from_table.front()).epsilon(1e-6)); +} + +// Without -S the cell is taken as primitive, which for a centred standard is NOT the same list: the +// face-centred absences are what move the first ring from 100 out to 111. The point of the test is that +// the difference is real, so that "assumed primitive" in the log is a warning worth reading. +TEST_CASE("Calibration_PrimitiveAssumptionDiffersForACentredCell", "[DetGeomCalib]") { + const UnitCell ceo2(5.4115, 5.4115, 5.4115, 90, 90, 90); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("F m -3 m"); + REQUIRE(sg != nullptr); + const auto centred = CalculateXtalRings(ceo2, *sg); + const auto primitive = CalculateXtalRings(ceo2); + REQUIRE(!centred.empty()); + REQUIRE(!primitive.empty()); + CHECK(primitive.size() > centred.size()); + CHECK(centred.front() > primitive.front()); +}