From f4e281b2f5908f0551865619a8bb99124a1f8b45 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 9 Aug 2026 20:16:39 +0200 Subject: [PATCH] Rotation: give prediction and partiality the energy bandwidth The rotation predictor and RotationPartiality used the mosaicity alone. Energy bandwidth broadens a reflection's rocking curve as (dlambda/lambda)*tan(theta_B), resolution-dependent and negligible at low angle, so on a large-bandwidth beam the modelled reflecting range was too narrow exactly where the crystal still diffracts: 0.064 deg of broadening against a fitted 0.083, i.e. 26% at the detector edge. The stills predictor has carried the term since it was written; only rotation was missing it. Add it in the three places that have to agree. The predictor widens both its acceptance window and the partiality it hands to integration; the merge widens the partiality it recomputes from the smoothed mosaicity; and the per-image mosaicity fit subtracts the same term before fitting, so what it returns is the intrinsic mosaicity rather than the mosaicity plus the beam. Without that last part the bandwidth would be counted twice. The term goes in without the 1/zeta of the usual expression: the erf already divides by zeta, so adding a per-reflection width that itself carries 1/zeta would divide by it twice - up to 20x at the minimum zeta. dphi = delta*tan(theta_B), and the zeta stays where it was. The rotation identity dtheta/dphi = zeta was checked against a numerical solve of the diffraction condition at four resolutions and three orientations. Monochromatic data is untouched by construction - the term is guarded on a non-zero bandwidth and is an assignment, not arithmetic, when there is none. Verified: 246456 reflections byte-identical through the predictor, 4.7 million rocking-fraction evaluations with no bitwise difference, and identical merge tables end to end. The bandwidth is read from the file (incident_wavelength_spread) or from --bandwidth, and is absent from every dataset in the rotation battery. On the bandwidth dataset the fitted mosaicity becomes resolution-independent (0.0745 -> 0.0719 deg), the prediction window widens, frames per rocking event go 4.6 -> 5.3, per-image correlation to the merge rises 0.710 -> 0.725, and R_meas improves 0.1-0.7 pp in every shell while CC1/2 falls 0.6-0.9 pp in the outer two. Merged quality is net neutral: the combine normalises by sum(partiality), so a uniform widening largely cancels. Co-Authored-By: Claude Opus 5 (1M context) --- image_analysis/scale_merge/RotationScaleMerge.cpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) 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);