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 <noreply@anthropic.com>

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<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));
+    // The solid angle of a flat pixel depends on the incidence angle to the detector
-    float cos_2theta_3 = cos_2theta * cos_2theta * cos_2theta;
+    // normal, cos(alpha) = det_distance / |detector-frame position|. This is evaluated
-    return cos_2theta_3;
+    // 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

@@ -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

@@ -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

@@ -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):
+        // PyFAI convention: poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1).
-        // poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1)
+        // detector_rot[0] = rot1, detector_rot[1] = rot2 are refined; rot3 is fixed
-        // detector_rot[0] = rot1, detector_rot[1] = rot2 (rot3 = 0 assumed)
+        // (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 t2_x = t1_x;
-        const T y = c2 * t1_y + s2 * t1_z;
+        const T t2_y = c2 * t1_y + s2 * t1_z;
-        const T z = -s2 * t1_y + c2 * 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

@@ -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);
+    // At the beam centre the correction is 1
-    REQUIRE(geom.CalcAzIntSolidAngleCorr(2 * M_PI) == Catch::Approx(0.5f * 0.5f * 0.5f));
-
-    // theta = 30 deg
-    // cos (2 * theta) = 1/2
     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();