diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp
index 7e952786..2af4a9d6 100644
--- a/image_analysis/geom_refinement/PostRefine.cpp
+++ b/image_analysis/geom_refinement/PostRefine.cpp
@@ -8,6 +8,7 @@
 #include <array>
 #include <cmath>
 #include <limits>
+#include <memory>
 #include <numeric>
 
 #include "../../common/JFJochMath.h"   // PI
@@ -18,13 +19,15 @@
 
 namespace {
 
-// One integrated partial, flattened across all images.
+// One integrated partial, flattened across all images. Kept as narrow as the sort and the event
+// split allow: on a large cell this array is gigabytes, and the scatter and every level of the
+// per-bucket sort move all of it. The goniometer angle is not stored - it is a function of the
+// image number alone, and is rebuilt from it where it is needed.
 struct Partial {
     int h, k, l;
     float img;
-    double I, sigma;
-    double angle_rad;   // frame mid-exposure goniometer angle
-    double obs_x, obs_y; // observed spot centroid (pixels); NAN if the box sum found no centroid
+    float I, sigma;
+    float obs_x, obs_y;  // observed spot centroid (pixels); NAN if the box sum found no centroid
 };
 
 // A rocking event and its precomputed reference reciprocal vector (phi=0 frame, from the indexed lattice).
@@ -108,46 +111,82 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
         const Coord ax = axis.GetAxis();
         const Coord Astar = reference_latt.Astar(), Bstar = reference_latt.Bstar(), Cstar = reference_latt.Cstar();
 
+        // The goniometer angle is a function of the image number alone, so it is recomputed where it
+        // is used rather than carried through the array below: eight bytes per partial cost more in
+        // the fill, the scatter and every level of the sort than the multiply-add that rebuilds them.
+        const auto angle_rad = [&](float img) { return (axis.GetAngle_deg(img) + wedge_half) * PI / 180.0; };
+
         // Count first, then fill. Growing one vector by push_back over tens of millions of
         // reflections copies the whole thing every time it doubles - several gigabytes of pure
         // copying - and the counts are cheap to take. Each outcome then owns a slice, so the fill
-        // runs on all threads and lands in the order the serial loop produced.
+        // runs on all threads and lands in the order the serial loop produced. The h range comes out
+        // of the same sweep: the sort below buckets by h and needs to know how many buckets that is,
+        // and this pass already reads every reflection.
         const size_t nthreads = std::max(1, settings.num_threads);
         const int n_out = static_cast<int>(outcomes.size());
         std::vector<size_t> pts_offset(n_out + 1, 0);
+        std::vector<int> h_lo_of(n_out), h_hi_of(n_out);
         ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
             for (int o = lo; o < hi; o++) {
                 size_t keep = 0;
+                int lmin = std::numeric_limits<int>::max(), lmax = std::numeric_limits<int>::min();
                 for (const auto &r : outcomes[o].reflections)
-                    if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f)
+                    if (std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0.0f) {
                         keep++;
+                        lmin = std::min(lmin, r.h);
+                        lmax = std::max(lmax, r.h);
+                    }
                 pts_offset[o + 1] = keep;
+                h_lo_of[o] = lmin;
+                h_hi_of[o] = lmax;
             }
         });
-        for (int o = 0; o < n_out; o++)
+        int h_lo = std::numeric_limits<int>::max(), h_hi = std::numeric_limits<int>::min();
+        for (int o = 0; o < n_out; o++) {
             pts_offset[o + 1] += pts_offset[o];
+            h_lo = std::min(h_lo, h_lo_of[o]);
+            h_hi = std::max(h_hi, h_hi_of[o]);
+        }
+        const size_t n_pts = pts_offset[n_out];
+        const int H = (h_lo <= h_hi) ? (h_hi - h_lo + 1) : 1;
 
-        std::vector<Partial> pts(pts_offset[n_out]);
+        // A vector of n partials VALUE-initialises them: on a large cell that is gigabytes of zeroing
+        // on one thread, and it is that one thread which first touches every page - which on a
+        // multi-socket machine leaves the whole array on its node, so every pass that follows runs at
+        // one node's memory bandwidth. new[] leaves the partials untouched, so the parallel fill is
+        // the first touch and each page lands on the node of the thread that filled it.
+        std::unique_ptr<Partial[]> pts(new Partial[n_pts]);
+        // The bucket histogram the sort needs is taken here rather than in a pass of its own, since
+        // the fill already has h in hand. Its chunks are the outcome chunks ParallelChunks makes, so
+        // the scatter below has to be split the same way.
+        const int nt = static_cast<int>(std::clamp<size_t>(nthreads, 1, std::max(1, n_out)));
+        const int chunk = (n_out + nt - 1) / nt;
+        std::vector<std::vector<int32_t>> hist(nt, std::vector<int32_t>(H, 0));
         ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
+            std::vector<int32_t> &h_count = hist[lo / chunk];
             for (int o = lo; o < hi; o++) {
                 size_t at = pts_offset[o];
                 for (const auto &r : outcomes[o].reflections) {
                     if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) continue;
-                    const double mid_deg = axis.GetAngle_deg(r.image_number) + wedge_half;
-                    const double ox = std::isfinite(r.observed_x) ? r.observed_x : NAN;
-                    const double oy = std::isfinite(r.observed_y) ? r.observed_y : NAN;
-                    pts[at++] = Partial{r.h, r.k, r.l, r.image_number, r.I, r.sigma,
-                                        mid_deg * PI / 180.0, ox, oy};
+                    const float ox = std::isfinite(r.observed_x) ? r.observed_x : NAN;
+                    const float oy = std::isfinite(r.observed_y) ? r.observed_y : NAN;
+                    pts[at++] = Partial{r.h, r.k, r.l, r.image_number, r.I, r.sigma, ox, oy};
+                    h_count[r.h - h_lo]++;
                 }
             }
         });
-        logger.Info("Post-refine: {} partials gathered", pts.size());
-        if (pts.size() < static_cast<size_t>(settings.min_events)) return result;
+        logger.Info("Post-refine: {} partials gathered", n_pts);
+        if (n_pts < static_cast<size_t>(settings.min_events)) return result;
+
+        // Where each bucket starts, and the buckets largest first: the sort lays the array out this
+        // way and the event split below walks the same buckets.
+        std::vector<int32_t> bstart(H + 1, 0);
+        std::vector<int> order(H);
 
         // Bucket by h, then sort the buckets. h is the leading key, so the sorted array is the
         // buckets laid end to end, and each bucket sorts on its own thread. Sorting the whole thing
-        // in one pass moved 56 bytes per element through every level of a comparison sort, on one
-        // thread, over tens of millions of reflections.
+        // in one pass moved every partial through every level of a comparison sort, on one thread,
+        // over tens of millions of reflections.
         {
             // Not a total order: two partials of one reflection on one image still tie, as they did
             // before this was bucketed. What makes the result reproducible is the scatter below
@@ -161,21 +200,6 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
                 if (a.l != b.l) return a.l < b.l;
                 return a.img < b.img;
             };
-            const int n = static_cast<int>(pts.size());
-            int h_lo = std::numeric_limits<int>::max(), h_hi = std::numeric_limits<int>::min();
-            for (const auto &q : pts) { h_lo = std::min(h_lo, q.h); h_hi = std::max(h_hi, q.h); }
-            const int H = (h_lo <= h_hi) ? (h_hi - h_lo + 1) : 1;
-            const int nt = static_cast<int>(std::clamp<size_t>(nthreads, 1, std::max(1, n)));
-            const int chunk = (n + nt - 1) / nt;
-
-            std::vector<std::vector<int32_t>> hist(nt, std::vector<int32_t>(H, 0));
-            ParallelChunks(nt, nthreads, [&](int tlo, int thi) {
-                for (int t = tlo; t < thi; ++t) {
-                    const int lo = t * chunk, hi = std::min(n, lo + chunk);
-                    for (int i = lo; i < hi; ++i) hist[t][pts[i].h - h_lo]++;
-                }
-            });
-            std::vector<int32_t> bstart(H + 1, 0);
             int32_t acc = 0;
             for (int b = 0; b < H; ++b) {
                 bstart[b] = acc;
@@ -183,21 +207,18 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             }
             bstart[H] = acc;
 
-            std::vector<Partial> sorted(pts.size());
-            ParallelChunks(nt, nthreads, [&](int tlo, int thi) {
-                for (int t = tlo; t < thi; ++t) {
-                    std::vector<int32_t> fill = hist[t];
-                    const int lo = t * chunk, hi = std::min(n, lo + chunk);
-                    for (int i = lo; i < hi; ++i) sorted[fill[pts[i].h - h_lo]++] = pts[i];
-                }
+            std::unique_ptr<Partial[]> sorted(new Partial[n_pts]);
+            ParallelChunks(n_out, nthreads, [&](int lo, int hi) {
+                std::vector<int32_t> fill = hist[lo / chunk];
+                for (size_t i = pts_offset[lo]; i < pts_offset[hi]; ++i)
+                    sorted[fill[pts[i].h - h_lo]++] = pts[i];
             });
-            std::vector<int> order(H);
             std::iota(order.begin(), order.end(), 0);
             std::sort(order.begin(), order.end(),
                       [&](int a, int b) { return (bstart[a + 1] - bstart[a]) > (bstart[b + 1] - bstart[b]); });
             ParallelFor(H, nthreads, [&](int oi) {
                 const int b = order[oi];
-                std::sort(sorted.begin() + bstart[b], sorted.begin() + bstart[b + 1], part_less);
+                std::sort(sorted.get() + bstart[b], sorted.get() + bstart[b + 1], part_less);
             });
             pts.swap(sorted);
         }
@@ -205,50 +226,88 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
         // Split into rocking events (same raw hkl, adjacent frames). Only >=2-frame events carry an
         // unbiased phi_obs (a single-frame centroid is just the frame centre); precompute e_ref per event.
         constexpr float MAX_FRAME_GAP = 2.0f;
-        std::vector<Event> events;
-        size_t i = 0, event_frames = 0;
-        while (i < pts.size()) {
+        const auto run_end = [&](size_t i, size_t end) {
             size_t j = i + 1;
-            while (j < pts.size() && pts[j].h == pts[i].h && pts[j].k == pts[i].k && pts[j].l == pts[i].l
+            while (j < end && pts[j].h == pts[i].h && pts[j].k == pts[i].k && pts[j].l == pts[i].l
                    && pts[j].img - pts[j - 1].img <= MAX_FRAME_GAP)
                 ++j;
-            if (j - i >= 2) {
-                double sumI = 0, sumIphi = 0, sumSig = 0;
-                for (size_t m = i; m < j; ++m) {
-                    const double Ipos = std::max(0.0, pts[m].I);
-                    sumI += Ipos; sumIphi += Ipos * pts[m].angle_rad; sumSig += pts[m].sigma;
-                }
-                if (sumI > 0.0 && sumSig > 0.0) {
-                    const double phi = sumIphi / sumI;
-                    event_frames += j - i;
-                    const Coord e = Astar * static_cast<float>(pts[i].h) + Bstar * static_cast<float>(pts[i].k)
-                                    + Cstar * static_cast<float>(pts[i].l);
-                    events.push_back(Event{phi, std::sqrt(sumI / sumSig), {e.x, e.y, e.z},
-                                           pts[i].h, pts[i].k, pts[i].l});
-                }
+            return j;
+        };
+        // The event the partials [i, j) make, or false where their intensities cannot place a
+        // centroid. Counting the events and writing them both walk the buckets, and both build the
+        // event this way.
+        const auto make_event = [&](size_t i, size_t j, Event &out) {
+            double sumI = 0, sumIphi = 0, sumSig = 0;
+            for (size_t m = i; m < j; ++m) {
+                const double Ipos = std::max(0.0, static_cast<double>(pts[m].I));
+                sumI += Ipos; sumIphi += Ipos * angle_rad(pts[m].img); sumSig += pts[m].sigma;
             }
-            i = j;
-        }
+            if (!(sumI > 0.0 && sumSig > 0.0)) return false;
+            const Coord e = Astar * static_cast<float>(pts[i].h) + Bstar * static_cast<float>(pts[i].k)
+                            + Cstar * static_cast<float>(pts[i].l);
+            out = Event{sumIphi / sumI, std::sqrt(sumI / sumSig), {e.x, e.y, e.z},
+                        pts[i].h, pts[i].k, pts[i].l};
+            return true;
+        };
+        // An event never crosses an h boundary - h is the leading sort key - so the buckets can be
+        // walked independently, and laying their events out in bucket order gives exactly the order
+        // the serial walk produced. Counting first also sizes the array in one go, in place of a
+        // push_back that grew a gigabyte by doubling.
+        std::vector<int32_t> ev_count(H, 0);
+        std::vector<size_t> ev_frames(H, 0);
+        ParallelFor(H, nthreads, [&](int oi) {
+            const int b = order[oi];
+            const size_t end = bstart[b + 1];
+            int c = 0;
+            Event ev;
+            for (size_t i = bstart[b]; i < end; ) {
+                const size_t j = run_end(i, end);
+                if (j - i >= 2 && make_event(i, j, ev)) ++c;
+                i = j;
+            }
+            ev_count[b] = c;
+        });
+        std::vector<int32_t> ev_start(H + 1, 0);
+        for (int b = 0; b < H; ++b) ev_start[b + 1] = ev_start[b] + ev_count[b];
+        const size_t n_events = ev_start[H];
+        std::unique_ptr<Event[]> events(new Event[n_events]);
+        ParallelFor(H, nthreads, [&](int oi) {
+            const int b = order[oi];
+            const size_t end = bstart[b + 1];
+            int at = ev_start[b];
+            size_t frames = 0;
+            for (size_t i = bstart[b]; i < end; ) {
+                const size_t j = run_end(i, end);
+                if (j - i >= 2 && make_event(i, j, events[at])) { frames += j - i; ++at; }
+                i = j;
+            }
+            ev_frames[b] = frames;
+        });
+        size_t event_frames = 0;
+        for (int b = 0; b < H; ++b) event_frames += ev_frames[b];
         // Frames per event is the phi_obs sampling: near 2 the reflections barely rock, so the angle
         // this refinement is fitted to is under-determined. It is a geometry count, so unlike an
         // intensity-weighted width it cannot be inflated by noise.
-        logger.Info("Post-refine: {} multi-frame rocking events ({:.1f} frames per event)", events.size(),
-                    events.empty() ? 0.0 : static_cast<double>(event_frames) / events.size());
-        if (static_cast<int>(events.size()) < settings.min_events) return result;
+        logger.Info("Post-refine: {} multi-frame rocking events ({:.1f} frames per event)", n_events,
+                    n_events == 0 ? 0.0 : static_cast<double>(event_frames) / n_events);
+        if (static_cast<int>(n_events) < 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<Event> scale_events = events;
-
+        // has to see them. Both sets come out of the one array by selecting on indices instead: with
+        // the weights in the same places nth_element takes the same decisions it would take on the
+        // events themselves, so the selection is the same one in the same order and the whole list no
+        // longer has to be duplicated to survive it.
         constexpr size_t MAX_EVENTS = 20000;
-        if (events.size() > MAX_EVENTS) {
-            std::nth_element(events.begin(), events.begin() + MAX_EVENTS, events.end(),
-                             [](const Event &a, const Event &b) { return a.weight > b.weight; });
-            events.resize(MAX_EVENTS);
+        std::vector<int32_t> selected(n_events);
+        std::iota(selected.begin(), selected.end(), 0);
+        if (selected.size() > MAX_EVENTS) {
+            std::nth_element(selected.begin(), selected.begin() + MAX_EVENTS, selected.end(),
+                             [&](int32_t a, int32_t b) { return events[a].weight > events[b].weight; });
+            selected.resize(MAX_EVENTS);
         }
 
         // ---- GEOMETRY REFINEMENT: the XtalOptimizer-equivalent, done as TWO SEPARATE
@@ -288,7 +347,8 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             // ===== Step A: cell scale s + rotation axis from phi_obs =====
             auto excit_cost = [&](Subset s, double sc, const double axv[3]) {
                 double c = 0.0; int n = 0;
-                for (const auto &ev : events) {
+                for (const int32_t i : selected) {
+                    const Event &ev = events[i];
                     if (!in(ev.h, ev.k, ev.l, s)) continue;
                     ScaleAxisExcitationResidual r(lambda_l, ev.phi_obs, 1.0, ev.e_ref);
                     double sd = sc, av[3] = {axv[0], axv[1], axv[2]}, resid = 0.0;
@@ -299,7 +359,8 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             auto solve_scale_axis = [&](Subset s, double &s_out, double ax_out[3]) {
                 double sc = 1.0, axv[3] = {ax0[0], ax0[1], ax0[2]};
                 ceres::Problem p;
-                for (const auto &ev : events) {
+                for (const int32_t i : selected) {
+                    const Event &ev = events[i];
                     if (!in(ev.h, ev.k, ev.l, s)) continue;
                     p.AddResidualBlock(new ceres::AutoDiffCostFunction<ScaleAxisExcitationResidual, 1, 1, 3>(
                         new ScaleAxisExcitationResidual(lambda_l, ev.phi_obs, settings.excitation_weight, ev.e_ref)),
@@ -338,13 +399,13 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             // 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);
+            double phi_c = 0.0, phi_lo = events[0].phi_obs, phi_hi = events[0].phi_obs;
+            for (size_t e = 0; e < n_events; ++e) {
+                phi_c += events[e].phi_obs;
+                phi_lo = std::min(phi_lo, events[e].phi_obs);
+                phi_hi = std::max(phi_hi, events[e].phi_obs);
             }
-            phi_c /= static_cast<double>(scale_events.size());
+            phi_c /= static_cast<double>(n_events);
             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
@@ -366,13 +427,13 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             // Coefficients are computed in double and stored narrowed: their rounding perturbs the
             // minimiser by ~1e-10, and k is carried downstream as a float.
             struct ScaleTerm { float a, C, R, psi; };
-            std::vector<ScaleTerm> terms(scale_events.size());
-            std::vector<int> fifth_of(scale_events.size());
+            std::vector<ScaleTerm> terms(n_events);
+            std::vector<int> fifth_of(n_events);
             const double aa_c[3] = {-phi_c * u[0], -phi_c * u[1], -phi_c * u[2]};
-            ParallelChunks(static_cast<int>(scale_events.size()), nthreads, [&](int lo, int hi) {
+            ParallelChunks(static_cast<int>(n_events), nthreads, [&](int lo, int hi) {
                 for (int e = lo; e < hi; ++e) {
-                    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};
+                    const double p[3] = {events[e].e_ref[0] / s, events[e].e_ref[1] / s,
+                                         events[e].e_ref[2] / s};
                     double em[3];
                     ceres::AngleAxisRotatePoint(aa_c, p, em);
                     const double ue = u[0] * em[0] + u[1] * em[1] + u[2] * em[2];
@@ -380,11 +441,11 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
                     const double C = 0.5 * lambda_l * e2 + u[2] * ue;
                     const double A = em[2] - u[2] * ue;
                     const double B = u[0] * em[1] - u[1] * em[0];
-                    terms[e] = ScaleTerm{static_cast<float>(scale_events[e].phi_obs - phi_c),
+                    terms[e] = ScaleTerm{static_cast<float>(events[e].phi_obs - phi_c),
                                          static_cast<float>(C), static_cast<float>(std::hypot(A, B)),
                                          static_cast<float>(std::atan2(B, A))};
                     fifth_of[e] = std::clamp(static_cast<int>(
-                            5.0 * (scale_events[e].phi_obs - phi_lo) / std::max(1e-9, phi_hi - phi_lo)), 0, 4);
+                            5.0 * (events[e].phi_obs - phi_lo) / std::max(1e-9, phi_hi - phi_lo)), 0, 4);
                 }
             });
 
@@ -492,7 +553,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             // 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<int>(scale_events.size()) >= MIN_SCALE_EVENTS
+            const bool enough_data = static_cast<int>(n_events) >= 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;
@@ -503,7 +564,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             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,
+                        k_fit, sweep_deg, phi_c * 180.0 / PI, n_events, 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)");
@@ -541,13 +602,14 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
 
             // ===== Step B: detector distance + beam from the observed positions, cell fixed =====
             std::vector<const Partial *> obs;
-            for (const auto &pp : pts)
-                if (std::isfinite(pp.obs_x) && std::isfinite(pp.obs_y)) obs.push_back(&pp);
+            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]);
             constexpr size_t MAX_OBS = 20000;
             if (obs.size() > MAX_OBS) {
                 std::nth_element(obs.begin(), obs.begin() + MAX_OBS, obs.end(),
                                  [](const Partial *a, const Partial *b) {
-                                     return a->I / std::max(1e-9, a->sigma) > b->I / std::max(1e-9, b->sigma); });
+                                     return a->I / std::max(1e-9, static_cast<double>(a->sigma))
+                                            > b->I / std::max(1e-9, static_cast<double>(b->sigma)); });
                 obs.resize(MAX_OBS);
             }
             result.obs_used = static_cast<int>(obs.size());
@@ -557,7 +619,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
                 double c = 0.0; int n = 0;
                 for (const Partial *pp : obs) {
                     if (!in(pp->h, pp->k, pp->l, s)) continue;
-                    XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, pp->angle_rad,
+                    XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img),
                                    pp->h, pp->k, pp->l, sys);
                     double resid[3] = {0, 0, 0};
                     r(beam, dist, det_rot, rot_vec, p0, p1, p2, resid);
@@ -571,7 +633,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
                 for (const Partial *pp : obs) {
                     if (!in(pp->h, pp->k, pp->l, s)) continue;
                     p.AddResidualBlock(new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
-                        new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, pp->angle_rad,
+                        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<double *>(det_rot), rot_vec, p0, p1, p2);
@@ -625,7 +687,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
             result.distance_after_mm = dist;
             result.beam_x_before_px = beam_x0; result.beam_x_after_px = beam[0];
             result.beam_y_before_px = beam_y0; result.beam_y_after_px = beam[1];
-            result.events_used = static_cast<int>(events.size());
+            result.events_used = static_cast<int>(selected.size());
             result.ok = result.cell_refined || result.detector_refined;
             if (!result.ok)
                 logger.Info("Post-refine GEOM: neither step passed cross-validation - geometry left at nominal");