PixelRefine: Improvements to accept more reasonable count of reflections

image_analysis/pixel_refinement/PixelRefine.cpp

@@ -551,6 +551,35 @@ void PixelRefine::BuildParameterBlocks(const PixelRefineData &data,
     }
 }
 
+BraggPredictionSettings PixelRefine::BuildPredictionSettings(const PixelRefineData &data) const {
+    // Radial Ewald-acceptance band: predict exactly the reflections the merge's
+    // partiality floor would still keep, and no tighter. A reflection's partiality
+    // is exp(-s^2 / R0_eff^2) (s = excitation error), so it survives min_partiality
+    // when |s| <= R0 * sqrt(-ln(min_partiality)). Using that as the cutoff keeps
+    // prediction and the partiality cut consistent. The struct default (0.0005) is
+    // far narrower than R0 (~0.005), so previously most keepable reflections were
+    // never predicted and per-reflection multiplicity collapsed ~4x.
+    const double min_part = std::clamp(experiment.GetScalingSettings().GetMinPartiality(), 1e-6, 0.999);
+    const double r0 = std::max(data.R[0], 1e-4);
+    const double cutoff = r0 * std::sqrt(-std::log(min_part));
+
+    // Relative bandwidth (sigma of dlambda/lambda): the explicit PixelRefine model
+    // value if set, else the experiment's nominal bandwidth (FWHM -> sigma). >0
+    // thickens the band radially at high resolution (the 1/d^2 pink-beam smear),
+    // matching the integrator and preventing the outer shells from being clipped.
+    float bw_sigma = static_cast<float>(data.bandwidth);
+    if (bw_sigma <= 0.0f)
+        bw_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f;
+
+    return BraggPredictionSettings{
+        .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
+        .ewald_dist_cutoff = static_cast<float>(cutoff),
+        .max_hkl = 100,
+        .centering = data.centering,
+        .bandwidth_sigma = bw_sigma
+    };
+}
+
 template<class T>
 void PixelRefine::Run(const T *image,
                       BraggPrediction &prediction,
@@ -561,12 +590,7 @@ void PixelRefine::Run(const T *image,
     const double lambda = data.geom.GetWavelength_A();
     const double pixel_size = data.geom.GetPixelSize_mm();
 
-    const BraggPredictionSettings settings_prediction{
-        .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
-        .max_hkl = 100,
-        .centering = data.centering,
-        .bandwidth_sigma = static_cast<float>(data.bandwidth)  // relative Δλ/λ sigma
-    };
+    const BraggPredictionSettings settings_prediction = BuildPredictionSettings(data);
 
     const int radius = data.shoebox_radius;
     const int bkg_outer_radius = std::max(radius + 1, data.bkg_outer_radius);
@@ -1045,12 +1069,7 @@ std::vector<float> PixelRefine::PredictImage(const T *image,
             .PoniRot1_rad(data.geom.GetPoniRot1_rad())
             .PoniRot2_rad(data.geom.GetPoniRot2_rad());
 
-    const BraggPredictionSettings settings_prediction{
-        .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
-        .max_hkl = 100,
-        .centering = data.centering,
-        .bandwidth_sigma = static_cast<float>(data.bandwidth)  // relative Δλ/λ sigma
-    };
+    const BraggPredictionSettings settings_prediction = BuildPredictionSettings(data);
     const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
     const auto &predicted = prediction.GetReflections();
     const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);
@@ -1139,12 +1158,7 @@ std::vector<float> PixelRefine::ChiSquaredImage(const T *image,
             .PoniRot1_rad(data.geom.GetPoniRot1_rad())
             .PoniRot2_rad(data.geom.GetPoniRot2_rad());
 
-    const BraggPredictionSettings settings_prediction{
-        .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
-        .max_hkl = 100,
-        .centering = data.centering,
-        .bandwidth_sigma = static_cast<float>(data.bandwidth)
-    };
+    const BraggPredictionSettings settings_prediction = BuildPredictionSettings(data);
     const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction);
     const auto &predicted = prediction.GetReflections();
     const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius);