From acd20256760f647585e458c430f6917dba032ec5 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 11 Jul 2026 11:18:19 +0200 Subject: [PATCH] azint: tilt-correct solid-angle correction; honour imported rot3 in refinement 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 --- common/DiffractionGeometry.cpp | 19 ++++++------ common/DiffractionGeometry.h | 1 - docs/CHANGELOG.md | 1 + .../geom_refinement/XtalOptimizer.cpp | 27 +++++++++++----- tests/DiffractionGeometryTest.cpp | 31 +++++++++++++++---- 5 files changed, 55 insertions(+), 24 deletions(-) diff --git a/common/DiffractionGeometry.cpp b/common/DiffractionGeometry.cpp index 1234ffa8..9f3ed32c 100644 --- a/common/DiffractionGeometry.cpp +++ b/common/DiffractionGeometry.cpp @@ -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(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 { diff --git a/common/DiffractionGeometry.h b/common/DiffractionGeometry.h index 171977b2..bf1924f2 100644 --- a/common/DiffractionGeometry.h +++ b/common/DiffractionGeometry.h @@ -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 ResPhiToPxl(float d_A, float phi_rad) const; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 177887fe..d261ced7 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -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). diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index a4676450..368f1e67 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -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)), diff --git a/tests/DiffractionGeometryTest.cpp b/tests/DiffractionGeometryTest.cpp index 8b60f95f..306304b5 100644 --- a/tests/DiffractionGeometryTest.cpp +++ b/tests/DiffractionGeometryTest.cpp @@ -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();