diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 2317b2592..cd54fe056 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -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 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 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); }