From 38bb511d233d17f303472e4070ee3e46b9649b7a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 23 Sep 2026 16:36:51 +0200 Subject: [PATCH] Make a first-pass detector tilt no mounting can have prove itself on the spots The first-pass rotation fit refines the spindle-perpendicular detector tilt freely, and on a sweep whose seed spots reach only a few degrees of 2theta the keystone that would determine it is a pixel or two. The fit then commits whatever the centroids' own systematics prefer - measured 2.5 deg on one sweep seeded to 11.7 A - and that tilt mispredicts the detector corners by 20-60 px against a 6 px integration disc, collapsing the run from 2.8 A to 6.35 A and biasing the metric-symmetry arbiter to a = b on the way. On synthetic data the same runaway is reproducible: a 1.5 deg detector error on 6 A reflections walks the free fit to 3.9 deg, on 8 A reflections to 37. Three data-driven discriminators were tried first and refuted on the corpus: a blanket restraint (0.13-0.16 A and up to 28% of ISa lost on the crystals whose tilt is largest and real, and in-house runs pinned at a header tilt known to be wrong), the distance's held-out excitation criterion (the rocking angles move by a tenth of their noise whether the tilt is real or not), and a keystone comparison of the fitted and header tilts with the beam free in both arms (236-set battery: ~30 sets moved, one collapsed from P1 to C2, one lost a screw axis; differences of 0.007-0.6 sigma held the header on right and wrong cases alike). For a tilt of a few tenths of a degree the keystone is a fraction of a pixel and nothing in the spots says whether it is real. What separated the cases was the tilt's size: every set the keystone arm moved carries 0.08-0.46 deg, and a survey of 211 corpus sets leaves the file's tilt by more than 0.56 deg on exactly one - a 2theta arm swung out 12.8 deg that its file records as square, which the free fit recovers and the held arm refuses by 90 sigma. So the hardware prior decides whether to ask, and the spots decide. Below 1 deg of walk from the tilt the pass started at the fit is trusted as before, by construction. Beyond it the winning candidate is refined again from its start with the tilt held there - the beam centre taking the shift the tilt is equivalent to, so the header tilt is never paired with a beam fitted beside a refused tilt - and the two are judged as the rounds of one chain already are, on the spots each indexes inside the wide gate: the walked tilt stands only when it leads by more than the count's own noise. A real tilt of degrees has a keystone of tens of pixels over the seed spots and wins outright; an artefact has none and loses on a tie. Held rather than bounded, because a box the fit lands on is the same wrong answer at a smaller size. Decided at the fit, so every pass and arm of a run handles itself and nothing is carried between passes; the unconstrained alternative is re-solved beside a refused tilt too. The chain that drives a candidate to its fixed point becomes a lambda so it can be run on the held candidate. The run logs the walk, what it would have moved the far corner by, the shift the beam took instead and both spot counts; where it was refused the report's REFINED_DETECTOR_TILT is the starting tilt and REFUSED_DETECTOR_TILT what the fit had walked to. Tests: a half-degree tilt is found, kept and not reported as a walk; a walk past the prior is made on synthetic data, and the result's verdict, counts and geometry agree with each other. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- docs/CPU_DATA_ANALYSIS_INDEXING.md | 2 + .../rotation_indexer/RotationIndexer.cpp | 104 +++++++++++++--- .../rotation_indexer/RotationIndexer.h | 13 ++ rugnux/ResultReport.cpp | 12 ++ rugnux/Rugnux.cpp | 36 ++++++ rugnux/Rugnux.h | 4 + tests/RotationIndexerTest.cpp | 115 ++++++++++++++++++ 7 files changed, 271 insertions(+), 15 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_INDEXING.md b/docs/CPU_DATA_ANALYSIS_INDEXING.md index 712c252aa..8f7e5cc7f 100644 --- a/docs/CPU_DATA_ANALYSIS_INDEXING.md +++ b/docs/CPU_DATA_ANALYSIS_INDEXING.md @@ -309,6 +309,8 @@ The loose first stage necessarily admits some spots that are not reflections of So every round is scored on the **widest** gate — the population the first pass selects on and the last pass does not fit, which makes it the one a converged solve is not optimising — and the best-scoring round is committed. Two conditions keep that from acting on noise. The score is a *count* of spots, so a lead of fewer than $\sqrt{\text{count}}$ of them leaves the last round standing. And the round taken has to be the less distorted lattice as well as the better-fitting one: the lattice search (§6) is re-asked every round to **measure** how far the cell sits from the ideal metric of the class it matches (imposing that class is measured fatal — the snap puts almost everything outside the refinement's own gate), and an earlier round is taken only when it matched the *same* class and sits closer to it. Same class is a precondition and not a precaution: the deviation is a fraction of whichever class's tolerance admitted it, so two classes' deviations are not the same quantity, and a round that matched no class reports zero, which means "nothing was asserted" rather than "undistorted". A chain that has settled scores its rounds within a spot or two of each other and a symmetry-constrained solve holds its distortion at zero throughout, so the rule fires on neither: measured over 914 chains, an earlier round scores higher on 44 % of them and the committed round changes on 1 dataset in 54. +**A tilt no mounting can have has to prove itself.** The detector tilt is refined freely because on a sweep whose spots reach far enough in $2\theta$ it is a measurement, and restraining it costs those crystals resolution (measured: 0.13–0.16 Å and up to a quarter of ISa on the crystals whose fitted tilt is largest). What makes it a measurement is the *keystone* — a tilted plane puts one side of the detector nearer and the other further, so the spots move by an amount that grows with their distance from the beam — and a first pass made of spots reaching a few degrees of $2\theta$ sees a keystone of a pixel or two at most. To such a fit a tilt is a whole-pattern shift the beam centre imitates exactly, its size is whatever the centroids' own systematics happen to prefer, and the value it commits then mispredicts the far corner of the detector by tens of pixels against an integration disc of a few: a first pass seeded to $2\theta = 5°$ committed 2.5° and the run collapsed from 2.8 Å to 6.4 Å. For a tilt of a few tenths of a degree nothing in the spots says whether it is real — the held-out positional residual, the rocking angles and a re-fit at the header tilt were all measured unable to, at the same insignificance on crystals whose tilt is real and on the one whose tilt was the artefact — so below what a mounting can be off square by the fit is trusted as before: a detector is mounted square to the beam to a fraction of a degree, and over 211 datasets the fitted tilt left the file's by more than 0.56° on one. A chain that has walked more than 1° from the tilt it started at has either measured nothing or found a detector the file misdescribes (that one: a $2\theta$ arm swung out 12.8° that the file records as square), and at that size the spots *do* tell the two apart, because a real tilt of degrees has a keystone of tens of pixels over the spots the fit is made of and an artefact has none. So the candidate is refined again from where it started with the tilt held there — the beam centre takes the shift the tilt is equivalent to — and the two are judged as the rounds of one chain are, on the spots each indexes inside the wide gate: the walked tilt stands only when it leads by more than the count's own noise. Held, not bounded — a box the fit lands on is the same wrong answer at a smaller size. The log says what the walked tilt would have moved the far corner by and how the two counts came out; where it was refused, the report's `REFINED_DETECTOR_TILT` is the tilt the pass started at and `REFUSED_DETECTOR_TILT` the tilt the fit had walked to. + ### 7.5 Rotation geometry post-refinement (two-pass) The refinement above (§7.2) runs per image against that image's spots. For rotation data an additional **post-refinement** (on by default; `--rotation-no-postrefine` disables it) improves the detector distance, beam centre and crystal cell/axis using **all** frames at once, then re-integrates: diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 568f21297..676a2d799 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -142,6 +142,7 @@ void RotationIndexer::RunIndexing() { } const auto indexer_result = indexer_.Run(experiment, coords); indexer_error_ = indexer_result.error; + tilt_walk_.reset(); if (!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) { DiffractionExperiment experiment_copy(experiment); @@ -460,41 +461,112 @@ void RotationIndexer::RunIndexing() { // towards a symmetry the data do not support. constexpr int ROT_REFINE_OUTER_ROUNDS = 20; constexpr double ROT_REFINE_TILT_SETTLED_RAD = 1.0e-5; // ~0.6 mdeg - if (have_best && !real_time) { - auto wide_count = [&](const XtalOptimizerData &d) { - const auto c = accumulate(d.geom, d.axis); - return IndexedFraction(d.latt, c, XTAL_OPTIMIZER_WIDE_TOLERANCE) * static_cast(c.size()); - }; + auto wide_count = [&](const XtalOptimizerData &d) { + const auto c = accumulate(d.geom, d.axis); + return IndexedFraction(d.latt, c, XTAL_OPTIMIZER_WIDE_TOLERANCE) * static_cast(c.size()); + }; + auto drive_to_fixed_point = [&](XtalOptimizerData best) { // True when a matched the same Bravais class as b and sits closer to its ideal metric. auto less_distorted = [](const LatticeSearchResult &a, const LatticeSearchResult &b) { return a.system == b.system && a.centering == b.centering && MetricDeviation(a) < MetricDeviation(b); }; - XtalOptimizerData kept = best_data; + XtalOptimizerData kept = best; float kept_count = -1.0f; float last_count = 0.0f; LatticeSearchResult kept_class, last_class; for (int r = 0; r < ROT_REFINE_OUTER_ROUNDS; ++r) { - XtalOptimizerData d = best_data; + XtalOptimizerData d = best; if (!XtalOptimizer(d, v_, kCeresThreads)) break; - const double d1 = std::abs(d.geom.GetPoniRot1_rad() - best_data.geom.GetPoniRot1_rad()); - const double d2 = std::abs(d.geom.GetPoniRot2_rad() - best_data.geom.GetPoniRot2_rad()); - best_data = std::move(d); - last_count = wide_count(best_data); - last_class = LatticeSearch(best_data.latt); + const double d1 = std::abs(d.geom.GetPoniRot1_rad() - best.geom.GetPoniRot1_rad()); + const double d2 = std::abs(d.geom.GetPoniRot2_rad() - best.geom.GetPoniRot2_rad()); + best = std::move(d); + last_count = wide_count(best); + last_class = LatticeSearch(best.latt); if (last_count >= kept_count) { kept_count = last_count; kept_class = last_class; - kept = best_data; + kept = best; } if (std::max(d1, d2) < ROT_REFINE_TILT_SETTLED_RAD) break; } if (kept_count > last_count + std::sqrt(last_count) && less_distorted(kept_class, last_class)) - best_data = std::move(kept); - best_frac = IndexedFraction(best_data.latt, accumulate(best_data.geom, best_data.axis), index_tol); + best = std::move(kept); + return best; + }; + if (have_best && !real_time) + best_data = drive_to_fixed_point(best_data); + + // The tilt is refined freely above because on a sweep whose spots reach far enough in 2theta + // it is a measurement, and restraining it costs those crystals resolution (measured: 0.13 to + // 0.16 A and up to a quarter of ISa on the crystals whose fitted tilt is largest). What makes + // it a measurement is the KEYSTONE - a tilted plane puts one side of the detector nearer and + // the other further, so the spots move by an amount that grows with their distance from the + // beam - and a first pass made of spots reaching a few degrees of 2theta sees a keystone of a + // pixel or two at most. To such a fit a tilt is a whole-pattern shift the beam centre + // imitates exactly, its size is whatever the centroids' own systematics happen to prefer, and + // the value it commits then mispredicts the far corner of the detector by tens of pixels + // against an integration disc of a few. Measured: a first pass seeded to 2theta 5 deg + // committed 2.5 deg and the run collapsed from 2.8 A to 6.4 A. + // + // Whether a fitted tilt is real cannot be read off the spots for a tilt of a few tenths of a + // degree: there the keystone is a fraction of a pixel, and the held-out residual, the rocking + // angles and a re-fit at the header tilt were all measured unable to tell the fitted tilt + // from the header's - at the same insignificance on crystals whose tilt is real and on the + // one whose tilt was the artefact. So below what a mounting can be off square by the fit is + // trusted, as it always was: a detector is mounted square to the beam to a fraction of a + // degree, and over 211 corpus datasets the fitted tilt leaves the file's by more than 0.56 deg + // on one. A fit that has walked further than that has either measured nothing or found a + // detector the file misdescribes (the one: a 2theta arm swung out 12.8 deg that the file + // records as square) - and at that size the two ARE told apart by the spots, because a real + // tilt of degrees has a keystone of tens of pixels over the spots the fit is made of and an + // artefact has none. So the candidate is refined again from where it started with the tilt + // held there, the beam centre taking the shift the tilt stood in for, and the two are judged + // as the rounds of one chain are: on the spots each indexes inside the wide gate, the walked + // tilt standing only when it leads by more than the count's own noise. Held, not bounded - a + // box the fit lands on is the same wrong answer at a smaller size. + constexpr double ROT_TILT_PRIOR_DEG = 1.0; + if (have_best && !real_time) { + const double walk_deg = std::hypot(best_data.geom.GetPoniRot1_rad() - geom_.GetPoniRot1_rad(), + best_data.geom.GetPoniRot2_rad() - geom_.GetPoniRot2_rad()) + * 180.0 / PI; + if (walk_deg > ROT_TILT_PRIOR_DEG) { + const bool tri_won = work[best_ci].has_tri + && best_sr.system == gemmi::CrystalSystem::Triclinic; + XtalOptimizerData held = tri_won ? work[best_ci].tri : work[best_ci].constrained; + held.refine_detector_angles = false; + if (XtalOptimizer(held, v_, kCeresThreads)) { + held = drive_to_fixed_point(std::move(held)); + const float spots_walked = wide_count(best_data); + const float spots_held = wide_count(held); + const bool refused = !(spots_walked > spots_held + std::sqrt(spots_held)); + tilt_walk_ = RotationIndexerResult::TiltWalk{ + .tilt_rad = {best_data.geom.GetPoniRot1_rad(), best_data.geom.GetPoniRot2_rad()}, + .spots_walked = spots_walked, + .spots_held = spots_held, + .refused = refused}; + if (refused) { + best_data = std::move(held); + // The unconstrained alternative was refined beside the refused tilt too. + if (best_alt) { + XtalOptimizerData alt = work[best_ci].tri; + alt.refine_detector_angles = false; + if (XtalOptimizer(alt, v_, kCeresThreads)) + best_alt = std::make_shared(RotationIndexerResult{ + .lattice = alt.latt, + .search_result = best_alt->search_result, + .geom = alt.geom, + .axis = alt.axis, + }); + } + } + } + } } + if (have_best && !real_time) + best_frac = IndexedFraction(best_data.latt, accumulate(best_data.geom, best_data.axis), index_tol); if (have_best) { search_result_ = best_sr; @@ -656,6 +728,7 @@ std::optional RotationIndexer::GetLattice() const { .geom = updated_geom_, .axis = axis_, .unconstrained = unconstrained_, + .tilt_walk = tilt_walk_, }; } @@ -672,6 +745,7 @@ void RotationIndexer::ForceResult(const RotationIndexerResult &result) { updated_geom_ = result.geom; axis_ = result.axis; unconstrained_ = result.unconstrained; + tilt_walk_ = result.tilt_walk; } bool RotationIndexer::AccumulationFull() const { diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index de7f6c264..29dac0ab1 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -24,6 +24,18 @@ struct RotationIndexerResult { // constraint then snaps a real angle to the ideal one; here the caller can settle that on a // statistic the accumulated-spot fraction is too blunt for. Null on the alternative itself. std::shared_ptr unconstrained; + // Set where the free fit walked the detector tilt further from the tilt the indexer was handed + // than a mounted detector is off square by (ROT_TILT_PRIOR_DEG in RotationIndexer.cpp), and the + // lattice was therefore refined again with the tilt held where it started and the two compared + // on the spots they index: the tilt the fit walked to, the two counts, and which won. Refused, + // geom carries the tilt the indexer started at; confirmed, the walked tilt. + struct TiltWalk { + std::array tilt_rad; + float spots_walked = 0.0f; + float spots_held = 0.0f; + bool refused = false; + }; + std::optional tilt_walk; }; class RotationIndexer { @@ -52,6 +64,7 @@ class RotationIndexer { LatticeSearchResult search_result_; std::vector extra_lattices_; std::shared_ptr unconstrained_; + std::optional tilt_walk_; IndexerThreadPool &indexer_; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index e760ecec3..7e787b1fe 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -363,6 +363,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyText("REFINED_DETECTOR_TILT", fmt::format("{:.4f} {:.4f}", (*result.refined_detector_tilt_deg)[0], (*result.refined_detector_tilt_deg)[1]))); + if (result.refused_detector_tilt_deg) + Add(s, KeyText("REFUSED_DETECTOR_TILT", + fmt::format("{:.4f} {:.4f}", (*result.refused_detector_tilt_deg)[0], + (*result.refused_detector_tilt_deg)[1]))); // What the beam flew through, and what assuming it was worth. NOTHING in any file rugnux // reads states the medium, so this is an ASSUMPTION the run made on the user's behalf and @@ -461,6 +465,14 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " does not replace the calibration.\n" " rot3 is omitted because a rotation about the beam is an exact null of this experiment\n" " and the fit cannot move it.", true)); + if (result.refused_detector_tilt_deg) + Add(s, Prose(" REFUSED_DETECTOR_TILT is where this pass's free fit of the tilt had walked to - more than\n" + " a degree from the tilt the pass started at, further than a mounted detector is off\n" + " square by. The lattice was refined again with the tilt held where the pass started,\n" + " and the spots indexed no fewer at that tilt than at the walked one, so the walk measured\n" + " nothing: REFINED_DETECTOR_TILT above is that starting tilt, not a fit. It happens on a\n" + " sweep whose spots reach only a few degrees of 2theta, where a tilt is a whole-pattern\n" + " shift the beam centre imitates exactly and nothing sets its size.", true)); if (result.pass_count > 1) Add(s, Prose(fmt::format( " A rotation run integrates twice: once at the geometry in the input file, then again at\n" diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 54cb0c99f..c4ce5ed0e 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -4638,6 +4638,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b "{}/{} validation spots = {:.1f}% against {:.1f}% at a wrong spindle angle)", best.name, best.score, static_cast(validation.size()), evidence.on_lattice, evidence.spots, 100.0 * pooled, 100.0 * pooled_chance); + if (best.result->tilt_walk) { + // What the walked tilt does to the prediction: the keystone at the far corner of the + // detector, against the disc the integration sums over, and the whole-pattern shift + // it is equivalent to, which the beam centre carries instead where it is refused. + const auto &tw = *best.result->tilt_walk; + const auto &g = best.result->geom; + const double start_rot1 = tw.refused ? g.GetPoniRot1_rad() : experiment_.GetPoniRot1_rad(); + const double start_rot2 = tw.refused ? g.GetPoniRot2_rad() : experiment_.GetPoniRot2_rad(); + const double walk_rad = std::hypot(tw.tilt_rad[0] - start_rot1, tw.tilt_rad[1] - start_rot2); + const double corner_px = std::hypot( + std::max(g.GetBeamX_pxl(), experiment_.GetXPixelsNumConv() - g.GetBeamX_pxl()), + std::max(g.GetBeamY_pxl(), experiment_.GetYPixelsNumConv() - g.GetBeamY_pxl())); + const double lever_px = g.GetDetectorDistance_mm() / g.GetPixelSize_mm(); + const std::string what = fmt::format( + "Detector tilt: the free fit walked to {:.4f},{:.4f} deg, {:.2f} deg from the " + "{:.4f},{:.4f} deg this pass started at - further than a mounted detector is off " + "square by - which moves the far corner of the detector by {:.0f} px against an " + "integration disc of {:.1f} px, so the lattice was refined again at the starting " + "tilt with the beam centre free to take the {:.0f} px shift the tilt is equivalent " + "to, and the two were judged on the spots they index: {:.0f} at the walked tilt " + "against {:.0f} at the starting one", + tw.tilt_rad[0] * 180.0 / PI, tw.tilt_rad[1] * 180.0 / PI, walk_rad * 180.0 / PI, + start_rot1 * 180.0 / PI, start_rot2 * 180.0 / PI, + walk_rad * corner_px * corner_px / lever_px, + experiment_.GetBraggIntegrationSettings().GetR1(), walk_rad * lever_px, + tw.spots_walked, tw.spots_held); + if (tw.refused) + logger.Warning("{} - the spots cannot tell the two apart, the walked tilt measured " + "nothing, and the pass runs at the tilt it started at", what); + else + logger.Warning("{} - the spots confirm the walked tilt, and the pass runs at it; the " + "file misdescribes this detector's tilt", what); + } // The spots in hand were already cut to the budget in force, so the measurement can only // shorten a budget, never lengthen one - which is why RunAllPasses hands the second pass the @@ -5057,6 +5090,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.refined_detector_tilt_deg = std::array{ rot->geom.GetPoniRot1_rad() * 180.0 / PI, rot->geom.GetPoniRot2_rad() * 180.0 / PI}; + if (rot->tilt_walk && rot->tilt_walk->refused) + result.refused_detector_tilt_deg = std::array{ + rot->tilt_walk->tilt_rad[0] * 180.0 / PI, rot->tilt_walk->tilt_rad[1] * 180.0 / PI}; if (rot->axis) end_msg.refined_rotation_axis = rot->axis->GetAxis(); end_msg.rotation_lattice_type = LatticeMessage{ diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 8c43eb14e..1ed2bee04 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -216,6 +216,10 @@ struct ProcessResult { // can be checked against a powder calibration or handed to another program; empty on stills and // when no rotation lattice was finalized. std::optional> refined_detector_tilt_deg; + // The tilt the free fit walked to where rotation indexing REFUSED it - further from the tilt the + // pass started at than a mounted detector can be off square by (see RotationIndexer) - so that + // refined_detector_tilt_deg above is the tilt the pass started at, not a fit. Empty otherwise. + std::optional> refused_detector_tilt_deg; // The tilt this pass integrated at, and where the beam actually lands under it - taken from the // same experiment_ as used_beam_*/used_distance_mm, so the four describe ONE geometry. Reading any // of them off the caller's pre-run copy instead would mix a post-refined centre with the file's diff --git a/tests/RotationIndexerTest.cpp b/tests/RotationIndexerTest.cpp index 7fdf2b2a5..73cfa44e4 100644 --- a/tests/RotationIndexerTest.cpp +++ b/tests/RotationIndexerTest.cpp @@ -180,3 +180,118 @@ TEST_CASE("RotationIndexer::RefineConstrained puts a free metric back on its cla CHECK(refit->indexed_fraction > refit->indexed_fraction_before); CHECK(refit->indexed_fraction > 0.5f); } + +// Index a synthetic sweep recorded on a detector tilted by true_tilt_deg beyond the tilt the indexer +// is handed, with reflections to res_A. +static std::optional IndexOnTiltedDetector(double true_tilt_deg, float res_A) { + constexpr double header_rot2_rad = 0.02; + DiffractionExperiment exp_header; + exp_header.IncidentEnergy_keV(WVL_1A_IN_KEV) + .BeamX_pxl(1000) + .BeamY_pxl(1000) + .PoniRot1_rad(0.01) + .PoniRot2_rad(header_rot2_rad) + .DetectorDistance_mm(200) + .ImagesPerTrigger(50); + + IndexingSettings settings; +#ifdef JFJOCH_USE_CUDA + settings.Algorithm(IndexingAlgorithmEnum::FFT); +#elif JFJOCH_USE_FFTW + settings.Algorithm(IndexingAlgorithmEnum::FFTW); +#else + return {}; +#endif + settings.RotationIndexing(true).RotationIndexingAngularStride_deg(1.0).RotationIndexingMinAngularRange_deg(30.0); + exp_header.ImportIndexingSettings(settings); + + GoniometerAxis axis("omega", 0.0f, 1.0f, Coord(1, 0, 0), std::nullopt); + exp_header.Goniometer(axis); + + // The detector the spots were actually recorded on. + DiffractionExperiment exp_true = exp_header; + exp_true.PoniRot2_rad(header_rot2_rad + true_tilt_deg * PI / 180.0); + + const CrystalLattice latt_base = + CrystalLattice(40, 50, 80, 90, 90, 90).Multiply(RotMatrix(2.0, Coord(sqrt(3)/3, sqrt(3)/3, sqrt(3)/3))); + + BraggPredictionSettings prediction_settings{ .high_res_A = res_A, .ewald_dist_cutoff = 0.002 }; + IndexerThreadPool indexer_thread_pool(exp_header.GetIndexingSettings()); + RotationIndexer indexer(exp_header, indexer_thread_pool); + BraggPrediction prediction; + + for (int img = 0; img < 50; ++img) { + std::vector spots; + const float angle_deg = axis.GetAngle_deg(img) + axis.GetWedge_deg() / 2.0f; + const CrystalLattice latt_img = latt_base.Multiply(axis.GetTransformationAngle(angle_deg).transpose()); + const auto n = prediction.Calc(exp_true, 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 = img; + s.intensity = 1.0f; + s.phi = angle_deg; + s.ice_ring = false; + s.indexed = true; + spots.push_back(s); + } + indexer.ProcessImage(img, spots); + if (img == 30) + indexer.RunIndexing(); + } + // Round-trip through ForceResult - how a canonical pass takes over the result of the scheme + // indexer that found the lattice, and what the report then reads - so what comes back is what a + // run sees, the tilt walk included. + const auto found = indexer.GetLattice(); + if (!found) + return {}; + RotationIndexer forced(exp_header, indexer_thread_pool); + forced.ForceResult(*found); + return forced.GetLattice(); +} + +// The detector tilt is refined freely, and a fit that walks it further from where it started than a +// mounted detector can be off square by (ROT_TILT_PRIOR_DEG) is refused and made again with the tilt +// held. The prior must not touch a tilt a mounting can have: on a detector tilted half a degree +// beyond the tilt the indexer is handed, with reflections to 2.5 A, the fit finds it, keeps it and +// reports no refusal. +TEST_CASE("RotationIndexer keeps a tilt a mounting can have") { + const auto ret = IndexOnTiltedDetector(0.5, 2.5f); + REQUIRE(ret.has_value()); + CHECK_FALSE(ret->tilt_walk.has_value()); + CHECK((ret->geom.GetPoniRot2_rad() - 0.02) * 180.0 / PI == Catch::Approx(0.5).margin(0.05)); + CHECK(ret->geom.GetPoniRot1_rad() == Catch::Approx(0.01).margin(1e-3)); + const auto uc = ret->lattice.GetUnitCell(); + CHECK(uc.a == Catch::Approx(40.0).margin(0.3)); + CHECK(uc.b == Catch::Approx(50.0).margin(0.3)); + CHECK(uc.c == Catch::Approx(80.0).margin(0.5)); +} + +// The other side of the prior, and the unidentifiability that makes it necessary: the same detector +// tilted 1.5 deg beyond the handed tilt, but with reflections only to 6 A, where the keystone the +// fit could read the tilt off is a fraction of a pixel. The free fit does not find 1.5 deg - it runs +// away to about 4 deg (measured 3.9; at 8 A it reaches 37), because at that 2theta reach the tilt is +// a whole-pattern shift the beam centre imitates and nothing pins its size. That walk is past the +// prior, so the lattice is refined again with the tilt held and the two are judged on the spots they +// index; the result records the walk, both counts, the verdict, and a geometry that matches it. The +// lattice is not checked: a 1.5 deg detector error on 6 A data already puts the FFT on a different +// cell before any fit. +TEST_CASE("RotationIndexer judges a tilt no mounting can have on the spots") { + const auto ret = IndexOnTiltedDetector(1.5, 6.0f); + REQUIRE(ret.has_value()); + REQUIRE(ret->tilt_walk.has_value()); + const auto &tw = *ret->tilt_walk; + CHECK(std::hypot(tw.tilt_rad[0] - 0.01, tw.tilt_rad[1] - 0.02) * 180.0 / PI > 1.0); + CHECK(tw.spots_walked > 0.0f); + CHECK(tw.spots_held > 0.0f); + CHECK(tw.refused == !(tw.spots_walked > tw.spots_held + std::sqrt(tw.spots_held))); + if (tw.refused) { + CHECK(ret->geom.GetPoniRot1_rad() == Catch::Approx(0.01).margin(1e-7)); + CHECK(ret->geom.GetPoniRot2_rad() == Catch::Approx(0.02).margin(1e-7)); + } else { + CHECK(ret->geom.GetPoniRot1_rad() == Catch::Approx(tw.tilt_rad[0]).margin(1e-7)); + CHECK(ret->geom.GetPoniRot2_rad() == Catch::Approx(tw.tilt_rad[1]).margin(1e-7)); + } +}