leonarski_fandClaude Opus 5
5b8ce26c83
A reflection arriving at an angle to the detector normal crosses D/cos(alpha) of whatever lies between the sample and the sensor, not D, so it is attenuated more than one arriving head-on and reads low. That is the same geometry as the sensor crossing already corrected here and the opposite sign, and it was missing. The factor is exp(D/L*(1/cos(alpha)-1)) from the NIST attenuation coefficient of the medium, the stated distance and the stated wavelength. Nothing in it is fitted, and it is not justified by any measured amplitude: the flight path and the sensor crossing are collinear to better than 0.998 over the angular range any single experiment samples, so no fit of one can be evidence for the other. It is the tabulated absorption of a known thickness of a known material over a known path. The medium cannot be detected. No field of the NXmx application definition describes it, none of the masters this program reads carries one, and it cannot be inferred from the implied transmission either - in this corpus a station confirmed to use helium sits at 51% implied air transmission and one confirmed to use air at 63%, so any rule separating them is a threshold fitted between two points. It is therefore assumed, stated, and overridable: --flight-path air|helium|vacuum, defaulting to air. Helium is its own material rather than an alias for vacuum, attenuating about a six hundredth of air rather than nothing. On an untilted detector the correction is a function of resolution alone, so its entire effect on merged data is a shift in the Wilson B - which is what the report now prints beside the assumption, accurate to better than a tenth of an angstrom squared against measurement from 0.05 up to 28. Where that shift is large the report warns, because a wrong medium is then the largest number in the run: applied to data from the confirmed helium station it returns a B of 14 A^2 at 3.0 A resolution, which is not a value a crystal can have. The corpus contains its own control. One crystal, one station, three collections a quarter of an hour apart at falling energy through the same air: corrected, the Wilson B rises monotonically with the dose, as it must. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
579 lines
23 KiB
C++
579 lines
23 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 "../image_analysis/bragg_prediction/BraggPrediction.h"
|
|
#include <iostream>
|
|
#include "../image_analysis/SensorAbsorption.h"
|
|
|
|
// 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);
|
|
}
|
|
}
|
|
}
|
|
|
|
#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
|