diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 7f3ce011..25a6427e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -11,6 +11,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Spot finding: `--spot-sigma` now defaults to 4.0 (was 3.0) and `--max-spots` to 1000 spots per image (was 250). * Resolution limits: Bragg integration, azimuthal integration and spot finding all default to **as far as the detector reaches**; `--integration-high-resolution` and `--spot-high-resolution` still set one by hand, and 0 means "no limit" at either end. * Bragg prediction: How far the predictor walks the lattice is a setting (`bragg_integration_settings.max_hkl`, `--max-hkl`) instead of a fixed 100, derived per crystal offline and held fixed online. +* Rotation data: a set X-ray bandwidth (`--bandwidth`) now widens each reflection's rocking curve by Δλ/λ·tan(θ) in prediction and partiality, and is deconvolved out of the fitted mosaicity; monochromatic data are unaffected. * Bragg integration: The local background ring is made robust with a **high-side sigma clip** (`--background-clip `, default 4) instead of the symmetric trimmed mean; the trim stays reachable with `--background-trim `. * Bragg integration: The **uncertainty of the background estimate** is now propagated into `sigma`, which both engines previously omitted. * Bragg integration: Profile-fit `sigma` is **no longer inflated on weak reflections** - the fit weights take the signal estimate as it is instead of clamping it at zero, and the per-pixel variance floor is 0.01 counts instead of 1/12. diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index 1c634ce3..b70344f4 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -25,9 +25,15 @@ struct BraggPredictionSettings { float min_zeta = 0.05; float mosaicity_multiplier = 4.0; // Relative X-ray bandwidth Δλ/λ expressed as a Gaussian sigma (0 = monochromatic). - // When > 0 the Ewald-shell acceptance is thickened radially per reflection by + // Stills: the Ewald-shell acceptance is thickened radially per reflection by // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no // longer clips high-resolution reflections. + // Rotation: differentiating Bragg's law at fixed d gives an extra rocking width + // Δθ = bandwidth_sigma·tan(θ_B), a spread in the same glancing angle the mosaic spread + // smears, so it adds to σ_M in quadrature. It is NOT divided by ζ: rotating the crystal + // by Δφ changes θ by ζ·Δφ, so the 1/ζ that turns an angular width into a rotation width is + // already the one the partiality applies to σ_M. CalcMosaicityXDS deconvolves the same term + // out of the fitted σ_M, so it is not counted twice. float bandwidth_sigma = 0.0f; }; diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index d7e6b0c2..3b202019 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -50,6 +50,16 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys const float mos_angle_rad = settings.mosaicity_deg * static_cast(PI) / 180.f; const float half_wedge_angle_rad = settings.wedge_deg * static_cast(PI) / 180.f / 2.0f ; + // Energy bandwidth widens the rocking curve. Differentiating Bragg's law at fixed d gives + // dtheta = (dlambda/lambda) tan(theta_B), a spread in the same glancing angle the mosaic spread + // smears, so it adds to sigma_M in quadrature. It is NOT divided by zeta here: rotating the + // crystal by dphi changes theta by zeta*dphi, so the 1/zeta that turns an angular width into a + // rotation width is already the one c1 (and the epsilon3 cutoff) applies to sigma_M. The fitted + // sigma_M has this term deconvolved out (CalcMosaicityXDS), so it is not counted twice. + // sin(theta_B) = lambda/(2d) = lambda*|p0|/2. Zero bandwidth leaves every reflection untouched. + const float bandwidth_sigma = settings.bandwidth_sigma; + const float half_wavelength_A = geom.GetWavelength_A() / 2.0f; + for (int h = -settings.max_h; h <= settings.max_h; h++) { // Precompute A* h contribution @@ -85,6 +95,15 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys if ((rho_sq < p_m3 * p_m3) || (p0_sq > 4 * S0 * S0)) continue; + // Effective rocking width for this reflection: mosaicity broadened by the bandwidth + // term. sin(theta_B) <= 1 is guaranteed by the p0_sq test just above. + float mos_eff_rad = mos_angle_rad; + if (bandwidth_sigma > 0.0f) { + const float sin_theta = half_wavelength_A * std::sqrt(p0_sq); + const float dphi_bw = bandwidth_sigma * sin_theta / std::sqrt(1.0f - sin_theta * sin_theta); + mos_eff_rad = std::sqrt(mos_angle_rad * mos_angle_rad + dphi_bw * dphi_bw); + } + for (const auto& p_m1 : p_m1_opt) { if (i >= max_reflections) continue; @@ -105,7 +124,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys float epsilon3 = std::fabs(phi * zeta_abs); - if (epsilon3 > settings.mosaicity_multiplier * mos_angle_rad) + if (epsilon3 > settings.mosaicity_multiplier * mos_eff_rad) continue; // Reciprocal Lorentz (Kabsch 2010): L^-1 = |m2 . (S x S0)| / (|S| |S0|) = @@ -114,7 +133,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys // corrupts the absolute/Wilson scale (it cancels within a resolution shell, so // CC1/2 / CCref / R-meas are neutral). const float lorentz_reciprocal = std::fabs(m2 * (S % S0)) / (S.Length() * S0.Length()); - const float c1 = zeta_abs / (std::sqrt(2.0f) * mos_angle_rad); + const float c1 = zeta_abs / (std::sqrt(2.0f) * mos_eff_rad); const float partiality = (std::erf((phi + half_wedge_angle_rad) * c1) - std::erf((phi - half_wedge_angle_rad) * c1)) / 2.0f; diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index f4d2188e..0c624751 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -86,6 +86,15 @@ namespace { float p_m1_pos = sqrtf(rho_sq - p_m3 * p_m3); float p_m1_arr[2] = {p_m1_pos, -p_m1_pos}; + // Effective rocking width: mosaicity broadened by the bandwidth term (consistent with CPU), + // dtheta = (dlambda/lambda) tan(theta_B) in quadrature with sigma_M, zeta-free. + float mos_eff_rad = C.mos_angle_rad; + if (C.bandwidth_sigma > 0.0f) { + float sin_theta = C.half_wavelength_A * sqrtf(p0_sq); + float dphi_bw = C.bandwidth_sigma * sin_theta / sqrtf(1.0f - sin_theta * sin_theta); + mos_eff_rad = sqrtf(C.mos_angle_rad * C.mos_angle_rad + dphi_bw * dphi_bw); + } + int count = 0; for (int idx = 0; idx < 2; ++idx) { float p_m1 = p_m1_arr[idx]; @@ -117,7 +126,7 @@ namespace { // epsilon3 cutoff check (consistent with CPU, Kabsch formulation) float epsilon3 = fabsf(phi * zeta_abs); - if (epsilon3 > C.mosaicity_multiplier * C.mos_angle_rad) + if (epsilon3 > C.mosaicity_multiplier * mos_eff_rad) continue; float cx, cy, cz; @@ -131,7 +140,7 @@ namespace { // Partiality calculation (Kabsch formulation) // c1 = sqrt(2) * sigma / zeta, where sigma = mosaicity - float c1 = zeta_abs / (sqrtf(2.0f) * C.mos_angle_rad); + float c1 = zeta_abs / (sqrtf(2.0f) * mos_eff_rad); float half_wedge = C.wedge_angle_rad / 2.0f; float partiality = (erff((phi + half_wedge) * c1) - erff((phi - half_wedge) * c1)) / 2.0f; @@ -218,6 +227,8 @@ namespace { kc.wedge_angle_rad = settings.wedge_deg * static_cast(PI) / 180.0f; kc.min_zeta = settings.min_zeta; kc.mosaicity_multiplier = settings.mosaicity_multiplier; + kc.bandwidth_sigma = settings.bandwidth_sigma; + kc.half_wavelength_A = geom.GetWavelength_A() / 2.0f; kc.Astar = lattice.Astar(); kc.Bstar = lattice.Bstar(); diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.h b/image_analysis/bragg_prediction/BraggPredictionRotGPU.h index f1429900..f58ec0ab 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.h @@ -20,6 +20,8 @@ struct KernelConstsRot { float wedge_angle_rad; float min_zeta; float mosaicity_multiplier; + float bandwidth_sigma; // relative dlambda/lambda as a sigma; 0 = monochromatic + float half_wavelength_A; Coord Astar, Bstar, Cstar, S0; Coord m1, m2, m3; float m2_S0; diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index e1befb77..23d82a9e 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -54,11 +54,17 @@ namespace { // delta_phi = oscillation range (radians) // sigma_M = mosaicity (radians) // zeta = |m2 · e1| Lorentz factor component - inline float R_fraction(float tau, float delta_phi, float sigma_M, float zeta) { + // sigma_bw = this reflection's energy-bandwidth rocking width (radians, zeta-free). It broadens the + // observed tau spread on top of the mosaicity, so carrying it here makes the fitted sigma_M the + // INTRINSIC width that prediction and scaling re-broaden per reflection - the same deconvolution + // FitProfileRadius does for the profile radius. 0 = monochromatic. + inline float R_fraction(float tau, float delta_phi, float sigma_M, float zeta, float sigma_bw) { if (zeta < 1e-6f || sigma_M < 1e-9f) return 0.0f; - const float sigma_eff = sigma_M / zeta; + const float sigma_total = sigma_bw > 0.0f + ? std::sqrt(sigma_M * sigma_M + sigma_bw * sigma_bw) : sigma_M; + const float sigma_eff = sigma_total / zeta; const float sqrt2_sigma = std::sqrt(2.0f) * sigma_eff; if (sqrt2_sigma < 1e-12f) @@ -74,11 +80,12 @@ namespace { // Returns sum of log(R) for all reflections inline double log_likelihood(const std::vector& tau_values, const std::vector& zeta_values, + const std::vector& sigma_bw_values, float delta_phi, float sigma_M) { double ll = 0.0; for (size_t i = 0; i < tau_values.size(); ++i) { - float R = R_fraction(tau_values[i], delta_phi, sigma_M, zeta_values[i]); + float R = R_fraction(tau_values[i], delta_phi, sigma_M, zeta_values[i], sigma_bw_values[i]); if (std::isfinite(R) && R > 1e-30f) { ll += std::log(static_cast(R)); } else { @@ -91,6 +98,7 @@ namespace { // Golden section search for maximum likelihood sigma_M inline float find_sigma_M_mle(const std::vector& tau_values, const std::vector& zeta_values, + const std::vector& sigma_bw_values, float delta_phi, float sigma_min_deg = 0.001f, float sigma_max_deg = 2.0f) { @@ -105,8 +113,8 @@ namespace { const float tol = 1e-6f; int iter = 0; while (std::fabs(b - a) > tol && iter++ < 100) { - double fc = log_likelihood(tau_values, zeta_values, delta_phi, c); - double fd = log_likelihood(tau_values, zeta_values, delta_phi, d); + double fc = log_likelihood(tau_values, zeta_values, sigma_bw_values, delta_phi, c); + double fd = log_likelihood(tau_values, zeta_values, sigma_bw_values, delta_phi, d); if (fc > fd) { b = d; @@ -237,6 +245,12 @@ namespace { const Coord S0 = experiment.GetScatteringVector(); const float delta_phi_rad = deg_to_rad(axis.GetWedge_deg()); + // Energy bandwidth adds sigma_bw = (dlambda/lambda)*tan(theta_B) to each reflection's rocking + // width; deconvolving it here leaves sigma_M the intrinsic mosaicity, so prediction and scaling + // can re-add it per reflection without counting it twice. sin(theta_B) = lambda*|pstar|/2. + const float bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f; + const float half_wavelength_A = experiment.GetWavelength_A() / 2.0f; + // Fit from the strongest spots only, never the whole list. A spot is detected when // I_full * R(tau) clears the finder's threshold, so the deeper the spot list reaches the more // large-|tau| partially-recorded spots it holds - and sigma_M is fitted from exactly that tau @@ -253,8 +267,10 @@ namespace { std::vector tau_values; std::vector zeta_values; + std::vector sigma_bw_values; tau_values.reserve(n_fit); zeta_values.reserve(n_fit); + sigma_bw_values.reserve(n_fit); for (size_t si = 0; si < n_fit; ++si) { const auto &s = spots[si]; @@ -301,15 +317,22 @@ namespace { if (!std::isfinite(zeta) || !std::isfinite(tau_rad) || zeta < 0.1f) continue; + float sigma_bw = 0.0f; + if (bandwidth_sigma > 0.0f) { + const float sin_theta = half_wavelength_A * pstar.Length(); + sigma_bw = bandwidth_sigma * sin_theta / std::sqrt(1.0f - sin_theta * sin_theta); + } + tau_values.push_back(tau_rad); zeta_values.push_back(zeta); + sigma_bw_values.push_back(sigma_bw); } if (tau_values.size() < 10) return std::nullopt; // Find sigma_M by maximizing log-likelihood - float sigma_M_rad = find_sigma_M_mle(tau_values, zeta_values, delta_phi_rad); + float sigma_M_rad = find_sigma_M_mle(tau_values, zeta_values, sigma_bw_values, delta_phi_rad); return rad_to_deg(sigma_M_rad); }