mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 80 Wiki Activity
Files
e9da892790c192a7deff4229106da0db6f2a14bc
Jungfraujoch/tests/CalcBraggPredictionTest.cpp
T
leonarski_fandClaude Opus 5.5 e7ef72ba40 Rotation prediction: place each partial at its slice of the rocking curve 
A rotation frame records the part of a reflection's rocking curve inside its own oscillation, and
while the crystal turns through the curve the spot walks along its Debye ring. The predictor put
every partial at the exact diffracting condition, so a partial recorded on the curve's flank was
integrated pixels away from where its flux landed. The walk is largest where the reflection moves
nearly tangent to the Ewald sphere (low |zeta|), whose curves are widest.

The prediction now turns S about the beam by the rotation's component along the ring times the
flux-weighted centre of the frame's slice of the curve (a truncated-normal mean, RockingSlice.h),
on the CPU and the GPU predictor alike. 2theta is unchanged; a frame that straddles the condition
symmetrically gets no shift.

Measured (observed r1 centroid minus prediction, tangential, I/sig > 10): against the slice model
correlation 0.93-0.98 before, 0.0 after; residual rms 1.0-1.8 px -> 0.35-0.48 px at |zeta| < 0.3.
Emulated capture of low-|zeta| partials at high angle 0.50 -> 0.94 of a 12 px aperture. CPU runs
on four rotation sets (400-800 frames): R_meas -0.03..-0.09 pp, ISa +1.2..+1.8, rejected
observations roughly halved, space groups unchanged.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
2026-09-24 14:06:24 +02:00

655 lines
28 KiB
C++
Raw Blame History

// 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 "../image_analysis/bragg_prediction/BraggPrediction.h"
#include <iostream>
#include "../image_analysis/SensorAbsorption.h"
#include "../image_analysis/bragg_prediction/BraggPredictionRot.h"
#include "../image_analysis/bragg_prediction/RockingSlice.h"
#include <map>
// The flight-path term is a zero-parameter prediction: a NIST attenuation coefficient, the stated
// sample-to-detector distance, and Beer-Lambert. Nothing about it is fitted, so it can be checked
// against the tables it is computed from rather than against any measurement.
TEST_CASE("BraggPrediction_FlightPathAttenuation", "[air_path]") {
using sensor_absorption::FlightPathAttenuation;
constexpr double D_mm = 160.0;
auto lambda_of = [](double keV) { return 12.39842 / keV; };
auto air = [&](double keV) {
return FlightPathAttenuation::Build(FlightPathMedium::Air, D_mm, lambda_of(keV));
};
// Attenuation length of dry air, as D/L at a known distance. These follow from the NIST dry-air
// mass attenuation coefficients and 1.205e-3 g/cm^3, and span four decades over the corpus's
// energy range: L is 8.2 m at 18 keV but only 8.9 cm at 3.8 keV.
CHECK(air(18.0).d_over_L == Catch::Approx(0.01959).epsilon(0.01));
CHECK(air(12.4).d_over_L == Catch::Approx(0.05351).epsilon(0.01));
CHECK(air(3.76).d_over_L == Catch::Approx(1.7972).epsilon(0.01));
// The argon K edge at 3.2029 keV is in the table and is a real step: argon is 1.28% of air by
// mass but dominates its absorption here, so interpolating straight across it would understate
// the correction in exactly the regime where the correction is largest.
CHECK(air(3.21).d_over_L > air(3.19).d_over_L * 1.05f);
const auto a18 = air(18.0);
// Normalised at normal incidence: head-on is untouched, by construction rather than by rounding.
CHECK(a18.Factor(1.0f) == 1.0f);
// Always >= 1 and monotone in the obliquity, because a reflection that arrives at a larger angle
// crossed strictly more of the medium.
float prev = 1.0f;
for (double alpha_deg : {10.0, 20.0, 30.0, 40.0, 50.0, 55.0}) {
const float f = a18.Factor(static_cast<float>(std::cos(alpha_deg * M_PI / 180.0)));
CHECK(f > prev);
prev = f;
}
const float cos55 = static_cast<float>(std::cos(55.0 * M_PI / 180.0));
// 18 keV over 160 mm at 55 degrees: +1.5%.
CHECK(a18.Factor(cos55) == Catch::Approx(1.0147).epsilon(0.002));
// It runs the other way to the sensor term at the same angle: the sensor makes an oblique
// reflection read high, the medium makes it read low, and neither is the other's undoing.
const auto qe = sensor_absorption::SensorQE::Build("Si", 450.0, lambda_of(18.0));
CHECK(qe.Factor(cos55) < 1.0f);
CHECK(a18.Factor(cos55) > 1.0f);
// HELIUM is why a long-wavelength station is usable at all: at 3.8 keV it attenuates some three
// orders of magnitude less than air, so the same geometry that costs a factor of several in air
// costs a fraction of a per cent in helium. It is NOT vacuum, and is not modelled as one.
const auto he = FlightPathAttenuation::Build(FlightPathMedium::Helium, D_mm, lambda_of(3.76));
CHECK(he.d_over_L > 0.0f);
CHECK(he.d_over_L < air(3.76).d_over_L / 100.0f);
CHECK(he.Factor(cos55) > 1.0f);
CHECK(he.Factor(cos55) < 1.01f);
// VACUUM leaves every intensity exactly untouched, at every angle - bit-identical to applying
// no correction at all.
const auto vac = FlightPathAttenuation::Build(FlightPathMedium::Vacuum, D_mm, lambda_of(3.76));
CHECK(vac.d_over_L == 0.0f);
CHECK(vac.Factor(1.0f) == 1.0f);
CHECK(vac.Factor(0.5f) == 1.0f);
CHECK(vac.Factor(0.1f) == 1.0f);
// Same for a file that states no distance or no wavelength.
CHECK(FlightPathAttenuation::Build(FlightPathMedium::Air, 0.0, lambda_of(12.4)).d_over_L == 0.0f);
CHECK(FlightPathAttenuation::Build(FlightPathMedium::Air, D_mm, 0.0).d_over_L == 0.0f);
// What the report quotes as the worth of the assumption. On an untilted detector the correction
// is a pure function of resolution, so its whole effect on merged data is a Wilson-B shift; the
// sign is negative because lifting the high-angle data flattens the fall-off. Measured on real
// data at three geometries, this predicts the observed shift to within 8%.
const double dB_hard = air(13.0).WilsonBShift_A2(1.28, 50.0, lambda_of(13.0));
const double dB_soft = air(3.76).WilsonBShift_A2(3.02, 50.0, lambda_of(3.76));
CHECK(dB_hard < 0.0);
CHECK(std::fabs(dB_hard) < 1.0); // hard X-rays: nothing a user could read off the data
CHECK(dB_soft < -5.0); // long wavelength: the dominant correction in the run
CHECK(std::fabs(dB_soft) > std::fabs(dB_hard) * 10.0);
// Vacuum is worth exactly nothing, which is what makes the report line meaningful.
CHECK(vac.WilsonBShift_A2(3.02, 50.0, lambda_of(3.76)) == 0.0);
}
TEST_CASE("BraggPrediction_11keV") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(11.0);
DiffractionGeometry geom = experiment.GetDiffractionGeometry();
CrystalLattice lattice(Coord{20, 10, 0}, Coord{-20, 40, 0},
Coord{0, 0, 200});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.001,
.max_h = 40, .max_k = 40, .max_l = 40
};
BraggPrediction prediction;
int count = prediction.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (int i = 0; i < count; i++) {
auto r = prediction.GetReflections().at(i);
auto recip = r.h * lattice.Astar() + r.k * lattice.Bstar() + r.l * lattice.Cstar();
REQUIRE(std::abs(r.h) < settings.max_h );
REQUIRE(std::abs(r.k) < settings.max_k );
REQUIRE(std::abs(r.l) < settings.max_l );
REQUIRE(r.d >= settings.high_res_A);
REQUIRE(r.d == Catch::Approx(1/std::sqrt(recip * recip)).margin(0.01f));
REQUIRE(r.dist_ewald == Catch::Approx(std::abs(geom.DistFromEwaldSphere(recip))).epsilon(1e-4));
auto [x,y] = geom.RecipToDetector(recip);
REQUIRE(r.predicted_x == Catch::Approx(x).margin(0.01));
REQUIRE(r.predicted_y == Catch::Approx(y).margin(0.01));
}
}
TEST_CASE("BraggPrediction_15keV") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(15.0);
DiffractionGeometry geom = experiment.GetDiffractionGeometry();
CrystalLattice lattice(Coord{20, 10, 0}, Coord{-20, 40, 0},
Coord{0, 0, 200});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.001,
.max_h = 40, .max_k = 40, .max_l = 40
};
BraggPrediction prediction;
int count = prediction.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (int i = 0; i < count; i++) {
auto r = prediction.GetReflections().at(i);
auto recip = r.h * lattice.Astar() + r.k * lattice.Bstar() + r.l * lattice.Cstar();
REQUIRE(std::abs(r.h) < settings.max_h );
REQUIRE(std::abs(r.k) < settings.max_k );
REQUIRE(std::abs(r.l) < settings.max_l );
REQUIRE(r.d >= settings.high_res_A);
REQUIRE(r.d == Catch::Approx(1/std::sqrt(recip * recip)).margin(0.01f));
REQUIRE(r.dist_ewald == Catch::Approx(std::abs(geom.DistFromEwaldSphere(recip))).epsilon(1e-3));
auto [x,y] = geom.RecipToDetector(recip);
REQUIRE(r.predicted_x == Catch::Approx(x).margin(0.01));
REQUIRE(r.predicted_y == Catch::Approx(y).margin(0.01));
}
}
TEST_CASE("BraggPrediction_Rot1_Rot2") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.PoniRot1_rad(2.0/180.0 * M_PI).PoniRot2_rad(3.0/180.0 * M_PI)
.IncidentEnergy_keV(11.0);
DiffractionGeometry geom = experiment.GetDiffractionGeometry();
CrystalLattice lattice(Coord{20, 10, 0}, Coord{-20, 40, 0},
Coord{0, 0, 200});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.001,
.max_h = 40, .max_k = 40, .max_l = 40
};
BraggPrediction prediction;
int count = prediction.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (int i = 0; i < count; i++) {
auto r = prediction.GetReflections().at(i);
auto recip = r.h * lattice.Astar() + r.k * lattice.Bstar() + r.l * lattice.Cstar();
REQUIRE(std::abs(r.h) < settings.max_h );
REQUIRE(std::abs(r.k) < settings.max_k );
REQUIRE(std::abs(r.l) < settings.max_l );
REQUIRE(r.d >= settings.high_res_A);
REQUIRE(r.d == Catch::Approx(1/std::sqrt(recip * recip)).margin(0.01f));
REQUIRE(r.dist_ewald == Catch::Approx(std::abs(geom.DistFromEwaldSphere(recip))).epsilon(1e-4));
auto [x,y] = geom.RecipToDetector(recip);
REQUIRE(r.predicted_x == Catch::Approx(x).margin(0.01));
REQUIRE(r.predicted_y == Catch::Approx(y).margin(0.01));
}
}
TEST_CASE("BraggPrediction_backscattering") {
DiffractionExperiment experiment(DetJF9M());
experiment.DetectorDistance_mm(120.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(3.0/WVL_1A_IN_KEV);
// Orthogonal basis, sufficient to test parity rules
CrystalLattice lattice(
Coord{40, 0, 0},
Coord{0, 50, 0},
Coord{0, 0, 60}
);
BraggPredictionSettings settings{
.high_res_A = 3.0f,
.ewald_dist_cutoff = 0.1f, // Very large cutoff, to be able to see as many reflections as possible
.max_h = 50, .max_k = 50, .max_l = 50
};
BraggPrediction pred;
int count = pred.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
REQUIRE(r.d > 3.0 / sqrt(2.0));
}
}
TEST_CASE("BraggPrediction_systematic_absences") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(120.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(12.0);
// Orthogonal basis, sufficient to test parity rules
CrystalLattice lattice(
Coord{40, 0, 0},
Coord{0, 50, 0},
Coord{0, 0, 60}
);
BraggPredictionSettings settings{
.high_res_A = 3.0f,
.ewald_dist_cutoff = 0.1f, // Very large cutoff, to be able to see as many reflections as possible
.max_h = 50, .max_k = 50, .max_l = 50
};
BraggPrediction pred;
SECTION("I centering") {
// 1) Body-centered I: reflections with h+k+l odd must be absent
settings.centering = 'I';
int count_I = pred.Calc(experiment, lattice, settings);
REQUIRE(count_I > 0);
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break; // ignore unfilled tail if any
REQUIRE(((r.h + r.k + r.l) % 2) == 0);
}
}
SECTION ("F centering") {
// 2) Face-centered F: h,k,l all even or all odd
settings.centering = 'F';
int count_F = pred.Calc(experiment, lattice, settings);
REQUIRE(count_F > 0);
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
const bool he = (r.h & 1) == 0, ke = (r.k & 1) == 0, le = (r.l & 1) == 0;
const bool all_even = he && ke && le;
const bool all_odd = (!he) && (!ke) && (!le);
REQUIRE((all_even || all_odd));
}
}
SECTION("R centering") {
settings.centering = 'R';
int count_R = pred.Calc(experiment, lattice, settings);
REQUIRE(count_R > 0);
// R (hexagonal setting): -h + k + l = 3n
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
int cond = (-r.h + r.k + r.l) % 3;
if (cond < 0) cond += 3;
REQUIRE(cond == 0);
}
}
}
TEST_CASE("RockingSliceCentroid_TruncatedNormalMean", "[rocking_slice]") {
// The closed form against a direct quadrature of the slice, on a frame that holds the centre, one on
// the curve's flank and one on its far tail.
const float sigma = 0.004f, half_wedge = 0.0015f, c1 = 1.0f / (std::sqrt(2.0f) * sigma);
for (float phi : {0.0f, 0.001f, -0.006f, 0.012f}) {
double sw = 0.0, stw = 0.0;
for (int i = 0; i <= 20000; ++i) {
const double t = phi - half_wedge + 2.0 * half_wedge * i / 20000.0;
const double w = std::exp(-t * t / (2.0 * sigma * sigma));
sw += w; stw += t * w;
}
const float partiality = (std::erf((phi + half_wedge) * c1) - std::erf((phi - half_wedge) * c1)) / 2.0f;
INFO("phi " << phi);
CHECK(RockingSliceCentroid_rad(phi, half_wedge, c1, partiality) == Catch::Approx(stw / sw).margin(2e-6));
}
CHECK(RockingSliceCentroid_rad(0.0f, half_wedge, c1, 0.3f) == 0.0f);
}
TEST_CASE("BraggPredictionRot_PartialWalksAlongItsRing", "[rocking_slice]") {
// Each rotation frame predicts a partial where the frame's slice of its rocking curve puts it: the
// exact-condition position walked along the Debye ring. Over the frames of one reflection the
// predicted positions therefore stay at one distance from the beam and move monotonically along
// the ring; a prediction at the exact condition would be the same point on every frame.
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1100.0).BeamY_pxl(1000.0).IncidentEnergy_keV(12.4);
const GoniometerAxis axis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {});
experiment.Goniometer(axis);
const CrystalLattice lattice(Coord{40, 0, 0}, Coord{0, 50, 0}, Coord{0, 0, 60});
BraggPredictionSettings settings{.high_res_A = 2.5, .ewald_dist_cutoff = 0.0015,
.max_h = 20, .max_k = 25, .max_l = 30,
.wedge_deg = 0.1f, .mosaicity_deg = 0.1f};
const float bx = 1100.0f, by = 1000.0f;
std::map<std::tuple<int, int, int>, std::vector<std::pair<float, float>>> track; // (radius, azimuth)
BraggPredictionRot pred;
for (int frame = 0; frame < 300; ++frame) {
const auto latt = lattice.Multiply(axis.GetTransformationAngle(frame * 0.1f));
const int n = pred.Calc(experiment, latt, settings);
for (int i = 0; i < n; ++i) {
const auto &r = pred.GetReflections().at(i);
if (r.zeta > 0.3f) continue;
const float dx = r.predicted_x - bx, dy = r.predicted_y - by;
track[{r.h, r.k, r.l}].emplace_back(std::hypot(dx, dy), std::atan2(dy, dx));
}
}
int tested = 0;
for (const auto &[hkl, pts] : track) {
if (pts.size() < 20) continue;
float rmin = pts[0].first, rmax = pts[0].first;
float along_first = 0.0f, along_last = 0.0f;
int sign_changes = 0;
float prev_step = 0.0f;
for (size_t j = 0; j < pts.size(); ++j) {
rmin = std::min(rmin, pts[j].first);
rmax = std::max(rmax, pts[j].first);
const float along = (pts[j].second - pts[0].second) * pts[0].first; // px along the ring
if (j == 0) along_first = along;
along_last = along;
if (j > 0) {
const float step = (pts[j].second - pts[j - 1].second);
if (step * prev_step < 0.0f) ++sign_changes;
if (step != 0.0f) prev_step = step;
}
}
INFO("hkl " << std::get<0>(hkl) << " " << std::get<1>(hkl) << " " << std::get<2>(hkl) << " frames " << pts.size());
CHECK(rmax - rmin < 0.05f); // stays on its ring
CHECK(std::fabs(along_last - along_first) > 0.1f); // but walks along it
CHECK(sign_changes == 0); // in one direction
++tested;
}
CHECK(tested > 20);
}
#ifdef JFJOCH_USE_CUDA
#include "../image_analysis/bragg_prediction/BraggPredictionGPU.h"
#include "../image_analysis/bragg_prediction/BraggPredictionRotGPU.h"
#include "../image_analysis/bragg_prediction/BraggPredictionRot.h"
#include <map>
#include <algorithm>
TEST_CASE("BraggPredictionGPU") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(13.0);
DiffractionGeometry geom = experiment.GetDiffractionGeometry();
CrystalLattice lattice(Coord{20, 10, 0}, Coord{-20, 40, 0},
Coord{0, 0, 200});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.001,
.max_h = 40, .max_k = 40, .max_l = 40
};
BraggPredictionGPU prediction;
int count = prediction.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (int i = 0; i < count; i++) {
auto r = prediction.GetReflections().at(i);
auto recip = r.h * lattice.Astar() + r.k * lattice.Bstar() + r.l * lattice.Cstar();
REQUIRE(std::abs(r.h) < settings.max_h );
REQUIRE(std::abs(r.k) < settings.max_k );
REQUIRE(std::abs(r.l) < settings.max_l );
REQUIRE(r.d >= settings.high_res_A);
REQUIRE(r.d == Catch::Approx(1/std::sqrt(recip * recip)).margin(0.01f));
REQUIRE(r.dist_ewald == Catch::Approx(std::abs(geom.DistFromEwaldSphere(recip))).epsilon(2e-2));
auto [x,y] = geom.RecipToDetector(recip);
REQUIRE(r.predicted_x == Catch::Approx(x).margin(0.05));
REQUIRE(r.predicted_y == Catch::Approx(y).margin(0.05));
}
}
TEST_CASE("BraggPredictionGPU_Rot1_Rot2") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.PoniRot1_rad(2.0/180.0 * M_PI).PoniRot2_rad(3.0/180.0 * M_PI)
.IncidentEnergy_keV(11.0);
DiffractionGeometry geom = experiment.GetDiffractionGeometry();
CrystalLattice lattice(Coord{20, 10, 0}, Coord{-20, 40, 0},
Coord{0, 0, 200});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.001,
.max_h = 40, .max_k = 40, .max_l = 40
};
BraggPredictionGPU prediction;
int count = prediction.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (int i = 0; i < count; i++) {
auto r = prediction.GetReflections().at(i);
auto recip = r.h * lattice.Astar() + r.k * lattice.Bstar() + r.l * lattice.Cstar();
REQUIRE(std::abs(r.h) < settings.max_h );
REQUIRE(std::abs(r.k) < settings.max_k );
REQUIRE(std::abs(r.l) < settings.max_l );
REQUIRE(r.d >= settings.high_res_A);
REQUIRE(r.d == Catch::Approx(1/std::sqrt(recip * recip)).margin(0.01f));
REQUIRE(r.dist_ewald == Catch::Approx(std::abs(geom.DistFromEwaldSphere(recip))).epsilon(1e-3));
auto [x,y] = geom.RecipToDetector(recip);
REQUIRE(r.predicted_x == Catch::Approx(x).margin(0.05));
REQUIRE(r.predicted_y == Catch::Approx(y).margin(0.05));
}
}
TEST_CASE("BraggPredictionGPU_systematic_absences") {
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(120.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(12.0);
CrystalLattice lattice(
Coord{40, 0, 0},
Coord{0, 50, 0},
Coord{0, 0, 60}
);
BraggPredictionSettings settings{
.high_res_A = 3.0f,
.ewald_dist_cutoff = 0.1f,
.max_h = 50, .max_k = 50, .max_l = 50
};
BraggPredictionGPU pred;
SECTION ("I centering") {
// 1) Body-centered I
settings.centering = 'I';
int count_I = pred.Calc(experiment, lattice, settings);
REQUIRE(count_I > 0);
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
REQUIRE(((r.h + r.k + r.l) % 2) == 0);
}
}
SECTION ("F centering") {
// 2) Face-centered F
settings.centering = 'F';
int count_F = pred.Calc(experiment, lattice, settings);
REQUIRE(count_F > 0);
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
const bool he = (r.h & 1) == 0, ke = (r.k & 1) == 0, le = (r.l & 1) == 0;
const bool all_even = he && ke && le;
const bool all_odd = (!he) && (!ke) && (!le);
REQUIRE((all_even || all_odd));
}
}
SECTION("R centering") {
settings.centering = 'R';
int count_R = pred.Calc(experiment, lattice, settings);
REQUIRE(count_R > 0);
// R (hexagonal setting): -h + k + l = 3n
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
int cond = (-r.h + r.k + r.l) % 3;
if (cond < 0) cond += 3;
REQUIRE(cond == 0);
}
}
}
TEST_CASE("BraggPredictionGPU_backscattering") {
DiffractionExperiment experiment(DetJF9M());
experiment.DetectorDistance_mm(120.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.IncidentEnergy_keV(3.0/WVL_1A_IN_KEV);
// Orthogonal basis, sufficient to test parity rules
CrystalLattice lattice(
Coord{40, 0, 0},
Coord{0, 50, 0},
Coord{0, 0, 60}
);
BraggPredictionSettings settings{
.high_res_A = 3.0f,
.ewald_dist_cutoff = 0.1f, // Very large cutoff, to be able to see as many reflections as possible
.max_h = 50, .max_k = 50, .max_l = 50
};
BraggPredictionGPU pred;
int count = pred.Calc(experiment, lattice, settings);
REQUIRE(count > 0);
for (const auto& r : pred.GetReflections()) {
if (r.d == 0) break;
REQUIRE(r.d > 3.0 / sqrt(2.0));
}
}
TEST_CASE("BraggPrediction_CPU_GPU_consistency_tilted") {
// Verify CPU and GPU implementations produce identical results with tilted detector
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.PoniRot1_rad(0.04).PoniRot2_rad(-0.025)
.IncidentEnergy_keV(12.0);
CrystalLattice lattice(Coord{30, 10, 0}, Coord{-15, 45, 0}, Coord{0, 0, 150});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.0015,
.max_h = 30, .max_k = 30, .max_l = 30
};
BraggPrediction cpu_pred;
BraggPredictionGPU gpu_pred;
int cpu_count = cpu_pred.Calc(experiment, lattice, settings);
int gpu_count = gpu_pred.Calc(experiment, lattice, settings);
REQUIRE(cpu_count > 0);
REQUIRE(gpu_count > 0);
// Build map of GPU reflections by hkl
std::map<std::tuple<int,int,int>, const Reflection*> gpu_refl_map;
for (int i = 0; i < gpu_count; ++i) {
const auto& r = gpu_pred.GetReflections().at(i);
gpu_refl_map[{r.h, r.k, r.l}] = &r;
}
// Check that each CPU reflection has a matching GPU reflection
int matched = 0;
float min_corr = 1.0f;
float max_flight = 1.0f;
for (int i = 0; i < cpu_count; ++i) {
const auto& cpu_r = cpu_pred.GetReflections().at(i);
auto key = std::make_tuple(cpu_r.h, cpu_r.k, cpu_r.l);
auto it = gpu_refl_map.find(key);
if (it != gpu_refl_map.end()) {
const auto& gpu_r = *it->second;
CHECK(cpu_r.predicted_x == Catch::Approx(gpu_r.predicted_x).margin(0.1));
CHECK(cpu_r.predicted_y == Catch::Approx(gpu_r.predicted_y).margin(0.1));
CHECK(cpu_r.d == Catch::Approx(gpu_r.d).margin(0.01));
// Both halves of the correction are part of the prediction, not a downstream product:
// the sensor-efficiency term lived on the CPU path alone for a while because nothing here
// compared it, and it is now a field of its own - so it is compared as a field of its own.
CHECK(cpu_r.prescaling_corr == Catch::Approx(gpu_r.prescaling_corr).epsilon(1e-4));
CHECK(cpu_r.qe_corr == Catch::Approx(gpu_r.qe_corr).epsilon(1e-4));
CHECK(cpu_r.flight_corr == Catch::Approx(gpu_r.flight_corr).epsilon(1e-4));
CHECK(cpu_r.image_scale_corr == Catch::Approx(gpu_r.image_scale_corr).epsilon(1e-4));
min_corr = std::min(min_corr, cpu_r.qe_corr);
max_flight = std::max(max_flight, cpu_r.flight_corr);
matched++;
}
}
// Most reflections should match (allow for some numerical differences at boundaries)
CHECK(matched > cpu_count * 0.95);
// ... and the comparison above must not be vacuous: on this geometry (320 um Si at 12 keV,
// reflections out to 2 A) the sensor correction is several per cent, so a qe_corr that stayed
// at 1 on both sides would mean the correction had been dropped from BOTH paths. It is checked
// on qe_corr and not on prescaling_corr, which no longer carries it.
CHECK(min_corr < 0.99f);
// The same for the air term, which runs the other way: at 12 keV over this distance the air
// correction reaches a few tenths of a per cent at the detector corner, so an flight_corr pinned at
// 1 on both sides would mean it had been dropped from BOTH paths rather than agreeing.
CHECK(max_flight > 1.0f);
}
TEST_CASE("BraggPredictionRot_CPU_GPU_consistency_tilted") {
// The rotation counterpart of the test above. Same purpose: every per-reflection quantity the
// two implementations both produce has to agree - the Lorentz-polarization factor and the sensor
// efficiency each in their own field, so neither can hide behind the other in a product.
DiffractionExperiment experiment(DetJF4M());
experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0)
.PoniRot1_rad(0.04).PoniRot2_rad(-0.025)
.IncidentEnergy_keV(12.0)
.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
CrystalLattice lattice(Coord{30, 10, 0}, Coord{-15, 45, 0}, Coord{0, 0, 150});
BraggPredictionSettings settings{
.high_res_A = 2.0,
.ewald_dist_cutoff = 0.0015,
.max_h = 30, .max_k = 30, .max_l = 30
};
BraggPredictionRot cpu_pred;
BraggPredictionRotGPU gpu_pred;
int cpu_count = cpu_pred.Calc(experiment, lattice, settings);
int gpu_count = gpu_pred.Calc(experiment, lattice, settings);
REQUIRE(cpu_count > 0);
REQUIRE(gpu_count > 0);
std::map<std::tuple<int,int,int>, const Reflection*> gpu_refl_map;
for (int i = 0; i < gpu_count; ++i) {
const auto& r = gpu_pred.GetReflections().at(i);
gpu_refl_map[{r.h, r.k, r.l}] = &r;
}
int matched = 0;
float min_corr = 1.0f;
float max_flight = 1.0f;
for (int i = 0; i < cpu_count; ++i) {
const auto& cpu_r = cpu_pred.GetReflections().at(i);
auto it = gpu_refl_map.find(std::make_tuple(cpu_r.h, cpu_r.k, cpu_r.l));
if (it != gpu_refl_map.end()) {
const auto& gpu_r = *it->second;
CHECK(cpu_r.predicted_x == Catch::Approx(gpu_r.predicted_x).margin(0.1));
CHECK(cpu_r.predicted_y == Catch::Approx(gpu_r.predicted_y).margin(0.1));
CHECK(cpu_r.d == Catch::Approx(gpu_r.d).margin(0.01));
CHECK(cpu_r.prescaling_corr == Catch::Approx(gpu_r.prescaling_corr).epsilon(1e-3));
CHECK(cpu_r.qe_corr == Catch::Approx(gpu_r.qe_corr).epsilon(1e-3));
CHECK(cpu_r.flight_corr == Catch::Approx(gpu_r.flight_corr).epsilon(1e-3));
min_corr = std::min(min_corr, cpu_r.qe_corr);
max_flight = std::max(max_flight, cpu_r.flight_corr);
matched++;
}
}
CHECK(matched > cpu_count * 0.95);
CHECK(min_corr < 0.99f);
CHECK(max_flight > 1.0f);
}
#endif
Reference in New Issue View Git Blame Copy Permalink