calibration: -C overrides the calibrant, with absences from -S

The five named standards are a convenience, not the limit of what a powder
calibration can be run against. A unit cell given with -C now IS the standard
in --mode calibration, and its rings are enumerated from that cell.

A cell alone does not give a ring list, though: the centring and any glide
decide which hkl the lattice actually diffracts into, and the fit pairs the
innermost OBSERVED ring with the innermost LISTED one - so a list opening with
a reflection that is not there scales the whole calibration by the ratio
between them. Where -S is given, the absences come from the space group itself
via gemmi, which covers centring, glides and screws in one mechanism rather
than the three hand-written conditions the built-in table uses. The two agree
exactly on LaB6 and CeO2, which is the cross-check that says the gemmi route is
safe to hand a user's cell.

They do NOT agree on silicon, and the test now pins that: Fd-3m's symmetry
absences are only the F centring, while silicon's 222 and its relatives are
extinguished by its two-atom basis - a structure-factor absence, not a symmetry
one, so no symmetry handler can know it. gemmi offers 24 rings where the
diamond condition gives 18. The extra rings do not move the first one, so the
distance is not scaled, but they are rings with no intensity offered to the
matcher - which is why --calibrant si stays, and why the usage text says -S
supplies symmetry absences only.

Without -S the cell is taken as primitive and the log says "assumed primitive".

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27

image_analysis/geom_refinement/AssignSpotsToRings.cpp

@@ -186,7 +186,12 @@ namespace {
     }
 }
 
-std::vector<float> 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 <typename Present>
+std::vector<float> XtalRings(const UnitCell &cell, int hkl_max, Present present) {
     CrystalLattice latt(cell);
 
     Coord Astar = latt.Astar();
@@ -201,7 +206,7 @@ std::vector<float> 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<float> CalculateXtalRings(const UnitCell &cell, ReflectionCondition
 
     return u;
 }
+}   // namespace
+
+std::vector<float> 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<float> 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<float> CalculateCubicXtalRings(float a, int hkl_max) {

image_analysis/geom_refinement/AssignSpotsToRings.h

@@ -4,6 +4,8 @@
 #pragma once
 
 #include <vector>
+#include <gemmi/symmetry.hpp>
+
 #include "../../common/UnitCell.h"
 #include "../../common/SpotToSave.h"
 #include "RingOptimizer.h"
@@ -45,6 +47,14 @@ enum class ReflectionCondition {
 std::vector<float> 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<float> CalculateXtalRings(const UnitCell &cell, const gemmi::SpaceGroup &sg,
+                                      int hkl_max = 6);
 std::vector<float> CalculateCubicXtalRings( float a, int hkl_max = 4);
 float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_pxl, float d_A);
 

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 <name>               Powder standard: lab6|agbh|ceo2|si|ice (default: lab6, case-insensitive)" << std::endl;
+    std::cout << "   --calibrant <name>               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 <txt>              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");

tests/CalibrationTest.cpp

@@ -7,6 +7,7 @@
 #include <fstream>
 #include <map>
 #include <sstream>
+#include <algorithm>
 
 #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<Standard> 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());
+}