rugnux: index the measured beam centre too, and say when the two answers differ

The estimate the previous commit makes is free, and a first pass re-uses spots it has already found,
so running the first pass a second time at the measured centre costs about what one rung of
--beam-center-search costs. That buys the one comparison an on-failure trigger structurally cannot
make. A centre wrong ACROSS the spindle announces itself - the indexed fraction collapses - but a
centre wrong ALONG it does not: the run holds 96-100 % of frames indexed and quietly returns a 2x,
3x or sqrt(3) axis harmonic. Nothing fails, so nothing fires. Indexing both centres and comparing
the two lattices is what can see that at all.

What must not arbitrate the two is the indexed frame count. Acceptance is a fractional-Miller test,
so a cell twice as long has to place every spot twice as accurately to score the same; measured on
real data, two centres 0.12 px apart gave the deposited cell at 99.23 % and a halved axis at
100.00 %, and the wrong answer indexed better. A rule of the form "take the centre that indexes
more" picks wrong in exactly the case the comparison exists for.

So only what needs no arbiter is decided. The file's centre indexing nothing where the measured one
indexes a majority is not a comparison, it is a run that produced nothing and now does: take the
measured centre. The two agreeing on the lattice is reported and nothing else - it is a free
statement that the header is good enough for this crystal, which is most runs and is worth saying.
A disagreement between two passes that both worked is reported with both cells, their primitive
volumes and the ratio, flagged when that ratio is an axis harmonic, and left undecided.

Runs before the blind ladder, so a measured hypothesis is tried before any grid, and only where the
move exceeds three times the fit's own sigma - below that the two centres are one measurement twice.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f
This commit is contained in:
2026-08-30 09:00:55 +02:00
co-authored by Claude Opus 5
parent 0644404286
commit 60f6a46b25
2 changed files with 134 additions and 7 deletions
+127 -5
View File
@@ -896,11 +896,15 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
logger.Info("Beam centre check: that is under three times the fit's own sigma, so the "
"file's centre is as good as this measurement can tell");
else if (across && *across > need)
logger.Warning("Beam centre check: {:.2f} px of the difference is ACROSS the spindle, "
"more than the {:.2f} px this geometry absorbs. A centre error is fixed "
"in the lab frame, so it smears the accumulated reciprocal-space cloud "
"and the first pass loses peaks or takes an axis harmonic - expect the "
"cell below to be wrong", *across, need);
// Info and not a warning: `need` is where the FFT amplitude at a 200 A axis has fallen
// to 0.70, which is a floor under what matters and not a threshold for failure - it is
// exceeded on two of three ordinary in-house sweeps whose lattice comes out right. What
// is worth a warning is what the second pass below then finds, not this.
logger.Info("Beam centre check: {:.2f} px of the difference is ACROSS the spindle, more "
"than the {:.2f} px this geometry absorbs, and a centre error is fixed in "
"the lab frame - it smears the accumulated reciprocal-space cloud, so the "
"first pass can lose peaks or take an axis harmonic. The second first pass "
"below indexes both centres and reports whether it did", *across, need);
else
logger.Info("Beam centre check: that is a real difference, and {} - the second first "
"pass below indexes both centres and compares what they give",
@@ -2175,6 +2179,124 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
}
}
// The background centre, tried as a SECOND FIRST PASS.
//
// The estimate itself is free - it is read off the projection --detect-beam-stop has already
// built - and a first pass re-uses the spots it has already found, so indexing the run a
// second time at the measured centre costs about what one rung of the ladder below costs.
// That buys the comparison the on-failure ladder structurally cannot make: a centre wrong
// ALONG the spindle does not fail. It keeps 96-100 % of frames indexed and quietly returns a
// 2x, 3x or sqrt(3) axis harmonic, so no failure trigger ever fires on it. Indexing both
// centres and comparing the two ANSWERS is what can see it.
//
// What must NOT arbitrate the two is the frame count. Acceptance is a fractional-Miller test,
// so a cell twice as long has to place every spot twice as accurately to score the same:
// measured on real data, two centres 0.12 px apart gave the deposited cell at 99.23 % and a
// halved axis at 100.00 %, and the WRONG answer indexed better. A rule of the form "take the
// centre that indexes more" would therefore pick wrong in exactly the case this comparison
// exists for. So only what needs no arbiter is decided here:
//
// * the file's centre indexes nothing and the measured one indexes a majority -> take it.
// There is no comparison to get wrong, only a run that produced nothing and now does;
// * the two agree on the lattice -> keep the file's centre, and say they agree. This is the
// overwhelming majority, and saying it is worth the pass: it is a free statement that the
// header is good enough for this crystal, which nothing else in the run makes;
// * they disagree -> say so, with both cells and their primitive-volume ratio, and change
// nothing. Naming a disagreement is useful; resolving it needs an arbiter nobody has.
if (!cancelled_ && config_.beam_center_check && background_center_) {
// Consumed, not re-read: a two-pass rotation run comes back through here with the geometry
// the first pass post-refined, and a pre-scan estimate of the FILE's centre says nothing
// about that one. The pre-scan does not re-measure it either - by then the beam-stop shadow
// is already in the mask, so the projection is not rebuilt.
const BeamCenterEstimate measured = *background_center_;
background_center_.reset();
const float header_x = experiment_.GetBeamX_pxl(), header_y = experiment_.GetBeamY_pxl();
const float measured_x = measured.beam_x_pxl;
const float measured_y = measured.beam_y_pxl;
const float moved = std::hypot(measured_x - header_x, measured_y - header_y);
// A move the fit cannot resolve is not a different hypothesis, only the same one twice.
if (moved > 3.0f * measured.sigma_pxl) {
const int majority = static_cast<int>(validation.size()) / 2;
experiment_.BeamX_pxl(measured_x).BeamY_pxl(measured_y);
const FirstPass alt = pick_best(*indexer_pool, *indexer);
const bool header_indexes = best.result.has_value() && best.score > majority;
const bool measured_indexes = alt.result.has_value() && alt.score > majority;
if (measured_indexes && !header_indexes) {
logger.Warning("Beam centre from the background: the file's ({:.2f},{:.2f}) indexes "
"{}/{} validation frames and the measured ({:.2f},{:.2f}) indexes "
"{}/{} - the beam centre in this file is wrong, and the run adopts "
"the measured one", header_x, header_y, best.score,
static_cast<int>(validation.size()), measured_x, measured_y,
alt.score, static_cast<int>(validation.size()));
best = alt;
} else {
experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y);
if (header_indexes && measured_indexes) {
// Primitive volumes, as everywhere else here: two passes can find the same
// lattice in different settings, and a centred setting's cell is an exact
// integer multiple of its primitive one.
const double ratio = alt.vol < best.vol ? best.vol / alt.vol : alt.vol / best.vol;
const auto &hc = best.result->search_result.conventional.GetUnitCell();
const auto &mc = alt.result->search_result.conventional.GetUnitCell();
const bool same_class =
best.result->search_result.system == alt.result->search_result.system
&& best.result->search_result.centering == alt.result->search_result.centering;
if (same_class && ratio < 1.02) {
logger.Info("Beam centre check: indexing at the measured centre "
"({:.2f},{:.2f}) returns the same lattice as the file's - "
"{}-centred {}, cells agreeing to {:.2f} % in volume. The {:.2f} px "
"difference does not decide this crystal's cell",
measured_x, measured_y, alt.result->search_result.centering,
gemmi::crystal_system_str(alt.result->search_result.system),
100.0 * (ratio - 1.0), moved);
} else {
const double nearest = std::round(ratio);
const bool harmonic = (nearest >= 2.0 && std::abs(ratio - nearest) < 0.15)
|| std::abs(ratio - std::sqrt(3.0)) < 0.15;
logger.Warning("Beam centre check: the two centres DISAGREE about the "
"lattice. The file's ({:.2f},{:.2f}) gives {}-centred {}, "
"{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} (primitive volume "
"{:.0f} A^3, {}/{} frames); the measured ({:.2f},{:.2f}) gives "
"{}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} "
"(primitive volume {:.0f} A^3, {}/{} frames){}",
header_x, header_y, best.result->search_result.centering,
gemmi::crystal_system_str(best.result->search_result.system),
hc.a, hc.b, hc.c, hc.alpha, hc.beta, hc.gamma, best.vol,
best.score, static_cast<int>(validation.size()),
measured_x, measured_y, alt.result->search_result.centering,
gemmi::crystal_system_str(alt.result->search_result.system),
mc.a, mc.b, mc.c, mc.alpha, mc.beta, mc.gamma, alt.vol,
alt.score, static_cast<int>(validation.size()),
harmonic ? fmt::format(" - the volumes differ by {:.2f}x, so one "
"is an axis harmonic of the other and the "
"beam centre is what settles which",
ratio)
: "");
// Deliberately no decision. The frame counts above are the one statistic
// that must not choose: on the measured case the harmonic indexed MORE.
logger.Warning("Beam centre check: the run keeps the file's centre. The "
"frame counts cannot arbitrate this - a cell twice as long "
"must place every spot twice as accurately to score the "
"same, so the harmonic is often the one that indexes more. "
"Re-run with --estimate-beam-center to measure the centre "
"from the spot symmetry, which is exact geometry");
}
} else if (header_indexes) {
logger.Info("Beam centre check: the measured centre ({:.2f},{:.2f}) indexes "
"{}/{} validation frames against the file's {}/{} - the file's "
"centre stands", measured_x, measured_y, alt.score,
static_cast<int>(validation.size()), best.score,
static_cast<int>(validation.size()));
} else {
logger.Info("Beam centre check: neither the file's centre ({}/{} frames) nor "
"the measured one ({}/{}) indexes this run", best.score,
static_cast<int>(validation.size()), alt.score,
static_cast<int>(validation.size()));
}
}
}
}
// Beam centre as an INDEXING HYPOTHESIS. A beam-centre error is not repairable downstream: it
// 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
+7 -2
View File
@@ -111,7 +111,10 @@ struct ProcessConfig {
// Beam centre from the scattered background, measured on EVERY run and reported (--beam-center-check,
// on by default). The fit reads the projection --detect-beam-stop has already built and no frame of
// its own, so it is close to free, and its result is worth having on a run that comes out right:
// it is what lets a job reprocessing a database say a deposited header is 73 px out.
// it is what lets a job reprocessing a database say a deposited header is 73 px out. It also
// indexes a SECOND first pass at that centre, where the move is significant, and compares the
// two answers - the one comparison an on-failure trigger cannot make, because a centre wrong
// along the spindle does not fail.
bool beam_center_check = true;
// Adaptive integration radius (--adaptive-integration-radius). When set, the pre-scan measures how
@@ -338,7 +341,9 @@ class Rugnux {
// cost of the failure.
bool beam_center_searched_ = false;
// What the scattered background makes the beam centre, measured in the pre-scan (--beam-center-check)
// and never committed there: the run keeps the centre it was given.
// and never committed there: the run keeps the centre it was given. RunPipeline consumes it for
// the second first pass, and clears it, so the two-pass loop does not re-ask on a post-refined
// geometry the pre-scan estimate says nothing about.
std::optional<BeamCenterEstimate> background_center_;
// The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the