leonarski_fandClaude Opus 5
659ba0d5b9
Nothing constrained the tilted geometry. Every existing test is either self-consistent or moves one angle at a time, and every file CI writes has zero tilt - where a swapped axis, a reversed composition order and a wrong pivot all give exactly the same answer. Two real bugs lived in that gap. Two checks, because there are two things to guard. DiffractionGeometry_Tilted_vs_PyFAI_and_DIALS pins the model itself: all three PONI angles non-zero and of mixed sign, compared per pixel against reference positions from pyFAI (an independent implementation of the convention) and from DIALS, to 2 um. Because the references are quoted in their own frames and the test applies the documented mappings - pyFAI (t1,t2,t3) -> our (x,y,z), and imgCIF = ours turned 180 degrees about x - it pins those relations too, not just the arithmetic. The comment gives the snippets to regenerate both sets. tests/nxmx_geometry_dials_test.py guards the writer, which is where the bugs actually were and which the unit test cannot reach. CI writes a master, patches the geometry in and asks DIALS where the panel is. Only the angle VALUES are patched; the axis vectors, the depends_on chain and the pivot stay as the writer emitted them, so they remain under test. Verified to fail on the pre-fix encoding: 7.3 mm, exit 1, naming the chain as the thing to look at. Recorded in both, because it cost an hour: compare via get_origin() and the fast/slow axes, NOT get_pixel_lab_coord(), which applies a parallax correction from the sensor thickness that Jungfraujoch does not model - about 0.1 mm at the detector edge, easily mistaken for a geometry error. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
634 lines
24 KiB
C++
634 lines
24 KiB
C++
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include <iostream>
|
|
#include "../common/DiffractionGeometry.h"
|
|
#include "../common/DiffractionExperiment.h"
|
|
|
|
TEST_CASE("RecipToDetector_1", "[LinearAlgebra][Coord]") {
|
|
DiffractionExperiment x(DetJF(8, 2));
|
|
x.BeamX_pxl(1024).BeamY_pxl(1024).DetectorDistance_mm(120);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
float pos_x = 512, pos_y = 512;
|
|
|
|
auto recip = geom.DetectorToRecip(pos_x, pos_y);
|
|
auto [proj_x, proj_y] = geom.RecipToDetector(recip);
|
|
REQUIRE(proj_x == Catch::Approx(pos_x));
|
|
REQUIRE(proj_y == Catch::Approx(pos_y));
|
|
REQUIRE((recip - geom.DetectorToRecip(proj_x, proj_y)).Length() < 0.00000001f);
|
|
REQUIRE(std::fabs(geom.DistFromEwaldSphere(recip)) < 4e-4);
|
|
}
|
|
|
|
TEST_CASE("RecipToDetector_2", "[LinearAlgebra][Coord]") {
|
|
DiffractionExperiment x(DetJF(8, 2));
|
|
x.BeamX_pxl(1024).BeamY_pxl(1024).DetectorDistance_mm(120);
|
|
float pos_x = 1023, pos_y = 1023;
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
auto recip = geom.DetectorToRecip(pos_x, pos_y);
|
|
auto [proj_x, proj_y] = geom.RecipToDetector(recip);
|
|
REQUIRE(proj_x == Catch::Approx(pos_x));
|
|
REQUIRE(proj_y == Catch::Approx(pos_y));
|
|
REQUIRE((recip - geom.DetectorToRecip(proj_x, proj_y)).Length() < 0.00000001f);
|
|
REQUIRE(std::fabs(geom.DistFromEwaldSphere(recip)) < 4e-4);
|
|
}
|
|
|
|
TEST_CASE("RecipToDetector_3", "[LinearAlgebra][Coord]") {
|
|
DiffractionExperiment x(DetJF(8, 2));
|
|
x.BeamX_pxl(1024).BeamY_pxl(1024).DetectorDistance_mm(120);
|
|
float pos_x = 30, pos_y = 30;
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
auto recip = geom.DetectorToRecip(pos_x, pos_y);
|
|
auto [proj_x, proj_y] = geom.RecipToDetector(recip);
|
|
REQUIRE(proj_x == Catch::Approx(pos_x));
|
|
REQUIRE(proj_y == Catch::Approx(pos_y));
|
|
REQUIRE((recip - geom.DetectorToRecip(proj_x, proj_y)).Length() < 0.00000001f);
|
|
REQUIRE(std::fabs(geom.DistFromEwaldSphere(recip)) < 4e-4);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_Phi","") {
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
x.BeamX_pxl(1000).BeamY_pxl(1000);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
CHECK(geom.Phi_rad(2000, 1000) * (180.0 / M_PI) == Catch::Approx(0.0));
|
|
CHECK(geom.Phi_rad(2000, 0) * (180.0 / M_PI) == Catch::Approx(315.0f));
|
|
CHECK(geom.Phi_rad(1000, 0) * (180.0 / M_PI) == Catch::Approx(270.0f));
|
|
CHECK(geom.Phi_rad(0, 0) * (180.0 / M_PI) == Catch::Approx(225.0f));
|
|
CHECK(geom.Phi_rad(0, 1000) * (180.0 / M_PI) == Catch::Approx(180.0f));
|
|
CHECK(geom.Phi_rad(1000, 2000) * (180.0 / M_PI) == Catch::Approx(90.f));
|
|
CHECK(geom.Phi_rad(2000, 2000) * (180.0 / M_PI) == Catch::Approx(45.0f));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_Cos2Theta","") {
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
x.BeamX_pxl(1000).BeamY_pxl(1000);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
// det distance == 1000 pixel
|
|
// theta = 30 deg
|
|
// tan(2 * theta) = sqrt(3)
|
|
REQUIRE(cosf(geom.TwoTheta_rad(1000, 1000 * (1.0 + sqrt(3)))) == Catch::Approx(0.5f));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PxlToRes","") {
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
// sin(theta) = 1/2
|
|
// theta = 30 deg
|
|
// tan(2 * theta) = sqrt(3)
|
|
REQUIRE(geom.PxlToRes( 0, 1000 * sqrt(3)) == Catch::Approx(1.0));
|
|
|
|
// sin(theta) = 1/4
|
|
// theta = 14.47 deg
|
|
// tan(2 * theta) = 0.55328333517
|
|
|
|
REQUIRE(geom.PxlToRes(1000 * 0.55328333517 * cosf(1), 1000 * 0.55328333517 * sinf(1)) == Catch::Approx(2.0));
|
|
}
|
|
|
|
|
|
TEST_CASE("DiffractionGeometry_ResToPxl","") {
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
// sin(theta) = 1/2
|
|
// theta = 30 deg
|
|
// tan(2 * theta) = sqrt(3)
|
|
REQUIRE(geom.ResToPxl(1.0) == Catch::Approx(1000 * sqrt(3)));
|
|
|
|
// sin(theta) = 1/4
|
|
// theta = 14.47 deg
|
|
// tan(2 * theta) = 0.55328333517
|
|
REQUIRE(geom.ResToPxl(2.0) == Catch::Approx(1000 * 0.55328333517));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_SolidAngleCorrection","") {
|
|
DiffractionExperiment x;
|
|
x.IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
x.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(75);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
// At the beam centre the correction is 1
|
|
REQUIRE(geom.CalcAzIntSolidAngleCorr(1000, 1000) == 1.0f);
|
|
|
|
// 2 * theta = 60 deg -> cos(2 * theta) = 1/2 -> correction = (1/2)^3
|
|
REQUIRE(geom.CalcAzIntSolidAngleCorr(1000 * (1.0 + sqrt(3)), 1000) == Catch::Approx(0.5f * 0.5f * 0.5f));
|
|
REQUIRE(geom.CalcAzIntSolidAngleCorr(1000, 1000 * (1.0 + sqrt(3))) == Catch::Approx(0.5f * 0.5f * 0.5f));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_SolidAngleCorrection_TiltInvariant","") {
|
|
// The solid-angle correction depends on the incidence angle to the detector
|
|
// normal, so for a given pixel it must be invariant under a rigid detector tilt
|
|
// (rot1/rot2/rot3) -- the same behaviour as PyFAI solidAngleArray.
|
|
DiffractionExperiment x;
|
|
x.IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
x.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(75);
|
|
DiffractionGeometry flat = x.GetDiffractionGeometry();
|
|
|
|
x.PoniRot1_rad(0.2).PoniRot2_rad(-0.1).PoniRot3_rad(0.5);
|
|
DiffractionGeometry tilted = x.GetDiffractionGeometry();
|
|
|
|
CHECK(tilted.CalcAzIntSolidAngleCorr(100, 100) == Catch::Approx(flat.CalcAzIntSolidAngleCorr(100, 100)));
|
|
CHECK(tilted.CalcAzIntSolidAngleCorr(1500, 400) == Catch::Approx(flat.CalcAzIntSolidAngleCorr(1500, 400)));
|
|
CHECK(tilted.CalcAzIntSolidAngleCorr(800, 1900) == Catch::Approx(flat.CalcAzIntSolidAngleCorr(800, 1900)));
|
|
CHECK(tilted.CalcAzIntSolidAngleCorr(1000, 1000) == Catch::Approx(flat.CalcAzIntSolidAngleCorr(1000, 1000)));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PolarizationCorrection","") {
|
|
DiffractionExperiment x;
|
|
x.IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
x.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(75);
|
|
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
// Circular polarization 0.5*(1+cos(2theta)^2)
|
|
x.PolarizationFactor(0);
|
|
REQUIRE(geom.CalcAzIntPolarizationCorr(1000 * (1.0 + sqrt(3)), 1000, 0) == Catch::Approx(0.5f * (1 + 0.5f * 0.5f)));
|
|
REQUIRE(geom.CalcAzIntPolarizationCorr(1000, 1000 * (1.0 + sqrt(3)), 0) == Catch::Approx(0.5f * (1 + 0.5f * 0.5f)));
|
|
|
|
// Horizontal polarization
|
|
x.PolarizationFactor(1);
|
|
|
|
// No correction in vertical direction
|
|
REQUIRE(geom.CalcAzIntPolarizationCorr(1000, 1000 * (1.0 + sqrt(3)), 1) == Catch::Approx(1.0f));
|
|
REQUIRE(geom.CalcAzIntPolarizationCorr(1000, 1000 * (1.0 - sqrt(3)), 1) == Catch::Approx(1.0f));
|
|
|
|
// cos(2*theta)^2 in horizontal direction
|
|
REQUIRE(geom.CalcAzIntPolarizationCorr(1000 * (1.0 + sqrt(3)), 1000, 1) == Catch::Approx(0.5f * 0.5f));
|
|
REQUIRE(geom.CalcAzIntPolarizationCorr(1000 * (1.0 - sqrt(3)), 1000, 1) == Catch::Approx(0.5f * 0.5f));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_AngleFromEwaldSphere") {
|
|
DiffractionGeometry geom;
|
|
geom.Wavelength_A(1.0);
|
|
|
|
// Center of Ewald sphere == (0,0,-1)
|
|
|
|
// Points on Ewald sphere
|
|
REQUIRE(geom.AngleFromEwaldSphere_deg(Coord(1, 0, -1)) == 0.0f);
|
|
REQUIRE(geom.AngleFromEwaldSphere_deg(Coord(1.0f / sqrtf(2.0f), 1.0f / sqrtf(2.0f), -1)) == 0.0f);
|
|
|
|
REQUIRE(geom.AngleFromEwaldSphere_deg(Coord(1, 0, 1)) == Catch::Approx(90.0f));
|
|
REQUIRE(geom.AngleFromEwaldSphere_deg(Coord(-sqrtf(2.0f), 0, 0)) == Catch::Approx(45.0f));
|
|
REQUIRE(geom.AngleFromEwaldSphere_deg(Coord(-sqrtf(3.0f), 0, 0)) == Catch::Approx(60.0f));
|
|
|
|
float cos_1deg = cosf(1.0f * M_PI / 180.0f);
|
|
float sin_1deg = sinf(1.0f * M_PI / 180.0f);
|
|
|
|
REQUIRE(fabsf(geom.AngleFromEwaldSphere_deg((Coord(cos_1deg - sin_1deg, 0, -(cos_1deg + sin_1deg)))) - 1.0f) < 0.0005);
|
|
|
|
// Cannot be rotated to fit into the Ewald sphere
|
|
REQUIRE(isnanf(geom.AngleFromEwaldSphere_deg(Coord(0, 0, 1))));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_AngleFromEwaldSphere_Wvl2A") {
|
|
DiffractionGeometry geom;
|
|
geom.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(100).Wavelength_A(2.0);
|
|
|
|
CHECK(geom.AngleFromEwaldSphere_deg(geom.DetectorToRecip(300,300)) < 0.05f);
|
|
CHECK(geom.AngleFromEwaldSphere_deg(geom.DetectorToRecip(200,1700)) < 0.05f);
|
|
CHECK(geom.AngleFromEwaldSphere_deg(geom.DetectorToRecip(1200,1800)) < 0.05f);
|
|
CHECK(geom.AngleFromEwaldSphere_deg(geom.DetectorToRecip(1500,100)) < 0.05f);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_ProjectToEwaldSphere") {
|
|
DiffractionGeometry geom;
|
|
geom.BeamX_pxl(1000).BeamY_pxl(437).DetectorDistance_mm(100).Wavelength_A(2.0);
|
|
Coord p0 = geom.DetectorToRecip(300,300);
|
|
Coord p1 = geom.ProjectToEwaldSphere(p0);
|
|
REQUIRE(p0.x == Catch::Approx(p1.x));
|
|
REQUIRE(p0.y == Catch::Approx(p1.y));
|
|
REQUIRE(p0.z == Catch::Approx(p1.z));
|
|
|
|
Coord p2 = Coord(1,0,0);
|
|
|
|
REQUIRE(std::fabs(geom.DistFromEwaldSphere(p2) > 0.01));
|
|
REQUIRE(std::fabs(geom.DistFromEwaldSphere(geom.ProjectToEwaldSphere(p2))) < 0.0001);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_DirectBeam") {
|
|
DiffractionGeometry geom;
|
|
geom.Wavelength_A(1.0);
|
|
geom.BeamX_pxl(1230).BeamY_pxl(1450);
|
|
|
|
auto [x, y] = geom.GetDirectBeam_pxl();
|
|
REQUIRE(x == Catch::Approx(1230.0f));
|
|
REQUIRE(y == Catch::Approx(1450.0f));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_DirectBeam_RotZ") {
|
|
DiffractionGeometry geom;
|
|
geom.Wavelength_A(1.0);
|
|
geom.BeamX_pxl(1230).BeamY_pxl(1450);
|
|
geom.PoniRot3_rad(-M_PI_2);
|
|
|
|
auto [x, y] = geom.GetDirectBeam_pxl();
|
|
REQUIRE(x == Catch::Approx(1230.0f));
|
|
REQUIRE(y == Catch::Approx(1450.0f));
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_DirectBeam_RotY") {
|
|
DiffractionGeometry geom;
|
|
geom.Wavelength_A(1.0);
|
|
geom.DetectorDistance_mm(100);
|
|
geom.PixelSize_mm(1.0);
|
|
geom.BeamX_pxl(1230).BeamY_pxl(1450);
|
|
geom.PoniRot2_rad(-M_PI_4); // 45 deg rotation
|
|
|
|
auto [x, y] = geom.GetDirectBeam_pxl();
|
|
CHECK(x == Catch::Approx(1230.0f)); // no Change for X
|
|
CHECK(y >1450.0f);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 1.0
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.0
|
|
Rot2: 0.0
|
|
Rot3: 0.0
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI uses nm^-1 for Q?
|
|
// The beam centre is Poni/pixel_size - 0.5 in every PONI test here: our coordinates are pixel-centred
|
|
// (0.0 is the centre of the first pixel) while pyFAI measures from the edge of the sensor and puts the
|
|
// centre of pixel i at (i + 0.5) * pixel size - see docs/DETECTOR_GEOMETRY.md. So 0.150 m / 75 um gives
|
|
// 1999.5, not 2000. With the half pixel the reference values below are reproduced to float precision;
|
|
// without it every one of them is out by 2.6e-3 nm^-1, which the old 1e-2 tolerance hid.
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
float diff_800_400 = fabs(geom.PxlToQ( 800,400)*10.0 - 6.295358803860941);
|
|
float diff_400_800 = fabs(geom.PxlToQ( 400,800)*10.0 - 7.554628215027982);
|
|
float diff_1300_2000 = fabs(geom.PxlToQ( 1300,2000)*10.0 - 5.73479724964891);
|
|
|
|
REQUIRE(diff_800_400 < 1e-4);
|
|
REQUIRE(diff_400_800 < 1e-4);
|
|
REQUIRE(diff_1300_2000 < 1e-4);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI_phi","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 1.0
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.0
|
|
Rot2: 0.0
|
|
Rot3: 0.0
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI uses nm^-1 for Q?
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
float phi_2000_0 = fabs(geom.Phi_rad(2000,0) - 2 * M_PI + 1.5702959937284997);
|
|
float phi_2000_2000 = fabs(geom.Phi_rad(2000,2000) - 1.5702964938446844);
|
|
float phi_0_1000 = fabs(geom.Phi_rad(0,1000) - 3.1413425992666903);
|
|
float phi_2000_1300 = fabs(geom.Phi_rad(1300,2000) - 2.1809518509415025);
|
|
|
|
CHECK(phi_2000_0 < 1e-4);
|
|
CHECK(phi_2000_2000 < 1e-4);
|
|
CHECK(phi_0_1000 < 1e-4);
|
|
CHECK(phi_2000_1300 < 1e-4);
|
|
}
|
|
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI_phi_rot3","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 1.0
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.0
|
|
Rot2: 0.0
|
|
Rot3: 0.5
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI uses nm^-1 for Q?
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.PoniRot3_rad(0.5);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
REQUIRE(geom.GetPoniRot3_rad() == Catch::Approx(0.5f));
|
|
|
|
float phi_800_400 = fabs(geom.Phi_rad(800,400) - 3.105073518019684);
|
|
float phi_2000_1300 = fabs(geom.Phi_rad(1300,2000) - 1.6809518509415027);
|
|
|
|
CHECK(phi_800_400 < 1e-4);
|
|
CHECK(phi_2000_1300 < 1e-4);
|
|
}
|
|
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI_phi_rot1_rot2_rot3","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 1.0
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.2
|
|
Rot2: 0.1
|
|
Rot3: 0.5
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI uses nm^-1 for Q?
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV)
|
|
.PoniRot1_rad(0.2).PoniRot2_rad(-0.1).PoniRot3_rad(0.5);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
REQUIRE(geom.GetPoniRot1_rad() == Catch::Approx(0.2f));
|
|
REQUIRE(geom.GetPoniRot2_rad() == Catch::Approx(-0.1f));
|
|
REQUIRE(geom.GetPoniRot3_rad() == Catch::Approx(0.5f));
|
|
|
|
float phi_800_400 = fabs(geom.Phi_rad(800,400) - 2 * M_PI + 1.4175001633470816);
|
|
float phi_2000_1300 = fabs(geom.Phi_rad(1300,2000) - 2 * M_PI + 0.6630282166663707);
|
|
|
|
CHECK(phi_800_400 < 1e-4);
|
|
CHECK(phi_2000_1300 < 1e-4);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI_rot1","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 1.0
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.2
|
|
Rot2: 0.0
|
|
Rot3: 0.0
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI uses nm^-1 for Q?
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
geom.PoniRot1_rad(0.2);
|
|
|
|
float diff_800_400 = fabs(geom.PxlToQ( 800,400)*10.0 - 7.471276390173706);
|
|
float diff_400_800 = fabs(geom.PxlToQ( 400,800)*10.0 - 5.148411999405654);
|
|
float diff_1300_2000 = fabs(geom.PxlToQ( 1300,2000)*10.0 - 10.37635963741911);
|
|
|
|
CHECK(diff_800_400 < 1e-4);
|
|
CHECK(diff_400_800 < 1e-4);
|
|
CHECK(diff_1300_2000 < 1e-4);
|
|
}
|
|
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI_rot1_rot2","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 1.0
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.2
|
|
Rot2: 0.1
|
|
Rot3: 0.0
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI uses nm^-1 for Q?
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(1000).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
geom.PoniRot1_rad(0.2).PoniRot2_rad(-0.1);
|
|
|
|
float diff_800_400 = fabs(geom.PxlToQ( 800,400)*10.0 - 11.412737079654118);
|
|
float diff_400_800 = fabs(geom.PxlToQ( 400,800)*10.0 - 8.805012278158177);
|
|
float diff_1300_2000 = fabs(geom.PxlToQ( 1300,2000)*10.0 - 9.363455481328781);
|
|
|
|
CHECK(diff_800_400 < 1e-4);
|
|
CHECK(diff_400_800 < 1e-4);
|
|
CHECK(diff_1300_2000 < 1e-4);
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PyFAI_Solid_angle","") {
|
|
/*
|
|
poni_version: 2
|
|
Detector: Eiger4M
|
|
Detector_config: {}
|
|
Distance: 0.2
|
|
Poni1: 0.075
|
|
Poni2: 0.150
|
|
Rot1: 0.0
|
|
Rot2: 0.0
|
|
Rot3: 0.0
|
|
Wavelength: 1e-10
|
|
*/
|
|
|
|
// PyFAI solidAngleArray is computed from the incidence angle to the detector normal,
|
|
// so it is independent of the poni rotation (tilt). CalcAzIntSolidAngleCorr matches this;
|
|
// the invariance is checked in DiffractionGeometry_SolidAngleCorrection_TiltInvariant.
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(200).BeamX_pxl(1999.5).BeamY_pxl(999.5).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
float diff_100_100 = fabs(geom.CalcAzIntSolidAngleCorr( 100,100) - 0.4844596502755233);
|
|
CHECK(diff_100_100 < 1e-5);
|
|
|
|
float diff_400_800 = fabs(geom.CalcAzIntSolidAngleCorr( 400,800)- 0.6267921080721112);
|
|
CHECK(diff_400_800 < 1e-5);
|
|
}
|
|
|
|
TEST_CASE("ResPhiToPxl") {
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
|
|
auto out = geom.ResPhiToPxl(1.0, 0);
|
|
CHECK(geom.PxlToRes(out.first, out.second) == Catch::Approx(1.0));
|
|
CHECK(fabs(geom.Phi_rad(out.first, out.second)) < 0.001 );
|
|
|
|
out = geom.ResPhiToPxl(1.0, M_PI);
|
|
CHECK(geom.PxlToRes(out.first, out.second) == Catch::Approx(1.0));
|
|
CHECK(fabs(geom.Phi_rad(out.first, out.second) - M_PI) < 0.001 );
|
|
|
|
out = geom.ResPhiToPxl(2.0, 0.7567);
|
|
CHECK(geom.PxlToRes(out.first, out.second) == Catch::Approx(2.0));
|
|
CHECK(fabs(geom.Phi_rad(out.first, out.second) - 0.7567) < 0.001 );
|
|
}
|
|
|
|
TEST_CASE("ResPhiToPxl_poni_rot") {
|
|
DiffractionExperiment x(DetJF4M());
|
|
x.DetectorDistance_mm(75).IncidentEnergy_keV(WVL_1A_IN_KEV);
|
|
DiffractionGeometry geom = x.GetDiffractionGeometry();
|
|
geom.PoniRot3_rad(0.5).PoniRot2_rad(-0.1).PoniRot2_rad(0.3);
|
|
|
|
auto out = geom.ResPhiToPxl(1.0, 0);
|
|
CHECK(geom.PxlToRes(out.first, out.second) == Catch::Approx(1.0));
|
|
CHECK(fabs(geom.Phi_rad(out.first, out.second)) < 0.001 );
|
|
|
|
out = geom.ResPhiToPxl(1.0, M_PI);
|
|
CHECK(geom.PxlToRes(out.first, out.second) == Catch::Approx(1.0));
|
|
CHECK(fabs(geom.Phi_rad(out.first, out.second) - M_PI) < 0.001 );
|
|
|
|
out = geom.ResPhiToPxl(2.0, 0.7567);
|
|
CHECK(geom.PxlToRes(out.first, out.second) == Catch::Approx(2.0));
|
|
CHECK(fabs(geom.Phi_rad(out.first, out.second) - 0.7567) < 0.001 );
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_DetectorToRecip_RecipToDetector_tilted") {
|
|
// Verify roundtrip consistency with non-zero rot1/rot2
|
|
DiffractionGeometry geom;
|
|
geom.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150)
|
|
.PixelSize_mm(0.075).Wavelength_A(1.0)
|
|
.PoniRot1_rad(0.05).PoniRot2_rad(-0.03);
|
|
|
|
// Test multiple points across the detector
|
|
std::vector<std::pair<float, float>> test_points = {
|
|
{500, 500}, {1500, 500}, {500, 1500}, {1500, 1500},
|
|
{800, 1200}, {1200, 800}, {300, 1700}, {1700, 300}
|
|
};
|
|
|
|
for (const auto& [x, y] : test_points) {
|
|
Coord recip = geom.DetectorToRecip(x, y);
|
|
auto [proj_x, proj_y] = geom.RecipToDetector(recip);
|
|
|
|
CHECK(proj_x == Catch::Approx(x).margin(0.001));
|
|
CHECK(proj_y == Catch::Approx(y).margin(0.001));
|
|
}
|
|
}
|
|
|
|
TEST_CASE("DiffractionGeometry_PONI_matrix_consistency") {
|
|
// Verify that the PONI rotation matrix gives consistent results
|
|
// when used for both forward and inverse transformations
|
|
DiffractionGeometry geom;
|
|
geom.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(100)
|
|
.PixelSize_mm(0.075).Wavelength_A(1.0)
|
|
.PoniRot1_rad(0.04).PoniRot2_rad(-0.025);
|
|
|
|
const auto& poni_rot = geom.GetPoniRotMatrix();
|
|
const auto poni_rot_T = poni_rot.transpose();
|
|
|
|
// Test: poni_rot * poni_rot^T should be identity (orthogonal matrix)
|
|
for (int i = 0; i < 3; ++i) {
|
|
for (int j = 0; j < 3; ++j) {
|
|
Coord ei, ej;
|
|
ei[i] = 1.0f;
|
|
ej[j] = 1.0f;
|
|
float expected = (i == j) ? 1.0f : 0.0f;
|
|
CHECK((poni_rot * (poni_rot_T * ej))[i] == Catch::Approx(expected).margin(1e-6));
|
|
}
|
|
}
|
|
|
|
// Test: S0 vector transformation
|
|
Coord S0 = geom.GetScatteringVector();
|
|
// For beam along z, S0 = (0, 0, 1/λ)
|
|
CHECK(S0.x == Catch::Approx(0.0f));
|
|
CHECK(S0.y == Catch::Approx(0.0f));
|
|
CHECK(S0.z == Catch::Approx(1.0f));
|
|
}
|
|
// Cross-check of a TILTED detector against two independent implementations, pyFAI and DIALS/dxtbx.
|
|
//
|
|
// Every other geometry test here is either self-consistent (round trips) or exercises one angle at a
|
|
// time. This one pins all three PONI angles at once, non-zero and of mixed sign, against reference
|
|
// positions computed outside Jungfraujoch. That matters because the errors this guards against are
|
|
// second order: a wrong composition order or a swapped axis is invisible unless two angles are
|
|
// non-zero simultaneously, and a wrong pivot is invisible to anything that only checks directions.
|
|
//
|
|
// HOW TO REGENERATE THE NUMBERS
|
|
//
|
|
// pyFAI (`pip install pyFAI`), which is an independent implementation of the PONI convention:
|
|
//
|
|
// from pyFAI.geometry import Geometry
|
|
// from pyFAI.detectors import Detector
|
|
// px = 75e-6
|
|
// det = Detector(pixel1=px, pixel2=px, max_shape=(2164, 2030))
|
|
// g = Geometry(dist=0.150, poni1=1275*px + px/2, poni2=1000*px + px/2,
|
|
// rot1=0.05, rot2=+0.03, rot3=0.02, # (+rot1, -rot2, +rot3); see WritePoniFile
|
|
// detector=det, wavelength=1e-10)
|
|
// t3, t1, t2 = g.calc_pos_zyx(d1=[y], d2=[x]) # metres, pyFAI's own axes
|
|
// lab_mm = (t2*1e3, t1*1e3, t3*1e3) # pyFAI (t1,t2,t3) -> our (x,y,z)
|
|
//
|
|
// The half pixel in poni1/poni2 is the origin convention (docs/DETECTOR_GEOMETRY.md): our beam
|
|
// centre is pixel-centred, pyFAI measures from the sensor edge.
|
|
//
|
|
// DIALS: write a master, patch this geometry into it, and read the panel back.
|
|
//
|
|
// source /opt/dials-v3-27-0/dials_env.sh
|
|
// build/tools/jfjoch_hdf5_test <input.h5> -n1 -S -o g # writes g_master.h5
|
|
// # with h5py, set /entry/instrument/detector/{beam_center_x,beam_center_y,distance},
|
|
// # transformations/{rot1,rot2,rot3}, and recompute transformations/translation - both its
|
|
// # magnitude and its @vector - as {bx*px, by*px, distance} normalised, since the writer
|
|
// # derives it from the beam centre and distance.
|
|
// p = ExperimentListFactory.from_filenames(['g_master.h5'])[0].detector[0]
|
|
// lab = p.get_origin() + x*px_mm*p.get_fast_axis() + y*px_mm*p.get_slow_axis()
|
|
//
|
|
// Use get_origin()/get_fast_axis()/get_slow_axis() as above, NOT get_pixel_lab_coord(): that applies
|
|
// a parallax correction from the sensor thickness and material which DiffractionGeometry does not
|
|
// model, and it costs ~0.1 mm at the detector edge - enough to look like a geometry error.
|
|
//
|
|
// DIALS reports in the imgCIF frame, which is ours turned 180 degrees about x (diag(1,-1,-1)) - a
|
|
// proper rotation, not a mirror. The test applies that mapping, so it pins the frame relation too.
|
|
TEST_CASE("DiffractionGeometry_Tilted_vs_PyFAI_and_DIALS", "[DiffractionGeometry]") {
|
|
DiffractionGeometry geom;
|
|
geom.BeamX_pxl(1000.0f).BeamY_pxl(1275.0f).DetectorDistance_mm(150.0f)
|
|
.PixelSize_mm(0.075f).Wavelength_A(1.0f)
|
|
.PoniRot1_rad(0.05f).PoniRot2_rad(-0.03f).PoniRot3_rad(0.02f);
|
|
|
|
struct Reference {
|
|
int x, y;
|
|
float pyfai[3]; // our frame: x, y, z [mm]
|
|
float dials[3]; // imgCIF frame: x, y, z [mm]
|
|
};
|
|
|
|
const std::vector<Reference> reference = {
|
|
{ 0, 0, {-69.399541334f, -98.819975327f, 150.623559791f},
|
|
{-69.399544982f, 98.819979800f, -150.623559832f}},
|
|
{ 2029, 2163, { 85.802269836f, 60.488193357f, 147.887421982f},
|
|
{ 85.802273560f, -60.488196450f, -147.887421893f}},
|
|
{ 1000, 1275, { 7.405508016f, -4.642730847f, 149.745128473f},
|
|
{ 7.405508016f, 4.642730847f, -149.745128473f}},
|
|
{ 300, 1800, {-44.232874953f, 35.676608756f, 153.548927011f},
|
|
{-44.232877406f, -35.676610671f, -153.548927191f}},
|
|
{ 1700, 400, { 58.519162201f, -71.195011374f, 145.153948055f},
|
|
{ 58.519164629f, 71.195014535f, -145.153947836f}},
|
|
};
|
|
|
|
// 2 um, i.e. 1/37 of a pixel. The references agree with each other to ~4e-6 mm; the margin is
|
|
// set by float32 rounding in DiffractionGeometry and in the HDF5 file the DIALS values came from.
|
|
const double margin = 2e-3;
|
|
|
|
for (const auto &r: reference) {
|
|
const Coord lab = geom.LabCoord(static_cast<float>(r.x), static_cast<float>(r.y));
|
|
|
|
CHECK(lab.x == Catch::Approx(r.pyfai[0]).margin(margin));
|
|
CHECK(lab.y == Catch::Approx(r.pyfai[1]).margin(margin));
|
|
CHECK(lab.z == Catch::Approx(r.pyfai[2]).margin(margin));
|
|
|
|
CHECK(lab.x == Catch::Approx( r.dials[0]).margin(margin));
|
|
CHECK(lab.y == Catch::Approx(-r.dials[1]).margin(margin));
|
|
CHECK(lab.z == Catch::Approx(-r.dials[2]).margin(margin));
|
|
}
|
|
}
|