diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 12540b3c..98f0a7af 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/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(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(cutoff), + .max_hkl = 100, + .centering = data.centering, + .bandwidth_sigma = bw_sigma + }; +} + template 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(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 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(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 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(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); diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 1a2535c1..ffcc08bd 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -163,6 +163,13 @@ class PixelRefine { double beam[2], double &dist_mm, double detector_rot[2], double latt_vec0[3], double latt_vec1[3], double latt_vec2[3]) const; + + // Bragg-prediction settings shared by Run/PredictImage/ChiSquaredImage so all + // three predict an identical reflection set. Critically it sets the radial + // Ewald-acceptance band (ewald_dist_cutoff) and the bandwidth broadening - see + // the definition for why these must track the partiality model, not the + // BraggPredictionSettings struct defaults. + BraggPredictionSettings BuildPredictionSettings(const PixelRefineData &data) const; public: PixelRefine(const DiffractionExperiment &experiment, const std::vector &reference);