leonarski_fandClaude Opus 5.5
698d06a0e2
zeta >= min_zeta, so the rocking-curve test rejects every solution with min_zeta * (|phi| - half wedge) > multiplier * mosaicity. Tested right after phi, with a 1e-5 rad margin (float rounding of the full test is below 1e-6), it skips the cross product, normalisation and zeta of the 85-95% of solutions no frame keeps, and cannot reject one the full test would pass. CPU build: on the order of 10 s of a 1.2 A sweep's prediction. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
241 lines
13 KiB
C++
241 lines
13 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "../../common/JFJochMath.h"
|
|
#include "BraggPredictionRot.h"
|
|
#include "../SensorAbsorption.h"
|
|
#include "../bragg_integration/SystematicAbsence.h"
|
|
#include "RockingSlice.h"
|
|
|
|
|
|
int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const CrystalLattice &lattice,
|
|
const BraggPredictionSettings &settings) {
|
|
|
|
const auto geom = experiment.GetDiffractionGeometry();
|
|
const auto det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
|
|
const auto det_height_pxl = static_cast<float>(experiment.GetYPixelsNum());
|
|
|
|
const float one_over_dmax = 1.0f / settings.high_res_A;
|
|
const float one_over_dmax_sq = one_over_dmax * one_over_dmax;
|
|
|
|
float one_over_wavelength = 1.0f / geom.GetWavelength_A();
|
|
|
|
const Coord Astar = lattice.Astar();
|
|
const Coord Bstar = lattice.Bstar();
|
|
const Coord Cstar = lattice.Cstar();
|
|
const Coord S0 = geom.GetScatteringVector();
|
|
|
|
std::vector<float> rot = geom.GetDetectorMatrix().transpose().arr();
|
|
|
|
// Precompute detector geometry constants
|
|
float beam_x = geom.GetBeamX_pxl();
|
|
float beam_y = geom.GetBeamY_pxl();
|
|
float det_distance = geom.GetDetectorDistance_mm();
|
|
float pixel_size = geom.GetPixelSize_mm();
|
|
float F = det_distance / pixel_size;
|
|
|
|
const auto gon_opt = experiment.GetGoniometer();
|
|
if (!gon_opt.has_value())
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"BraggPredictionRotationCPU requires a goniometer axis");
|
|
const GoniometerAxis& gon = *gon_opt;
|
|
|
|
const Coord m2 = gon.GetAxis().Normalize();
|
|
const Coord m1 = (m2 % S0).Normalize();
|
|
const Coord m3 = (m1 % m2).Normalize();
|
|
|
|
const float m2_S0 = m2 * S0;
|
|
const float m3_S0 = m3 * S0;
|
|
|
|
int i = 0;
|
|
const Coord beam_dir = S0.Normalize();
|
|
|
|
const float mos_angle_rad = settings.mosaicity_deg * static_cast<float>(PI) / 180.f;
|
|
const float half_wedge_angle_rad = settings.wedge_deg * static_cast<float>(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;
|
|
|
|
// Angle-dependent sensor efficiency. Per-dataset constants (thickness, material, wavelength)
|
|
// collapse to two numbers here; the per-reflection part is one exponential below. Inert - and
|
|
// bit-identical to not applying it - wherever the sensor is opaque, which is every long
|
|
// wavelength, so it needs no flag and no threshold anyone has to choose.
|
|
const auto &det = experiment.GetDetectorSetup();
|
|
const auto sensor_qe = sensor_absorption::SensorQE::Build(
|
|
det.GetSensorMaterial(), det.GetSensorThickness_um(), geom.GetWavelength_A());
|
|
// The air in the sample-to-pixel flight path carries the same cos(alpha) dependence, with the
|
|
// opposite sign. See sensor_absorption::FlightPathAttenuation; off leaves d_over_L at 0, which is
|
|
// bit-identical to not applying it.
|
|
const auto air = sensor_absorption::FlightPathAttenuation::Build(
|
|
experiment.GetBraggIntegrationSettings().GetFlightPath(), geom.GetDetectorDistance_mm(),
|
|
geom.GetWavelength_A());
|
|
|
|
for (int h = -settings.max_h; h <= settings.max_h; h++) {
|
|
// Precompute A* h contribution
|
|
|
|
for (int k = -settings.max_k; k <= settings.max_k; k++) {
|
|
// Accumulate B* k contribution
|
|
|
|
for (int l = -settings.max_l; l <= settings.max_l; l++) {
|
|
if (systematic_absence(h, k, l, settings.centering))
|
|
continue;
|
|
|
|
if (i >= max_reflections)
|
|
continue;
|
|
Coord p0 = Astar * h + Bstar * k + Cstar * l;
|
|
|
|
float p0_sq = p0 * p0;
|
|
if (p0_sq <= 0.0f || p0_sq > one_over_dmax_sq)
|
|
continue;
|
|
|
|
const float p0_m1 = p0 * m1;
|
|
const float p0_m2 = p0 * m2;
|
|
const float p0_m3 = p0 * m3;
|
|
|
|
const float rho_sq = p0_sq - (p0_m2 * p0_m2);
|
|
|
|
const float p_m3 = (- p0_sq / 2 - p0_m2 * m2_S0) / m3_S0;
|
|
const float p_m2 = p0_m2;
|
|
const float p_m1_opt[2] = {
|
|
std::sqrt(rho_sq - p_m3 * p_m3),
|
|
-std::sqrt(rho_sq - p_m3 * p_m3)
|
|
};
|
|
|
|
// No solution for Laue equations
|
|
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;
|
|
|
|
const float cosphi = (p_m1 * p0_m1 + p_m3 * p0_m3) / rho_sq;
|
|
const float sinphi = (p_m1 * p0_m3 - p_m3 * p0_m1) / rho_sq;
|
|
float phi = -1.0f * std::atan2(sinphi, cosphi);
|
|
|
|
// Most solutions are far from this frame. zeta >= min_zeta, so the rocking-curve test
|
|
// below (epsilon3) rejects every solution with min_zeta * (|phi| - half wedge) above the
|
|
// limit; testing that first, with a margin well above float rounding, skips the rest
|
|
// of the arithmetic for them and changes nothing about the ones that pass.
|
|
if (settings.min_zeta > 0.0f
|
|
&& settings.min_zeta * (std::fabs(static_cast<double>(phi)) - half_wedge_angle_rad)
|
|
> settings.mosaicity_multiplier * static_cast<double>(mos_eff_rad) + 1e-5)
|
|
continue;
|
|
|
|
Coord p = m1 * p_m1 + m2 * p_m2 + m3 * p_m3; // p0 vector "rotated" to diffracting condition
|
|
Coord S = S0 + p;
|
|
|
|
const Coord e1 = (S % S0).Normalize();
|
|
|
|
const float zeta_abs = std::fabs(m2 * e1);
|
|
|
|
if (zeta_abs < settings.min_zeta)
|
|
continue;
|
|
|
|
// Is any of this reflection's rocking curve inside THIS image's oscillation range?
|
|
// phi is the offset from the frame's mid-exposure angle to the exact diffracting
|
|
// condition, so the curve - a Gaussian of width mos_eff/zeta in phi - has to be
|
|
// measured against the frame's EDGE, not its centre. Testing |phi|*zeta alone
|
|
// asks whether the frame centre is within the curve, which is a different and
|
|
// stricter question: it rejects a reflection whose curve overlaps the exposure
|
|
// but whose exact condition falls outside it. Since consecutive frame centres are
|
|
// one wedge apart, the nearest centre can be half a wedge away, so a reflection is
|
|
// rejected on EVERY frame - lost entirely, not merely clipped - once
|
|
// mos_eff < zeta * wedge / (2 * multiplier). That is coarse slicing on a sharp
|
|
// crystal at high zeta, where the reflection is fully recorded and cleanest.
|
|
float epsilon3 = std::fabs(phi * zeta_abs) - half_wedge_angle_rad * zeta_abs;
|
|
|
|
if (epsilon3 > settings.mosaicity_multiplier * mos_eff_rad)
|
|
continue;
|
|
|
|
// Reciprocal Lorentz (Kabsch 2010): L^-1 = |m2 . (S x S0)| / (|S| |S0|) =
|
|
// |zeta * sin angle(S,S0)|. The original divided by the scalar product
|
|
// S.S0 = |S||S0|cos(2theta), adding a spurious 1/cos(2theta) (1.8x at 1 A) that
|
|
// 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_eff_rad);
|
|
|
|
const float partiality = (std::erf((phi + half_wedge_angle_rad) * c1)
|
|
- std::erf((phi - half_wedge_angle_rad) * c1)) / 2.0f;
|
|
// Place the partial where this frame recorded it: at the centre of the frame's slice of
|
|
// the rocking curve, which is the exact-condition position walked along its Debye ring by
|
|
// the rotation between the two (RockingSlice.h). Walking along the ring is turning S
|
|
// about the beam, by the rotation's component along the ring over the ring's radius.
|
|
const float c_slice = RockingSliceCentroid_rad(phi, half_wedge_angle_rad, c1, partiality);
|
|
const Coord S_par = beam_dir * (beam_dir * S), S_perp = S - S_par;
|
|
const Coord S_turn = beam_dir % S; // |S_turn| = |S_perp|
|
|
const float psi = c_slice * ((m2 % p) * S_turn) / (S_perp * S_perp);
|
|
S = S_par + S_perp * std::cos(psi) + S_turn * std::sin(psi);
|
|
|
|
// Inlined RecipToDetector: the full transposed detector matrix, tilt and discrete orientation
|
|
// Apply rotation matrix transpose
|
|
float S_rot_x = rot[0] * S.x + rot[1] * S.y + rot[2] * S.z;
|
|
float S_rot_y = rot[3] * S.x + rot[4] * S.y + rot[5] * S.z;
|
|
float S_rot_z = rot[6] * S.x + rot[7] * S.y + rot[8] * S.z;
|
|
|
|
if (S_rot_z <= 0)
|
|
continue;
|
|
|
|
float x = beam_x + F * S_rot_x / S_rot_z;
|
|
float y = beam_y + F * S_rot_y / S_rot_z;
|
|
|
|
if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl))
|
|
continue;
|
|
|
|
float dist_ewald_sphere = std::fabs(S.Length() - one_over_wavelength);
|
|
|
|
// Sensor quantum efficiency at this reflection's own angle of incidence on the
|
|
// detector. The angle is taken against the DETECTOR NORMAL - S_rot is the
|
|
// diffracted direction in the detector's own frame, so its z component over its
|
|
// length is that cosine already, at no cost. Taking it here rather than from the
|
|
// resolution is what makes it right on a tilted detector, where the incidence
|
|
// angle stops being a function of resolution and the correction stops
|
|
// cancelling within a resolution shell.
|
|
const float cos_alpha = S_rot_z / S.Length();
|
|
const float qe_corr = sensor_qe.Factor(cos_alpha);
|
|
const float flight_corr = air.Factor(cos_alpha);
|
|
|
|
float d = 1.0f / sqrtf(p0_sq);
|
|
reflections[i] = Reflection{
|
|
.h = h,
|
|
.k = k,
|
|
.l = l,
|
|
.delta_phi_deg = phi * 180.0f / static_cast<float>(PI),
|
|
.predicted_x = x,
|
|
.predicted_y = y,
|
|
.observed_x = NAN,
|
|
.observed_y = NAN,
|
|
.d = d,
|
|
.dist_ewald = dist_ewald_sphere,
|
|
.prescaling_corr = lorentz_reciprocal,
|
|
.qe_corr = qe_corr,
|
|
.flight_corr = flight_corr,
|
|
.partiality = partiality,
|
|
.zeta = zeta_abs,
|
|
.image_scale_corr = lorentz_reciprocal * qe_corr * flight_corr / partiality,
|
|
};
|
|
i++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return TruncateToOutput(i);
|
|
}
|