diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 75e6275de..f8f06a4cf 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -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(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(validation.size()), measured_x, measured_y, + alt.score, static_cast(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(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(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(validation.size()), best.score, + static_cast(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(validation.size()), alt.score, + static_cast(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 diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 7d1b25771..d45dcabf2 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -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 background_center_; // The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the