Jungfraujoch/tests/CalibrationTest.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <catch2/catch_all.hpp>

#include <cstdio>
#include <fstream>
#include <map>
#include <sstream>
#include <algorithm>

#include "../common/Definitions.h"
#include "../common/JFJochMath.h"
#include <nlohmann/json.hpp>

#include "Calibration_output.h"
#include "Dataset_settings.h"
#include "../image_analysis/geom_refinement/AssignSpotsToRings.h"
#include "../image_analysis/geom_refinement/Calibrants.h"
#include "../image_analysis/geom_refinement/PowderCalibration.h"

TEST_CASE("Calibrants_LookupIsCaseInsensitive", "[DetGeomCalib]") {
    CHECK(CalibrantRings("LaB6") == CalibrantRings("lab6"));
    CHECK(CalibrantRings("AgBh") == CalibrantRings("agbh"));
    CHECK(CalibrantRings("nonsense").empty());
}

// GuessInitialGeometry pairs the innermost OBSERVED ring with the first entry of this list to fix the
// detector distance, so the first entry has to be a reflection that is really there. Only the primitive
// standard starts at (100): both face-centred ones extinguish it and start at (111), i.e. sqrt(3) times
// further out. Listing a forbidden ring first would scale every calibration by that ratio.
TEST_CASE("Calibrants_FirstRingIsThePresentOne", "[DetGeomCalib]") {
    struct Standard { std::string name; double a_A; double first_hkl_norm; };
    // sqrt(h^2+k^2+l^2) of the innermost present reflection: LaB6 Pm-3m -> 100, CeO2 Fm-3m and Si Fd-3m -> 111.
    const std::vector<Standard> standards = {
        {"lab6", LAB6_CELL_A, 1.0},
        {"ceo2", 5.4115,      std::sqrt(3.0)},
        {"si",   5.43102,     std::sqrt(3.0)}
    };
    for (const auto &s : standards) {
        const auto q = CalibrantRings(s.name);
        REQUIRE(!q.empty());
        CHECK(q.front() == Catch::Approx(2.0 * PI * s.first_hkl_norm / s.a_A).epsilon(1e-5));
    }
}

// The centring conditions themselves, checked ring by ring rather than only on the first one: an
// extinct reflection anywhere in the list mis-assigns the observed rings around it.
TEST_CASE("Calibrants_CentredStandardsOmitTheExtinctRings", "[DetGeomCalib]") {
    auto has_ring = [](const std::vector<float> &q, double a_A, double hkl_norm) {
        const auto want = static_cast<float>(2.0 * PI * hkl_norm / a_A);
        return std::any_of(q.begin(), q.end(),
                           [&](float v) { return std::fabs(v - want) < 1e-3f; });
    };

    const auto ceo2 = CalibrantRings("ceo2");
    CHECK(has_ring(ceo2, 5.4115, std::sqrt(3.0)));    // 111 - all odd
    CHECK(has_ring(ceo2, 5.4115, std::sqrt(4.0)));    // 200 - all even
    CHECK(has_ring(ceo2, 5.4115, std::sqrt(12.0)));   // 222 - all even, present without a glide plane
    CHECK_FALSE(has_ring(ceo2, 5.4115, 1.0));         // 100 - mixed parity
    CHECK_FALSE(has_ring(ceo2, 5.4115, std::sqrt(2.0)));  // 110 - mixed parity

    const auto si = CalibrantRings("si");
    CHECK(has_ring(si, 5.43102, std::sqrt(3.0)));     // 111 - all odd
    CHECK(has_ring(si, 5.43102, std::sqrt(8.0)));     // 220 - all even, h+k+l = 4n
    CHECK_FALSE(has_ring(si, 5.43102, 1.0));          // 100 - mixed parity
    CHECK_FALSE(has_ring(si, 5.43102, std::sqrt(4.0)));   // 200 - all even, h+k+l = 2
    CHECK_FALSE(has_ring(si, 5.43102, std::sqrt(12.0)));  // 222 - the diamond glide takes it out
}

// Ice is the reason the calibrant abstraction is a ring list and not a UnitCell: its entries are
// measured ring positions, and enumerating hkl from the hexagonal cell would add rings that are
// systematically absent in P6_3/mmc.
TEST_CASE("Calibrants_IceIsTheRingList", "[DetGeomCalib]") {
    const auto q = CalibrantRings("ice");
    REQUIRE(q.size() == ICE_RING_RES_A.size());
    CHECK(std::is_sorted(q.begin(), q.end()));
    CHECK(q.front() == Catch::Approx(2.0 * PI / ICE_RING_RES_A[0]).epsilon(1e-5));   // 3.895 A, the widest
}

// pyFAI's Poni1 is the SLOW axis (rows, our y) and Poni2 the FAST axis (columns, our x), both in
// metres. Transposing them produces a file that is silently wrong, so pin the mapping with a geometry
// whose two axes differ.
TEST_CASE("Calibration_PoniFileAxisConvention", "[DetGeomCalib]") {
    DiffractionExperiment x(DetJF4M());
    x.BeamX_pxl(1000.0f).BeamY_pxl(1275.0f).DetectorDistance_mm(150.0f);

    DiffractionGeometry geom = x.GetDiffractionGeometry();
    geom.PoniRot1_rad(0.01f).PoniRot2_rad(-0.02f).PoniRot3_rad(0.03f);

    const std::string path = "poni_test.poni";
    WritePoniFile(path, x, geom);

    std::map<std::string, std::string> keys;
    std::ifstream f(path);
    std::string line;
    while (std::getline(f, line)) {
        const auto colon = line.find(':');
        if (line.empty() || line[0] == '#' || colon == std::string::npos)
            continue;
        keys[line.substr(0, colon)] = line.substr(colon + 2);
    }
    f.close();
    std::remove(path.c_str());

    const double pixel_m = geom.GetPixelSize_mm() * 1e-3;
    CHECK(keys["poni_version"] == "2.1");
    // orientation 2 = "top left seen from the sample", the MX convention we assemble to. Without it
    // pyFAI applies its own default (3, bottom left) and gets the azimuth sense backwards.
    CHECK(keys["Detector_config"].find("\"orientation\": 2") != std::string::npos);
    // The half pixel is the origin convention (docs/DETECTOR_GEOMETRY.md): our beam centre is
    // pixel-centred, pyFAI measures from the edge of the sensor and puts the centre of pixel i at
    // (i + 0.5) * pixel size.
    // Declaring orientation 2 anchors Poni1 at the top edge, so the same physical point is
    // (height - 1 - beam_y) rows down from it.
    CHECK(std::stod(keys["Poni1"])
          == Catch::Approx((x.GetYPixelsNumConv() - 1 - 1275 + 0.5) * pixel_m));   // slow axis = y
    CHECK(std::stod(keys["Poni2"]) == Catch::Approx(1000.5 * pixel_m));   // fast axis = x
    CHECK(std::stod(keys["Distance"]) == Catch::Approx(0.150));
    // With orientation declared, (Rot1, Rot2, Rot3) = (+rot1, +rot2, -rot3 + pi): a row flip is
    // improper, so it reverses rotations about x and about the beam and leaves the one about the
    // vertical, and the half turn sets the azimuthal reference - pyFAI's in-plane axes are the
    // negatives of ours, so without it every chi is 180 degrees out. Being a rotation about the
    // beam it leaves 2theta alone, which is why radial integration was right while the azimuth was
    // not. Pinned against pyFAI 2026.5.0 on a tilted detector, against the lab positions of the
    // NXmx chain: 2theta to 3.6e-15 deg and chi to 2.8e-14 deg. Do not "fix" these without
    // repeating that check - a powder-ring test cannot see rot3, which moves only the azimuth.
    CHECK(std::stod(keys["Rot1"]) == Catch::Approx(0.01));
    CHECK(std::stod(keys["Rot2"]) == Catch::Approx(-0.02));
    CHECK(std::stod(keys["Rot3"]) == Catch::Approx(-0.03 + PI));
    CHECK(std::stod(keys["Wavelength"]) == Catch::Approx(geom.GetWavelength_A() * 1e-10));
    // max_shape is [rows, cols] - the same slow-then-fast order as Poni1/Poni2.
    const std::string shape = "[" + std::to_string(x.GetYPixelsNumConv()) + ", "
                            + std::to_string(x.GetXPixelsNumConv()) + "]";
    CHECK(keys["Detector_config"].find(shape) != std::string::npos);
}

// The match window may never reach the neighbouring ring, for any calibrant. Where two rings are closer
// together than twice the nominal window, a fixed window takes in the neighbour's flank - which the
// rings path reads as this ring's background, and which the spots path (before it took the NEAREST ring)
// resolved by assigning both to the lower-q one.
TEST_CASE("Calibration_RingMatchWindowNeverReachesTheNeighbour", "[DetGeomCalib]") {
    for (const auto &c : Calibrants()) {
        const auto q = CalibrantRings(c.name);
        REQUIRE(q.size() > 1);
        for (size_t i = 0; i < q.size(); ++i) {
            const float w = RingMatchWindow(q, i, RING_MATCH_Q_RECIPA);
            CHECK(w <= RING_MATCH_Q_RECIPA);
            CHECK(w > 0.0f);
            if (i > 0)
                CHECK(q[i] - w >= 0.5f * (q[i] + q[i - 1]) - 1e-6f);
            if (i + 1 < q.size())
                CHECK(q[i] + w <= 0.5f * (q[i] + q[i + 1]) + 1e-6f);
        }
    }
}

// ...and the clamp is not a no-op. Silver behenate's orders sit about 0.108 1/A apart and hexagonal ice
// has rings inside 0.06, so both are narrowed below the nominal window - while LaB6, whose rings are
// well separated at low q, keeps it. Without a standard that actually crowds, the test above would pass
// on a clamp that never fired.
TEST_CASE("Calibration_CrowdedStandardsNarrowTheWindow", "[DetGeomCalib]") {
    auto narrowed = [](const std::string &name) {
        const auto q = CalibrantRings(name);
        size_t n = 0;
        for (size_t i = 0; i < q.size(); ++i)
            if (RingMatchWindow(q, i, RING_MATCH_Q_RECIPA) < RING_MATCH_Q_RECIPA) ++n;
        return n;
    };
    CHECK(narrowed("agbh") > 0);
    CHECK(narrowed("ice") > 0);
    // The innermost LaB6 rings are more than 0.2 1/A apart, so nothing narrows them.
    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());
}

// The JSON a calibration writes has to BE a dataset_settings body, not merely resemble one. Its
// "dataset_settings" member is fed straight into the model generated from broker/jfjoch_api.yaml, which
// is the only definition of that schema this project has - so if the spec grows a field, renames one or
// changes what it requires, this fails rather than a beamline discovering it at a POST.
TEST_CASE("Calibration_JsonIsADatasetSettingsBody", "[DetGeomCalib]") {
    DiffractionExperiment x(DetJF4M());
    x.IncidentEnergy_keV(12.4f).BeamX_pxl(1000.0f).BeamY_pxl(1050.0f).DetectorDistance_mm(150.0f);

    CalibrationResult result;
    result.geometry = x.GetDiffractionGeometry();
    result.geometry.BeamX_pxl(1001.25f).BeamY_pxl(1049.5f).DetectorDistance_mm(151.5f)
                   .PoniRot1_rad(0.01f).PoniRot2_rad(-0.02f);
    result.ring_points = 321;
    result.rms_radial_pxl = 0.42;
    result.tilt_refined = true;
    result.tilt_significance = 17.5f;
    result.header_distance_mm = 150.0f;

    const std::string path = "calibration_json_test.json";
    WriteCalibrationJson(path, x, result, "lab6", "rings");
    std::ifstream in(path);
    REQUIRE(in.good());
    nlohmann::json j;
    in >> j;
    in.close();
    std::remove(path.c_str());

    REQUIRE(j.contains("dataset_settings"));
    const auto &settings = j.at("dataset_settings");

    // The WHOLE file is a calibration_output, not just its geometry member - so the quality block a
    // reader needs in order to tell a calibration that worked from one that did not is part of the
    // published contract too, and a python client can decode the file without knowing anything else.
    org::openapitools::server::model::Calibration_output output;
    REQUIRE_NOTHROW(from_json(j, output));
    std::stringstream output_msg;
    CHECK(output.validate(output_msg));
    CHECK(output.getCalibration().getRingPoints() == 321);
    CHECK(output.getCalibration().getMethod() == "rings");
    CHECK(output.getCalibration().isTiltRefined());

    // Every key is a property the schema knows, and the four it requires are all there.
    org::openapitools::server::model::Dataset_settings model;
    REQUIRE_NOTHROW(from_json(settings, model));
    std::stringstream msg;
    CHECK(model.validate(msg));

    CHECK(model.getBeamXPxl() == Catch::Approx(1001.25));
    CHECK(model.getBeamYPxl() == Catch::Approx(1049.5));
    CHECK(model.getDetectorDistanceMm() == Catch::Approx(151.5));
    CHECK(model.getIncidentEnergyKeV() == Catch::Approx(12.4));
    // beam_x_pxl is the PONI, so it must be the fitted PONI and NOT the direct beam - those differ by
    // distance*tan(tilt)/pixel here, and writing the wrong one would move a beamline's geometry.
    const auto [direct_x, direct_y] = result.geometry.GetDirectBeam_pxl();
    CHECK(std::abs(direct_x - model.getBeamXPxl()) > 0.5f);
    CHECK(j.at("calibration").at("direct_beam_x_pxl").get<double>() == Catch::Approx(direct_x).margin(0.01));

    // A tilted geometry carries all three rotations, because a body without them states a FLAT
    // detector rather than an unstated one.
    CHECK(model.getPoniRot1Rad() == Catch::Approx(0.01).margin(1e-6));
    CHECK(model.getPoniRot2Rad() == Catch::Approx(-0.02).margin(1e-6));
    CHECK(settings.contains("poni_rot3_rad"));
}

// ...and an untilted result leaves the rotations out altogether, which means the same thing: the API's
// own default for each is 0.0. The test is here so the two branches cannot drift apart.
TEST_CASE("Calibration_JsonOmitsTheRotationsWhenTheyAreZero", "[DetGeomCalib]") {
    DiffractionExperiment x(DetJF4M());
    x.IncidentEnergy_keV(12.4f).DetectorDistance_mm(150.0f);

    CalibrationResult result;
    result.geometry = x.GetDiffractionGeometry();
    result.geometry.PoniRot1_rad(0.0f).PoniRot2_rad(0.0f).PoniRot3_rad(0.0f);
    result.tilt_refined = false;

    const std::string path = "calibration_json_flat_test.json";
    WriteCalibrationJson(path, x, result, "lab6", "rings");
    std::ifstream in(path);
    nlohmann::json j;
    in >> j;
    in.close();
    std::remove(path.c_str());

    const auto &settings = j.at("dataset_settings");
    CHECK_FALSE(settings.contains("poni_rot1_rad"));
    CHECK_FALSE(settings.contains("poni_rot2_rad"));
    CHECK_FALSE(settings.contains("poni_rot3_rad"));
    org::openapitools::server::model::Dataset_settings model;
    REQUIRE_NOTHROW(from_json(settings, model));
    std::stringstream msg;
    CHECK(model.validate(msg));
}