diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 379090b51..74ea92a32 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* rugnux: A goniometer that turned further than the angles stored in the file — which are the commanded ones — is now measured and **corrected in the second rotation pass**, or set by hand with `--rotation-scale `; the fitted correction is applied only when it exceeds 0.5 %, moves each end of the sweep by at least 0.5°, and comes out the same with any fifth of the sweep left out. * rugnux: Every `mx` and `scale` run now writes `_report.txt`, a results report modelled on XDS's `CORRECT.LP` — `KEY= value` lines, fixed-width tables and `WARNING:` sentences covering indexing, geometry post-refinement, the space-group decision, merging, twinning, radiation damage and the stretches of the sweep over which the crystal delivered much less than the rest of the run. * Space-group search: a screw axis is now claimed from how decisively its predicted-absent reflections are weaker than the rest of their own axial row, instead of from a minimum count of them, so a screw survives a sweep that recorded few axial reflections and is refused on a row too weak to decide either way. * Bragg integration: the profile fit's `background_variance` now takes the fitted intensity itself out of the fit variance instead of `max(0, I)`, so a reflection that fluctuated below zero no longer reports a background variance two to three times too small and is no longer weighted up for it. diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 4bc39ef81..fae7cf9f6 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -4,6 +4,7 @@ #include "PostRefine.h" #include +#include #include #include "../../common/JFJochMath.h" // PI @@ -54,6 +55,37 @@ struct ScaleAxisExcitationResidual { const double inv_lambda, angle_rad, weight, ex, ey, ez; }; +// GONIOMETER ROTATION SCALE k: the same Ewald excitation residual, but with the cell scale and the axis +// DIRECTION already committed by step A, so the single free quantity is how far the stage actually turned +// per unit of commanded angle. Two differences from step A matter: +// * the angle is measured from the CENTRE of the sweep, not from the goniometer's zero. The reference +// orientation is the one rotation indexing fitted against the commanded angles, so it has already +// absorbed the MEAN angle error; only the part that varies across the sweep is left to fit. Scaling the +// absolute angle instead - which is what reading k off the length of step A's axis vector does - asks +// the fit to also produce a constant offset it has no parameter for, and the least-squares compromise +// shrinks k towards 1 by var(phi) / (var(phi) + phi_centre^2): exactly a factor of four for the common +// case of a sweep starting at zero. +// * e_mid is the reference reciprocal vector already turned to the sweep centre and divided by the +// committed cell scale, so nothing but k is free. +struct RotationScaleResidual { + RotationScaleResidual(double lambda, double dangle_rad, const double u[3], const double e_mid[3]) + : inv_lambda(1.0 / lambda), dangle_rad(dangle_rad), + ux(u[0]), uy(u[1]), uz(u[2]), ex(e_mid[0]), ey(e_mid[1]), ez(e_mid[2]) {} + template + bool operator()(const T *const k, T *residual) const { + const T a = T(-dangle_rad) * k[0]; + const T aa[3] = {a * T(ux), a * T(uy), a * T(uz)}; + const T p_ref[3] = {T(ex), T(ey), T(ez)}; + T p_lab[3]; + ceres::AngleAxisRotatePoint(aa, p_ref, p_lab); + const T zeta = p_lab[0] * p_lab[0] + p_lab[1] * p_lab[1] + p_lab[2] * p_lab[2] + + T(2.0) * p_lab[2] * T(inv_lambda); + residual[0] = zeta * T(0.5) / T(inv_lambda); + return true; + } + const double inv_lambda, dangle_rad, ux, uy, uz, ex, ey, ez; +}; + } // namespace PostRefineResult PostRefineRotationGeometry(const std::vector &outcomes, @@ -128,6 +160,13 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(events.size()) < settings.min_events) return result; + // The rotation-scale fit further down is a single scalar whose whole point is how the residual + // varies ALONG the sweep, so it keeps every event. The cap below ranks by I/sigma, and on the + // crystals that have a stage fault the strong events sit in the middle of the sweep - the part + // that still indexes - so a capped set would leave the ends unrepresented in exactly the fit that + // has to see them. + const std::vector scale_events = events; + constexpr size_t MAX_EVENTS = 20000; if (events.size() > MAX_EVENTS) { std::nth_element(events.begin(), events.begin() + MAX_EVENTS, events.end(), @@ -215,29 +254,109 @@ PostRefineResult PostRefineRotationGeometry(const std::vector {:.3e} => {}", s, axdev, cvA_nom, cvA_ref, result.cell_refined ? "COMMIT" : "reject (kept nominal cell)"); - // Goniometer rotation SCALE, read off the length of the axis vector step A already fits. - // The residual rotates by -angle_rad * axis[], so if the stage turned k times the angle the - // file records, the fit drives |axis| -> k. Nothing applies it and no parameter is added - - // this length was always refined and then normalised away. It is reported only when the - // same cross-validation that gates the cell move says the fit is real, so a fold that - // merely soaked up noise cannot raise the flag. - const double k_fit = std::sqrt(ax_fit[0]*ax_fit[0] + ax_fit[1]*ax_fit[1] + ax_fit[2]*ax_fit[2]); - const bool k_cross_validated = convA && cvA_ref < 0.98 * cvA_nom; - result.rotation_scale = k_cross_validated ? k_fit : 1.0; - // 0.5 %: a well-calibrated stage sits inside it (measured: 36 of 37 rotation datasets peak - // at exactly 1.0000 on a direct scan of this factor), and the one real fault measured 1.2 %. + // ===== Goniometer rotation SCALE k, its own one-parameter fit on the same rocking events ===== + // The angles stored in the file are the COMMANDED ones, so a stage that turned k times as far + // is invisible in the header. Nothing else here can represent it: the cell scale, the axis + // direction, the distance and the beam are all orthogonal to a rotation MAGNITUDE error. Fitted + // after step A so the cell scale and the axis direction are fixed at their committed values and + // k is the only free quantity. + const double u[3] = {axv[0] / axlen, axv[1] / axlen, axv[2] / axlen}; + double phi_c = 0.0, phi_lo = scale_events[0].phi_obs, phi_hi = scale_events[0].phi_obs; + for (const auto &ev : scale_events) { + phi_c += ev.phi_obs; + phi_lo = std::min(phi_lo, ev.phi_obs); + phi_hi = std::max(phi_hi, ev.phi_obs); + } + phi_c /= static_cast(scale_events.size()); + const double sweep_deg = (phi_hi - phi_lo) * 180.0 / PI; + // The reference reciprocal vector turned to the sweep centre, at the committed cell scale. The + // angle then enters the fit measured FROM that centre. A constant crystal missetting about the + // spindle is k with a slope in phi, so measuring the angle from the goniometer's zero instead + // lets a missetting leak into k with gain / - which depends only on where the sweep + // happens to sit. On a short sweep starting near zero that gain is enormous: a 0.14 deg + // missetting on a 10 deg wedge fakes 1.4 % of k. Referred to the sweep centre the leak is + // identically zero at any width, and no parameter has to be added to get it. + std::vector> e_mid(scale_events.size()); + for (size_t e = 0; e < scale_events.size(); ++e) { + const double aa[3] = {-phi_c * u[0], -phi_c * u[1], -phi_c * u[2]}; + const double p[3] = {scale_events[e].e_ref[0] / s, scale_events[e].e_ref[1] / s, + scale_events[e].e_ref[2] / s}; + ceres::AngleAxisRotatePoint(aa, p, e_mid[e].data()); + } + // Robust-loss scale from the scatter the events actually have: it varies by more than a decade + // between datasets, so any fixed constant is either inert or throws away real data. + double rms = 0.0; + for (size_t e = 0; e < scale_events.size(); ++e) { + RotationScaleResidual r(lambda_l, scale_events[e].phi_obs - phi_c, u, e_mid[e].data()); + double one = 1.0, resid = 0.0; + r(&one, &resid); rms += resid * resid; + } + rms = std::sqrt(rms / static_cast(scale_events.size())); + // Fit k over the events whose phi lies outside the given fifth of the sweep (-1 = all of it). + auto solve_scale = [&](int drop_fifth) { + double kv = 1.0; + ceres::Problem p; + for (size_t e = 0; e < scale_events.size(); ++e) { + const int fifth = std::clamp(static_cast( + 5.0 * (scale_events[e].phi_obs - phi_lo) / std::max(1e-9, phi_hi - phi_lo)), 0, 4); + if (fifth == drop_fifth) continue; + p.AddResidualBlock(new ceres::AutoDiffCostFunction( + new RotationScaleResidual(lambda_l, scale_events[e].phi_obs - phi_c, u, e_mid[e].data())), + new ceres::HuberLoss(std::max(1e-12, 2.0 * rms)), &kv); + } + p.SetParameterLowerBound(&kv, 0, 0.95); p.SetParameterUpperBound(&kv, 0, 1.05); + ceres::Solver::Options o; o.linear_solver_type = ceres::DENSE_QR; o.max_num_iterations = 50; + o.num_threads = std::max(1, settings.num_threads); o.logging_type = ceres::LoggingType::SILENT; + ceres::Solver::Summary sum; ceres::Solve(o, &p, &sum); + return sum.IsSolutionUsable() ? kv : 1.0; + }; + const double k_fit = solve_scale(-1); + result.rotation_scale = k_fit; + + // ----- Whether to COMMIT it. A stage fault is rare - 36 of 37 rotation datasets sit at 1.0000 + // on a direct scan - and a 1 % angle correction applied to a healthy dataset would damage it + // silently, so every test below has to pass. + // Preconditions: below these the fit is reported but never acted on. Under ~30 deg of sweep k + // entangles with the axis direction and 10-20 deg truncations of a perfect dataset wander by + // +-0.6 %; a screening wedge must not trigger a correction. + constexpr int MIN_SCALE_EVENTS = 5000; + constexpr double MIN_SCALE_SWEEP_DEG = 30.0; + // T1 significance: 0.5 % is 18 sigma on the between-dataset scatter of healthy stages + // (robust sd 2.8e-4) and still 3.5x below the one measured fault. constexpr double ROTATION_SCALE_TOL = 0.005; - result.rotation_scale_suspect = k_cross_validated && std::fabs(k_fit - 1.0) > ROTATION_SCALE_TOL; + // T2 relevance: the misorientation the error produces at each end of the sweep. A large k over + // a short sweep moves nothing and is not worth correcting. + constexpr double MIN_SCALE_END_ERROR_DEG = 0.5; + // T3 uniformity: a stage error is a ramp present in EVERY part of the sweep, so dropping any + // fifth of it must leave the same k. A second lattice that dominates ONE END of the sweep - + // exactly what happens where the primary stops indexing - fakes a k indistinguishable from a + // real fault on T1 and T2, and is the reason this test is not optional. It replaces the + // hkl-hash split used elsewhere here, which cannot see it: both halves of that split sit at + // the same angles, so anything structured in phi survives in both folds. + constexpr double MIN_SCALE_JACKKNIFE_FRAC = 0.5; + const double end_error_deg = std::fabs(k_fit - 1.0) * sweep_deg / 2.0; + const bool enough_data = static_cast(scale_events.size()) >= MIN_SCALE_EVENTS + && sweep_deg >= MIN_SCALE_SWEEP_DEG; + const bool big_enough = enough_data && std::fabs(k_fit - 1.0) >= ROTATION_SCALE_TOL + && end_error_deg >= MIN_SCALE_END_ERROR_DEG; + double jackknife = 1.0; + if (big_enough) + for (int f = 0; f < 5; ++f) + jackknife = std::min(jackknife, (solve_scale(f) - 1.0) / (k_fit - 1.0)); + result.rotation_scale_suspect = big_enough && jackknife >= MIN_SCALE_JACKKNIFE_FRAC; + logger.Info("Post-refine rotation SCALE: k = {:.5f} over {:.0f} deg of sweep centred on {:.1f} " + "deg ({} events): end error {:.2f} deg, leave-a-fifth-out {:.2f} => {}", + k_fit, sweep_deg, phi_c * 180.0 / PI, scale_events.size(), end_error_deg, jackknife, + result.rotation_scale_suspect ? "COMMIT" + : !enough_data ? "report only (too little sweep or too few events)" + : "reject (kept the stored angles)"); if (result.rotation_scale_suspect) logger.Warning("Goniometer rotation scale looks off by {:+.2f} % (fitted {:.5f}): the stage " "appears to have turned {} than the angles stored in the file, which are the " - "COMMANDED values. This is a hardware calibration fault, not a data problem, " - "and nothing in the processing corrects it - expect inflated mosaicity, a " - "biased cell and lost high-resolution reflections", + "COMMANDED values. This is a hardware calibration fault, not a data problem - " + "left uncorrected it inflates mosaicity, biases the cell and loses " + "high-resolution reflections", 100.0 * (k_fit - 1.0), k_fit, k_fit > 1.0 ? "further" : "less far"); - else if (k_cross_validated) - logger.Info("Goniometer rotation scale {:.5f} (within {:.1f} % of the stored angles)", - k_fit, 100.0 * ROTATION_SCALE_TOL); // Cell (scale s, shape fixed) as the XtalResidual parameter blocks p0/p1/p2, held CONSTANT in step B. double p0[3] = {0, 0, 0}, p1[3] = {0, 0, 0}, p2[3] = {0, 0, 0}; diff --git a/image_analysis/geom_refinement/PostRefine.h b/image_analysis/geom_refinement/PostRefine.h index f9482bfa6..09799d8f6 100644 --- a/image_analysis/geom_refinement/PostRefine.h +++ b/image_analysis/geom_refinement/PostRefine.h @@ -34,13 +34,14 @@ struct PostRefineResult { double beam_y_before_px = 0.0, beam_y_after_px = 0.0; bool cell_refined = false; // GEOM step A (cell scale + axis) passed cross-validation bool detector_refined = false; // GEOM step B (distance + beam) passed cross-validation - // Implied GONIOMETER ROTATION SCALE: step A's axis vector is unnormalised, so the length it fits is - // 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). Reported only - - // nothing here applies it, and it adds no degree of freedom: the parameter was always refined, its - // length was simply normalised away and thrown out. 1.0 = header and stage agree. + // 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 step A as a single free parameter, with the cell scale and the axis + // direction held at their committed values. Always the fitted value; 1.0 = header and stage agree. double rotation_scale = 1.0; - bool rotation_scale_suspect = false; // |rotation_scale - 1| over the tolerance AND cross-validated + // Whether the fit passed every test needed to ACT on it: enough sweep and events, a significant and + // physically relevant size, and the same k from every fifth of the sweep. Only then is it applied. + bool rotation_scale_suspect = false; }; struct PostRefineSettings { diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index c541e3009..99bec2b6d 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -118,13 +118,17 @@ std::string RenderResultReport(const std::string &output_prefix, pr.beam_x_before_px, pr.beam_y_before_px, pr.beam_x_after_px, pr.beam_y_after_px)); Key(os, "GONIOMETER_ROTATION_SCALE", fmt::format("{:.5f}", pr.rotation_scale)); + Key(os, "GONIOMETER_ROTATION_SCALE_SUSPECT", pr.rotation_scale_suspect ? "TRUE" : "FALSE"); os << "\n GONIOMETER_ROTATION_SCALE is the factor by which the stage actually turned relative to\n" - << " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It is\n" - << " reported only; nothing here corrects for it.\n"; + << " the angles stored in the file (which are the commanded ones). 1.0 = they agree. It drives\n" + << " the second integration pass only when SUSPECT is TRUE - both cross-validated and outside\n" + << " the tolerance - since a stage that is in fact well calibrated must be left alone. A\n" + << " manual --rotation-scale replaces it and is applied to both passes.\n"; if (pr.rotation_scale_suspect) warnings.emplace_back(fmt::format( "The goniometer turned by a factor {:.5f} of the angles stored in the file - the " - "stage rotation looks mis-calibrated by {:+.2f}%", + "stage rotation looks mis-calibrated by {:+.2f}%. The correction was applied to this " + "run, but the fault is in the hardware and should be fixed there", pr.rotation_scale, 100.0 * (pr.rotation_scale - 1.0))); } diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 05a3dcf2f..efae37088 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -75,6 +75,21 @@ namespace { return ret; } + // Apply a goniometer rotation SCALE: the stage turned k times the angle the file records. What is + // wrong is the SWEEP, not where it began, so the per-frame increment and the per-frame oscillation + // width both scale by k while the starting angle is left alone - the stage reached that position + // before the sweep started. A constant offset in phi is in any case exactly degenerate with a + // rotation of the crystal orientation about the spindle, which indexing refines away, so anchoring + // the stretched sweep at image 0 costs nothing and keeps the first frame at the angle the file gives + // it. Every angle in the pipeline comes from this object, so scaling it is the whole correction. + GoniometerAxis ScaleRotation(const GoniometerAxis &g, float k) { + GoniometerAxis scaled(g.GetName(), g.GetStart_deg(), g.GetIncrement_deg() * k, + g.GetAxis(), g.GetHelicalStep()); + if (const auto wedge = g.GetScreeningWedge()) + scaled.ScreeningWedge(*wedge * k); + return scaled; + } + } Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment, @@ -84,6 +99,14 @@ Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment, // Bit 9 describes where THIS run found the beam stop, so a mask read back from a file that // already carries one starts clear; the user mask (bit 8) is left as it was loaded. pixel_mask_.ClearBeamStopMask(experiment_); + + // A manually asserted stage calibration applies to everything, before anything reads an angle. + if (config_.rotation_scale.has_value()) + if (const auto g = experiment_.GetGoniometer()) { + Logger("Rugnux").Info("Goniometer rotation scale {:.5f} applied from the command line", + *config_.rotation_scale); + experiment_.Goniometer(ScaleRotation(*g, *config_.rotation_scale)); + } } void Rugnux::FindBeamStop(int start_image, int images_to_process, int frame_count) { @@ -320,6 +343,7 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { const std::string base_prefix = config_.output_prefix; const auto gonio_snapshot = experiment_.GetGoniometer(); prepass_detector_geometry_.reset(); + prepass_rotation_scale_.reset(); prepass_result_.reset(); force_rotation_result_.reset(); prepass_sg_reindexed_ = false; @@ -345,6 +369,21 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { const auto &g = *prepass_detector_geometry_; experiment_.BeamX_pxl(g[0]).BeamY_pxl(g[1]).DetectorDistance_mm(g[2]); } + // ... and the goniometer rotation scale, on the same measure-then-re-integrate footing: the angles + // in the file are the commanded ones, so a stage that ran fast is a geometry error like any other. + if (prepass_rotation_scale_ && gonio_snapshot) { + experiment_.Goniometer(ScaleRotation(*gonio_snapshot, *prepass_rotation_scale_)); + // The pre-pass mosaicity is a width in degrees fitted against the angles the second pass has + // just stopped using, and the override can only ever raise the second pass's own estimate (it + // takes the larger of the two). Carrying it over would hold the second pass at the rocking + // width the uncorrected angles produced - the correction half-applied. Drop it and let the + // second pass fit its own. + prepass_mosaicity_.clear(); + logger.Info("Two-pass: goniometer rotation scale {:.5f} will drive the second integration pass " + "(oscillation {:.4f} -> {:.4f} deg per image)", *prepass_rotation_scale_, + gonio_snapshot->GetIncrement_deg(), + gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_); + } // Space group: the second pass RE-INDEXES DE NOVO (clear the group here) so the indexer's pseudo- // symmetry safeguards recover the true cell - reusing pass-1's group in the indexer forces build_sr // onto a wrong / doubled cell (a huge oblique cell collapses; a pseudo-centred cell doubles). Pass-1's @@ -383,6 +422,12 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { "spurious supercell; re-running the second pass with pass-1's result forced", v2, v2 / v1, v1); force_rotation_result_ = *prepass_result_; + // Pass 1's result carries pass 1's goniometer, and the per-image path takes its angles + // from there - forcing it whole would put the uncorrected angles back and silently undo + // the rotation scale for the re-run. + if (prepass_rotation_scale_ && force_rotation_result_->axis) + force_rotation_result_->axis = + ScaleRotation(*force_rotation_result_->axis, *prepass_rotation_scale_); pass2 = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); pass2.post_refine = pass1.post_refine; pass2.pass_number = 3; @@ -429,6 +474,8 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { // only costs a pass on a crystal that was going to be wrong otherwise. experiment_.BeamX_pxl(header_geometry[0]).BeamY_pxl(header_geometry[1]) .DetectorDistance_mm(header_geometry[2]); + if (prepass_rotation_scale_ && gonio_snapshot) + experiment_.Goniometer(*gonio_snapshot); // and the header rotation angles with it experiment_.SpaceGroupNumber(std::nullopt); config_.output_prefix = base_prefix; auto redo = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false); @@ -449,7 +496,7 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { "CC1/2 {:.3f} vs {:.3f})", compl2, completeness(pass1), cc2, cc1); } if (pass2.pass_decision.empty()) - pass2.pass_decision = prepass_detector_geometry_ + pass2.pass_decision = (prepass_detector_geometry_ || prepass_rotation_scale_) ? "post-refined geometry adopted" : "the post-refinement committed no geometry change, so this pass reproduces the first"; return pass2; @@ -1542,6 +1589,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Two-pass: geometry post-refine committed no detector change " "- 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. + if (pr.rotation_scale_suspect && !config_.rotation_scale.has_value()) + prepass_rotation_scale_ = static_cast(pr.rotation_scale); } } } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 3ea8cd43b..dfba2bce0 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -89,6 +89,12 @@ struct ProcessConfig { // CANONICAL "_*" output; the header-geometry first pass is kept alongside as "_01_*". bool rotation_postrefine_geometry = false; + // Manual goniometer rotation scale (--rotation-scale): the factor by which the stage really turned + // relative to the angles stored in the file, which are the COMMANDED ones. Applied to the goniometer + // before anything reads it, so BOTH passes use the corrected angles and the post-refine then reports + // whatever scale error is LEFT. Overrides the fitted correction. + std::optional rotation_scale; + // Scaling / merging (FullAnalysis). reference_data (from a reference MTZ) also drives the // per-image live scaling path when present. bool run_scaling = false; @@ -212,6 +218,10 @@ class Rugnux { // committed no change (the second pass then reproduces the first). std::optional> prepass_detector_geometry_; + // Two-pass geometry pre-pass: the goniometer rotation scale the post-refine fitted and flagged as a + // stage fault, applied by Run() to the second pass's goniometer. Empty when the fit found nothing. + std::optional prepass_rotation_scale_; + // Two-pass geometry pre-pass: pass-1's FULL indexing result (correct lattice + orientation + refined // search-result metadata + geometry + axis, found at the header geometry). Kept as a fallback seed: if the // second pass's de-novo re-index at the refined geometry collapses onto a spurious supercell (a bistable diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 5ddd94bd4..cebdf72c7 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -184,6 +184,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config, // disable flag only when the GUI turned it off on a rotation dataset (to reproduce that choice). if (config.rotation_indexing && !config.rotation_postrefine_geometry) args.emplace_back("--rotation-no-postrefine"); + if (config.rotation_scale.has_value()) + add("--rotation-scale", num(*config.rotation_scale)); // Merging is on by default; emit --no-merge only when it was turned off. if (config.run_scaling) { diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index f2b8022f2..3e0461917 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -118,6 +118,7 @@ void print_usage() { std::cout << " -R, --two-pass-rotation[=num] Two-pass offline rotation indexing (default for goniometer data; optional first-pass image count, default: 100)" << std::endl; std::cout << " --single-pass-rotation[=num] Use online-like single-pass rotation indexing (optional: min angular range deg)" << std::endl; std::cout << " --force-rotation-lattice Force rotation indexer with external lattice (in Angstrom) : \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, skips first pass)" << std::endl; + std::cout << " --rotation-scale Goniometer rotation scale: the stage turned k times the angle stored in the file (the commanded one). Applied to both passes; overrides the fitted correction" << std::endl; std::cout << " --rotation-no-postrefine Disable the (default-on) two-pass rotation post-refine (post-refine detector distance/beam + cell/axis, then re-integrate; the refined pass is the canonical _* output, the header-geometry pass is kept as _01_*)" << std::endl; std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; std::cout << " -S, --space-group Space group number (92) or symbol (P43212) - for indexing and scaling" << std::endl; @@ -210,6 +211,7 @@ enum { OPT_SINGLE_PASS_ROTATION, OPT_FORCE_ROTATION_LATTICE, OPT_ROTATION_NO_POSTREFINE, + OPT_ROTATION_SCALE, OPT_BACKGROUND_CLIP, OPT_BACKGROUND_RADIAL, OPT_REFINE_GEOMETRY, @@ -307,6 +309,7 @@ static option long_options[] = { {"polarization", required_argument, nullptr, OPT_POLARIZATION}, {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, {"rotation-no-postrefine", no_argument, nullptr, OPT_ROTATION_NO_POSTREFINE}, + {"rotation-scale", required_argument, nullptr, OPT_ROTATION_SCALE}, {"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY}, {"detect-beam-stop", optional_argument, nullptr, OPT_DETECT_BEAM_STOP}, @@ -573,6 +576,7 @@ static int RunRugnux(int argc, char **argv) { // file when it has them, so those options would otherwise not reach the pass that determines the // lattice - the setting would appear to do nothing at all on rotation data. bool rotation_postrefine_geometry = true; // default on; --rotation-no-postrefine disables it + std::optional rotation_scale; // --rotation-scale: asserted stage calibration int rotation_indexing_image_count = 100; std::optional rotation_indexing_range; bool run_scaling = true; // merge is on by default; --no-merge turns it off @@ -700,6 +704,9 @@ static int RunRugnux(int argc, char **argv) { case OPT_ROTATION_NO_POSTREFINE: rotation_postrefine_geometry = false; break; + case OPT_ROTATION_SCALE: + rotation_scale = parse_number_arg(optarg, "--rotation-scale", logger, 0.9f, 1.1f); + break; case OPT_DETECT_BEAM_STOP: // Frames projected to find the shadow. The default is what the detection was validated // on; fewer leaves the background too sparsely counted to tell a shadow from noise. @@ -2030,6 +2037,7 @@ static int RunRugnux(int argc, char **argv) { config.two_pass_rotation = two_pass_rotation; config.detect_beam_stop = detect_beam_stop; config.rotation_postrefine_geometry = rotation_postrefine_geometry; + config.rotation_scale = rotation_scale; config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; config.refine_geometry = refine_geometry;