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

// The innermost ring of a cubic standard is its (100), so the first q is 2*pi/a. This is what fixes
// the distance in GuessInitialGeometry, so a wrong table would put every calibration off by that scale.
TEST_CASE("Calibrants_CubicStandardsHaveTheirLatticeConstant", "[DetGeomCalib]") {
    const std::map<std::string, double> a_A = {{"lab6", LAB6_CELL_A}, {"ceo2", 5.4115}, {"si", 5.43102}};
    for (const auto &[name, a] : a_A) {
        const auto q = CalibrantRings(name);
        REQUIRE(!q.empty());
        CHECK(q.front() == Catch::Approx(2.0 * PI / a).epsilon(1e-5));
    }
}

// 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);

    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));
    // rot2 and rot3 are NEGATED into pyFAI's frame and rot1 is not: pyFAI's slow axis runs bottom to
    // top where the MX convention runs top to bottom, so the frames differ by a reflection in y. That
    // reverses the sense of a rotation about x or about the beam, while for a rotation about y itself
    // the axis reverses too and the two cancel. Cross-checked against pyFAI on a real LaB6 image - the
    // unflipped file integrates rings broader than a zero-tilt one. Do not "fix" these signs to match
    // the stored values without repeating that check.
    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.0));
    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);
}