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 "../common/Definitions.h"
#include "../common/JFJochMath.h"
#include "../image_analysis/geom_refinement/Calibrants.h"
#include "../rugnux/RugnuxCalibration.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_IceIsTheMeasuredRingList", "[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");
    // 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.
    CHECK(std::stod(keys["Poni1"]) == Catch::Approx(1275.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));
    // Only rot2 is NEGATED into pyFAI's frame; rot1 and rot3 are not. Pinned by feeding these exact
    // values to pyFAI and checking it reproduces the lab position DiffractionGeometry computes: the
    // mapping below agrees to 1.4e-17 m over the whole detector, with each angle tested separately,
    // while negating rot3 puts a pixel 25 mm out. Do not "fix" these signs without repeating that
    // check against pyFAI itself - its own documentation contradicts itself on the direction of its
    // axis 2, and a powder-ring check cannot test rot3, which is a rotation about the beam and 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));
    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);
}