BraggIntegrate: Lorentz-Polarization correction consistent (roughly) with XDS

common/Reflection.h

@@ -23,6 +23,7 @@ struct Reflection {
     float dist_ewald;
     float lp;
     bool observed = false;
+    float S_x, S_y, S_z;
 };
 
 #endif //JFJOCH_REFLECTION_H

image_analysis/bragg_integration/BraggIntegrate2D.cpp

@@ -88,6 +88,11 @@ std::vector<Reflection> IntegrateInternal(const DiffractionExperiment &experimen
     const float r_2_sq = settings.GetR2() * settings.GetR2();
     const float r_3_sq = settings.GetR3() * settings.GetR3();
 
+    const Coord S0 = geom.GetScatteringVector();
+    Coord m2 = Coord(0,0,0); // zero vector
+    if (experiment.GetGoniometer().has_value())
+        m2 = experiment.GetGoniometer()->GetAxis().Normalize();
+
     for (int i = 0; i < npredicted; i++) {
         auto r = predicted.at(i);
         IntegrateReflection(r, ptr, image.GetWidth(), image.GetHeight(), special_value, saturation,
@@ -97,11 +102,10 @@ std::vector<Reflection> IntegrateInternal(const DiffractionExperiment &experimen
             if (experiment.GetPolarizationFactor())
                 corr /= geom.CalcAzIntPolarizationCorr(r.predicted_x, r.predicted_y,
                                                            experiment.GetPolarizationFactor().value());
+            Coord S{r.S_x, r.S_y, r.S_z};
+            if (experiment.GetGoniometer().has_value())
+                corr *= std::fabs(m2 * (S % S0)) / (S * S0);
 
-            if (experiment.GetGoniometer().has_value()) {
-                float two_theta = geom.TwoTheta_rad(r.predicted_x, r.predicted_y);
-                corr *= std::sin(two_theta);
-            }
             r.lp = 1.0f / corr;
             r.I *= corr;
             r.bkg *= corr;

image_analysis/bragg_prediction/BraggPrediction.cpp

@@ -35,10 +35,11 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
     float beam_y = geom.GetBeamY_pxl();
     float det_distance = geom.GetDetectorDistance_mm();
     float pixel_size = geom.GetPixelSize_mm();
-    float coeff_const = det_distance / pixel_size;
+    float F = det_distance / pixel_size;
 
     int i = 0;
 
+
     for (int h = -settings.max_hkl; h <= settings.max_hkl; h++) {
         // Precompute A* h contribution
         const float Ah_x = Astar.x * h;
@@ -83,9 +84,9 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                         continue;
 
                     // Project to detector coordinates
-                    float coeff = coeff_const / S_rot_z;
-                    float x = beam_x + S_rot_x * coeff;
-                    float y = beam_y + S_rot_y * coeff;
+                    // Assume detector is along x,y,z coordinates after rotation
+                    float x = beam_x + F * S_rot_x / S_rot_z;
+                    float y = beam_y + F * S_rot_y / S_rot_z;
 
                     if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl))
                         continue;
@@ -98,7 +99,10 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                         .predicted_x = x,
                         .predicted_y = y,
                         .d = d,
-                        .dist_ewald = dist_ewald_sphere
+                        .dist_ewald = dist_ewald_sphere,
+                        .S_x = S_x,
+                        .S_y = S_y,
+                        .S_z = S_z
                     };
                     ++i;
                 }

image_analysis/bragg_prediction/BraggPredictionGPU.cu

@@ -81,6 +81,9 @@ namespace {
         out.predicted_y = y;
         out.d = 1.0f / sqrtf(recip_sq);
         out.dist_ewald = dist_ewald;
+        out.S_x = Sx;
+        out.S_y = Sy;
+        out.S_z = Sz;
         return true;
     }
 

tools/jfjoch_process.cpp

@@ -133,6 +133,7 @@ int main(int argc, char **argv) {
     experiment.Mode(DetectorMode::Standard); // Ensure full image analysis
     experiment.PixelSigned(true);
     experiment.OverwriteExistingFiles(true);
+    experiment.PolarizationFactor(0.0);
 
     // Configure Indexing
     IndexingSettings indexing_settings;