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

rugnux: decide the beam-centre axis harmonic in the one direction nothing else repairs

Browse Source 
The beam-centre check indexes a second first pass at the measured centre and,
where the two centres return cells related by an integer volume factor, says so
and decides nothing. That abstention is right in general - the validation frame
count leans to the sub-cell, since a cell twice as long must place every spot
twice as accurately to score the same - but it leaves one direction
unrecoverable. Where the MEASURED centre holds the LARGER cell, the run starts
on the sub-cell and stays there: the second pass re-finds the larger cell at the
refined geometry and its supercell guard refuses it for being larger than
pass 1's. The run then merges a halved axis at full indexing rate, with good
merging statistics and nothing to say it is wrong.

Decide that one direction, on the pooled spots rather than on the frames: the
excess each lattice puts on itself over its own wrong-spindle null, which is the
measurement the failed-pass branch of the same check already adopts a centre on,
and the one statistic that does not credit a denser lattice for accidental hits.
The opposite direction is deliberately left alone - a second pass is free to drop
to the smaller cell and does, on every corpus crystal that reaches this branch
that way.

Measured on a tetragonal crystal whose header centre is 3.6 px out across the
spindle, more than the geometry absorbs: the file's centre gives a halved a axis
(43/60 frames, 25.4% of the pooled spots on the lattice) and the measured one the
deposited cell (60/60, 43.3%). Before, the run merged the halved cell and the
deposited model refined to R-free 0.575; after, it merges the deposited lattice
to 1.62 A with R-free 0.254 against the deposition's 0.203.

Five other corpus crystals reach the same branch with the measured centre on the
SMALLER cell; all five are unchanged (space group, cell and R-free identical), as
are two crystals whose halved or doubled axis comes from elsewhere.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
This commit is contained in:
leonarski_f
2026-09-20 18:45:18 +02:00
co-authored by Claude Opus 5
parent 6dd9824eec
commit df4be47c84
3 changed files with 75 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/RUGNUX_ADVANCED.md
+1 -1
View File
@@ -235,7 +235,7 @@ Geometry:
| Option | Description |
| --- | --- |
| `--beam-center-check[=off]` | Measure the beam centre from the isotropy of the scattered background on **every** run, report how far the file's value is from it — against how right this particular geometry needs it to be — and index a **second first pass** at the measured centre to see whether the two centres give the same lattice. The fit reads the projection `--detect-beam-stop` already builds, so it costs no extra frames. **On by default**; `=off` disables. Nothing is committed on a run that indexes: the measured centre is adopted only where the file's centre indexes nothing and the measured one indexes a majority, and a disagreement between two passes that both worked is reported with both cells and left undecided |
| `--beam-center-check[=off]` | Measure the beam centre from the isotropy of the scattered background on **every** run, report how far the file's value is from it — against how right this particular geometry needs it to be — and index a **second first pass** at the measured centre to see whether the two centres give the same lattice. The fit reads the projection `--detect-beam-stop` already builds, so it costs no extra frames. **On by default**; `=off` disables. On a run that indexes, the measured centre is adopted in two cases only: where the file's centre indexes nothing and the measured one indexes a majority, and where the two cells are related by an integer volume factor, the measured centre's is the **larger** one, and it carries materially more of the pooled validation spots than the file's does (each cell scored against its own wrong-spindle null, so a denser lattice is not credited for accidental hits). That one direction is the axis harmonic a centre error along the spindle produces, and nothing downstream repairs it; every other disagreement is reported with both cells and left undecided |
| `--beam-center-search[=N\|off]` | After a first pass that indexes fewer than half the validation frames, step the centre a pixel at a time out to N px along **each** detector axis and keep the first rung that indexes a majority. **On by default** (12 px); `=off` disables. It runs only after a pass that has already failed, and spot finding is not repeated, so a run that indexes never pays for it. Both detector directions are searched: a centre error *across* the spindle collapses the indexed fraction and announces itself, while one *along* it holds the frame count up and quietly returns an axis harmonic |
| `--estimate-beam-center` | Measure the direct beam before indexing, from the symmetry of the spots where the sweep reaches at least half a turn and from the radial background profile where it does not; the value in the file is kept where neither can measure it. Off by default |
| `--no-fit-spindle` | With the above, keep the rotation axis given in the file instead of fitting its skew about the beam |
rugnux/Rugnux.cpp
+73 -19
View File
@@ -3144,6 +3144,20 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
return static_cast<double>(e.on_lattice - e.by_chance) > SPOT_BUDGET_SIGNIFICANCE_Z * sigma;
};
// The same quantity for COMPARING two lattices: the share of the pooled spots each one puts on
// itself over and above what a wrong spindle angle puts there. Each is scored against its own
// null, so a denser lattice is not credited for the spots it catches by accident - which is
// what makes the comparison fair between a cell and its axis harmonic.
auto excess = [](const PooledEvidence &e) {
return static_cast<double>(e.on_lattice - e.by_chance)
/ static_cast<double>(std::max<int64_t>(1, e.spots));
};
auto rate_var = [](const PooledEvidence &e) {
const double n = static_cast<double>(std::max<int64_t>(1, e.spots));
const double p = static_cast<double>(e.on_lattice) / n;
return p * (1.0 - p) / n;
};
// How deep into an image's intensity-ordered spot list this lattice is still being seen - the
// measured spot budget. Same frames, same per-image path as count_indexed above, but scoring the
// spots rather than the frames: each frame's spots are tallied by rank and the run keeps the
@@ -3738,9 +3752,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// * the file's centre indexes and 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;
// * they index and 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.
// * they index and disagree -> say so, with both cells and their primitive-volume ratio.
// Which of an integer-related pair is real is not decided on the frame counts, for the
// reason above; it is decided, in ONE direction only, on the pooled spots - see the
// branch itself. Naming the disagreement is useful even where nothing is decided.
//
// Where the file's centre indexes nothing and the measured one does not either, the axis-sign
// question is asked AGAIN at the measured centre before giving up. The two unknowns are
@@ -3888,15 +3903,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
restore_beam_center(header_x, header_y);
if (best.result.has_value())
at_header = pooled_evidence(*indexer, *best.result);
auto excess = [](const PooledEvidence &e) {
return static_cast<double>(e.on_lattice - e.by_chance)
/ static_cast<double>(std::max<int64_t>(1, e.spots));
};
auto rate_var = [](const PooledEvidence &e) {
const double n = static_cast<double>(std::max<int64_t>(1, e.spots));
const double p = static_cast<double>(e.on_lattice) / n;
return p * (1.0 - p) / n;
};
const double margin = excess(at_measured) - excess(at_header);
const double margin_sigma = std::sqrt(rate_var(at_measured) + rate_var(at_header));
if (alt.result.has_value() && beats_chance(at_measured)
@@ -4020,14 +4026,62 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
"beam centre is what settles which",
ratio)
: "");
// Deliberately no decision. The frame counts above are the one statistic
// Deliberately no decision on the frame counts. They 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");
//
// One direction is decided, and only one: where the MEASURED centre holds
// the LARGER of an integer-related pair. That is the case nothing else in
// the run repairs - the frame count leans to the sub-cell, and the second
// pass refuses a lattice larger than pass 1's (the supercell guard below),
// so a run that starts on the sub-cell stays on it. The opposite direction
// needs no help here: the second pass is free to drop to the smaller cell
// and does, on every corpus crystal that reaches this branch that way.
//
// Decided on the pooled excess over chance, the same measurement the
// failed-pass branch above adopts a centre on. Each lattice is scored
// against its own wrong-spindle null, so the larger cell is not credited
// for the spots a denser lattice catches by accident.
const double nearest_n = std::round(ratio);
const bool integer_harmonic = nearest_n >= 2.0 && nearest_n <= 4.0
&& std::abs(ratio - nearest_n) < 0.15;
bool adopted_measured = false;
if (integer_harmonic && alt.vol > best.vol) {
const PooledEvidence at_header = pooled_evidence(*indexer, *best.result);
try_beam_center(measured_x, measured_y);
const PooledEvidence at_measured = pooled_evidence(*indexer, *alt.result);
const double margin = excess(at_measured) - excess(at_header);
const double margin_sigma = std::sqrt(rate_var(at_measured)
+ rate_var(at_header));
adopted_measured = beats_chance(at_measured)
&& margin > SPOT_BUDGET_SIGNIFICANCE_Z * margin_sigma;
if (adopted_measured) {
logger.Warning("Beam centre check: the larger cell is the measured "
"centre's, and it is the one the spots are on - "
"{:.1f}% of the pooled validation spots lie on it "
"against {:.1f}% on the file's cell (chance {:.1f}% "
"and {:.1f}%) - so the {:.2f}x is a real axis and not "
"a harmonic, and the run adopts the measured centre",
100.0 * at_measured.on_lattice
/ std::max<int64_t>(1, at_measured.spots),
100.0 * at_header.on_lattice
/ std::max<int64_t>(1, at_header.spots),
100.0 * at_measured.by_chance
/ std::max<int64_t>(1, at_measured.spots),
100.0 * at_header.by_chance
/ std::max<int64_t>(1, at_header.spots),
ratio);
best = alt;
} else {
restore_beam_center(header_x, header_y);
}
}
if (!adopted_measured)
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 {
logger.Info("Beam centre check: the measured centre ({:.2f},{:.2f}) indexes "
rugnux/rugnux_cli.cpp
+1 -1
View File
@@ -112,7 +112,7 @@ void print_usage() {
std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl;
std::cout << " --estimate-beam-center Place the beam centre before anything is indexed, and use it in place of the header value when it is measured precisely enough. On a sweep that reaches half a turn it comes from the symmetry of the spot positions - the frames 180 deg apart are each other's mirror image, and every reflection is recorded twice - and where the sweep does not reach that far, from the isotropy of the scattered background, which needs only a few frames. It reads frames of its own, chosen as pairs half a turn apart, so it does not change the mask --detect-beam-stop finds. Ignored when --beam-x/--beam-y are given, and when a stills geometry refinement has already placed the centre from indexed spots" << std::endl;
std::cout << " --beam-center-search[=N|off] After a first pass that indexes fewer than half the validation frames, try the beam centre a pixel at a time out to N px along each detector axis and keep the first one that indexes a majority. A beam-centre error is fixed in the lab frame, so it smears the accumulated reciprocal-space cloud and the FFT takes an axis harmonic instead of the true axis; nothing downstream repairs that. ON by default (12 px); =off disables. It runs only after a pass that has already failed, so a run that indexes never pays for it" << std::endl;
std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. Nothing is committed on a run that indexes: the measured centre is adopted only where the file's indexes nothing and the measured one indexes a majority. ON by default; =off disables" << std::endl;
std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. On a run that indexes, the measured centre is adopted in two cases only: where the file's centre indexes nothing and the measured one indexes a majority, and where the two centres return cells related by an integer volume factor and the measured centre's is the LARGER one AND carries materially more of the pooled spots - the axis harmonic a centre error along the spindle produces, which nothing downstream repairs. ON by default; =off disables" << std::endl;
std::cout << " --no-fit-spindle Take the goniometer axis from the file rather than measuring the spindle's rotation about the beam from the spots. Measuring it is the DEFAULT and the file is never right: every master writes an exact lab axis and no goniometer is one. Both mirror lines of the beam-centre estimator turn with the spindle, so an axis a few tenths of a milliradian out smears the vote until a neighbouring tooth wins. Only has an effect with --estimate-beam-center" << std::endl;
std::cout << std::endl;