diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index de41f64e..9391fe29 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -468,7 +468,7 @@ The refinement jointly optimizes, depending on mode and constraints: - crystal orientation (a global rotation), - unit-cell parameters, with constraints determined by inferred crystal system. -By default only the beam center, unit cell and crystal orientation are refined; the detector distance, tilt angles and rotation-axis direction are held fixed unless explicitly enabled. A lighter **orientation-only** mode refines just the crystal orientation, for stills whose geometry is already trusted. It carries a weak small-rotation prior penalising the whole angle-axis vector (all three components, at a low weight); what it is there for is the poorly-determined out-of-plane component, which is the one the data barely constrain. +The detector distance is not refined against one crystal's spots at all: the positional residual leaves it degenerate with the cell scale, so it is fitted separately - by the rotation post-refinement, which holds the cell at the value its own cell/axis step settled on, and by the stills `--refine-geometry` bundle. Per image, the beam centre and the crystal orientation are refined, and the unit cell as well for stills. The first-pass rotation indexing refines the detector tilt and the rotation-axis direction too, against the spots accumulated across the sweep; everywhere else both are held fixed, because on a single crystal a tilt is absorbed almost exactly by the crystal orientation. A lighter **orientation-only** mode refines just the crystal orientation, for stills whose geometry is already trusted. It carries a weak small-rotation prior penalising the whole angle-axis vector (all three components, at a low weight); what it is there for is the poorly-determined out-of-plane component, which is the one the data barely constrain. For higher symmetries, constraints are enforced, e.g. - cubic: $a=b=c,\ \alpha=\beta=\gamma=90^\circ$, diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index d5a1034b..7e43d049 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -308,7 +308,6 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In .min_angle_deg = 30.0f, .max_angle_deg = 150.0f, .refine_beam_center = true, - .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = !experiment.IsRotationIndexing(), // The whole spot list is passed below, not the indexed subset, so on weak images most of what @@ -402,7 +401,6 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In .crystal_system = data.crystal_system, .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), .refine_beam_center = false, - .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = false, .refine_rotation_axis = false, diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index cedc68a8..11b438a3 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -475,32 +475,47 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(std::max(1, std::min(nthreads, static_cast(n_terms)))); const int cost_chunk = (n_terms + cost_nt - 1) / cost_nt; - const auto cost_by_fifth = [&](double k) { - std::vector> per_chunk(cost_nt); + // Several k are always wanted at once (the grid below asks for 101), and they all sweep the + // same event list, so sweep it ONCE and evaluate every k on each event while it is still in + // registers. The per-thread accumulator is one slot per (k, fifth) - 4 kB for the grid, small + // enough to stay in L1 - against re-reading the whole term array once per k. Each (k, fifth) + // still receives its events in the same order and the chunks are still summed in chunk order, + // so the sums are the ones a k-at-a-time loop produced, bit for bit. + const auto cost_grid = [&](const std::vector &ks) { + const int nk = static_cast(ks.size()); + std::vector>> per_chunk( + cost_nt, std::vector>(nk)); ParallelChunks(n_terms, nthreads, [&](int lo, int hi) { - std::array acc{}; + std::vector> acc(nk); for (int e = lo; e < hi; ++e) { - const double r = residual_at(terms[e], k); - const double s2 = r * r; - acc[fifth_of[e]] += (s2 <= huber_d2) ? s2 - : (2.0 * huber_delta * std::sqrt(s2) - huber_d2); + const ScaleTerm &term = terms[e]; + const int fifth = fifth_of[e]; + for (int g = 0; g < nk; ++g) { + const double r = residual_at(term, ks[g]); + const double s2 = r * r; + acc[g][fifth] += (s2 <= huber_d2) ? s2 + : (2.0 * huber_delta * std::sqrt(s2) - huber_d2); + } } - per_chunk[lo / cost_chunk] = acc; + per_chunk[lo / cost_chunk] = std::move(acc); }); - std::array total{}; + std::vector> total(nk); for (const auto &acc : per_chunk) - for (int j = 0; j < 5; ++j) total[j] += acc[j]; + for (int g = 0; g < nk; ++g) + for (int j = 0; j < 5; ++j) total[g][j] += acc[g][j]; return total; }; + const auto cost_by_fifth = [&](double k) { return cost_grid({k})[0]; }; // Scan the interval Ceres was bounded to, then close in. No event's phase can move by more than // a fraction of a period over an interval this narrow, so the objective has no structure the // grid could step over; the refinement is only there to place the minimum precisely. constexpr int SCALE_GRID = 101; constexpr double SCALE_K_LO = 0.95, SCALE_K_HI = 1.05; - std::vector> grid(SCALE_GRID); + std::vector grid_k(SCALE_GRID); for (int g = 0; g < SCALE_GRID; ++g) - grid[g] = cost_by_fifth(SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1)); + grid_k[g] = SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1); + const std::vector> grid = cost_grid(grid_k); auto solve_scale = [&](int drop_fifth) { const auto total = [&](const std::array &f) { double t = 0.0; @@ -601,23 +616,45 @@ PostRefineResult PostRefineRotationGeometry(const std::vector obs; - for (size_t i = 0; i < n_pts; ++i) - if (std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y)) obs.push_back(&pts[i]); + // Count first, then fill, exactly as the partial gather above does and for the same reason: + // this walks the same tens of millions of partials, and a pointer vector grown by push_back + // copies itself every time it doubles. new[] rather than a sized vector so the array is not + // zeroed on one thread before the parallel fill overwrites it. The fill lands in the order + // the serial loop produced, so the selection below sees the same sequence it always did. + const int n_obs_chunks = static_cast(std::clamp(nthreads, 1, + std::max(1, n_pts))); + const size_t obs_chunk = (n_pts + n_obs_chunks - 1) / n_obs_chunks; + std::vector obs_offset(n_obs_chunks + 1, 0); + const auto keep_obs = [&](size_t i) { + return std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y); + }; + ParallelChunks(static_cast(n_pts), nthreads, [&](int lo, int hi) { + size_t keep = 0; + for (int i = lo; i < hi; ++i) if (keep_obs(i)) keep++; + obs_offset[static_cast(lo) / obs_chunk + 1] = keep; + }); + for (int c = 0; c < n_obs_chunks; ++c) obs_offset[c + 1] += obs_offset[c]; + size_t n_obs = obs_offset[n_obs_chunks]; + std::unique_ptr obs(new const Partial *[n_obs]); + ParallelChunks(static_cast(n_pts), nthreads, [&](int lo, int hi) { + size_t at = obs_offset[static_cast(lo) / obs_chunk]; + for (int i = lo; i < hi; ++i) if (keep_obs(i)) obs[at++] = &pts[i]; + }); constexpr size_t MAX_OBS = 20000; - if (obs.size() > MAX_OBS) { - std::nth_element(obs.begin(), obs.begin() + MAX_OBS, obs.end(), + if (n_obs > MAX_OBS) { + std::nth_element(obs.get(), obs.get() + MAX_OBS, obs.get() + n_obs, [](const Partial *a, const Partial *b) { return a->I / std::max(1e-9, static_cast(a->sigma)) > b->I / std::max(1e-9, static_cast(b->sigma)); }); - obs.resize(MAX_OBS); + n_obs = MAX_OBS; } - result.obs_used = static_cast(obs.size()); + result.obs_used = static_cast(n_obs); const double beam_x0 = nominal_geom.GetBeamX_pxl(), beam_y0 = nominal_geom.GetBeamY_pxl(); const double dist0 = nominal_geom.GetDetectorDistance_mm(); auto pos_cost = [&](Subset s, const double beam[2], const double dist[1]) { double c = 0.0; int n = 0; - for (const Partial *pp : obs) { + for (size_t oi = 0; oi < n_obs; ++oi) { + const Partial *pp = obs[oi]; if (!in(pp->h, pp->k, pp->l, s)) continue; XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img), pp->h, pp->k, pp->l, sys); @@ -630,21 +667,20 @@ PostRefineResult PostRefineRotationGeometry(const std::vectorh, pp->k, pp->l, s)) continue; - p.AddResidualBlock(new ceres::AutoDiffCostFunction( - new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img), - pp->h, pp->k, pp->l, sys)), - new ceres::CauchyLoss(0.02), beam, dist, - const_cast(det_rot), rot_vec, p0, p1, p2); + p.AddResidualBlock(new ceres::AutoDiffCostFunction( + new XtalResidualBeamDistance( + XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img), + pp->h, pp->k, pp->l, sys), + det_rot, rot_vec, p0, p1, p2)), + new ceres::CauchyLoss(0.02), beam, dist); } if (p.NumResidualBlocks() == 0) { beam_out[0] = beam_x0; beam_out[1] = beam_y0; dist_out = dist0; return false; } - // Everything but the beam and the distance is held at its step-A value. Once per problem: - // the block only has to exist, and repeating it per observation was up to MAX_OBS times - // five calls for the same five blocks. - p.SetParameterBlockConstant(const_cast(det_rot)); - p.SetParameterBlockConstant(rot_vec); - p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2); + // Everything but the beam and the distance is held at its step-A value, and the residual + // above bakes those five blocks in rather than declaring them and freezing them, so there + // is nothing left to hold constant here. p.SetParameterLowerBound(dist, 0, dist0 * 0.95); p.SetParameterUpperBound(dist, 0, dist0 * 1.05); for (int j = 0; j < 2; ++j) { p.SetParameterLowerBound(beam, j, beam[j] - 15.0); p.SetParameterUpperBound(beam, j, beam[j] + 15.0); } @@ -655,7 +691,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector= static_cast(settings.min_events)) { + if (n_obs >= static_cast(settings.min_events)) { double beam_fit[2], dist_fit; const bool convB = solve_detector(FIT, beam_fit, dist_fit); const double b_nom[2] = {beam_x0, beam_y0}, d_nom[1] = {dist0}; @@ -679,7 +715,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector( - new XtalResidual(residual)), + new ceres::AutoDiffCostFunction( + new XtalResidualFixedDistance(residual, distance_mm)), loss, beam, - &distance_mm, detector_rot, rot_vec, latt_vec0, @@ -318,14 +316,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, // Distance, detector angles, rotation axis and cell are parameter blocks only in the general // seven-block residual; the reduced one bakes them in, so there is nothing left to configure. if (!beam_and_orientation_only) { - if (!data.refine_distance_mm) - problem.SetParameterBlockConstant(&distance_mm); - else { - const double dist_range = 0.1; - problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 - dist_range)); - problem.SetParameterUpperBound(&distance_mm, 0, distance_mm * (1.0 + dist_range)); - } - if (!data.refine_detector_angles) { problem.SetParameterBlockConstant(detector_rot); } else { @@ -374,7 +364,14 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, // Configure solver ceres::Solver::Options options; - options.linear_solver_type = ceres::DENSE_QR; + // Normal equations, not QR. The problem is very tall and thin - thousands of spots against at + // most 17 parameters - and that is the shape DENSE_QR handles worst: it copies the Jacobian out + // of Ceres' row-major storage into a column-major buffer on every solve, and Eigen's blocked + // Householder then degenerates to the unblocked path because its block size is min(48, columns). + // Accumulating J^T J reads the Jacobian once instead. Both solve the same damped system, so the + // step is the same to round-off; the column scaling Ceres applies by default and the LM diagonal + // keep the squared condition number in hand. + options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY; options.minimizer_progress_to_stdout = false; if (data.max_iterations > 0) options.max_num_iterations = data.max_iterations; @@ -400,8 +397,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]); } - if (data.refine_distance_mm) - data.geom.DetectorDistance_mm(distance_mm); if (data.refine_detector_angles) data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]); @@ -516,7 +511,7 @@ bool XtalOptimizerRotationOnly(XtalOptimizerData &data, } ceres::Solver::Options options; - options.linear_solver_type = ceres::DENSE_QR; + options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY; // tall and thin, as above options.minimizer_progress_to_stdout = false; if (data.max_iterations > 0) options.max_num_iterations = data.max_iterations; diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h index 35375be6..6d107a63 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.h +++ b/image_analysis/geom_refinement/XtalOptimizer.h @@ -23,7 +23,6 @@ struct XtalOptimizerData { float max_angle_deg = 120.0f; bool refine_beam_center = true; - bool refine_distance_mm = false; bool refine_detector_angles = false; bool refine_unit_cell = true; // This refines unit cell size + angles - orientation is always refined bool refine_rotation_axis = false; diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index 93f8dcc3..72c28145 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -228,32 +228,32 @@ struct XtalResidual { T col1_unrot[3] = {B(0, 1) * L1, B(1, 1) * L1, B(2, 1) * L1}; T col2_unrot[3] = {B(0, 2) * L2, B(1, 2) * L2, B(2, 2) * L2}; - T col0_rot[3], col1_rot[3], col2_rot[3]; - const AngleAxisRotator rot_p0(p0); - rot_p0.Rotate(col0_unrot, col0_rot); - rot_p0.Rotate(col1_unrot, col1_rot); - rot_p0.Rotate(col2_unrot, col2_rot); + // Build the reciprocal vector in the UNROTATED frame and turn it once at the end, rather than + // rotating all three direct columns first. A rotation commutes with the cross product and leaves + // the triple product alone - (Rb)x(Rc) = R(bxc), and (Ra).(R(bxc)) = a.(bxc) - so the reciprocal + // basis of the rotated cell is the rotation of the unrotated one, and h a* + k b* + l c* is the + // rotation of the same combination built from the unrotated columns. Three rotations become one, + // and the cross products now run on columns that are mostly zero for every system but triclinic. + const Eigen::Matrix a_unrot(col0_unrot[0], col0_unrot[1], col0_unrot[2]); + const Eigen::Matrix b_unrot(col1_unrot[0], col1_unrot[1], col1_unrot[2]); + const Eigen::Matrix c_unrot(col2_unrot[0], col2_unrot[1], col2_unrot[2]); - const Eigen::Matrix A(col0_rot[0], col0_rot[1], col0_rot[2]); - const Eigen::Matrix Bv(col1_rot[0], col1_rot[1], col1_rot[2]); - const Eigen::Matrix C(col2_rot[0], col2_rot[1], col2_rot[2]); + const Eigen::Matrix bxc = b_unrot.cross(c_unrot); + const Eigen::Matrix cxa = c_unrot.cross(a_unrot); + const Eigen::Matrix axb = a_unrot.cross(b_unrot); - const Eigen::Matrix BxC = Bv.cross(C); - const Eigen::Matrix CxA = C.cross(A); - const Eigen::Matrix AxB = A.cross(Bv); - - const T V = A.dot(BxC); - const T invV = T(1) / V; - - const Eigen::Matrix Astar = BxC * invV; - const Eigen::Matrix Bstar = CxA * invV; - const Eigen::Matrix Cstar = AxB * invV; + const T invV = T(1) / a_unrot.dot(bxc); const T h = T(exp_h); const T k = T(exp_k); const T l = T(exp_l); - const Eigen::Matrix e_pred_recip = Astar * h + Bstar * k + Cstar * l; + const Eigen::Matrix recip_unrot = (bxc * h + cxa * k + axb * l) * invV; + + T recip_in[3] = {recip_unrot[0], recip_unrot[1], recip_unrot[2]}; + T e_pred_recip[3]; + const AngleAxisRotator rot_p0(p0); + rot_p0.Rotate(recip_in, e_pred_recip); residual[0] = e_obs_recip[0] - e_pred_recip[0]; residual[1] = e_obs_recip[1] - e_pred_recip[1]; @@ -273,6 +273,69 @@ struct XtalResidual { gemmi::CrystalSystem symmetry; }; +// Same residual with the detector angles, the rotation axis and the whole crystal - orientation, cell +// lengths, cell angles - baked in as constants, so that only beam(2) and distance(1) remain parameter +// blocks. That is the rotation post-refinement's detector step, which deliberately holds the crystal at +// the value its cell/axis step settled on: the seven-block form differentiates 17 parameters in order to +// use 3, this one runs on Jet. It forwards to XtalResidual with the very same values, so the +// residual and the beam/distance columns of its Jacobian are unchanged. +struct XtalResidualBeamDistance { + XtalResidualBeamDistance(const XtalResidual &residual, + const double *detector_rot, + const double *rotation_axis, + const double *orientation, + const double *uc_len, + const double *uc_angle) + : residual(residual), + detector_rot{detector_rot[0], detector_rot[1]}, + rotation_axis{rotation_axis[0], rotation_axis[1], rotation_axis[2]}, + orientation{orientation[0], orientation[1], orientation[2]}, + uc_len{uc_len[0], uc_len[1], uc_len[2]}, + uc_angle{uc_angle[0], uc_angle[1], uc_angle[2]} { + } + + template + bool operator()(const T *const beam, const T *const distance, T *residual_out) const { + const T rot[2] = {T(detector_rot[0]), T(detector_rot[1])}; + const T axis[3] = {T(rotation_axis[0]), T(rotation_axis[1]), T(rotation_axis[2])}; + const T p0[3] = {T(orientation[0]), T(orientation[1]), T(orientation[2])}; + const T p1[3] = {T(uc_len[0]), T(uc_len[1]), T(uc_len[2])}; + const T p2[3] = {T(uc_angle[0]), T(uc_angle[1]), T(uc_angle[2])}; + return residual(beam, distance, rot, axis, p0, p1, p2, residual_out); + } + + const XtalResidual residual; + const double detector_rot[2]; + const double rotation_axis[3]; + const double orientation[3]; + const double uc_len[3]; + const double uc_angle[3]; +}; + +// Same residual with only the detector distance baked in, leaving the other six blocks free: 16 +// declared parameters instead of 17. Nothing in rugnux refines the distance from a single crystal - the +// rotation post-refinement fits it globally, with the cell held fixed so the two are no longer +// degenerate - so the block was only ever declared in order to be frozen. Dropping it is worth more +// than one parameter in seventeen suggests: Jet keeps its derivatives in an Eigen vector, and +// 17 doubles is one lane past the four-wide boundary, so it costs five packet operations and cannot be +// over-aligned where 16 costs four and can. The distance is forwarded unchanged, so the residual and +// every remaining Jacobian column are the same. +struct XtalResidualFixedDistance { + XtalResidualFixedDistance(const XtalResidual &residual, double distance_mm) + : residual(residual), distance_mm(distance_mm) { + } + + template + bool operator()(const T *const beam, const T *const detector_rot, const T *const rotation_axis, + const T *const p0, const T *const p1, const T *const p2, T *residual_out) const { + const T distance[1] = {T(distance_mm)}; + return residual(beam, distance, detector_rot, rotation_axis, p0, p1, p2, residual_out); + } + + const XtalResidual residual; + const double distance_mm; +}; + // Same residual with distance, detector angles, rotation axis and unit cell baked in as constants, so // that only beam(2) and orientation(3) remain parameter blocks. Ceres sizes its autodiff dual numbers // from the DECLARED blocks, not from which of them the caller then holds constant, so the seven-block diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 16965c50..8b2c62e3 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -212,7 +212,6 @@ void RotationIndexer::RunIndexing() { .min_angle_deg = 30.0f, .max_angle_deg = 150.0f, .refine_beam_center = true, - .refine_distance_mm = false, .refine_detector_angles = true, .refine_rotation_axis = true, .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), @@ -425,7 +424,6 @@ void RotationIndexer::RunIndexing() { .crystal_system = search_result_.system, .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), .refine_beam_center = false, - .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = false, .refine_rotation_axis = false, diff --git a/tests/XtalOptimizerTest.cpp b/tests/XtalOptimizerTest.cpp index f200693a..3e4a69d7 100644 --- a/tests/XtalOptimizerTest.cpp +++ b/tests/XtalOptimizerTest.cpp @@ -543,15 +543,17 @@ TEST_CASE("XtalOptimizer_rotation") { // direction that XtalOptimizer now deliberately restrains toward the header rather than refining // freely (see the BeamComponentPrior in XtalOptimizer.cpp). So only the perpendicular component // (Y) is seeded off-truth to exercise beam-centre recovery; X is seeded at its true value. + // The distance is seeded at its true value too: XtalOptimizer does not refine it - the rotation + // post-refinement fits it globally instead, where the cell is held fixed and the two are no longer + // degenerate - so a distance seeded off-truth here would simply be absorbed by the cell. XtalOptimizerData xtal_opt{}; xtal_opt.latt = CrystalLattice(39.7f, 50.6f, 79.6f, 90.0f, 94.5f, 90.5f); - xtal_opt.geom.BeamX_pxl(1000).BeamY_pxl(997).DetectorDistance_mm(203.0) + xtal_opt.geom.BeamX_pxl(1000).BeamY_pxl(997).DetectorDistance_mm(200.0) .PoniRot1_rad(0.01).PoniRot2_rad(0.02); xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic; xtal_opt.axis = axis; xtal_opt.min_spots = 200; xtal_opt.refine_beam_center = true; - xtal_opt.refine_distance_mm = true; xtal_opt.refine_detector_angles = false; xtal_opt.max_time = 30.0; @@ -568,9 +570,6 @@ TEST_CASE("XtalOptimizer_rotation") { CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.2f); CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.2f); - // 700 micron is a lot, but for now it is important to see that refinement improves in the right direction - CHECK(fabsf(xtal_opt.geom.GetDetectorDistance_mm() - exp_i.GetDetectorDistance_mm()) < 0.7f); - // Lattice checks CHECK(fabsf(uc_ref.a - uc_out.a) < 0.2f); CHECK(fabsf(uc_ref.b - uc_out.b) < 0.2f); @@ -641,7 +640,6 @@ TEST_CASE("XtalOptimizer_refine_rotation_axis") { std::nullopt); xtal_opt.min_spots = 200; xtal_opt.refine_beam_center = true; - xtal_opt.refine_distance_mm = false; xtal_opt.refine_detector_angles = false; xtal_opt.refine_rotation_axis = true; xtal_opt.max_time = 30.0;