indexing: deconvolve energy bandwidth from the profile radius

The profile radius (intrinsic excitation-error width = mosaicity + divergence)
was the plain RMS of dist_ewald over indexed spots. With a finite energy
bandwidth that spread is broadened by the bandwidth's radial smear
sigma_bw = bandwidth_sigma*lambda/(2 d^2), which prediction then re-applies per
reflection - so bandwidth was counted twice and the radius was inflated (most at
high resolution, sigma_bw ~ 1/d^2). Subtract the bandwidth variance from the
measured spread so the radius is the intrinsic width. bandwidth = 0
(monochromatic / rotation) is unchanged. Small for narrow bandwidths (~6% of the
variance, ~4% radius on the 1% jet); matters for wide-bandwidth / pink beam.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

image_analysis/indexing/AnalyzeIndexing.cpp

@@ -369,7 +369,9 @@ bool AnalyzeIndexing(DataMessage &message,
                 message.spots[i].indexed = indexed_spots[i];
                 message.spots[i].lattice = indexed_spots[i] ? 0 : -1;
             }
-            message.profile_radius = FitProfileRadius(message.spots);
+            message.profile_radius = FitProfileRadius(message.spots,
+                                                      experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f,
+                                                      experiment.GetWavelength_A());
             message.spot_count_indexed = nspots_indexed;
             message.indexing_lattice = latt;
             message.indexing_unit_cell = latt.GetUnitCell();

image_analysis/indexing/FitProfileRadius.cpp

@@ -31,21 +31,31 @@ std::optional<float> FitProfileRadius_MAD(const std::vector<SpotToSave>& xs) {
     return 1.4826f * med;
 }
 
-std::optional<float> FitProfileRadius(const std::vector<SpotToSave>& spots) {
-    float sum_squares = 0.0f;
+std::optional<float> FitProfileRadius(const std::vector<SpotToSave>& spots,
+                                      float bandwidth_sigma, float wavelength_A) {
+    double sum_squares = 0.0;     // measured excitation-error variance (sum dist_ewald^2)
+    double sum_bw_var = 0.0;      // energy-bandwidth contribution to subtract out
     int count = 0;
 
     for (const auto &s: spots) {
-        if (s.indexed) {
-            sum_squares += s.dist_ewald_sphere * s.dist_ewald_sphere;
-            count++;
+        if (!s.indexed)
+            continue;
+        sum_squares += static_cast<double>(s.dist_ewald_sphere) * s.dist_ewald_sphere;
+        // The energy bandwidth smears each reflection radially by sigma_bw = bandwidth_sigma*|recip_z|
+        // = bandwidth_sigma*lambda/(2 d^2) (the same term prediction re-adds per reflection, ~1/d^2 so
+        // largest at high resolution). Deconvolve it from the measured spread so the profile radius is
+        // the *intrinsic* mosaicity+divergence width and bandwidth is not double-counted at prediction.
+        if (bandwidth_sigma > 0.0f && s.d_A > 0.0f) {
+            const double sigma_bw = bandwidth_sigma * wavelength_A / (2.0 * static_cast<double>(s.d_A) * s.d_A);
+            sum_bw_var += sigma_bw * sigma_bw;
         }
+        count++;
     }
 
     if (count == 0)
         return std::nullopt;
 
-    auto std_dev = std::sqrt(sum_squares / count);
-    return std_dev;
+    const double variance = std::max(0.0, (sum_squares - sum_bw_var) / count);
+    return static_cast<float>(std::sqrt(variance));
 }
 

image_analysis/indexing/FitProfileRadius.h

@@ -9,5 +9,11 @@
 #include "../../common/SpotToSave.h"
 
 std::optional<float> FitProfileRadius_MAD(const std::vector<SpotToSave>& spots);
-std::optional<float> FitProfileRadius(const std::vector<SpotToSave>& spots);
+
+// Intrinsic excitation-error (mosaicity+divergence) width from the indexed-spot spread. When a finite
+// energy bandwidth is given, its radial smear (bandwidth_sigma*lambda/(2 d^2)) is deconvolved out, so
+// the result is the intrinsic width and bandwidth is not double-counted by prediction (which re-adds
+// it). bandwidth_sigma = 0 reproduces the plain RMS (monochromatic / rotation).
+std::optional<float> FitProfileRadius(const std::vector<SpotToSave>& spots,
+                                      float bandwidth_sigma = 0.0f, float wavelength_A = 0.0f);