diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 63971bf2..933edb40 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -29,6 +29,7 @@ struct ReflGroup { int h, k, l; double d; double Itrue; // reference intensity (held fixed) + double R_bw_sq; // bandwidth radial-width^2 contribution (0 = monochromatic) double predicted_x, predicted_y; std::vector pixels; }; @@ -47,10 +48,9 @@ double SafeInv(double x, double fallback) { // I_pred(pixel) = G * Itrue * B_term * P_radial * P_tangential + I_bkg // // B_term = exp(-B |q|^2 / 4) (Debye-Waller) -// P_radial = exp(-eps_r^2 / R0^2) (partiality: fraction of +// P_radial = exp(-eps_r^2 / R0_eff^2) (partiality: fraction of // the mosaic blob on the -// Ewald sphere; NOT -// normalized, <= 1) +// Ewald sphere; <= 1) // P_tangential = A_recip/(pi R1^2) * exp(-eps_t^2/R1^2)(spatial profile on the // detector, normalized so // that sum over pixels ~ 1) @@ -59,17 +59,28 @@ double SafeInv(double x, double fallback) { // I_pred - I_bkg over the shoebox reproduces G * Itrue * B_term * P_radial. // The 1/(pi R1^2) normalization is the missing piece that decouples the profile // width R1 from the overall scale G. +// +// X-ray bandwidth: a spread in lambda is a spread in the Ewald-sphere radius, +// i.e. a purely *radial* thickening of the shell. It adds (in quadrature) a +// resolution-dependent term to the radial width: +// R0_eff^2 = R0^2 + R_bw^2 , R_bw^2 = (b*lambda)^2 / (2 d^4) +// where b = relative bandwidth (sigma of dlambda/lambda). R_bw grows like 1/d^2, +// so bandwidth leaves low-resolution spots sharp and smears high-resolution ones +// radially - the pink-beam/DMM signature. R_bw_sq is a fixed per-reflection +// constant (b is known), so R0 keeps meaning "intrinsic" width (mosaic + +// divergence + beam). b = 0 makes R_bw = 0: a monochromatic no-op. // --------------------------------------------------------------------------- struct PixelResidual { PixelResidual(const PixelObs &obs, double Itrue, double lambda, double pixel_size, double exp_h, double exp_k, double exp_l, + double R_bw_sq, gemmi::CrystalSystem symmetry) : Itrue(Itrue), Iobs(obs.Iobs), Ibkg(obs.Ibkg), weight(obs.weight), A_recip(obs.A_recip), obs_x(obs.x), obs_y(obs.y), inv_lambda(1.0 / lambda), pixel_size(pixel_size), exp_h(exp_h), exp_k(exp_k), exp_l(exp_l), - angle_rad(obs.angle_rad), symmetry(symmetry) { + R_bw_sq(R_bw_sq), angle_rad(obs.angle_rad), symmetry(symmetry) { if (std::fabs(lambda) < 1e-6) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Lambda cannot be close to zero"); @@ -251,10 +262,11 @@ struct PixelResidual { // Caveat: a still samples the radial direction at a single offset, so the // sqrt(pi) R0 normalization makes g_r a density (1/A^-1) rather than a // dimensionless fraction. The leftover dimensional factor is absorbed by - // the free scale G; completing it physically needs an effective radial - // sampling width from the energy bandwidth + beam divergence (TODO). - const T g_radial = ceres::exp(-eps_radial * eps_radial / (R[0] * R[0])) - / (ceres::sqrt(T(M_PI)) * R[0]); + // the free scale G. The energy-bandwidth contribution to the radial width + // is folded in here via R_bw_sq (beam divergence is still TODO). + const T R0_eff_sq = R[0] * R[0] + T(R_bw_sq); + const T g_radial = ceres::exp(-eps_radial * eps_radial / R0_eff_sq) + / (ceres::sqrt(T(M_PI)) * ceres::sqrt(R0_eff_sq)); const T P_tang = T(A_recip) * ceres::exp(-eps_tang_sq / (R[1] * R[1])) / (T(M_PI) * R[1] * R[1]); @@ -289,6 +301,7 @@ struct PixelResidual { const double inv_lambda; const double pixel_size; const double exp_h, exp_k, exp_l; + const double R_bw_sq; // bandwidth radial-width^2 contribution (0 = monochromatic) const double angle_rad; gemmi::CrystalSystem symmetry; }; @@ -346,6 +359,17 @@ void PixelRefine::Run(const T *image, return (qx % qy).Length(); }; + // Bandwidth radial-width^2 (in the code's R = sqrt(2)*sigma convention): + // R_bw^2 = (b*lambda)^2 / (2 d^4), b = relative bandwidth (sigma). + // A fixed per-reflection constant; data.bandwidth == 0 -> monochromatic no-op. + const double bw = data.bandwidth; + auto bandwidth_radial_sq = [&](double d) -> double { + if (bw <= 0.0 || d <= 0.0) + return 0.0; + const double bl = bw * lambda; + return bl * bl / (2.0 * d * d * d * d); + }; + // Mutable experiment whose geometry is re-synced from the refined data.geom // before each prediction, so shoeboxes track the refined geometry/cell. DiffractionExperiment exp_iter = experiment; @@ -384,6 +408,7 @@ void PixelRefine::Run(const T *image, g.l = refl.l; g.d = refl.d; g.Itrue = reference_data[hkl]; + g.R_bw_sq = bandwidth_radial_sq(refl.d); g.predicted_x = refl.predicted_x; g.predicted_y = refl.predicted_y; @@ -489,7 +514,7 @@ void PixelRefine::Run(const T *image, auto *cost = new ceres::AutoDiffCostFunction< PixelResidual, 1, 2, 1, 2, 3, 3, 3, 3, 1, 1, 2>( new PixelResidual(obs, g.Itrue, lambda, pixel_size, - g.h, g.k, g.l, data.crystal_system)); + g.h, g.k, g.l, g.R_bw_sq, data.crystal_system)); problem.AddResidualBlock(cost, new ceres::HuberLoss(3.0), beam, &dist_mm, detector_rot, rot_vec, latt_vec0, latt_vec1, latt_vec2, @@ -625,7 +650,7 @@ void PixelRefine::Run(const T *image, size_t n = 0; for (const auto &obs : g.pixels) { - PixelResidual pr(obs, 1.0, lambda, pixel_size, g.h, g.k, g.l, data.crystal_system); + PixelResidual pr(obs, 1.0, lambda, pixel_size, g.h, g.k, g.l, g.R_bw_sq, data.crystal_system); double q_sq, eps_r, eps_t_sq; if (!pr.GeometryTerms(beam, &dist_mm, detector_rot, rot_vec, latt_vec0, latt_vec1, latt_vec2, q_sq, eps_r, eps_t_sq)) @@ -637,7 +662,9 @@ void PixelRefine::Run(const T *image, const double P_t = obs.A_recip * std::exp(-eps_t_sq / (data.R[1] * data.R[1])) / (M_PI * data.R[1] * data.R[1]); // Peak-normalized radial factor (the partiality), in (0,1]. - const double P_radial = std::exp(-eps_r * eps_r / (data.R[0] * data.R[0])); + // Bandwidth-broadened radial width, matching the model in Model(). + const double R0_eff_sq = data.R[0] * data.R[0] + g.R_bw_sq; + const double P_radial = std::exp(-eps_r * eps_r / R0_eff_sq); const double v = SafeInv(obs.weight * obs.weight, 1.0); // pixel variance const double signal = obs.Iobs - obs.Ibkg; diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 44ae21c8..07b280cf 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -70,6 +70,30 @@ // as refinement proceeds it diverges, so later predictions must use the *refined* // data.geom rather than the static experiment geometry. // +// On shoeboxes, tails, and gatekeeping +// ------------------------------------ +// We deliberately do NOT chase the full spot. A shoebox only needs to cover +// enough of a reflection for the fit to be meaningful; clipped tails are fine, +// because the partiality term already downweights whatever falls outside the +// well-modelled core. The premise - especially for serial crystallography - is +// that the problem is rarely *missing* information; it is *failing to gate out* +// information that is not meaningful. As long as the model knows a piece is +// missing (low partiality), it is safe to leave it missing. That flips the usual +// trade-off: rather than shrinking boxes to avoid contamination, we can grow them +// and let partiality decide, per pixel and per reflection, what actually carries +// signal. (For downstream integration this pixel-level gating is the point - keep +// only meaningful pixels, instead of a fixed geometric mask.) The bandwidth term +// below is part of the same idea: it tells the model where the radial tails are +// *expected* to be, so it can weight rather than blindly include them. +// +// X-ray bandwidth (optional) +// -------------------------- +// A finite bandwidth thickens the Ewald shell radially and smears spots along the +// radial direction, growing like 1/d^2 (the pink-beam/DMM signature). It enters +// as a fixed, resolution-dependent addition to the radial width R0 (see +// PixelRefineData::bandwidth and PixelRefine.cpp). It is OFF by default +// (bandwidth = 0, monochromatic); set it for DMM-type data, leave it for Si. +// // Status: experimental prototype. The forward model (esp. the still-image // Lorentz/partiality normalization) is deliberately simple and expected to // evolve. See PixelRefine.cpp for the physics conventions and known caveats. @@ -86,6 +110,11 @@ struct PixelRefineData { double scale_factor = 1.0; // overall scale G double R[2] = {0.005, 0.005}; // R[0] = radial (partiality) width, R[1] = tangential (profile) width (A^-1) + // Relative X-ray bandwidth (sigma of dlambda/lambda), e.g. ~0.004 for a 1% + // FWHM DMM, ~1e-4 for Si(111). Adds a resolution-dependent radial broadening + // to R[0]. 0 = monochromatic (the term switches off entirely). + double bandwidth = 0.0; + // Goniometer: for a still image keep angle_deg = 0. For a wedge, pass the // angle (deg) of the slice centre relative to the reference (phi=0) frame. double angle_deg = 0.0;