diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 52e4ba05..14748a65 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -594,10 +594,21 @@ void RotationScaleMerge::SmoothMosaicityAndPartiality() { // Recompute each partial's partiality from the smoothed mosaicity (same wedge the predictor used). // Frames without a mosaicity keep the stored partiality. const double wedge = gon ? std::fabs(gon->GetWedge_deg()) : 0.0; + // Energy bandwidth broadens each reflection's rocking curve by (dlambda/lambda)*tan(theta_B) on top + // of the mosaicity - the same term the predictor adds and the mosaicity fit takes out, so the + // partiality recomputed here matches the one integration was based on. sin(theta_B) = lambda/(2d). + const double bandwidth_sigma = x.GetBandwidthFWHM().value_or(0.0f) / 2.3548; + const double half_wavelength_A = x.GetWavelength_A() / 2.0; ParallelChunks(static_cast(partials.size()), nthreads, [&](int lo, int hi) { for (int i = lo; i < hi; ++i) { auto &o = partials[i]; - const float mos = mos_smooth[o.frame]; + float mos = mos_smooth[o.frame]; + if (bandwidth_sigma > 0.0 && o.d > 0.0f) { + const double sin_theta = half_wavelength_A / o.d; + const double sigma_bw = (180.0 / PI) * bandwidth_sigma * sin_theta + / std::sqrt(1.0 - sin_theta * sin_theta); + mos = static_cast(std::sqrt(mos * mos + sigma_bw * sigma_bw)); + } if (std::isfinite(mos) && mos > 1e-6f && std::isfinite(o.zeta) && o.zeta > 0.0f && std::isfinite(o.delta_phi)) o.partiality = RotationPartiality(o.delta_phi, o.zeta, mos, wedge);