mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 66 Wiki Activity

azint: tilt-correct solid-angle correction; honour imported rot3 in refinement

Browse Source 
The azimuthal-integration solid-angle correction used cos^3(2*theta), where
2*theta is the true scattering angle (from LabCoord, including detector tilt).
The solid angle of a flat pixel actually depends on the incidence angle to the
detector normal, cos(alpha) = det_distance / |detector-frame position|, which is
invariant under detector tilt (rot1/rot2/rot3). Only for an untilted detector do
the two agree. Switch CalcAzIntSolidAngleCorr(x,y) to the tilt-invariant form,
matching PyFAI solidAngleArray and MAX IV azint. Drop the q-only overload (it can
only ever be the untilted approximation and was used only in tests) and move its
test onto the (x,y) form; add a tilt-invariance test.

XtalOptimizer's residual reconstructed each spot's lab position from rot1/rot2
only, hardcoding rot3 = 0, while the rest of the pipeline (and its own spot
selection) used the full PONI rotation. An imported non-zero rot3 was therefore
silently dropped during refinement. Bake rot3 into the residual as a fixed
Rz(-rot3) so refinement stays consistent (no-op when rot3 == 0).

Polarization and azimuthal binning already honoured rot3 via the full PONI
rotation (Phi_rad), validated against PyFAI chi() by the existing rot3 phi tests.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-11 11:34:11 +02:00
co-authored by Claude Opus 4.8
parent e7fbeb527f
commit acd2025676
5 changed files with 55 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/DiffractionGeometry.cpp
+9 -10
View File
@@ -87,17 +87,16 @@ float DiffractionGeometry::DistFromEwaldSphere(const Coord &recip) const {
return S.Length() - (1.0f/wavelength_A);
}
float DiffractionGeometry::CalcAzIntSolidAngleCorr(float q) const {
float sin_theta = q * wavelength_A / (4 * static_cast<float>(PI));
float cos_2theta = 1.0f - 2.0f * sin_theta * sin_theta; // cos(2*alpha) = 1 - 2 * sin(alpha)^2
float cos_2theta_3 = cos_2theta * cos_2theta * cos_2theta;
return cos_2theta_3;
}
float DiffractionGeometry::CalcAzIntSolidAngleCorr(float x, float y) const {
float cos_2theta = cosf(TwoTheta_rad(x, y));
float cos_2theta_3 = cos_2theta * cos_2theta * cos_2theta;
return cos_2theta_3;
// The solid angle of a flat pixel depends on the incidence angle to the detector
// normal, cos(alpha) = det_distance / |detector-frame position|. This is evaluated
// in the detector's own frame, so it is invariant under detector tilt (rot1/rot2/rot3),
// matching PyFAI solidAngleArray and MAX IV azint. It reduces to cos^3(2*theta) only
// for an untilted detector.
float u = (x - beam_x_pxl) * pixel_size_mm;
float v = (y - beam_y_pxl) * pixel_size_mm;
float cos_alpha = det_distance_mm / sqrtf(u * u + v * v + det_distance_mm * det_distance_mm);
return cos_alpha * cos_alpha * cos_alpha;
}
float DiffractionGeometry::CalcAzIntPolarizationCorr(float x, float y, float coeff) const {
common/DiffractionGeometry.h
-1
View File
@@ -54,7 +54,6 @@ public:
[[nodiscard]] float ResToPxl(float d_A) const;
[[nodiscard]] Coord ResToPxl(float d_A, float phi) const;
[[nodiscard]] float DistFromEwaldSphere(const Coord& recip) const;
[[nodiscard]] float CalcAzIntSolidAngleCorr(float q) const;
[[nodiscard]] float CalcAzIntSolidAngleCorr(float x, float y) const;
[[nodiscard]] float CalcAzIntPolarizationCorr(float x, float y, float coeff) const;
[[nodiscard]] std::pair<float, float> ResPhiToPxl(float d_A, float phi_rad) const;
docs/CHANGELOG.md
+1
View File
@@ -3,6 +3,7 @@
### 1.0.0-rc.158
This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use.
* Analysis: The azimuthal-integration solid-angle correction now follows the incidence angle to the detector normal (`cos^3` of that angle) instead of `cos^3(2*theta)`, so it is correct for a tilted detector and matches PyFAI `solidAngleArray` and MAX IV azint (unchanged for an untilted detector). Crystal geometry refinement (`XtalOptimizer`) no longer silently ignores an imported PONI `rot3` (rotation about the beam): it is applied as a fixed rotation in the residual so refinement stays consistent with the rest of the pipeline. Polarization and azimuthal binning already honoured `rot3` through the full PONI rotation.
* jfjoch_viewer: Open datasets on the WSL2/UNC filesystem (paths starting `\\`); write processing outputs next to the input file, with a Browse button and independent `_process.h5` / merged `.mtz`/`.cif` toggles; and show the determined space group in the merge-statistics window.
* rugnux: Accept an absolute `-o` output prefix in offline processing.
* Packaging: The self-contained Linux viewer `.tgz` now bundles cuFFT, so it runs without a system CUDA toolkit (`.deb`/`.rpm` are unchanged, distro-managed).
image_analysis/geom_refinement/XtalOptimizer.cpp
+20 -7
View File
@@ -13,6 +13,7 @@ struct XtalResidual {
XtalResidual(double x, double y,
double lambda,
double pixel_size,
double rot3,
double angle_rad,
double exp_h, double exp_k,
double exp_l,
@@ -20,6 +21,7 @@ struct XtalResidual {
: obs_x(x), obs_y(y),
inv_lambda(1.0/lambda),
pixel_size(pixel_size),
rot3(rot3),
exp_h(exp_h),
exp_k(exp_k),
exp_l(exp_l),
@@ -39,10 +41,10 @@ struct XtalResidual {
const T *const p1,
const T *const p2,
T *residual) const {
// PyFAI convention (left-handed for rot1/rot2):
// poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1)
// detector_rot[0] = rot1, detector_rot[1] = rot2 (rot3 = 0 assumed)
// PyFAI convention: poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1).
// detector_rot[0] = rot1, detector_rot[1] = rot2 are refined; rot3 is fixed
// (e.g. from a PONI import) and baked in here as a constant so that a non-zero
// rot3 is not silently dropped during refinement.
const T rot1 = detector_rot[0];
const T rot2 = detector_rot[1];
@@ -54,6 +56,10 @@ struct XtalResidual {
const T c2 = ceres::cos(rot2);
const T s2 = ceres::sin(rot2);
// Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0
const T c3 = T(cos(rot3));
const T s3 = T(sin(rot3));
// Detector coordinates in mm
const T det_x = (T(obs_x) - beam[0]) * T(pixel_size);
const T det_y = (T(obs_y) - beam[1]) * T(pixel_size);
@@ -65,9 +71,14 @@ struct XtalResidual {
const T t1_z = -s1 * det_x + c1 * det_z;
// Then apply Rx(-rot2): rotate around X
const T x = t1_x;
const T y = c2 * t1_y + s2 * t1_z;
const T z = -s2 * t1_y + c2 * t1_z;
const T t2_x = t1_x;
const T t2_y = c2 * t1_y + s2 * t1_z;
const T t2_z = -s2 * t1_y + c2 * t1_z;
// Then apply Rz(-rot3): rotate around Z (beam)
const T x = c3 * t2_x + s3 * t2_y;
const T y = -s3 * t2_x + c3 * t2_y;
const T z = t2_z;
// convert to recip space
const T lab_norm = ceres::sqrt(x * x + y * y + z * z);
@@ -183,6 +194,7 @@ struct XtalResidual {
const double obs_x, obs_y;
const double inv_lambda;
const double pixel_size;
const double rot3;
const double exp_h;
const double exp_k;
const double exp_l;
@@ -349,6 +361,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
new XtalResidual(pt.x, pt.y,
data.geom.GetWavelength_A(),
data.geom.GetPixelSize_mm(),
data.geom.GetPoniRot3_rad(),
angle_rad,
h, k, l,
data.crystal_system)),
tests/DiffractionGeometryTest.cpp
+25 -6
View File
@@ -117,13 +117,30 @@ TEST_CASE("DiffractionGeometry_SolidAngleCorrection","") {
x.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(75);
DiffractionGeometry geom = x.GetDiffractionGeometry();
REQUIRE(geom.CalcAzIntSolidAngleCorr(0.0) == 1.0f);
REQUIRE(geom.CalcAzIntSolidAngleCorr(2 * M_PI) == Catch::Approx(0.5f * 0.5f * 0.5f));
// theta = 30 deg
// cos (2 * theta) = 1/2
// 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","") {
@@ -424,7 +441,9 @@ Rot3: 0.0
Wavelength: 1e-10
*/
// Not sure why, but PyFAI solidAngleArray doesn't take into account poni rotation (???)
// 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(2000).BeamY_pxl(1000).IncidentEnergy_keV(WVL_1A_IN_KEV);
DiffractionGeometry geom = x.GetDiffractionGeometry();