Rotation: land the rest of the bandwidth term

f4e281b2f described this change in full but committed only one of its six files.
What went in was RotationScaleMerge.cpp - the merge widening the partiality it
recomputes from the smoothed mosaicity. That is precisely the part which is unsafe on
its own, by the original message's own argument: without the mosaicity fit subtracting
the term before fitting, the bandwidth is counted twice, and without the predictor
widening its acceptance window, the partiality the merge recomputes no longer matches
the one integration measured.

Add the five files that were left behind: the rotation predictor and its GPU twin
widen the acceptance window and the partiality handed to integration, the settings
struct carries the term, and CalcMosaicityXDS deconvolves it before fitting so what it
returns is the intrinsic mosaicity rather than the mosaicity plus the beam.

Monochromatic data is untouched by construction - every hunk is guarded on a non-zero
bandwidth, which is read from incident_wavelength_spread or --bandwidth and is absent
from every dataset in the rotation battery. Verified on the one dataset that has a
bandwidth: at --bandwidth 0, the merge table is identical to the branch tip; with the
bandwidth set, the fitted mosaicity drops 0.0718 -> 0.0694 deg as the deconvolution
takes effect and CC1/2 in the outermost shell recovers 30.3 -> 31.4%.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

docs/CHANGELOG.md

@@ -11,6 +11,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as
 * Spot finding: `--spot-sigma` now defaults to 4.0 (was 3.0) and `--max-spots` to 1000 spots per image (was 250).
 * Resolution limits: Bragg integration, azimuthal integration and spot finding all default to **as far as the detector reaches**; `--integration-high-resolution` and `--spot-high-resolution` still set one by hand, and 0 means "no limit" at either end.
 * Bragg prediction: How far the predictor walks the lattice is a setting (`bragg_integration_settings.max_hkl`, `--max-hkl`) instead of a fixed 100, derived per crystal offline and held fixed online.
+* Rotation data: a set X-ray bandwidth (`--bandwidth`) now widens each reflection's rocking curve by Δλ/λ·tan(θ) in prediction and partiality, and is deconvolved out of the fitted mosaicity; monochromatic data are unaffected.
 * Bragg integration: The local background ring is made robust with a **high-side sigma clip** (`--background-clip <n>`, default 4) instead of the symmetric trimmed mean; the trim stays reachable with `--background-trim <f>`.
 * Bragg integration: The **uncertainty of the background estimate** is now propagated into `sigma`, which both engines previously omitted.
 * Bragg integration: Profile-fit `sigma` is **no longer inflated on weak reflections** - the fit weights take the signal estimate as it is instead of clamping it at zero, and the per-pixel variance floor is 0.01 counts instead of 1/12.

image_analysis/bragg_prediction/BraggPrediction.h

@@ -25,9 +25,15 @@ struct BraggPredictionSettings {
     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
+    // Stills: 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.
+    // Rotation: differentiating Bragg's law at fixed d gives an extra rocking width
+    // Δθ = bandwidth_sigma·tan(θ_B), a spread in the same glancing angle the mosaic spread
+    // smears, so it adds to σ_M in quadrature. It is NOT divided by ζ: rotating the crystal
+    // by Δφ changes θ by ζ·Δφ, so the 1/ζ that turns an angular width into a rotation width is
+    // already the one the partiality applies to σ_M. CalcMosaicityXDS deconvolves the same term
+    // out of the fitted σ_M, so it is not counted twice.
     float bandwidth_sigma = 0.0f;
 };
 

image_analysis/bragg_prediction/BraggPredictionRot.cpp

@@ -50,6 +50,16 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
     const float mos_angle_rad = settings.mosaicity_deg * static_cast<float>(PI) / 180.f;
     const float half_wedge_angle_rad = settings.wedge_deg * static_cast<float>(PI) / 180.f / 2.0f ;
 
+    // Energy bandwidth widens the rocking curve. Differentiating Bragg's law at fixed d gives
+    // dtheta = (dlambda/lambda) tan(theta_B), a spread in the same glancing angle the mosaic spread
+    // smears, so it adds to sigma_M in quadrature. It is NOT divided by zeta here: rotating the
+    // crystal by dphi changes theta by zeta*dphi, so the 1/zeta that turns an angular width into a
+    // rotation width is already the one c1 (and the epsilon3 cutoff) applies to sigma_M. The fitted
+    // sigma_M has this term deconvolved out (CalcMosaicityXDS), so it is not counted twice.
+    // sin(theta_B) = lambda/(2d) = lambda*|p0|/2. Zero bandwidth leaves every reflection untouched.
+    const float bandwidth_sigma = settings.bandwidth_sigma;
+    const float half_wavelength_A = geom.GetWavelength_A() / 2.0f;
+
     for (int h = -settings.max_h; h <= settings.max_h; h++) {
         // Precompute A* h contribution
 
@@ -85,6 +95,15 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
                 if ((rho_sq < p_m3 * p_m3) || (p0_sq > 4 * S0 * S0))
                     continue;
 
+                // Effective rocking width for this reflection: mosaicity broadened by the bandwidth
+                // term. sin(theta_B) <= 1 is guaranteed by the p0_sq test just above.
+                float mos_eff_rad = mos_angle_rad;
+                if (bandwidth_sigma > 0.0f) {
+                    const float sin_theta = half_wavelength_A * std::sqrt(p0_sq);
+                    const float dphi_bw = bandwidth_sigma * sin_theta / std::sqrt(1.0f - sin_theta * sin_theta);
+                    mos_eff_rad = std::sqrt(mos_angle_rad * mos_angle_rad + dphi_bw * dphi_bw);
+                }
+
                 for (const auto& p_m1 : p_m1_opt) {
                     if (i >= max_reflections)
                         continue;
@@ -105,7 +124,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
 
                     float epsilon3 = std::fabs(phi * zeta_abs);
 
-                    if (epsilon3 > settings.mosaicity_multiplier * mos_angle_rad)
+                    if (epsilon3 > settings.mosaicity_multiplier * mos_eff_rad)
                         continue;
 
                     // Reciprocal Lorentz (Kabsch 2010): L^-1 = |m2 . (S x S0)| / (|S| |S0|) =
@@ -114,7 +133,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys
                     // corrupts the absolute/Wilson scale (it cancels within a resolution shell, so
                     // CC1/2 / CCref / R-meas are neutral).
                     const float lorentz_reciprocal = std::fabs(m2 * (S % S0)) / (S.Length() * S0.Length());
-                    const float c1 = zeta_abs / (std::sqrt(2.0f) * mos_angle_rad);
+                    const float c1 = zeta_abs / (std::sqrt(2.0f) * mos_eff_rad);
 
                     const float partiality = (std::erf((phi + half_wedge_angle_rad) * c1)
                         - std::erf((phi - half_wedge_angle_rad) * c1)) / 2.0f;

image_analysis/bragg_prediction/BraggPredictionRotGPU.cu

@@ -86,6 +86,15 @@ namespace {
         float p_m1_pos = sqrtf(rho_sq - p_m3 * p_m3);
         float p_m1_arr[2] = {p_m1_pos, -p_m1_pos};
 
+        // Effective rocking width: mosaicity broadened by the bandwidth term (consistent with CPU),
+        // dtheta = (dlambda/lambda) tan(theta_B) in quadrature with sigma_M, zeta-free.
+        float mos_eff_rad = C.mos_angle_rad;
+        if (C.bandwidth_sigma > 0.0f) {
+            float sin_theta = C.half_wavelength_A * sqrtf(p0_sq);
+            float dphi_bw = C.bandwidth_sigma * sin_theta / sqrtf(1.0f - sin_theta * sin_theta);
+            mos_eff_rad = sqrtf(C.mos_angle_rad * C.mos_angle_rad + dphi_bw * dphi_bw);
+        }
+
         int count = 0;
         for (int idx = 0; idx < 2; ++idx) {
             float p_m1 = p_m1_arr[idx];
@@ -117,7 +126,7 @@ namespace {
 
             // epsilon3 cutoff check (consistent with CPU, Kabsch formulation)
             float epsilon3 = fabsf(phi * zeta_abs);
-            if (epsilon3 > C.mosaicity_multiplier * C.mos_angle_rad)
+            if (epsilon3 > C.mosaicity_multiplier * mos_eff_rad)
                 continue;
 
             float cx, cy, cz;
@@ -131,7 +140,7 @@ namespace {
 
             // Partiality calculation (Kabsch formulation)
             // c1 = sqrt(2) * sigma / zeta, where sigma = mosaicity
-            float c1 = zeta_abs / (sqrtf(2.0f) * C.mos_angle_rad);
+            float c1 = zeta_abs / (sqrtf(2.0f) * mos_eff_rad);
             float half_wedge = C.wedge_angle_rad / 2.0f;
             float partiality = (erff((phi + half_wedge) * c1)
                                 - erff((phi - half_wedge) * c1)) / 2.0f;
@@ -218,6 +227,8 @@ namespace {
          kc.wedge_angle_rad = settings.wedge_deg * static_cast<float>(PI) / 180.0f;
          kc.min_zeta = settings.min_zeta;
          kc.mosaicity_multiplier = settings.mosaicity_multiplier;
+         kc.bandwidth_sigma = settings.bandwidth_sigma;
+         kc.half_wavelength_A = geom.GetWavelength_A() / 2.0f;
 
          kc.Astar = lattice.Astar();
          kc.Bstar = lattice.Bstar();

image_analysis/bragg_prediction/BraggPredictionRotGPU.h

@@ -20,6 +20,8 @@ struct KernelConstsRot {
     float wedge_angle_rad;
     float min_zeta;
     float mosaicity_multiplier;
+    float bandwidth_sigma;      // relative dlambda/lambda as a sigma; 0 = monochromatic
+    float half_wavelength_A;
     Coord Astar, Bstar, Cstar, S0;
     Coord m1, m2, m3;
     float m2_S0;

image_analysis/indexing/AnalyzeIndexing.cpp

@@ -54,11 +54,17 @@ namespace {
     // delta_phi = oscillation range (radians)
     // sigma_M = mosaicity (radians)
     // zeta = |m2 · e1| Lorentz factor component
-    inline float R_fraction(float tau, float delta_phi, float sigma_M, float zeta) {
+    // sigma_bw = this reflection's energy-bandwidth rocking width (radians, zeta-free). It broadens the
+    //   observed tau spread on top of the mosaicity, so carrying it here makes the fitted sigma_M the
+    //   INTRINSIC width that prediction and scaling re-broaden per reflection - the same deconvolution
+    //   FitProfileRadius does for the profile radius. 0 = monochromatic.
+    inline float R_fraction(float tau, float delta_phi, float sigma_M, float zeta, float sigma_bw) {
         if (zeta < 1e-6f || sigma_M < 1e-9f)
             return 0.0f;
 
-        const float sigma_eff = sigma_M / zeta;
+        const float sigma_total = sigma_bw > 0.0f
+                                  ? std::sqrt(sigma_M * sigma_M + sigma_bw * sigma_bw) : sigma_M;
+        const float sigma_eff = sigma_total / zeta;
         const float sqrt2_sigma = std::sqrt(2.0f) * sigma_eff;
 
         if (sqrt2_sigma < 1e-12f)
@@ -74,11 +80,12 @@ namespace {
     // Returns sum of log(R) for all reflections
     inline double log_likelihood(const std::vector<float>& tau_values,
                                  const std::vector<float>& zeta_values,
+                                 const std::vector<float>& sigma_bw_values,
                                  float delta_phi,
                                  float sigma_M) {
         double ll = 0.0;
         for (size_t i = 0; i < tau_values.size(); ++i) {
-            float R = R_fraction(tau_values[i], delta_phi, sigma_M, zeta_values[i]);
+            float R = R_fraction(tau_values[i], delta_phi, sigma_M, zeta_values[i], sigma_bw_values[i]);
             if (std::isfinite(R) && R > 1e-30f) {
                 ll += std::log(static_cast<double>(R));
             } else {
@@ -91,6 +98,7 @@ namespace {
     // Golden section search for maximum likelihood sigma_M
     inline float find_sigma_M_mle(const std::vector<float>& tau_values,
                                   const std::vector<float>& zeta_values,
+                                  const std::vector<float>& sigma_bw_values,
                                   float delta_phi,
                                   float sigma_min_deg = 0.001f,
                                   float sigma_max_deg = 2.0f) {
@@ -105,8 +113,8 @@ namespace {
         const float tol = 1e-6f;
         int iter = 0;
         while (std::fabs(b - a) > tol && iter++ < 100) {
-            double fc = log_likelihood(tau_values, zeta_values, delta_phi, c);
-            double fd = log_likelihood(tau_values, zeta_values, delta_phi, d);
+            double fc = log_likelihood(tau_values, zeta_values, sigma_bw_values, delta_phi, c);
+            double fd = log_likelihood(tau_values, zeta_values, sigma_bw_values, delta_phi, d);
 
             if (fc > fd) {
                 b = d;
@@ -237,6 +245,12 @@ namespace {
         const Coord S0 = experiment.GetScatteringVector();
         const float delta_phi_rad = deg_to_rad(axis.GetWedge_deg());
 
+        // Energy bandwidth adds sigma_bw = (dlambda/lambda)*tan(theta_B) to each reflection's rocking
+        // width; deconvolving it here leaves sigma_M the intrinsic mosaicity, so prediction and scaling
+        // can re-add it per reflection without counting it twice. sin(theta_B) = lambda*|pstar|/2.
+        const float bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f;
+        const float half_wavelength_A = experiment.GetWavelength_A() / 2.0f;
+
         // Fit from the strongest spots only, never the whole list. A spot is detected when
         // I_full * R(tau) clears the finder's threshold, so the deeper the spot list reaches the more
         // large-|tau| partially-recorded spots it holds - and sigma_M is fitted from exactly that tau
@@ -253,8 +267,10 @@ namespace {
 
         std::vector<float> tau_values;
         std::vector<float> zeta_values;
+        std::vector<float> sigma_bw_values;
         tau_values.reserve(n_fit);
         zeta_values.reserve(n_fit);
+        sigma_bw_values.reserve(n_fit);
 
         for (size_t si = 0; si < n_fit; ++si) {
             const auto &s = spots[si];
@@ -301,15 +317,22 @@ namespace {
             if (!std::isfinite(zeta) || !std::isfinite(tau_rad) || zeta < 0.1f)
                 continue;
 
+            float sigma_bw = 0.0f;
+            if (bandwidth_sigma > 0.0f) {
+                const float sin_theta = half_wavelength_A * pstar.Length();
+                sigma_bw = bandwidth_sigma * sin_theta / std::sqrt(1.0f - sin_theta * sin_theta);
+            }
+
             tau_values.push_back(tau_rad);
             zeta_values.push_back(zeta);
+            sigma_bw_values.push_back(sigma_bw);
         }
 
         if (tau_values.size() < 10)
             return std::nullopt;
 
         // Find sigma_M by maximizing log-likelihood
-        float sigma_M_rad = find_sigma_M_mle(tau_values, zeta_values, delta_phi_rad);
+        float sigma_M_rad = find_sigma_M_mle(tau_values, zeta_values, sigma_bw_values, delta_phi_rad);
 
         return rad_to_deg(sigma_M_rad);
     }