diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 859ca9929..328b6a989 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,6 +8,7 @@ * `rugnux` places a detector swung out on a 2theta arm where the file says it stands. * `rugnux` writes the unmerged MTZ by default, with a P1 merge beside it, so a wrong space group can be re-merged without reprocessing. * `rugnux` handles symmetry better: the lattice, the point group, the setting and the systematic absences. +* `rugnux` fits the direction of the goniometer axis and not its length, so the cell chosen by the first pass is the one its own refinement scored. * `rugnux` reports the detector geometry it measured - the direct beam, the tilt and the beam centre - and what a single sweep can and cannot determine. * `rugnux` reports the resolution the CC1/2 fit reached and the strong-direction diffraction limit. * `rugnux` reports twinning measured before and after the space group was decided, and no longer reports it when the L-test contradicts it. diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 58c357762..55d1756bb 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -455,6 +455,18 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, if (!data.refine_rotation_axis) { problem.SetParameterBlockConstant(rot_vec); + } else { + // Only the DIRECTION of the goniometer axis is a parameter. The residual applies + // angle_rad * |rot_vec|, so a free three-vector also fits a rotation SCALE - which + // GoniometerAxis::Axis() then normalises away, leaving the candidate scored by + // RotationIndexer::accumulate() under a rotation model the fit did not use. Measured + // over the corpus, that length reached 1.2 % and the fit/score disagreement a whole + // degree of goniometer angle. It is not a usable measurement either: on synthetic + // data it recovers 54 % of a known scale error, repeated first passes on one dataset + // disagree with each other in SIGN, and on the one dataset with a real 1.3 % stage + // fault it comes out negative. The rotation scale is measured properly, once, with + // four gates and a jackknife, in PostRefine. + problem.SetManifold(rot_vec, new ceres::SphereManifold<3>); } if (!data.refine_unit_cell) { diff --git a/tests/XtalOptimizerTest.cpp b/tests/XtalOptimizerTest.cpp index 3e4a69d7f..70660c80b 100644 --- a/tests/XtalOptimizerTest.cpp +++ b/tests/XtalOptimizerTest.cpp @@ -669,5 +669,75 @@ TEST_CASE("XtalOptimizer_refine_rotation_axis") { CHECK(fabsf(xtal_opt.axis->GetAxis().z) < 0.01f); } +// The rotation residual applies `angle_rad * |rot_vec|`, so the goniometer axis block must be +// constrained to unit length or the solve also fits a rotation SCALE that GoniometerAxis::Axis() +// then normalises away. Generate spots whose true rotation is k times the angles the optimizer is +// handed - k = 1 for a healthy goniometer, k = 1.01 for a stage that turned 1 % further than it was +// commanded - and check the axis DIRECTION and the cell come back either way. (Before the axis block +// was put on the unit sphere, the free length absorbed part of that 1 %: on this noise-free data it +// reached only 1.00538 of the 1.01 it would have had to reach to be a measurement.) +TEST_CASE("XtalOptimizer_rotation_axis_direction_only") { + DiffractionExperiment exp_i; + exp_i.IncidentEnergy_keV(WVL_1A_IN_KEV) + .BeamX_pxl(1000).BeamY_pxl(1000) + .DetectorDistance_mm(200); + + const CrystalLattice latt_base(40, 50, 80, 90, 95, 90); + const auto uc_ref = latt_base.GetUnitCell(); + + const double k_true = GENERATE(1.000, 1.010); + + // The angles the FILE records. The stage actually turned k_true times as far. + GoniometerAxis nominal("omega", 0.0f, 3.0f, Coord(1, 0, 0), std::nullopt); + + BraggPredictionSettings prediction_settings{.high_res_A = 1.5, .ewald_dist_cutoff = 0.002}; + BraggPrediction prediction; + + const size_t nimages = 40; + std::vector> spots(nimages); + for (size_t img = 0; img < nimages; ++img) { + const float nominal_deg = nominal.GetAngle_deg(img) + nominal.GetWedge_deg() / 2.0f; + const RotMatrix rot = nominal.GetTransformationAngle(static_cast(nominal_deg * k_true)); + const CrystalLattice latt_img = latt_base.Multiply(rot.transpose()); + const auto n = prediction.Calc(exp_i, latt_img, prediction_settings); + for (int i = 0; i < n; ++i) { + const auto &r = prediction.GetReflections().at(i); + SpotToSave s{}; + s.x = r.predicted_x; + s.y = r.predicted_y; + s.image = static_cast(img); + s.intensity = 1.0f; + s.phi = nominal_deg; + s.indexed = true; + spots.at(img).push_back(s); + } + } + + XtalOptimizerData xtal_opt{}; + xtal_opt.latt = latt_base; + xtal_opt.geom.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(200.0); + xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic; + xtal_opt.axis = GoniometerAxis("omega", 0.0f, 3.0f, Coord(0.999, 0.02, 0.02).Normalize(), + std::nullopt); + xtal_opt.min_spots = 200; + xtal_opt.refine_beam_center = true; + xtal_opt.refine_detector_angles = true; + xtal_opt.refine_rotation_axis = true; + xtal_opt.max_iterations = 200; + + REQUIRE(XtalOptimizer(xtal_opt, spots)); + + const Coord axis_out = xtal_opt.axis->GetAxis(); + CHECK(std::fabs(axis_out.Length() - 1.0f) < 1e-5f); + CHECK(std::fabs(axis_out.x - 1.0f) < 0.01f); + CHECK(std::fabs(axis_out.y) < 0.01f); + CHECK(std::fabs(axis_out.z) < 0.01f); + + const auto uc_out = xtal_opt.latt.GetUnitCell(); + CHECK(std::fabs(uc_ref.a - uc_out.a) < 0.2f); + CHECK(std::fabs(uc_ref.b - uc_out.b) < 0.2f); + CHECK(std::fabs(uc_ref.c - uc_out.c) < 0.3f); +} + // --- helpers for lattice sanity tests --- #include