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

beam centre: gate on what the geometry needs, and on the move being significant

Browse Source 
The accept bound was one flat 1.0 px constant applied at three points. It was calibrated on a
244 mm / 0.95 A / 0.15 mm geometry, and the displacement the J0 smearing law actually allows runs
from 0.6 px to 9 px over the in-house and non-SLS corpora - a fifteen-fold spread - so on a loose
geometry the constant is nine times too tight and throws away answers that are perfectly usable.
The bound is now the LARGER of the constant and what the geometry asks, which is one-sided by
construction: nothing today's gate accepts can be lost.

What the geometry asks is printed and never tested against. It says how wrong the header may be;
the estimator's sigma says how well the estimator knows its own answer, and gating one on the other
rejects a centre correct to 0.03 px because its sigma was 0.98. It is not a relevance floor either
- measured on real data, a 0.12 px change of centre, 0.03x of what the law asks across the spindle,
is the difference between the deposited cell and a halved axis - so the printed line says so.

Added alongside: the move has to be worth making. Under three times the estimator's own sigma it is
not a measurement of anything, and a centre that is not moved cannot push the two-pass loop off its
fixed point. On the 39 rotation regression crystals this adopts 38 and keeps the header on one.

Two things that were happening silently now say so: the move resolved across the spindle, which is
the only component the law is about, and the fall-through itself. Below about 220 deg of sweep the
spot symmetry never clears the bound and the background answers instead - correctly, to within
0.23 px of the full-sweep centre at every span from 120 to 220 deg - and until now the log did not
say the spot arm had even been tried.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f
This commit is contained in:
leonarski_f
2026-08-30 08:37:55 +02:00
co-authored by Claude Opus 5
parent d1818b195e
commit 52e9e9da2b
1 changed files with 90 additions and 7 deletions
Download Patch File Download Diff File Expand all files Collapse all files
rugnux/Rugnux.cpp
+90 -7
View File
@@ -85,6 +85,54 @@ namespace {
// out, and on none of the 38 as they are.
constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f;
// The longest cell axis assumed where no cell is known yet. It has to be LARGE: a long axis makes
// the tolerance below small, and the ceiling is the larger of it and the flat bound above, so a
// large default is the one that changes nothing.
constexpr float UNKNOWN_CELL_A_MAX_A = 200.0f;
// How far the beam centre may be wrong before the FIRST PASS stops indexing. A centre error is
// fixed in the LAB frame, so accumulating a sweep smears every reciprocal-lattice point around a
// circle and the FFT amplitude at an axis of length a is multiplied by J0(2 pi delta p a/(D
// lambda)); this is the displacement at which that has fallen to 0.70, and the first zero is
// 2.09x further out. Over the in-house and non-SLS corpora it runs from 0.6 px to 9 px - a
// fifteen-fold spread, which is why the flat bound above is only ever used as a FLOOR under it.
//
// It is printed, and it is deliberately not a relevance test. It says how wrong the header may
// be; the estimator's sigma says how well the estimator knows its own answer, and gating one on
// the other rejects a centre correct to 0.03 px because its sigma was 0.98. Nor is a move that is
// small against it a move that can be ignored: measured on real data, a 0.12 px change of centre
// - 0.03x of this - is the difference between the deposited cell and a halved axis.
float BeamCenterNeed_pxl(const DiffractionExperiment &experiment) {
const auto cell = experiment.GetUnitCell();
const float a_max = cell ? std::max({cell->a, cell->b, cell->c}) : UNKNOWN_CELL_A_MAX_A;
return 0.183f * experiment.GetDetectorDistance_mm() * experiment.GetWavelength_A()
/ (experiment.GetPixelSize_mm() * a_max);
}
// A detector-plane displacement resolved ACROSS the spindle - the only component the tolerance
// above is about. Along the spindle the error translates the derotated cloud rigidly and the FFT
// amplitude cannot see it at all; across it the peaks smear and the amplitude falls. (The lattice
// fit that follows the FFT is not translation-invariant, so "along the spindle" is free for the
// transform and not for the run - which is why this is reported and not used to decide anything.)
// The detector direction wanted is the one whose lab image is perpendicular to the spindle n:
// (ux e1 + uy e2).n = 0, so u is proportional to (-(e2.n), e1.n).
std::optional<float> AcrossSpindle_pxl(const DiffractionExperiment &experiment, float dx, float dy) {
const auto goniometer = experiment.GetGoniometer();
if (!goniometer || !goniometer->IsScanning())
return {};
const auto geom = experiment.GetDiffractionGeometry();
const float beam_x = geom.GetBeamX_pxl(), beam_y = geom.GetBeamY_pxl();
const Coord n = goniometer->GetAxis().Normalize();
const Coord origin = geom.LabCoord(beam_x, beam_y);
const Coord e1 = geom.LabCoord(beam_x + 1.0f, beam_y) - origin;
const Coord e2 = geom.LabCoord(beam_x, beam_y + 1.0f) - origin;
const float ux = -(e2 * n), uy = e1 * n;
const float length = std::hypot(ux, uy);
if (length < 1e-3f)
return {};
return std::abs(dx * ux + dy * uy) / length;
}
// Images the projection is built from when the beam centre is wanted but the beam-stop pre-pass
// is off; with it on, that pre-pass's image count is used and one projection serves both.
constexpr int BEAM_CENTER_PROJECTION_IMAGES = 60;
@@ -815,6 +863,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
std::string source = "spot symmetry";
std::optional<BeamCenterEstimate> estimate;
SpindleEstimate spindle;
// The bound an estimate has to meet. One-sided by construction: the flat constant was calibrated
// on a 244 mm / 0.95 A / 0.15 mm geometry and is far too tight for a loose one, so where the
// geometry asks less of the centre the bound is relaxed to what it asks and never below the
// constant. Nothing it accepts today is lost.
const float need = BeamCenterNeed_pxl(experiment_);
const float ceiling = std::max(MAX_BEAM_CENTER_SIGMA_PXL, need);
if (want_spot_symmetry)
estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, beam_center_spots,
config_.fit_spindle ? &spindle : nullptr);
@@ -827,7 +881,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
// second: this runs before the fall-through, and only what survives it falls through. The extra
// frames go to the spot finder alone, so the beam-stop projection keeps the images it was
// validated on.
if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL)) {
if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > ceiling)) {
// What has already been read FOR THE SPOTS, which is not the whole sample: the frames the
// shadow was built from carry none of them, and testing against the union would silently
// skip the ones this pass exists to read.
@@ -899,7 +953,18 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
"(Friedel vote {:.0f} -> {:.0f})",
1e3 * spindle.azimuth_rad, 1e3 * spindle.tip_rad,
spindle.vote_excess_nominal, spindle.vote_excess);
if (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL) {
if (!estimate || estimate->sigma_pxl > ceiling) {
// Say so. Below about 220 deg of sweep the spot symmetry never clears this bound and the
// background answers instead - measured within 0.23 px of the full-sweep centre at every span
// from 120 to 220 deg - so the fall-through is what covers short sweeps, and until this line
// existed nothing in the log said the spot arm had even been tried.
if (want_spot_symmetry && estimate)
logger.Info("Beam centre: the spot symmetry answers ({:.2f},{:.2f}) at sigma {:.2f} px, "
"over the {:.2f} px ceiling - falling through to the background",
estimate->beam_x_pxl, estimate->beam_y_pxl, estimate->sigma_pxl, ceiling);
else if (want_spot_symmetry)
logger.Info("Beam centre: the spot symmetry does not come out on this sweep - falling "
"through to the background");
source = "background";
estimate = FindBeamCenterFromBackground(experiment_, pixel_mask_, finder.GetMeanProjection());
}
@@ -908,14 +973,32 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
return;
}
const float moved = std::hypot(estimate->beam_x_pxl - experiment_.GetBeamX_pxl(),
estimate->beam_y_pxl - experiment_.GetBeamY_pxl());
const bool commit = estimate->sigma_pxl <= MAX_BEAM_CENTER_SIGMA_PXL;
const float dx = estimate->beam_x_pxl - experiment_.GetBeamX_pxl();
const float dy = estimate->beam_y_pxl - experiment_.GetBeamY_pxl();
const float moved = std::hypot(dx, dy);
const auto across = AcrossSpindle_pxl(experiment_, dx, dy);
// Two conditions, and the second is the new one. A move under three times the estimator's own
// sigma is not a measurement of anything, and a centre that is not moved cannot move the
// two-pass loop off its fixed point either. On the 39 rotation regression crystals this adopts
// 38 and keeps the header on one, a 0.74 px move at sigma 0.32.
const bool commit = estimate->sigma_pxl <= ceiling && moved > 3.0f * estimate->sigma_pxl;
const char *verdict = "COMMIT";
if (!commit)
verdict = estimate->sigma_pxl > ceiling ? "reject: sigma over the ceiling (kept header)"
: "reject: the move is under 3 sigma (kept header)";
logger.Info("Beam centre from {}: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved {:.2f} px, "
"sigma {:.2f} px => {}",
"sigma {:.2f} px against a {:.2f} px ceiling => {}",
source, experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
estimate->beam_x_pxl, estimate->beam_y_pxl, moved, estimate->sigma_pxl,
commit ? "COMMIT" : "reject (kept header)");
ceiling, verdict);
if (across)
logger.Info("Beam centre: {:.2f} px of that move is across the spindle, where this geometry "
"asks for {:.2f} px. That is a floor and not a tolerance - a move 0.03x of it has "
"decided a halved axis", *across, need);
if (across && *across > need)
logger.Warning("Beam centre: the file's centre is {:.2f} px out across the spindle, more than "
"the {:.2f} px this geometry absorbs - expect the first pass to lose peaks or "
"to take an axis harmonic", *across, need);
if (commit)
experiment_.BeamX_pxl(estimate->beam_x_pxl).BeamY_pxl(estimate->beam_y_pxl);
}