mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 74 Wiki Activity
Files
7c10d62daba739ed6e06507fc795bc07df22b5fd
Jungfraujoch/tests/CalcBraggPredictionTest.cpp
T
leonarski_fandClaude Opus 5 7c10d62dab integration: the sensor efficiency is carried as its own quantity, not folded into the Lorentz-polarization factor 
It was multiplied into the per-reflection factor at prediction, so that factor held
Lorentz, polarization and efficiency at once and the two spellings that reach a file
- the wire key and the reflection dataset - meant something different from what
they had meant the day before. The unmerged MTZ had to divide the two apart again
at write time to fill its own columns, which is a good sign the wrong thing was
being carried.

Carry them separately. The prescaling factor is Lorentz and polarization again, what
its name and both reference implementations mean by it, and the efficiency is its own
field through prediction, integration, serialization and storage. Fifteen sites that
want the total now multiply the two - once per reflection, not once per pixel.

The efficiency is stored rather than recomputed on read, because the writer has no
geometry to recompute it from, and because a file written before the correction
existed would have had a radial trend invented for it. Sixty stored files were
checked for the one combination that would be ambiguous - the old meaning of the
factor beside a stored efficiency - and none carries it.

Output does not move. Re-scaling a file written before the efficiency existed is
byte-identical, which is a proof rather than a sample, since the stored factor is
exactly one there. Where the efficiency is live, one product is reassociated -
(L*Q)/P becomes (L/P)*Q - and about a third of the values differ in the last bit or
two: every structural column is identical, so no reflection is gained, lost or
reindexed, and no intensity in 1.4 million observations moves by as much as 1e-4 of
its own sigma.

The parity tests now compare the efficiency as well, and their non-vacuity guard
watches it rather than the factor it left - which is the same guard that went blind
when the efficiency was added to a field it was not watching.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
2026-09-05 16:28:15 +02:00

489 lines
18 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>
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;
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.image_scale_corr == Catch::Approx(gpu_r.image_scale_corr).epsilon(1e-4));
min_corr = std::min(min_corr, cpu_r.qe_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);
}
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;
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));
min_corr = std::min(min_corr, cpu_r.qe_corr);
matched++;
}
}
CHECK(matched > cpu_count * 0.95);
CHECK(min_corr < 0.99f);
}
#endif
Reference in New Issue View Git Blame Copy Permalink