BraggPrediction: Include X-ray bandwidth

image_analysis/IndexAndRefine.cpp

@@ -289,7 +289,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
         .max_hkl = 100,
         .centering = outcome.symmetry.centering,
         .wedge_deg = std::fabs(wedge_deg),
-        .mosaicity_deg = std::fabs(mos_deg)
+        .mosaicity_deg = std::fabs(mos_deg),
+        // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged.
+        .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f
     };
 
     // Select the integration path: classical 2D integration, or the experimental

image_analysis/bragg_prediction/BraggPrediction.cpp

@@ -4,6 +4,13 @@
 #include "BraggPrediction.h"
 #include "../bragg_integration/SystematicAbsence.h"
 
+namespace {
+    // Number of bandwidth sigmas included in the (radially thickened) Ewald-shell
+    // acceptance window. 3σ captures essentially the whole pink-beam smear; matches
+    // the conservative end of the mosaicity cutoff used by callers.
+    constexpr float kBandwidthCutoffSigmas = 3.0f;
+}
+
 BraggPrediction::BraggPrediction(int max_reflections)
     : max_reflections(max_reflections), reflections(max_reflections) {}
 
@@ -76,7 +83,18 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                 float S_len = sqrtf(S_x * S_x + S_y * S_y + S_z * S_z);
                 float dist_ewald_sphere = std::fabs(S_len - one_over_wavelength);
 
-                if (dist_ewald_sphere <= settings.ewald_dist_cutoff ) {
+                // Energy bandwidth thickens the Ewald shell radially: at the
+                // diffraction condition |S|-1/λ shifts by recip_z·(Δλ/λ), i.e.
+                // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d², identical to
+                // PixelRefine's R_bw). Broaden the acceptance window in quadrature so
+                // high-resolution shells (smeared most, ∝1/d²) are not clipped.
+                float radial_cutoff = settings.ewald_dist_cutoff;
+                if (settings.bandwidth_sigma > 0.0f) {
+                    const float bw_tol = kBandwidthCutoffSigmas * settings.bandwidth_sigma * std::fabs(recip_z);
+                    radial_cutoff = std::sqrt(radial_cutoff * radial_cutoff + bw_tol * bw_tol);
+                }
+
+                if (dist_ewald_sphere <= radial_cutoff ) {
                     const float s0_p0 = S0.x * recip_x + S0.y * recip_y + S0.z * recip_z;
                     const float val = s0_sq * recip_sq - s0_p0 * s0_p0;
 

image_analysis/bragg_prediction/BraggPrediction.h

@@ -18,6 +18,12 @@ struct BraggPredictionSettings {
     float mosaicity_deg = 0.2f;
     float min_zeta = 0.05;
     float mosaicity_multiplier = 4.0;
+    // Relative X-ray bandwidth Δλ/λ expressed as a Gaussian sigma (0 = monochromatic).
+    // When > 0 the Ewald-shell acceptance is thickened radially per reflection by
+    // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no
+    // longer clips high-resolution reflections. Must stay the last member so existing
+    // designated initializers remain valid.
+    float bandwidth_sigma = 0.0f;
 };
 
 class BraggPrediction {

image_analysis/bragg_prediction/BraggPredictionGPU.cu

@@ -8,6 +8,10 @@
 #include <cuda_runtime.h>
 
 namespace {
+    // Number of bandwidth sigmas included in the (radially thickened) Ewald-shell
+    // acceptance window. Mirrors the CPU BraggPrediction path.
+    constexpr float kBandwidthCutoffSigmas = 3.0f;
+
     __device__ inline bool is_odd(int v) { return (v & 1) != 0; }
 
     __device__ inline float angle_from_ewald_sphere_deg(const Coord &S0, float recip_x, float recip_y, float recip_z, float recip_sq) {
@@ -95,7 +99,14 @@ namespace {
         float Sz = recip_z + C.S0.z;
         float S_len = sqrtf(Sx * Sx + Sy * Sy + Sz * Sz);
         float dist_ewald = fabsf(S_len - C.one_over_wavelength);
-        if (dist_ewald > C.ewald_cutoff) return false;
+        // Energy bandwidth thickens the Ewald shell radially: σ_bw = |recip_z|·(Δλ/λ)
+        // (= bλ/2d²). Broaden the acceptance window in quadrature (see CPU path).
+        float radial_cutoff = C.ewald_cutoff;
+        if (C.bandwidth_sigma > 0.0f) {
+            const float bw_tol = kBandwidthCutoffSigmas * C.bandwidth_sigma * fabsf(recip_z);
+            radial_cutoff = sqrtf(radial_cutoff * radial_cutoff + bw_tol * bw_tol);
+        }
+        if (dist_ewald > radial_cutoff) return false;
         float Srx = C.rot[0] * Sx + C.rot[1] * Sy + C.rot[2] * Sz;
         float Sry = C.rot[3] * Sx + C.rot[4] * Sy + C.rot[5] * Sz;
         float Srz = C.rot[6] * Sx + C.rot[7] * Sy + C.rot[8] * Sz;
@@ -145,7 +156,8 @@ namespace {
                                           const CrystalLattice &lattice,
                                           float high_res_A,
                                           float ewald_dist_cutoff,
-                                          char centering) {
+                                          char centering,
+                                          float bandwidth_sigma) {
         KernelConsts kc{};
         auto geom = experiment.GetDiffractionGeometry();
         kc.det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
@@ -157,6 +169,7 @@ namespace {
         kc.one_over_dmax_sq = one_over_dmax * one_over_dmax;
         kc.one_over_wavelength = 1.0f / geom.GetWavelength_A();
         kc.ewald_cutoff = ewald_dist_cutoff;
+        kc.bandwidth_sigma = bandwidth_sigma;
         kc.Astar = lattice.Astar();
         kc.Bstar = lattice.Bstar();
         kc.Cstar = lattice.Cstar();
@@ -178,7 +191,8 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment,
                              const CrystalLattice &lattice,
                              const BraggPredictionSettings &settings) {
     // Build constants on host
-    KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, settings.centering);
+    KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff,
+                                        settings.centering, settings.bandwidth_sigma);
     cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream);
     cudaMemsetAsync(d_count, 0, sizeof(int), stream);
 

image_analysis/bragg_prediction/BraggPredictionGPU.h

@@ -18,6 +18,7 @@ struct KernelConsts {
     float one_over_wavelength;
     float one_over_dmax_sq;
     float ewald_cutoff;
+    float bandwidth_sigma;   // relative Δλ/λ (sigma); 0 = monochromatic
     Coord Astar, Bstar, Cstar, S0;
     float rot[9];
     char centering;

image_analysis/pixel_refinement/PixelRefine.cpp

@@ -475,7 +475,8 @@ void PixelRefine::Run(const T *image,
     const BraggPredictionSettings settings_prediction{
         .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
         .max_hkl = 100,
-        .centering = data.centering
+        .centering = data.centering,
+        .bandwidth_sigma = static_cast<float>(data.bandwidth)  // relative Δλ/λ sigma
     };
 
     const auto azim_result = profile.GetResult();
@@ -942,7 +943,8 @@ std::vector<float> PixelRefine::PredictImage(const AzimuthalIntegrationProfile &
     const BraggPredictionSettings settings_prediction{
         .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
         .max_hkl = 100,
-        .centering = data.centering
+        .centering = data.centering,
+        .bandwidth_sigma = static_cast<float>(data.bandwidth)  // relative Δλ/λ sigma
     };
     const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
     const auto &predicted = prediction.GetReflections();