diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index c297de04..98b1aed6 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -289,7 +289,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .max_hkl = 100, .centering = outcome.symmetry.centering, .wedge_deg = std::fabs(wedge_deg), - .mosaicity_deg = std::fabs(mos_deg) + .mosaicity_deg = std::fabs(mos_deg), + // FWHM -> sigma; 0 when monochromatic, leaving the prediction unchanged. + .bandwidth_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f }; // Select the integration path: classical 2D integration, or the experimental diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 962ae45d..4c134e38 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -4,6 +4,13 @@ #include "BraggPrediction.h" #include "../bragg_integration/SystematicAbsence.h" +namespace { + // Number of bandwidth sigmas included in the (radially thickened) Ewald-shell + // acceptance window. 3σ captures essentially the whole pink-beam smear; matches + // the conservative end of the mosaicity cutoff used by callers. + constexpr float kBandwidthCutoffSigmas = 3.0f; +} + BraggPrediction::BraggPrediction(int max_reflections) : max_reflections(max_reflections), reflections(max_reflections) {} @@ -76,7 +83,18 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal float S_len = sqrtf(S_x * S_x + S_y * S_y + S_z * S_z); float dist_ewald_sphere = std::fabs(S_len - one_over_wavelength); - if (dist_ewald_sphere <= settings.ewald_dist_cutoff ) { + // Energy bandwidth thickens the Ewald shell radially: at the + // diffraction condition |S|-1/λ shifts by recip_z·(Δλ/λ), i.e. + // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d², identical to + // PixelRefine's R_bw). Broaden the acceptance window in quadrature so + // high-resolution shells (smeared most, ∝1/d²) are not clipped. + float radial_cutoff = settings.ewald_dist_cutoff; + if (settings.bandwidth_sigma > 0.0f) { + const float bw_tol = kBandwidthCutoffSigmas * settings.bandwidth_sigma * std::fabs(recip_z); + radial_cutoff = std::sqrt(radial_cutoff * radial_cutoff + bw_tol * bw_tol); + } + + if (dist_ewald_sphere <= radial_cutoff ) { const float s0_p0 = S0.x * recip_x + S0.y * recip_y + S0.z * recip_z; const float val = s0_sq * recip_sq - s0_p0 * s0_p0; diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index 6aaf4fc5..13d1c7ed 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -18,6 +18,12 @@ struct BraggPredictionSettings { float mosaicity_deg = 0.2f; 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 + // σ_bw = |recip_z|·bandwidth_sigma (= bλ/2d²), so the 1/d² pink-beam smear no + // longer clips high-resolution reflections. Must stay the last member so existing + // designated initializers remain valid. + float bandwidth_sigma = 0.0f; }; class BraggPrediction { diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index b15c34bc..cf8f8cc6 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -8,6 +8,10 @@ #include namespace { + // Number of bandwidth sigmas included in the (radially thickened) Ewald-shell + // acceptance window. Mirrors the CPU BraggPrediction path. + constexpr float kBandwidthCutoffSigmas = 3.0f; + __device__ inline bool is_odd(int v) { return (v & 1) != 0; } __device__ inline float angle_from_ewald_sphere_deg(const Coord &S0, float recip_x, float recip_y, float recip_z, float recip_sq) { @@ -95,7 +99,14 @@ namespace { float Sz = recip_z + C.S0.z; float S_len = sqrtf(Sx * Sx + Sy * Sy + Sz * Sz); float dist_ewald = fabsf(S_len - C.one_over_wavelength); - if (dist_ewald > C.ewald_cutoff) return false; + // Energy bandwidth thickens the Ewald shell radially: σ_bw = |recip_z|·(Δλ/λ) + // (= bλ/2d²). Broaden the acceptance window in quadrature (see CPU path). + float radial_cutoff = C.ewald_cutoff; + if (C.bandwidth_sigma > 0.0f) { + const float bw_tol = kBandwidthCutoffSigmas * C.bandwidth_sigma * fabsf(recip_z); + radial_cutoff = sqrtf(radial_cutoff * radial_cutoff + bw_tol * bw_tol); + } + if (dist_ewald > radial_cutoff) return false; float Srx = C.rot[0] * Sx + C.rot[1] * Sy + C.rot[2] * Sz; float Sry = C.rot[3] * Sx + C.rot[4] * Sy + C.rot[5] * Sz; float Srz = C.rot[6] * Sx + C.rot[7] * Sy + C.rot[8] * Sz; @@ -145,7 +156,8 @@ namespace { const CrystalLattice &lattice, float high_res_A, float ewald_dist_cutoff, - char centering) { + char centering, + float bandwidth_sigma) { KernelConsts kc{}; auto geom = experiment.GetDiffractionGeometry(); kc.det_width_pxl = static_cast(experiment.GetXPixelsNum()); @@ -157,6 +169,7 @@ namespace { kc.one_over_dmax_sq = one_over_dmax * one_over_dmax; kc.one_over_wavelength = 1.0f / geom.GetWavelength_A(); kc.ewald_cutoff = ewald_dist_cutoff; + kc.bandwidth_sigma = bandwidth_sigma; kc.Astar = lattice.Astar(); kc.Bstar = lattice.Bstar(); kc.Cstar = lattice.Cstar(); @@ -178,7 +191,8 @@ int BraggPredictionGPU::Calc(const DiffractionExperiment &experiment, const CrystalLattice &lattice, const BraggPredictionSettings &settings) { // Build constants on host - KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, settings.centering); + KernelConsts hK = BuildKernelConsts(experiment, lattice, settings.high_res_A, settings.ewald_dist_cutoff, + settings.centering, settings.bandwidth_sigma); cudaMemcpyAsync(dK, &hK, sizeof(KernelConsts), cudaMemcpyHostToDevice, stream); cudaMemsetAsync(d_count, 0, sizeof(int), stream); diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.h b/image_analysis/bragg_prediction/BraggPredictionGPU.h index 2bf7a27b..8c6a8486 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.h @@ -18,6 +18,7 @@ struct KernelConsts { float one_over_wavelength; float one_over_dmax_sq; float ewald_cutoff; + float bandwidth_sigma; // relative Δλ/λ (sigma); 0 = monochromatic Coord Astar, Bstar, Cstar, S0; float rot[9]; char centering; diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index b2d2578f..8f78e16c 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -475,7 +475,8 @@ void PixelRefine::Run(const T *image, const BraggPredictionSettings settings_prediction{ .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), .max_hkl = 100, - .centering = data.centering + .centering = data.centering, + .bandwidth_sigma = static_cast(data.bandwidth) // relative Δλ/λ sigma }; const auto azim_result = profile.GetResult(); @@ -942,7 +943,8 @@ std::vector PixelRefine::PredictImage(const AzimuthalIntegrationProfile & const BraggPredictionSettings settings_prediction{ .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), .max_hkl = 100, - .centering = data.centering + .centering = data.centering, + .bandwidth_sigma = static_cast(data.bandwidth) // relative Δλ/λ sigma }; const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction); const auto &predicted = prediction.GetReflections();