diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index a6da5bbb2..8e40c0399 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -5006,6 +5006,99 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const auto &rot_ss = experiment_.GetScalingSettings(); const bool is_rotation = experiment_.IsRotationIndexing(); // rotation indexing -> rotation scaling/merge std::optional rsm; + std::optional prepass_postrefine_obs; + // The rotation geometry post-refinement (see its call sites below for what it is for). + const auto post_refine_geometry = [&] { + const auto &outcomes = indexer->GetIntegrationOutcome(); + 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) { + 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(config_.nthreads); + prs.refine_geometry = true; + // The centre the pre-scan measured, where the fit was precise enough to be used as + // a centre at all (the same ceiling --estimate-beam-center adopts one on). It does + // not move the geometry - it only tells the bound on step B's beam move where else + // the centre is independently believed to be. + if (measured_beam_center_ + && measured_beam_center_->sigma_pxl <= std::max(MAX_BEAM_CENTER_SIGMA_PXL, + BeamCenterNeed_pxl(experiment_))) + prs.measured_beam_px = std::array{ + static_cast(measured_beam_center_->beam_x_pxl), + static_cast(measured_beam_center_->beam_y_pxl)}; + // The first pass's own geometry, not the header's. The outcomes being post-refined + // were integrated at it, and step B holds the tilt fixed while it fits the beam and + // the distance - so handing it the header tilt fits a beam that has to absorb a tilt + // error the images do not have, and the answer describes neither geometry. + // + // On the pre-pass the observations were gathered right after Ingest (see there) and + // the scaling engine is finished - nothing after this point reads it on that pass + // (the search and the reference re-merge run only on the canonical one) - so it is + // freed BEFORE the fit: the engine holds about twice the partial payload again, + // and with it gone the fit no longer sits on the pass's memory peak. + if (prepass_postrefine_obs) + rsm.reset(); + auto pr_obs = prepass_postrefine_obs + ? std::move(*prepass_postrefine_obs) + : GatherPostRefineObservations(indexer->GetIntegrationOutcome(), + config_.nthreads, + /*release_reflections=*/false, logger); + prepass_postrefine_obs.reset(); + const auto pr = PostRefineRotationGeometry(std::move(pr_obs), *rot->axis, + rot->geom, rot->lattice, prs, logger); + result.post_refine = pr; + if (pr.ok && pr.detector_refined) { + // Store, don't apply: this pass continues to write its own (header-geometry) output; + // Run() applies the refined geometry before the second pass. + prepass_detector_geometry_ = std::array{ + static_cast(pr.beam_x_after_px), + static_cast(pr.beam_y_after_px), + static_cast(pr.distance_after_mm), + 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) - {}", + 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, + 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{}", + 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 + // exactly the ones whose cell and detector steps do NOT pass, because the angle error + // is what their residual is made of. A calibration fault is rare (36 of 37 rotation + // datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset + // 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. + // 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(pr.rotation_scale); + } + } + }; + // On a CANONICAL pass the post-refinement only measures (see below), and it reads nothing the + // scaling engine produces: the reflections' hkl, I, sigma and positions are final once the + // images are integrated, and the engine writes back per-frame fields only. So measure it + // before the engine is built rather than beside it - the gathered observations and their + // rocking events are gigabytes on a long axis, and beside the engine's own arrays they set + // the pass's memory high-water mark. + if (postrefine_probe_ && !geometry_prepass && is_rotation) + post_refine_geometry(); if (is_rotation) { if (experiment_.GetRefineRotationWedgeInScaling() || rot_ss.GetRotationWedgeForScaling().has_value()) @@ -5029,7 +5122,6 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // it is read: on a fine-sliced long axis the vectors are gigabytes, held to the end of the // pass for those 32 bytes a partial. A run with a reference MTZ keeps them: resolving the // indexing ambiguity re-reads them (and re-ingests). - std::optional prepass_postrefine_obs; if (geometry_prepass && rsm && config_.reference_data.empty()) prepass_postrefine_obs = GatherPostRefineObservations( indexer->GetIntegrationOutcome(), config_.nthreads, @@ -5211,89 +5303,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // 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(); - 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) { - 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(config_.nthreads); - prs.refine_geometry = true; - // The centre the pre-scan measured, where the fit was precise enough to be used as - // a centre at all (the same ceiling --estimate-beam-center adopts one on). It does - // not move the geometry - it only tells the bound on step B's beam move where else - // the centre is independently believed to be. - if (measured_beam_center_ - && measured_beam_center_->sigma_pxl <= std::max(MAX_BEAM_CENTER_SIGMA_PXL, - BeamCenterNeed_pxl(experiment_))) - prs.measured_beam_px = std::array{ - static_cast(measured_beam_center_->beam_x_pxl), - static_cast(measured_beam_center_->beam_y_pxl)}; - // The first pass's own geometry, not the header's. The outcomes being post-refined - // were integrated at it, and step B holds the tilt fixed while it fits the beam and - // the distance - so handing it the header tilt fits a beam that has to absorb a tilt - // error the images do not have, and the answer describes neither geometry. - // - // On the pre-pass the observations were gathered right after Ingest (see there) and - // the scaling engine is finished - nothing after this point reads it on that pass - // (the search and the reference re-merge run only on the canonical one) - so it is - // freed BEFORE the fit: the engine holds about twice the partial payload again, - // and with it gone the fit no longer sits on the pass's memory peak. - if (prepass_postrefine_obs) - rsm.reset(); - auto pr_obs = prepass_postrefine_obs - ? std::move(*prepass_postrefine_obs) - : GatherPostRefineObservations(indexer->GetIntegrationOutcome(), - config_.nthreads, - /*release_reflections=*/false, logger); - prepass_postrefine_obs.reset(); - const auto pr = PostRefineRotationGeometry(std::move(pr_obs), *rot->axis, - rot->geom, rot->lattice, prs, logger); - result.post_refine = pr; - if (pr.ok && pr.detector_refined) { - // Store, don't apply: this pass continues to write its own (header-geometry) output; - // Run() applies the refined geometry before the second pass. - prepass_detector_geometry_ = std::array{ - static_cast(pr.beam_x_after_px), - static_cast(pr.beam_y_after_px), - static_cast(pr.distance_after_mm), - 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) - {}", - 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, - 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{}", - 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 - // exactly the ones whose cell and detector steps do NOT pass, because the angle error - // is what their residual is made of. A calibration fault is rare (36 of 37 rotation - // datasets sit at 1.0000) and applying a 1 % angle correction to a healthy dataset - // 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. - // 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(pr.rotation_scale); - } - } - } + if (geometry_prepass && rsm && experiment_.IsRotationIndexing()) + post_refine_geometry(); // For a de-novo search, gate the P1 scaling+merge (the pass that feeds the space-group search) at // >= 1 over thin resolution shells: noise-dominated high-res shells otherwise mislead