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

post-refinement: a one per cent bound on a step, not on the answer

Browse Source 
The fit was refused outright if it wanted to move the detector distance or a cell
edge by one per cent or more. The bound exists to stop a runaway, but it was
applied to the whole move rather than to a step, so a geometry that was genuinely
several per cent out could never be reached - on one dataset the refused fit was
310.000 -> 305.692 mm with a cell within 0.06 per cent of the deposited one, and it
was thrown away for being too large a correction.

A move of one per cent or more is now re-fitted as a walk of one per cent steps,
each seeded where the last arrived and each required to lower the held-out
split-half residual, and a walk of more than one step is ratified by re-indexing
at where it arrived. That is what the bound was really for: a second lattice does
not index better at the new geometry, and a real distance error does - on the
dataset this was found with, pass-1 frames go 37/60 to 60/60, R_meas 0.74 to
0.081, ISa 0.98 to 15.0, and the walk stops at a minimum bracketed on both sides
without ever being told what the answer should be. A walk that uses every step it
is allowed has not settled and is refused; that test is what stops the runaway
this replaces the veto with.

The first solve is untouched. An earlier version that also widened its box was
measurably not inert - two microns on one control moved a 694-reflection merge by
three points of R_meas - so a move inside the bound still commits exactly as before.

Refusals now say so in the report rather than passing silently. Four of ninety-nine
open datasets carry a uniform cell-scale error; ninety-eight of a hundred and three
are byte-identical or differ only by the new line, and no run on forty-eight
datasets with calibrated headers moves at all.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011GxZqDiFP3KqriBhNdcR56
This commit is contained in:
leonarski_fandClaude Opus 5 committed 2026-09-14 00:18:32 +02:00
1 parent 33c0124833
commit 5d643498c6
5 files changed
+245 -55
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/geom_refinement/PostRefine.cpp
+140 -41
View File
@@ -43,6 +43,21 @@ constexpr double BEAM_BOUND_PXL = 15.0;
// +-2.86 deg box the solver works in, so the box reaches everywhere the gate accepts.
constexpr double ANGLE_BOUND_DEG = 1.0;
// How far ONE step of the joint fit may move the detector distance or a cell length. This was a cap
// on the whole fit, and as a cap it was the opposite of its job: the header is most worth correcting
// when it is most wrong, and it threw away fits that were right (measured on a crystal whose refused
// 1.39 % move landed on the deposited cell to 0.06 %). What it is still for - a fit corrupted by
// something every cross-validation fold shares, a second lattice most often, which the relative "the
// held-out residual improved" gate cannot see - is kept by making it a TRUST REGION on the step: a
// solve reaches at most this far, and a fit that wants more takes another step from there. A move of
// more than one step is then ratified by re-indexing at it, which is what actually tells a second
// lattice from a wrong header (see PostRefineResult::large_move).
constexpr double STEP_BOUND = 0.01;
// ... and the walk is bounded, so a fit that never settles still terminates. Five steps is the +-5 %
// box the solver was always given on the distance, i.e. as far as this fit has ever been trusted.
constexpr int MAX_STEPS = 5;
// One integrated partial, flattened across all images. Kept as narrow as the sort and the event
// split allow: on a large cell this array is gigabytes, and the scatter and every level of the
// per-bucket sort move all of it. The goniometer angle is not stored - it is a function of the
@@ -491,14 +506,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
return (np + ne) ? (cp + ce) / static_cast<double>(np + ne) : 0.0;
};
// One problem, both residual families, every block seeded at nominal. The detector tilt is
// declared and held constant so that a non-zero rot1/rot2 still acts on the observed side.
auto solve_joint = [&](Subset s, double bm[2], double ds[1], double rv[3],
// One problem, both residual families, every block seeded at the geometry it is handed -
// nominal for the first step of the walk below, the last step's answer after that. The
// detector tilt is declared and held constant so that a non-zero rot1/rot2 still acts on
// the observed side.
auto solve_joint = [&](Subset s, bool wide, double bm[2], double ds[1], double rv[3],
double q0[3], double q1[3], double q2[3]) {
bm[0] = beam_x0; bm[1] = beam_y0; ds[0] = dist0;
for (int j = 0; j < 3; ++j) {
rv[j] = ax0[j]; q0[j] = p0_0[j]; q1[j] = p1_0[j]; q2[j] = p2_0[j];
}
const double ds_seed = ds[0];
const double q1_seed[3] = {q1[0], q1[1], q1[2]};
ceres::Problem p;
for (size_t oi = 0; oi < n_obs; ++oi) {
const Partial *pp = obs[oi];
@@ -519,20 +534,30 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
if (p.NumResidualBlocks() == 0) return false;
p.SetParameterBlockConstant(det_rot);
if (!p2_free) p.SetParameterBlockConstant(q2);
p.SetParameterLowerBound(ds, 0, dist0 * 0.95); p.SetParameterUpperBound(ds, 0, dist0 * 1.05);
// The box has to reach everywhere the gate below would accept, or the gate is never the
// thing that decides: a fit pinned at a box face lands exactly ON the bound and is then
// refused for being there.
// Two boxes on the distance and the cell lengths. WIDE is the one this fit has always
// had, around nominal: it is what the first solve is asked in, so a fit that settles
// inside the bound is the fit this code has always committed. Otherwise the box is one
// STEP around wherever this solve starts - the trust region the walk below steps in.
// The beam, the axis and the cell angles keep their own absolute bounds around NOMINAL
// either way, since a box that travelled with the seed would let them drift a whole
// bound per step. Each of those has to reach everywhere the gate below would accept, or
// the gate is never the thing that decides - a fit pinned at a box face lands exactly ON
// the bound and is then refused for being there.
p.SetParameterLowerBound(ds, 0, wide ? dist0 * 0.95 : ds_seed * (1.0 - STEP_BOUND));
p.SetParameterUpperBound(ds, 0, wide ? dist0 * 1.05 : ds_seed * (1.0 + STEP_BOUND));
for (int j = 0; j < 2; ++j) {
const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : bm[j];
p.SetParameterLowerBound(bm, j, std::min(bm[j], m) - BEAM_BOUND_PXL);
p.SetParameterUpperBound(bm, j, std::max(bm[j], m) + BEAM_BOUND_PXL);
const double b0 = (j == 0) ? beam_x0 : beam_y0;
const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : b0;
p.SetParameterLowerBound(bm, j, std::min(b0, m) - BEAM_BOUND_PXL);
p.SetParameterUpperBound(bm, j, std::max(b0, m) + BEAM_BOUND_PXL);
}
for (int j = 0; j < 3; ++j) {
p.SetParameterLowerBound(rv, j, ax0[j] - 0.05);
p.SetParameterUpperBound(rv, j, ax0[j] + 0.05);
p.SetParameterLowerBound(q1, j, 0.97 * p1_0[j]);
p.SetParameterUpperBound(q1, j, 1.03 * p1_0[j]);
p.SetParameterLowerBound(q1, j, wide ? 0.97 * p1_0[j]
: (1.0 - STEP_BOUND) * q1_seed[j]);
p.SetParameterUpperBound(q1, j, wide ? 1.03 * p1_0[j]
: (1.0 + STEP_BOUND) * q1_seed[j]);
}
if (p2_free)
for (int j = 0; j < 3; ++j) {
@@ -551,30 +576,78 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
if (n_obs >= static_cast<size_t>(settings.min_events)) {
double pos_nom = 0.0, exc_nom = 0.0, pos_ref = 0.0, exc_ref = 0.0;
const double cv_nom = joint_cost(VAL, beam, dist, axv, p0, p1, p2, pos_nom, exc_nom);
double bm_f[2], ds_f[1], rv_f[3], q0_f[3], q1_f[3], q2_f[3];
const bool convJ = solve_joint(FIT, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f);
const double cv_ref = joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, pos_ref, exc_ref);
// Commit only for a small, credible move: distance < 1 % (a calibrated header needs
// < ~0.6 %), each cell length < 1 %, each free cell angle < 1 deg (ANGLE_BOUND_DEG), and
// the beam inside the bound measured from whichever centre anything already believes is
// nearer. A larger move is the red flag for an unreliable fit - typically a second
// lattice whose spots bias every cross-validation fold identically, so the relative "it
// improved" gate is blind to it. The absolute size of the move discriminates a genuine
// header correction from that failure far better than the absolute residual, which real
// marginal (noisy / iced) data shares with it. Names the test it failed, or nullptr.
const auto out_of_bounds = [&](const double bm[2], const double ds[1],
const double q1c[3], const double q2c[3]) -> const char * {
if (std::fabs(ds[0] - dist0) >= 0.01 * dist0)
return "the distance moved more than 1 %";
double len_shift = 0.0, ang_shift = 0.0;
for (int j = 0; j < 3; ++j) {
// How far a geometry has moved from nominal, in units of one step.
const auto beyond_a_step = [&](const double ds[1], const double q1c[3]) {
double len_shift = 0.0;
for (int j = 0; j < 3; ++j)
len_shift = std::max(len_shift,
std::fabs(q1c[j] - p1_0[j]) / std::max(1e-9, p1_0[j]));
ang_shift = std::max(ang_shift, std::fabs(q2c[j] - p2_0[j]));
return std::fabs(ds[0] - dist0) >= STEP_BOUND * dist0 || len_shift >= STEP_BOUND;
};
// Ask the fit the way it has always been asked: ONE solve, seeded at nominal, in the
// wide box. A move that stays inside the bound is the fit this code has always
// committed and is committed unchanged - which is every crystal on the corpus but a
// handful, so nothing else about their runs moves either.
double bm_f[2] = {beam_x0, beam_y0}, ds_f[1] = {dist0};
double rv_f[3], q0_f[3], q1_f[3], q2_f[3];
for (int j = 0; j < 3; ++j) {
rv_f[j] = ax0[j]; q0_f[j] = p0_0[j]; q1_f[j] = p1_0[j]; q2_f[j] = p2_0[j];
}
bool convJ = solve_joint(FIT, /*wide=*/true, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f);
double cv_ref = convJ ? joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f,
pos_ref, exc_ref)
: cv_nom;
int steps = 1;
// A larger move is no longer refused for its size - but it is not taken as one solve
// asserted it either. It is re-fitted as a WALK: each step seeded where the last one
// arrived, each bounded by the trust region, and each required to lower the held-out
// residual, so the size is arrived at by a sequence of checked steps. The walk stops
// when a step stops paying, which is the same cross-validated quantity the commit gate
// decides on rather than a new rule, and a walk that never stops is refused below.
const bool walked = convJ && beyond_a_step(ds_f, q1_f);
if (walked) {
bm_f[0] = beam_x0; bm_f[1] = beam_y0; ds_f[0] = dist0;
for (int j = 0; j < 3; ++j) {
rv_f[j] = ax0[j]; q0_f[j] = p0_0[j]; q1_f[j] = p1_0[j]; q2_f[j] = p2_0[j];
}
if (len_shift >= 0.01)
return "a cell length moved more than 1 %";
convJ = false; cv_ref = cv_nom; steps = 0;
for (int step = 0; step < MAX_STEPS; ++step) {
double bm_s[2] = {bm_f[0], bm_f[1]}, ds_s[1] = {ds_f[0]};
double rv_s[3], q0_s[3], q1_s[3], q2_s[3];
for (int j = 0; j < 3; ++j) {
rv_s[j] = rv_f[j]; q0_s[j] = q0_f[j]; q1_s[j] = q1_f[j]; q2_s[j] = q2_f[j];
}
if (!solve_joint(FIT, /*wide=*/false, bm_s, ds_s, rv_s, q0_s, q1_s, q2_s))
break;
double pos_s = 0.0, exc_s = 0.0;
const double cv_s = joint_cost(VAL, bm_s, ds_s, rv_s, q0_s, q1_s, q2_s,
pos_s, exc_s);
if (step > 0 && !(cv_s < cv_ref))
break;
bm_f[0] = bm_s[0]; bm_f[1] = bm_s[1]; ds_f[0] = ds_s[0];
for (int j = 0; j < 3; ++j) {
rv_f[j] = rv_s[j]; q0_f[j] = q0_s[j]; q1_f[j] = q1_s[j]; q2_f[j] = q2_s[j];
}
cv_ref = cv_s; pos_ref = pos_s; exc_ref = exc_s;
convJ = true;
++steps;
}
}
// Commit only for a credible move: each free cell angle under 1 deg (ANGLE_BOUND_DEG)
// and the beam inside the bound measured from whichever centre anything already
// believes is nearer. The distance and the cell lengths are no longer bounded here -
// they are bounded one step at a time by the walk above, and a move of more than one
// step is ratified by re-indexing at where it arrived (result.large_move), which is the
// discriminator the size of the move was standing in for. Names the test it failed, or
// nullptr.
const auto out_of_bounds = [&](const double bm[2], const double q2c[3]) -> const char * {
double ang_shift = 0.0;
for (int j = 0; j < 3; ++j)
ang_shift = std::max(ang_shift, std::fabs(q2c[j] - p2_0[j]));
if (ang_shift >= ANGLE_BOUND_DEG * PI / 180.0)
return "a cell angle moved more than 1 deg";
const double from_nominal = std::hypot(bm[0] - beam_x0, bm[1] - beam_y0);
@@ -594,10 +667,17 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
// quantity the cell is committed for has got worse. Neither may degrade.
const char *refused =
!convJ ? "the fit did not converge"
// A walk that used every step it was given never settled: it stopped because it
// ran out of steps, not because it had arrived, and where that happens it is
// walking away rather than to something. This is the runaway the size of the
// move used to stand in for, tested where it can actually be seen - measured on
// a weak 7 A crystal whose fit asked for another 1 % at every one of five steps
// and 4.3 % in all, against the 1.2 % its cell is actually out by.
: steps >= MAX_STEPS ? "the fit used every step it was given and never settled"
: !(cv_ref < 0.98 * cv_nom) ? "the held-out residual did not improve enough"
: !(pos_ref < pos_nom) ? "the held-out positional residual did not improve"
: !(exc_ref < exc_nom) ? "the held-out excitation residual did not improve"
: out_of_bounds(bm_f, ds_f, q1_f, q2_f);
: out_of_bounds(bm_f, q2_f);
// What the half that was held out earned the right to fit is re-fitted on all of it, and
// that second solve moves the geometry - so the bounds are asked again of what will
@@ -607,10 +687,19 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
double q1_log[3], q2_log[3];
for (int j = 0; j < 3; ++j) { q1_log[j] = q1_f[j]; q2_log[j] = q2_f[j]; }
if (refused == nullptr) {
solve_joint(ALL, beam, dist, axv, p0, p1, p2); // commit: re-fit on all data
// A walk is re-fitted from where it arrived, one step wide, so the re-fit is its
// last step taken over all the data rather than the whole walk taken again on it.
// Everything else is re-fitted from nominal in the wide box, exactly as before.
if (walked) {
beam[0] = bm_f[0]; beam[1] = bm_f[1]; dist[0] = ds_f[0];
for (int j = 0; j < 3; ++j) {
axv[j] = rv_f[j]; p0[j] = q0_f[j]; p1[j] = q1_f[j]; p2[j] = q2_f[j];
}
}
solve_joint(ALL, !walked, beam, dist, axv, p0, p1, p2); // commit: re-fit on all data
ds_log[0] = dist[0]; bm_log[0] = beam[0]; bm_log[1] = beam[1];
for (int j = 0; j < 3; ++j) { q1_log[j] = p1[j]; q2_log[j] = p2[j]; }
refused = out_of_bounds(beam, dist, p1, p2);
refused = out_of_bounds(beam, p2);
if (refused != nullptr) {
beam[0] = beam_x0; beam[1] = beam_y0; dist[0] = dist0;
for (int j = 0; j < 3; ++j) {
@@ -623,8 +712,10 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
// comes out is the same in all of them, so a run that silently keeps its header geometry
// says nothing about whether the fit was bad, the improvement too small, or the move too
// large for the bound - which is the one case where the number worth reading is the one
// that was thrown away.
const std::string verdict = commit ? std::string("COMMIT")
// that was thrown away. The report carries the same sentence, so a refusal is visible
// without the log.
const std::string verdict = commit ? fmt::format("COMMIT ({} step{})", steps,
steps == 1 ? "" : "s")
: "reject (" + std::string(refused) + ")";
// The cell the log names is the EFFECTIVE one - what the residual's B matrix builds
// from the blocks - not the blocks themselves, whose unused components a high-symmetry
@@ -645,6 +736,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
len_fit[0], len_fit[1], len_fit[2],
ang_fit[0] * 180.0 / PI, ang_fit[1] * 180.0 / PI, ang_fit[2] * 180.0 / PI,
pos_nom, pos_ref, exc_nom, exc_ref, verdict);
if (!commit)
result.refused_reason = fmt::format(
"{} - it wanted {:.3f} mm and cell {:.3f} {:.3f} {:.3f}", refused,
ds_log[0], len_fit[0], len_fit[1], len_fit[2]);
// A move of more than one step is a walk, not a correction next to where it started,
// and it is the one the old bound refused outright. It is committed here and ratified
// where it can be: by re-indexing at it (see RunAllPasses).
result.large_move = commit && beyond_a_step(dist, p1);
} else {
logger.Info("Post-refine GEOM: only {} positional observations - the joint fit needs the "
"spot positions as well as the rocking angles, so nothing is refined", n_obs);
image_analysis/geom_refinement/PostRefine.h
+10
View File
@@ -5,6 +5,7 @@
#include <array>
#include <optional>
#include <string>
#include <vector>
#include "../../common/DiffractionGeometry.h"
@@ -40,6 +41,15 @@ struct PostRefineResult {
// together - so these two are always equal. Both are kept because the report names them separately.
bool cell_refined = false; // the joint fit passed cross-validation (crystal half)
bool detector_refined = false; // the joint fit passed cross-validation (detector half)
// The fit moved the distance or a cell length by more than one step of its trust region (see
// STEP_BOUND in PostRefine.cpp), so it walked rather than settling beside the geometry it started
// from. Such a move is not refused for its size - it is RATIFIED by re-integrating and re-indexing
// at it, which a second lattice does not survive and a wrong header does (see RunAllPasses).
bool large_move = false;
// Which test refused the fit and what it wanted, where it was refused; empty where it committed.
// Nothing else says it - the geometry that comes out of a refusal is the nominal one whatever
// refused it, so a silent refusal cannot be told from a fit that never wanted to move.
std::string refused_reason;
// GONIOMETER ROTATION SCALE: the factor by which the stage actually turned relative to the angle
// stored in the file (which is the COMMANDED value, hence a stage calibration error is invisible in
// the header). Fitted after the joint fit as a single free parameter, with the crystal and the axis
rugnux/ResultReport.cpp
+6
View File
@@ -457,6 +457,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
Add(s, KeyText("POSTREFINE_BEAM_CENTRE",
fmt::format("{:.2f} {:.2f} -> {:.2f} {:.2f}", pr.beam_x_before_px,
pr.beam_y_before_px, pr.beam_x_after_px, pr.beam_y_after_px)));
// Which test refused the fit, where one did. Every refusal leaves the same geometry - the
// header's - so without this the report cannot tell a fit that was thrown away from one that
// never wanted to move, and a refusal is the one case where the number worth reading is the
// one that was thrown away.
if (!pr.refused_reason.empty())
Add(s, KeyText("POSTREFINE_REFUSED", pr.refused_reason));
Add(s, KeyReal("GONIOMETER_ROTATION_SCALE", pr.rotation_scale, "{:.5f}"));
Add(s, KeyBool("GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect));
if (pr.rotation_scale_suspect)
rugnux/Rugnux.cpp
+84 -14
View File
@@ -1969,12 +1969,68 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) {
logger.Info("Rotation two-pass geometry post-refinement: second pass (refined geometry, canonical) -> {}_*",
base_prefix);
config_.output_prefix = base_prefix; // the refined pass is the canonical result (no _02 suffix)
// The canonical pass measures the post-refinement again where the pre-pass moved the geometry
// (it applies nothing), so that the walk below has something to decide on.
postrefine_probe_ = prepass_detector_geometry_.has_value();
auto pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
// Walk the geometry to the fit's fixed point. The post-refinement is fitted to reflections
// integrated at the very geometry it is correcting, so one fit reaches only as far as those
// reflections see: on a header whose distance is several percent long, the pre-pass absorbs the
// error into an inflated cell, and from that self-consistent pair the fit finds a fraction of it
// and stops. The canonical pass has just re-integrated and re-indexed at what it committed, so
// asking the same fit there is the next step of the same walk - and that re-indexing is also the
// RATIFICATION the old size bound stood in for: a second lattice does not index better at a
// geometry it did not come from, a wrong header does, and a pass that lost the crystal is thrown
// out by the quality guard below, which returns the run to the header geometry.
//
// What starts the walk is a move of more than one step: that is the fit saying the geometry is
// somewhere else entirely rather than a fraction of a percent away, and it leaves every run
// whose fit settles beside its header at the two passes it always had. Once the run IS walking
// it keeps walking while the fit still commits, because the steps do not get smaller in
// proportion to what is left: measured on a header 3.9 % long, the first round took a third of
// the error and each of the next took a third of the rest. The walk is bounded rather than
// stopped at a size, because there is no size that separates the last step of a real walk from
// the first step of a runaway - and it does not run away: started BELOW the optimum the same fit
// turns back up, so it settles around the answer instead of past it.
constexpr int MAX_GEOMETRY_ROUNDS = 4;
std::optional<std::array<float, 5>> adopted_geometry;
int geometry_rounds = 0;
while (!cancelled_ && geometry_rounds < MAX_GEOMETRY_ROUNDS && prepass_detector_geometry_
&& pass2.post_refine && pass2.post_refine->detector_refined
&& (pass2.post_refine->large_move || geometry_rounds > 0)) {
adopted_geometry = *prepass_detector_geometry_;
const auto &g = *adopted_geometry;
logger.Info("Two-pass: the post-refinement at the adopted geometry still moves it "
"(distance {:.3f} -> {:.3f} mm) - re-integrating and re-indexing at what it "
"asks for (round {} of at most {})", pass2.post_refine->distance_before_mm,
pass2.post_refine->distance_after_mm, geometry_rounds + 1,
MAX_GEOMETRY_ROUNDS);
experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2])
.PoniRot1_rad(g[3]).PoniRot2_rad(g[4]);
pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
++geometry_rounds;
}
postrefine_probe_ = false;
// The post-refinement is measured by pass 1 and consumed by pass 2, so carry it onto whichever
// result is returned - it is a result of the run, not of the pass that happened to fit it.
// result is returned - it is a result of the run, not of the pass that happened to fit it. Where
// the walk took rounds they are ONE geometry decision, so what the report names is where the
// geometry started and where the run left it, not the first step of the way there.
if (geometry_rounds > 0 && pass1.post_refine && adopted_geometry) {
pass1.post_refine->beam_x_after_px = (*adopted_geometry)[0];
pass1.post_refine->beam_y_after_px = (*adopted_geometry)[1];
pass1.post_refine->distance_after_mm = (*adopted_geometry)[2];
}
pass2.post_refine = pass1.post_refine;
pass2.pass_number = 2;
pass2.pass_count = 2;
pass2.pass_number = 2 + geometry_rounds;
pass2.pass_count = 2 + geometry_rounds;
if (geometry_rounds > 0)
pass2.pass_decision = fmt::format(
"post-refined geometry adopted over {} rounds: the fit was still walking at each of "
"them, and the detector distance ended at {:.3f} mm against the header's {:.3f}",
geometry_rounds + 1, (*adopted_geometry)[2],
pass1.post_refine ? pass1.post_refine->distance_before_mm : 0.0);
// r2 = r1 + 2 is also the inner edge of the background ring, so widening the signal disk pushes
// that ring out into the neighbours; a reflection whose ring is left with five or fewer clean
@@ -4901,14 +4957,23 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// apply the detector geometry from this integration. That is all the pre-pass is for: everything
// the second pass consumes from it is fitted here, and the space group is left to the second pass
// (see the search gate below).
if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) {
// postrefine_probe_ brings the post-refinement into a CANONICAL pass as well, where the
// pre-pass moved the geometry: the fit it committed was made on reflections integrated at the
// geometry BEFORE the move, and whether it has finished walking can only be asked of
// reflections integrated at the geometry after it. The probe measures - it applies nothing and
// leaves the smoothed mosaicity and the pre-pass lattice alone, which stay pass 1's for every
// pass that consumes them; RunAllPasses decides what to do with the answer.
if ((geometry_prepass || postrefine_probe_) && rsm && experiment_.IsRotationIndexing()) {
const auto &outcomes = indexer->GetIntegrationOutcome();
prepass_mosaicity_.assign(outcomes.size(), NAN);
for (size_t o = 0; o < outcomes.size(); ++o)
if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg;
if (geometry_prepass) {
prepass_mosaicity_.assign(outcomes.size(), NAN);
for (size_t o = 0; o < outcomes.size(); ++o)
if (outcomes[o].mosaicity_deg) prepass_mosaicity_[o] = *outcomes[o].mosaicity_deg;
}
if (result.consensus_cell.has_value()) {
if (const auto rot = indexer->FinalizeRotationIndexing(); rot && rot->axis) {
prepass_result_ = *rot; // full result (lattice+geom+axis), for the supercell fallback
if (geometry_prepass)
prepass_result_ = *rot; // full result (lattice+geom+axis), for the supercell fallback
PostRefineSettings prs;
prs.crystal_system = rot->search_result.system;
prs.num_threads = static_cast<int>(config_.nthreads);
@@ -4940,14 +5005,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
rot->geom.GetPoniRot1_rad(),
rot->geom.GetPoniRot2_rad()};
logger.Info("Two-pass: post-refined detector geometry (distance {:.3f} mm, "
"beam {:.2f},{:.2f} px, tilt {:.4f},{:.4f} deg) will drive the "
"second integration pass",
"beam {:.2f},{:.2f} px, tilt {:.4f},{:.4f} deg) - {}",
pr.distance_after_mm, pr.beam_x_after_px, pr.beam_y_after_px,
rot->geom.GetPoniRot1_rad() * 180.0 / PI,
rot->geom.GetPoniRot2_rad() * 180.0 / PI);
rot->geom.GetPoniRot2_rad() * 180.0 / PI,
geometry_prepass ? "it will drive the second integration pass"
: "what the fit still asks for at the geometry "
"this pass ran at");
} else {
logger.Info("Two-pass: geometry post-refine committed no detector change "
"- the second pass reproduces the first");
logger.Info("Two-pass: geometry post-refine committed no detector change{}",
geometry_prepass ? " - the second pass reproduces the first" : "");
}
// DECISION POINT for the goniometer rotation scale. It does not ride on pr.ok - the
// scale is its own cross-validated fit and the crystals that have a stage fault are
@@ -4957,7 +5024,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// would silently damage it, so the asymmetry is deliberate: committed only when the
// fit is both cross-validated and outside the tolerance. A manual --rotation-scale is
// already on the goniometer and is left alone.
if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value())
// Only the pre-pass fits the rotation scale. It is applied from the goniometer
// the file came with, so a scale measured again on angles that have already been
// corrected once is not a correction that can be applied on top of that one.
if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value() && geometry_prepass)
prepass_rotation_scale_ = static_cast<float>(pr.rotation_scale);
}
}
rugnux/Rugnux.h
+5
View File
@@ -501,6 +501,11 @@ class Rugnux {
// of itself as soon as its search merge is in, which is where it learns that its report and its
// files are about to be superseded.
std::optional<ProcessResult> quality_guard_pass1_;
// Whether a canonical pass also MEASURES the post-refinement (it never applies it). The fit the
// pre-pass committed was made on reflections integrated at the geometry before the move, so
// whether it has finished walking is a question only reflections integrated after the move can
// answer; RunAllPasses turns this on for the passes that answer it.
bool postrefine_probe_ = false;
// Pre-scan: read a spread sample of frames and take three things off them - the shadow of the beam
// stop and its holder, added to the pixel mask (config_.detect_beam_stop), the beam centre