mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases 67 Wiki Activity

PixelRefine: Add bandwidth contribution

Browse Source
This commit is contained in:
leonarski_f
2026-06-08 13:59:01 +02:00
parent 6f6098d508
commit 6f31159607
2 changed files with 67 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/pixel_refinement/PixelRefine.cpp
+38 -11
View File
@@ -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<PixelObs> 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;