diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 93f8bb98a..ec9c45588 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -12,6 +12,7 @@ * Systematic-absence calls rest on the evidence rather than on counts: a screw axis whose absent class the data show extinct is no longer refused because a handful of reflections in it read as present, and `SPACE_GROUP_ALTERNATIVES=` no longer drops a candidate that differs only on a zone the sweep never measured. * A reference correlation measured on too few reflections is refused instead of scored zero, so a run given a reference MTZ is no longer reindexed on an operator that mapped almost everything outside the reference's coverage. * A frame counts as indexed from 6 spots on its lattice rather than 9, so a weakly diffracting crystal whose frames cannot carry 9 is no longer refused the lattice it fits; `--min-indexed-spots` overrides it. +* Where the two-pass rotation indexer chooses between the file's geometry and its own refined one, a refined pass that measures far fewer of the low-order axial reflections is rejected unless it is also better overall, so a crystal no longer loses a screw axis to a pass that merged marginally better while missing the reflections a space group is read from. * `-C` accepts a known cell in any equivalent description - conventional or primitive, centred or not - instead of only the reduced primitive form, so a centred cell given the way it is published no longer makes the run report that it found no lattice. * Each reflection is corrected for the sensor's quantum efficiency at the angle it meets the detector (attenuation lengths from the NIST tables, which also fixes the spot-width parallax term on CdTe) and for the attenuation of the flight path between the sample and its pixel; `--flight-path air|helium|vacuum` declares the medium - default air, since no file states it - and the report says what was assumed and what it was worth. The unmerged MTZ records the factors in new `QE` and `FLIGHT` columns beside `LP`, so raw counts are `I / LP * QE * FLIGHT`, and `_process.h5` in new optional `qe` and `flight` datasets. * Rotation geometry post-refinement fits the crystal and the detector at once, against the observed spot positions and the observed rocking angles together, so the refined distance depends far less on how wrong the file's distance was. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 0f5bd6786..e47cddbe8 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1746,13 +1746,44 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { const double cc2 = pass2.search_merge_cc_half; constexpr double MAX_CREDIBLE_COMPLETENESS = 100.5; // rounding headroom, nothing more constexpr double MAX_CC_HALF_LOSS = 0.05; + + // Third arm: the axial rows. CC1/2 is an average over tens of thousands of reflections and + // the principal axial rows are a few dozen of them, so a pass can measure half the row and + // move CC1/2 by nothing - yet those are the reflections the systematic absences are read + // off, and losing them costs a screw axis and with it the space group. It happens because a + // reflection on an axial row crosses the Ewald sphere over a narrow range of the sweep, + // where a general reflection is spread over it: a pass that loses a wedge of the sweep loses + // whole axial rows while its bulk statistics improve, which is exactly what a wedge-losing + // pass looks like from CC1/2 - better. Measured over a dozen rotation crystals the two + // passes agree on this count to within 9% and usually exactly; the crystal this arm was + // written for measured 44%. + // + // It only decides where CC1/2 does not: a pass that wins on CC1/2 by more than the loss + // the arm above already tolerates has said something about the bulk of the data that a few + // dozen reflections cannot answer, and is kept. And it needs pass 1 to have measured enough + // of the rows for the comparison to mean anything - below eight there is no zone a screw + // could be claimed on either side, so the difference decides nothing. + constexpr double MIN_AXIAL_ROWS_KEPT = 0.75; + constexpr int64_t MIN_AXIAL_ROWS_TO_COMPARE = 8; + const int64_t axial1 = pass1.search_merge_axial_reflections; + const int64_t axial2 = pass2.search_merge_axial_reflections; + const bool lost_axial_rows = axial1 >= MIN_AXIAL_ROWS_TO_COMPARE + && axial2 < MIN_AXIAL_ROWS_KEPT * static_cast(axial1) + && cc2 < cc1 + MAX_CC_HALF_LOSS; + if (pass2.lattice_conflicts_with_fixed_sg || (compl_measured && compl2 > MAX_CREDIBLE_COMPLETENESS) - || cc2 < cc1 - MAX_CC_HALF_LOSS) { - logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass " - "({}, CC1/2 before corrections {:.3f} vs {:.3f}) - going back to the " - "header geometry. The refined geometry did not help this crystal.", - compl_text(), cc2, cc1); + || cc2 < cc1 - MAX_CC_HALF_LOSS || lost_axial_rows) { + const std::string worse = (cc2 < cc1 - MAX_CC_HALF_LOSS || !lost_axial_rows) + ? fmt::format("{}, CC1/2 before corrections {:.3f} vs {:.3f}", + compl_text(), cc2, cc1) + : fmt::format("it holds {} of the {} low-order axial reflections the " + "header-geometry pass measured, and its CC1/2 before corrections " + "{:.3f} vs {:.3f} does not say it is better", + axial2, axial1, cc2, cc1); + logger.Warning("Two-pass: the refined pass is worse than the header-geometry pass ({}) " + "- going back to the header geometry. The refined geometry did not help " + "this crystal.", worse); // The refined pass has already written the canonical files, so re-run at the header // geometry to replace them - the same remedy the supercell collapse above uses, and it // only costs a pass on a crystal that was going to be wrong otherwise. @@ -1774,9 +1805,7 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { redo.pass_number = pass2.pass_count + 1; redo.pass_count = pass2.pass_count + 1; redo.pass_decision = fmt::format( - "header geometry re-adopted: the post-refined pass was worse " - "({}, CC1/2 before corrections {:.3f} vs {:.3f})", - compl_text(), cc2, cc1); + "header geometry re-adopted: the post-refined pass was worse ({})", worse); return redo; } if (pass2.pass_decision.empty()) @@ -4026,6 +4055,18 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.search_merge_completeness = result.search_merge_completeness_measured ? 100.0 * o.unique_reflections / o.possible_unique_reflections : 0.0; result.search_merge_cc_half = sm.cc_half_before_corrections; + // ... and the reflections it holds on the principal axial rows, which the guard needs + // separately because CC1/2 cannot express them (see RunAllPasses). Only the low-order + // part is counted: a screw is read off the strong start of its row - the evidence gate + // wants four absences in one zone, and on a 70 A axis that is already at 7 A - while the + // high-order axial reflections are as numerous as they are uninformative and would dilute + // exactly the difference this number exists to show. + constexpr float AXIAL_ROW_LOW_ORDER_D_A = 4.0f; + result.search_merge_axial_reflections = std::count_if( + sm.merged.begin(), sm.merged.end(), [](const MergedReflection &r) { + return ((r.h == 0) + (r.k == 0) + (r.l == 0)) == 2 + && r.d > AXIAL_ROW_LOW_ORDER_D_A; + }); } // Rotation two-pass geometry pre-pass: this first scale/merge has now fitted a frame-order-smoothed diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 154e95102..de46a641d 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -256,6 +256,11 @@ struct ProcessResult { // not to test a number it does not have and the log says so instead of printing a zero. bool search_merge_completeness_measured = false; double search_merge_cc_half = NAN; + // How many reflections that merge holds on the three principal axial rows (h00, 0k0, 00l), at + // low order. These are the reflections whose systematic absences decide the screw axes, they + // number a few dozen against the merge's tens of thousands, and CC1/2 therefore cannot see them + // at all - which is why the two-pass quality guard counts them separately. + int64_t search_merge_axial_reflections = 0; // Per-reflection (I, sigma) of the final merged reflections, for the ISa diagnostic: I/sigma // plotted against I flattens off at the asymptote the error model reports as ISa, so the plot shows