Merge branch 'multisweep-joint' into ms-trial-merge

# Conflicts:
#	rugnux/rugnux_cli.cpp

common/Reflection.h

@@ -76,6 +76,13 @@ struct ScaledFull {
     float I = NAN;
     float sigma = NAN;
     float image_number = NAN;   // the frame of the event's peak partial, which tells two events of one index apart
+    // The full's counting variance as the merge rebuilds it at an intensity I, before its error model:
+    // var_bkg + var_per_I * I + (capture * I)^2, all at the intensity's scale (what a later merge rebuilds
+    // the full from). NaN where not exported.
+    float var_bkg = NAN;
+    float var_per_I = NAN;
+    float capture = NAN;
+    float d = NAN;              // resolution, A
 };
 
 struct MergedReflection {

common/ScalingSettings.cpp

@@ -315,3 +315,12 @@ ScalingSettings& ScalingSettings::ReportResolutionRange_A(std::optional<ReportRe
 std::optional<ReportResolutionRange> ScalingSettings::GetReportResolutionRange_A() const {
     return report_resolution_range;
 }
+
+ScalingSettings& ScalingSettings::SweepSegments(std::vector<SweepSegment> input) {
+    sweep_segments = std::move(input);
+    return *this;
+}
+
+const std::vector<SweepSegment> &ScalingSettings::GetSweepSegments() const {
+    return sweep_segments;
+}

common/ScalingSettings.h

@@ -3,7 +3,9 @@
 
 #pragma once
 
+#include <cstdint>
 #include <optional>
+#include <vector>
 #include "JFJochException.h"
 
 // How the high-resolution cutoff for the written reflections and the reported shell table is chosen
@@ -17,8 +19,20 @@ struct ReportResolutionRange {
     double d_max = 0.0;
 };
 
+// One sweep of several whose frames are laid end to end in one rotation merge (the joint merge of several
+// sweeps): frames [first_frame, first_frame + n_frames). Every observation of a segment's frames is a full
+// already combined in its own sweep (partiality 1, filed under the frame of its peak), so the merge
+// combines nothing again, recomputes nothing from the rocking geometry, and lets nothing that runs along
+// the frames - the per-frame scale's smoothness above all - reach across from one sweep into the next
+// (RotationScaleMerge).
+struct SweepSegment {
+    int32_t first_frame = 0;
+    int32_t n_frames = 0;
+};
+
 class ScalingSettings {
     bool refine_wedge = false;
+    std::vector<SweepSegment> sweep_segments;
 
     bool merge_friedel = true;
     std::optional<double> high_resolution_limit_A;
@@ -150,6 +164,8 @@ public:
     ScalingSettings& ResolutionCCTarget(double input);
     ScalingSettings& ReportShellCount(int input);
     ScalingSettings& ReportResolutionRange_A(std::optional<ReportResolutionRange> input);
+    // The sweeps a rotation merge's frames are made of; empty (the default): one sweep of partials.
+    ScalingSettings& SweepSegments(std::vector<SweepSegment> input);
 
     [[nodiscard]] bool GetRefineWedge() const;
 
@@ -189,4 +205,5 @@ public:
     [[nodiscard]] double GetResolutionCCTarget() const;
     [[nodiscard]] int GetReportShellCount() const;
     [[nodiscard]] std::optional<ReportResolutionRange> GetReportResolutionRange_A() const;
+    [[nodiscard]] const std::vector<SweepSegment> &GetSweepSegments() const;
 };

docs/ACKNOWLEDGEMENT.md

@@ -384,6 +384,15 @@ J. Beilsten-Edmands, G. Winter, R. Gildea et al., "Scaling diffraction data in t
 package: algorithms and new approaches for multi-crystal scaling" (2020), Acta Cryst. D76, 385-399
 [doi:10.1107/S2059798320003198](https://doi.org/10.1107/S2059798320003198).
 
+**Several sweeps merged as one** - each sweep's partials combined into fulls on their own, and the
+fulls of all the sweeps then scaled together, every sweep with its own frame scales and a relative B
+against the joint merge, is how XSCALE merges several data sets (W. Kabsch, "Integration, scaling,
+space-group assignment and post-refinement" (2010), above). Reading the radiation dose a sweep was taken
+at from its relative B, which orders the sweeps by dose from the data rather than from the order they
+were given in, follows K. Diederichs, "Some aspects of quantitative analysis and correction of radiation
+damage" (2006), Acta Cryst. D62, 96-101
+[doi:10.1107/S0907444905031537](https://doi.org/10.1107/S0907444905031537).
+
 **Uncertainty conventions** follow the IUCr Commission on Crystallographic Nomenclature:
 D. Schwarzenbach, S. C. Abrahams, H. D. Flack et al., "Statistical descriptors in crystallography:
 Report of the IUCr Subcommittee on Statistical Descriptors" (1989), Acta Cryst. A45, 63-75

docs/RUGNUX_CLI.md

@@ -30,9 +30,10 @@ rugnux -o myrun dataset_master.h5              # NXmx / EIGER / Jungfraujoch mas
 rugnux -o myrun sweep_00001.cbf                # CBF, marCCD or SMV sweep: name one frame
 rugnux -o myrun dataset_master.h5 model.cif    # with a model (or --model model.cif)
 rugnux -o myrun x_master.h5 http://host:5232/  # follow a broker's collection while it is written
+rugnux -o myrun a_master.h5 b_master.h5        # several sweeps of one crystal, merged as one
 ```
 
-Chosen when the file's goniometer axis turns. One sweep per run. A model is placed (rigid body, or
+Chosen when the file's goniometer axis turns. Several sweeps can be given (below, Inputs). A model is placed (rigid body, or
 molecular replacement), refined and validated after merging; it does not change the processing
 decisions ([Models and reference data](RUGNUX_MODEL_AND_REFERENCE.md)).
 
@@ -91,14 +92,34 @@ spot finding or integration ([Re-scaling and re-merging](RUGNUX_ADVANCED.md#re-s
 ## Inputs
 
 Inputs are files, never a directory, read in the order given (a shell glob gives them sorted by
-name). Several inputs are read together only as one stills dataset in `--mode mx`: all the files of
+name). Several inputs are one stills dataset in `--mode mx`, or several rotation sweeps (below): all the files of
 one dataset - `rugnux run*.h5` (SACLA), `*.cxi`, `*.xtc` (LCLS) - or several Jungfraujoch/NXmx
 masters of stills, read as one run: their images one after another, merged together, with the
 indexing rate reported per master as well. Such masters must share the detector, wavelength,
-distance and beam centre. Several inputs must be one format; a rotation sweep (name one frame of a
-CBF/marCCD/SMV sweep), any `--mode` other than `mx`, and following a broker take one. What each
-format needs is in [What Rugnux reads](RUGNUX_FORMATS.md); serial XFEL stills are read without a
-geometry file where their own metadata places the detector
+distance and beam centre. Several inputs must be one format; any `--mode` other than `mx`, and
+following a broker, take one.
+
+Several rotation sweeps are told apart from the files: one sweep per CBF/marCCD/SMV file-name
+template (the name with its frame number taken out - any frame of a sweep names the whole sweep, so
+`rugnux a_*.cbf` is the one sweep `a`), one per rotation master. Only the inputs given are processed;
+other sweeps beside them are named in the log. A mix of rotation and stills masters is refused, and
+so are sweeps at different wavelengths, `-s`/`-e` and `--force-still` with several sweeps.
+By default (`--sweeps joint`) the sweeps are one dataset, merged in two levels as XSCALE merges the
+data sets XDS made of each sweep: each sweep is processed on its own first - indexed, integrated, its
+partials combined and scaled, its space group determined - writing only its log `<prefix>_s<k>_log.txt`;
+the sweeps are then put on the axes of the one with the most reflections (the change of basis, an
+axis permutation or a twin law of the lattice, under which their merged intensities correlate best;
+a sweep correlating below 0.5 on every axes is left out and named in the log), and their fulls are
+scaled and merged together: a per-frame scale within each sweep, the sweeps' dose order read from
+their relative B factors (not from the input order or the header times), a decay along that order or
+one B per sweep, whichever predicts held-out frames better, then the space-group search, the
+resolution cut and everything after it on the joint data, into `<prefix>.*`. The detector-dependent
+corrections (absorption and detector surfaces) are each sweep's own; the joint merge writes no
+unmerged MTZ yet. `--sweeps separate` processes each sweep alone instead, as if it were the only
+input, into `<prefix>_s<k>_*` (k = 1, 2, ... in the order given).
+
+What each format needs is in [What Rugnux reads](RUGNUX_FORMATS.md); serial XFEL stills are read
+without a geometry file where their own metadata places the detector
 ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)).
 
 An atomic model and a reference reflection file can be given among the inputs, in any order,
@@ -130,6 +151,7 @@ With the input already on disk, the URL is ignored.
 | `-v, --verbose` | Show the whole log on the terminal in place of the steps; `<prefix>_log.txt` is written either way |
 | `--simulate-live <ms>` | Read a stored input as the collection it was while being written: image *i* becomes readable (*i*+1)·`ms` ms after the start (0: all at once). The run takes the live path a broker URL takes, with the same result; for testing the live workflow without a broker. `--mode mx` only |
 | `--fast-live` | With `--simulate-live` or a broker URL: run as the broker's rugnux live mode does — every decision stands where it was first taken, nothing is re-checked or run again once the sweep is in, the integration radius is fixed, and SIGINT cuts the run to what it has established. A difficult crystal can come out worse |
+| `--sweeps joint\|separate` | Several rotation sweeps (see Inputs): `joint` (the default) merges them as one dataset into `<prefix>.*`; `separate` processes each alone into `<prefix>_s<k>_*` |
 | `--no-gpu` | Run on the CPU only, as on a machine without a GPU, even where one is present: CUDA is not touched (nor, on macOS, the Apple GPU that Bragg integration, the model density and the molecular-replacement rotation function otherwise run on) and every GPU step takes its CPU path (`-X FFT` indexes with FFTW and `-X FFBIDX` with its CPU implementation; `-X cuFFT` is refused) |
 
 ### Mode

image_analysis/scale_merge/RotationScaleMerge.cpp

@@ -365,6 +365,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment,
     // Decay + absorption correction surfaces are one master toggle (on by default; both cross-validated,
     // so a no-op when their systematic is absent).
     refine_decay_b = s.GetCorrectionSurfaces();
+    segments = s.GetSweepSegments();
     absorption_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0;
     modulation_iter = s.GetCorrectionSurfaces() ? s.GetAbsorptionIter() : 0;
     relative_b_deg = s.GetRelativeBDegrees();   // per-batch relative-B batch width (0 = off)
@@ -531,11 +532,28 @@ void RotationScaleMerge::IngestHost() {
     if (frame_limit >= 0) n_frames = std::min(n_frames, frame_limit);
     partials_released = false;
     resident_ingest = false;
+    frame_segment.clear();
+    if (!segments.empty()) {
+        frame_segment.assign(n_frames, -1);
+        int32_t next = 0;
+        for (size_t k = 0; k < segments.size(); ++k) {
+            const SweepSegment &seg = segments[k];
+            if (seg.first_frame != next || seg.n_frames <= 0 || seg.first_frame + seg.n_frames > n_frames)
+                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                                      "RotationScaleMerge: the sweep segments do not cover the frames in order");
+            for (int32_t f = seg.first_frame; f < seg.first_frame + seg.n_frames; ++f)
+                frame_segment[f] = static_cast<int32_t>(k);
+            next = seg.first_frame + seg.n_frames;
+        }
+        if (next != n_frames)
+            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                                  "RotationScaleMerge: the sweep segments do not cover the frames in order");
+    }
 #ifdef JFJOCH_USE_CUDA
     // Probed before the observations are built: whether the build lays down the full Obs array or
     // only the ingest arrays depends on whether the device pipeline will be active.
     gpu_ = std::make_unique<RotationScaleMergeGPU>(gpu_device);
-    gpu_active_ = gpu_->Available() && !scale_on_cpu_;
+    gpu_active_ = gpu_->Available() && !scale_on_cpu_ && segments.empty();
     resident_ingest = gpu_active_ && !observation_dump;
 #endif
 
@@ -605,7 +623,12 @@ void RotationScaleMerge::IngestHost() {
         }
     });
 
-    MeasureIncidentFlux(mean_bkg);
+    // Fulls of several sweeps carry each sweep's own background: not one flux meter, so the per-frame
+    // scale alone takes the flux.
+    if (segments.empty())
+        MeasureIncidentFlux(mean_bkg);
+    else
+        frame_flux.assign(n_frames, 1.0);
 
     // The occupied range of h, which is the buckets the sort works on.
     int b_lo = H_ALL, b_hi = -1;
@@ -728,12 +751,19 @@ void RotationScaleMerge::IngestHost() {
             }
         });
     }
-    BuildInRangeObservations(keys, static_cast<int>(total), SmoothedFrameGeometry());
+    // Pre-combined fulls have no rocking geometry left to smooth or to recompute a partiality from.
+    BuildInRangeObservations(keys, static_cast<int>(total),
+                             segments.empty() ? SmoothedFrameGeometry() : SmoothedGeometry{});
     rawrun_group.assign(rawrun_start.size(), -1);
     logger->Info("RotationScaleMerge: ingested {} partial observations from {} frames ({} distinct hkl)",
                 n_partials_obs, n_frames, rawrun_start.size());
 
-    SmoothMosaicityAndPartiality();
+    if (segments.empty()) {
+        SmoothMosaicityAndPartiality();
+    } else {
+        ingest_zeta.reset();
+        ingest_d.reset();
+    }
     ParallelChunks(static_cast<int>(partials.size()), ThreadsForWork(partials.size(), nthreads),
                    [&](int lo, int hi) {
         for (int i = lo; i < hi; ++i)
@@ -2174,6 +2204,7 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop(
         has_fulls[f] = frame_obs_count[f] > 0 ? 1 : 0;
     }
     double best_step = INFINITY, lambda = 0.0;
+    std::vector<double> segment_lambda(segments.size(), 0.0);
     bool lambda_kept = false;
     std::vector<double> steps, g_best;
     AndersonMixing anderson(SCALING_ANDERSON_DEPTH);
@@ -2188,9 +2219,16 @@ RotationScaleMerge::ScalingLoopOutcome RotationScaleMerge::RunFullsScalingLoop(
             if (fitted)
                 y[f] += g_fit[f] > 0.0 ? std::max(std::log(g_fit[f] / g[f]), -1.0) : -1.0;
         }
-        if (!lambda_kept)
-            lambda = ChooseLogScaleSmoothing(y, J, std::max(1, rocking_event_frames_at_start), nthreads);
-        const std::vector<double> x = SmoothLogScale(y, J, lambda);
+        std::vector<double> x;
+        if (segments.empty()) {
+            if (!lambda_kept)
+                lambda = ChooseLogScaleSmoothing(y, J, std::max(1, rocking_event_frames_at_start), nthreads);
+            x = SmoothLogScale(y, J, lambda);
+        } else {
+            x = SmoothLogScaleBySegment(y, J, segment_lambda, !lambda_kept);
+            std::vector<double> l = segment_lambda;
+            lambda = median_of(l);
+        }
         for (int f = 0; f < n_frames; ++f)
             g_new[f] = std::exp(x[f]);
         const double g_ref = TypicalFrameScale(g_new, n_obs);
@@ -2342,8 +2380,14 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
         return o.group >= 0 && o.corr > 0.0f && std::isfinite(o.corr) && (o.d > 0.0f)
                && o.partiality >= min_partiality;
     };
+    // Where an observation sits on the dose axis: its frame, or with several sweeps its place in the dose
+    // order the data put them in (MeasureSegmentDose).
+    auto dose = [&](const auto &o) {
+        return segments.empty() ? static_cast<double>(o.image_number)
+                                : frame_dose[o.frame] + (o.image_number - static_cast<float>(o.frame));
+    };
     auto decay_factor = [&](const auto &o, double slope) {
-        return std::exp(slope * (o.image_number - fcenter) * s2_of(o.d));
+        return std::exp(slope * (dose(o) - fcenter) * s2_of(o.d));
     };
 
     // The fields the seven passes below read (three slope fits, four held-out scores), pulled out of the
@@ -2358,6 +2402,8 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
         BuildDecayTerms(n_groups, t);
     decay_terms_kept = {};
     decay_terms_kept_valid = false;
+    if (!segments.empty())
+        MeasureSegmentDose(t, n_groups);
     const int nt = static_cast<int>(ThreadsForWork(t.term.size(), nthreads));
     std::vector<double> sw(n_groups), swI(n_groups);
 
@@ -2385,7 +2431,7 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
             const double Is = static_cast<double>(o.I) * o.corr, sc = static_cast<double>(o.sigma) * o.corr;
             if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) return;
             const double w = (Is / sc) * (Is / sc);
-            const double x = (o.image_number - fcenter) * s2_of(o.d), y = std::log(Iref / Is);
+            const double x = (dose(o) - fcenter) * s2_of(o.d), y = std::log(Iref / Is);
             a[0] += w; a[1] += w * x; a[2] += w * y; a[3] += w * x * x; a[4] += w * x * y;
         }, [](Lin &a, const Lin &b) { for (int j = 0; j < 5; ++j) a[j] += b[j]; });
         const double Sw = S[0], Sx = S[1], Sy = S[2], Sxx = S[3], Sxy = S[4];
@@ -2435,6 +2481,125 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
     // sit at the same s^2 but different frames, so the correction spreads them) - net harmful. slope =
     // 2 dB/dframe, so the total relative-B change across the run is |slope/2 * n_frames|.
     constexpr double DECAY_MIN_DELTA_B = 2.0;   // A^2, minimum total relative-B over the run to engage
+    if (!segments.empty()) {
+        // Two hypotheses of how damage accrued over several sweeps, decided on held-out frames: a B rising
+        // linearly along the dose order (frame_dose), or one B per sweep - a step at each sweep boundary,
+        // flat within a sweep. Both are regressions of y = ln(I_ref / I_obs) on s^2 = 1/(4 d^2), the ramp
+        // on x = (dose - centre) s^2, the steps on x = s^2 with a slope per sweep, at the weights the slope
+        // fit above uses and within each frame (below). Each is fitted on the even frames and its
+        // weighted squared residual taken on the odd ones, against their own reference, and the other way
+        // round; the hypothesis that predicts the held-out frames better is applied, where it beats no
+        // decay at all by the margin a correction has to clear on one sweep (and the physical floor).
+        // The held-out frames' residual is read at weights fixed by the observations, so a correction
+        // cannot win by reshaping sigma.
+        using Sums = std::array<double, 5>;   // Sw, Sx, Sy, Sxx, Sxy
+        struct Line { double a = 0.0, b = 0.0; };
+        auto line_of = [](const Sums &S) {
+            Line l;
+            const double var = S[0] > 0.0 ? S[3] - S[1] * S[1] / S[0] : 0.0;
+            if (var > 0.0) l.b = (S[4] - S[1] * S[2] / S[0]) / var;
+            if (S[0] > 0.0) l.a = (S[2] - l.b * S[1]) / S[0];
+            return l;
+        };
+        auto reference = [&](int parity) {
+            std::fill(sw.begin(), sw.end(), 0.0);
+            std::fill(swI.begin(), swI.end(), 0.0);
+            for (const DecayTerm &o : t.term) {
+                if (parity >= 0 && (o.frame & 1) != parity) continue;
+                const double sc = static_cast<double>(o.sigma) * o.corr, w = 1.0 / (sc * sc);
+                sw[o.group] += w;
+                swI[o.group] += w * static_cast<double>(o.I) * o.corr;
+            }
+        };
+        // Every term of `parity` with its y and weight against the current reference.
+        auto for_terms = [&](int parity, const auto &fn) {
+            for (const DecayTerm &o : t.term) {
+                if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) continue;
+                const double Iref = swI[o.group] / sw[o.group];
+                const double Is = static_cast<double>(o.I) * o.corr, sc = static_cast<double>(o.sigma) * o.corr;
+                if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue;
+                fn(o, std::log(Iref / Is), (Is / sc) * (Is / sc));
+            }
+        };
+        auto x_ramp = [&](const DecayTerm &o) { return (dose(o) - fcenter) * s2_of(o.d); };
+        // The terms of `parity` with y and both x taken about their frame's weighted mean: each frame has a
+        // scale of its own (the per-frame G), which took up its part of the B averaged over its
+        // reflections, so only how y varies with resolution WITHIN a frame is decay.
+        struct Point { int32_t seg; double w, y, xr, xs; };
+        auto points = [&](int parity) {
+            reference(parity);
+            std::vector<Point> pts;
+            std::vector<int32_t> frame;
+            std::vector<std::array<double, 4>> mean(n_frames, std::array<double, 4>{});   // Sw, Sy, Sxr, Sxs
+            for_terms(parity, [&](const DecayTerm &o, double y, double w) {
+                pts.push_back({frame_segment[o.frame], w, y, x_ramp(o), s2_of(o.d)});
+                frame.push_back(o.frame);
+                auto &m = mean[o.frame];
+                m[0] += w; m[1] += w * y; m[2] += w * pts.back().xr; m[3] += w * pts.back().xs;
+            });
+            for (size_t k = 0; k < pts.size(); ++k) {
+                const auto &m = mean[frame[k]];
+                pts[k].y -= m[1] / m[0]; pts[k].xr -= m[2] / m[0]; pts[k].xs -= m[3] / m[0];
+            }
+            return pts;
+        };
+        struct Fit { double ramp = 0.0; std::vector<double> step; };
+        auto fit = [&](int parity) {
+            Sums r{};
+            std::vector<Sums> st(segments.size(), Sums{});
+            for (const Point &q : points(parity)) {
+                r[0] += q.w; r[1] += q.w * q.xr; r[2] += q.w * q.y; r[3] += q.w * q.xr * q.xr; r[4] += q.w * q.xr * q.y;
+                Sums &a = st[q.seg];
+                a[0] += q.w; a[1] += q.w * q.xs; a[2] += q.w * q.y; a[3] += q.w * q.xs * q.xs; a[4] += q.w * q.xs * q.y;
+            }
+            Fit out{line_of(r).b, {}};
+            for (const Sums &a : st) out.step.push_back(line_of(a).b);
+            return out;
+        };
+        // Held-out weighted squared residuals of no decay, the ramp and the steps.
+        auto score = [&](int test, const Fit &m, std::array<double, 3> &ssr) {
+            for (const Point &q : points(test)) {
+                const double e1 = q.y - m.ramp * q.xr, e2 = q.y - m.step[q.seg] * q.xs;
+                ssr[0] += q.w * q.y * q.y; ssr[1] += q.w * e1 * e1; ssr[2] += q.w * e2 * e2;
+            }
+        };
+        std::array<double, 3> ssr{};
+        const Fit even = fit(0), odd = fit(1);
+        score(1, even, ssr);
+        score(0, odd, ssr);
+        const Fit all = fit(-1);
+
+        const double slope = all.ramp;   // the slope decay_factor applies
+        std::vector<double> b(segments.size());
+        double b_mid = 0.0, b_lo = INFINITY, b_hi = -INFINITY;
+        for (size_t k = 0; k < segments.size(); ++k) {
+            b[k] = 0.5 * all.step[k];
+            b_mid += b[k] * segments[k].n_frames / static_cast<double>(n_frames);
+            b_lo = std::min(b_lo, b[k]);
+            b_hi = std::max(b_hi, b[k]);
+        }
+        const double delta_lin = std::fabs(0.5 * slope * n_frames), delta_step = b_hi - b_lo;
+        const bool lin_ok = delta_lin >= DECAY_MIN_DELTA_B && ssr[1] < (1.0 - CV_MIN_RELATIVE_GAIN) * ssr[0];
+        const bool step_ok = delta_step >= DECAY_MIN_DELTA_B && ssr[2] < (1.0 - CV_MIN_RELATIVE_GAIN) * ssr[0];
+        const bool use_step = step_ok && (!lin_ok || ssr[2] < ssr[1]);
+        const bool use_lin = lin_ok && !use_step;
+        segment_decay_model = use_step ? "per-sweep B" : use_lin ? "dose-linear" : "none";
+        const double rel = ssr[0] > 0.0 ? 100.0 / ssr[0] : 0.0;
+        logger->Info("Decay over {} sweeps, held-out residual against no decay: dose-linear (dB = {:.2f} A^2 over "
+                    "the dose order) {:.1f}%, one B per sweep (spread {:.2f} A^2) {:.1f}% => {}", segments.size(),
+                    delta_lin, ssr[1] * rel, delta_step, ssr[2] * rel, segment_decay_model);
+        if (use_step || use_lin)
+            ParallelChunks(static_cast<int>(fulls.size()), ThreadsForWork(fulls.size(), nthreads), [&](int lo, int hi) {
+                for (int i = lo; i < hi; ++i) {
+                    Obs &o = fulls[i];
+                    if (!usable(o)) continue;
+                    const double f = use_step ? std::exp(2.0 * (b[frame_segment[o.frame]] - b_mid) * s2_of(o.d))
+                                              : decay_factor(o, slope);
+                    o.corr = static_cast<float>(o.corr * f);
+                }
+            });
+        return;
+    }
     const double slope = fit_slope(-1);
     const double total_delta_B = std::fabs(0.5 * slope * n_frames);
     if (total_delta_B < DECAY_MIN_DELTA_B) {
@@ -2460,6 +2625,92 @@ void RotationScaleMerge::RefineDecay(int n_groups) {
                 0.5 * slope * n_frames, 0.5 * slope);
 }
 
+std::vector<double> RotationScaleMerge::SegmentB(const DecayTerms &t, int n_groups, int parity) const {
+    // Per segment, the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against
+    // s^2 = 1/(4 d^2), which is 2 B, with an intercept of its own (the scale is the per-frame G's), over
+    // the frames of `parity` (< 0: all) against a reference made of the same frames. NaN where a segment
+    // has nothing to fit.
+    const size_t ns = segments.size();
+    std::vector<double> sw(n_groups, 0.0), swI(n_groups, 0.0);
+    for (const DecayTerm &o : t.term) {
+        if (parity >= 0 && (o.frame & 1) != parity) continue;
+        const double sc = static_cast<double>(o.sigma) * o.corr, w = 1.0 / (sc * sc);
+        sw[o.group] += w;
+        swI[o.group] += w * static_cast<double>(o.I) * o.corr;
+    }
+    std::vector<std::array<double, 5>> S(ns, std::array<double, 5>{});   // Sw, Sx, Sy, Sxx, Sxy
+    for (const DecayTerm &o : t.term) {
+        if ((parity >= 0 && (o.frame & 1) != parity) || sw[o.group] <= 0.0) continue;
+        const double Iref = swI[o.group] / sw[o.group];
+        const double Is = static_cast<double>(o.I) * o.corr, sc = static_cast<double>(o.sigma) * o.corr;
+        if (!std::isfinite(Iref) || Iref <= 0.0 || !(Is > 0.0)) continue;
+        const double w = (Is / sc) * (Is / sc);
+        const double x = 1.0 / (4.0 * static_cast<double>(o.d) * o.d), y = std::log(Iref / Is);
+        auto &a = S[frame_segment[o.frame]];
+        a[0] += w; a[1] += w * x; a[2] += w * y; a[3] += w * x * x; a[4] += w * x * y;
+    }
+    std::vector<double> b(ns, NAN);
+    for (size_t k = 0; k < ns; ++k) {
+        const auto &a = S[k];
+        const double var = a[0] > 0.0 ? a[3] - a[1] * a[1] / a[0] : 0.0;
+        if (var > 0.0)
+            b[k] = 0.5 * (a[4] - a[1] * a[2] / a[0]) / var;
+    }
+    return b;
+}
+
+void RotationScaleMerge::MeasureSegmentDose(const DecayTerms &t, int n_groups) {
+    // A relative B per data set against the joint merge, as XSCALE fits it: Kabsch (2010) Acta Cryst.
+    // D66, 133-144; the dose read from it: Diederichs (2006) Acta Cryst. D62, 96-101.
+    // Radiation damage raises the relative B of a sweep collected later, so the B factors of the
+    // segments against their joint merge give the order the crystal took the dose in - read from the
+    // data, not from the order the sweeps were given in or the times their headers state. Per segment,
+    // the inverse-variance-of-log weighted slope of ln(I_ref / I_obs) against s^2 = 1/(4 d^2) is 2 B,
+    // with an intercept of its own (the scale is the per-frame G's). The least damaged segment is B = 0,
+    // and its frames start the dose axis.
+    const size_t ns = segments.size();
+    const std::vector<double> b = SegmentB(t, n_groups, -1);
+    // A segment with nothing to measure keeps its place in the input order, after the measured ones.
+    std::vector<int> order(ns);
+    std::iota(order.begin(), order.end(), 0);
+    std::stable_sort(order.begin(), order.end(), [&](int p, int q) {
+        if (std::isfinite(b[p]) != std::isfinite(b[q])) return std::isfinite(b[p]);
+        return std::isfinite(b[p]) && b[p] < b[q];
+    });
+    const double b0 = std::isfinite(b[order[0]]) ? b[order[0]] : 0.0;
+    frame_dose.assign(n_frames, 0.0);
+    double offset = 0.0;
+    std::string text;
+    for (size_t rank = 0; rank < ns; ++rank) {
+        const int k = order[rank];
+        const SweepSegment &seg = segments[k];
+        for (int32_t f = seg.first_frame; f < seg.first_frame + seg.n_frames; ++f)
+            frame_dose[f] = offset + (f - seg.first_frame);
+        offset += seg.n_frames;
+        segment_results[k].relative_b = b[k] - b0;
+        segment_results[k].dose_rank = static_cast<int>(rank);
+        text += fmt::format("{}segment {} (B {:+.2f} A^2)", rank == 0 ? "" : ", ", k + 1, b[k] - b0);
+    }
+    logger->Info("Dose order of the {} sweep segments, from their relative B against the joint merge: {}",
+                ns, text);
+}
+
+std::vector<double> RotationScaleMerge::SmoothLogScaleBySegment(const std::vector<double> &y,
+                                                                const std::vector<double> &J,
+                                                                std::vector<double> &lambda, bool choose) const {
+    std::vector<double> x(y.size());
+    for (size_t k = 0; k < segments.size(); ++k) {
+        const auto lo = y.begin() + segments[k].first_frame, hi = lo + segments[k].n_frames;
+        const auto jlo = J.begin() + segments[k].first_frame;
+        const std::vector<double> ys(lo, hi), Js(jlo, jlo + segments[k].n_frames);
+        if (choose)
+            lambda[k] = ChooseLogScaleSmoothing(ys, Js, std::max(1, rocking_event_frames_at_start), nthreads);
+        const std::vector<double> xs = SmoothLogScale(ys, Js, lambda[k]);
+        std::copy(xs.begin(), xs.end(), x.begin() + segments[k].first_frame);
+    }
+    return x;
+}
+
 // A^2, the bound the Gauss-Seidel solve clamps each batch's relative-B to. It is a numerical guard against
 // divergence, never an answer: a run of identical values in the solved curve means the solve reached the
 // bound, not that the crystal did anything. The CORRECTION keeps the tight bound it was tuned with; the
@@ -4716,12 +4967,15 @@ void RotationScaleMerge::Combine() {
         size_t i = 0;
         while (i < ev.size()) {
             size_t kk = i + 1;
-            float last_frame = partials[ev[i]].r->image_number;
-            while (kk < ev.size()) {
-                const float frame = partials[ev[kk]].r->image_number;
-                if (frame - last_frame > max_frame_gap) break;
-                last_frame = frame;
-                ++kk;
+            // A pre-combined full is a whole event already (SweepSegment).
+            if (segments.empty()) {
+                float last_frame = partials[ev[i]].r->image_number;
+                while (kk < ev.size()) {
+                    const float frame = partials[ev[kk]].r->image_number;
+                    if (frame - last_frame > max_frame_gap) break;
+                    last_frame = frame;
+                    ++kk;
+                }
             }
 
             double pooled_bkg = 0.0;
@@ -5860,18 +6114,31 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     // reflections are) - each with the sigma the merge weighted it by.
     if (export_scaled_fulls && full_stats && !for_search) {
         // Which fulls, in order, and then each one's record on its own - on all threads.
+        // The fulls of one sweep of several (uncut) also hand back the reflections the group extinguishes:
+        // a joint merge searches its own space group, and the screw axes are decided on exactly these.
+        const gemmi::GroupOps ops = Group().operations();
+        auto extinguished = [&](const Obs &o) {
+            return export_scaled_fulls_uncut && o.group < 0 && frame_cell_ok[o.frame] && o.corr > 0.0f
+                   && std::isfinite(o.corr) && std::isfinite(o.I) && o.d > 0.0f
+                   && ops.is_systematically_absent({{o.h, o.k, o.l}});
+        };
         const std::vector<int32_t> kept = KeptIndices(n_full, nthreads, [&](int i) {
             const int g = mf.group[i];
-            return g >= 0 && !rejected_obs[i] && std::isfinite(merged_I[g])
-                   && !(effective_d_min && acc[g].d < *effective_d_min);
+            if (g < 0)
+                return extinguished(fulls[i]);
+            return !rejected_obs[i] && std::isfinite(merged_I[g])
+                   && (export_scaled_fulls_uncut || !(effective_d_min && acc[g].d < *effective_d_min));
         });
         result.scaled_fulls.resize(kept.size());
         ParallelChunks(static_cast<int>(kept.size()), ThreadsForWork(kept.size(), nthreads), [&](int lo, int hi) {
             for (int j = lo; j < hi; ++j) {
                 const Obs &o = fulls[kept[j]];
                 const float I_corr = o.I * o.corr;
-                result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr),
-                                          o.image_number};
+                const float sigma = o.group >= 0 ? corrected_sigma(o, I_corr, o.sigma * o.corr)
+                                                 : model_sigma(o, I_corr, o.sigma * o.corr);
+                result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, sigma,
+                                          o.image_number, o.corr * o.corr * o.var_bkg, o.corr * o.var_per_I,
+                                          o.capture, o.d};
             }
         });
     }
@@ -6697,6 +6964,12 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search, bool full_st
                                                   const gemmi::SpaceGroup *group) {
     const ScopedGpuPin pin(gpu_device);   // see IngestHost
     run_group = group;
+    // Pre-combined fulls: there are no partials to fit a scale on, and every event is one frame.
+    if (!segments.empty()) {
+        scale_partials_choice = false;
+        rocking_event_frames_at_start = 1;
+        rocking_events_at_start = n_partials_obs;
+    }
     struct ClearGroup { const gemmi::SpaceGroup *&g; ~ClearGroup() { g = nullptr; } } clear_group{run_group};
     // Whether a frame's scale can be fitted on its partials depends on whether the partials tell a change
     // of scale from an error of the partiality model. The count of rocking events per frame is the prior
@@ -7047,7 +7320,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
     // kernels on it, and the restore re-uploads.
 
     // --- 2a. On the de-novo search pass only, drop observations whose Lorentz geometry is poor. ---
-    if (for_search && search_min_zeta > 0.0) {
+    if (for_search && search_min_zeta > 0.0 && segments.empty()) {   // a full has no rocking geometry
         int64_t n_dropped = 0;
         bool zeta_on_gpu = false;
 #ifdef JFJOCH_USE_CUDA
@@ -7357,20 +7630,27 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
     //         is a correction to the merged INTENSITIES, and a pass whose intensities are thrown away has
     //         nowhere to put one - the fitted gain lives on the fulls, which the next pass rebuilds. ---
     const bool refine_surfaces = !for_search && full_stats;
+    // What runs along the frames of ONE sweep - the damage monitor, the sweep ledger, the per-batch B and
+    // the surfaces, each read in one goniometer and one detector setting - has no meaning across the
+    // segments of several, and is left out there. The decay is fitted along the dose order instead.
+    const bool along_frames = segments.empty();
     const bool corrections =
-        refine_surfaces && (refine_decay_b || absorption_iter > 0 || modulation_iter > 0 || relative_b_deg > 0.0);
+        refine_surfaces && (refine_decay_b || (along_frames && (absorption_iter > 0 || modulation_iter > 0
+                                                                || relative_b_deg > 0.0)));
+    segment_results.assign(segments.size(), Result::Segment{});
+    segment_decay_model.clear();
     // Radiation-damage monitor: measure the per-batch relative-B on the scaled fulls BEFORE any correction
     // (report-only; captures the full damage signature, not a residual). Skipped on the P1 search pass.
     std::unique_ptr<TimingMark> sub_mark;
     auto sub = [&](const char *n) { sub_mark.reset(); sub_mark = std::make_unique<TimingMark>(*logger, std::string("rsm/s ") + n); };
     sub("damage B");
-    if (refine_surfaces)
+    if (refine_surfaces && along_frames)
         MeasureRadiationDamageB(n_groups);
     // Sweep-quality diagnostic, on the per-frame scale the partial scaling just fitted (the flux is
     // already out of it) and the per-frame CC computed above. It drops nothing itself; the ranges it
     // segments are the ledger 4c below decides on.
     sub("sweep quality");
-    if (refine_surfaces)
+    if (refine_surfaces && along_frames)
         MeasureSweepQuality(partial_scaled, cc, cc_n);
     // CC1/2 of these same fulls with the surfaces not yet folded in, for the caller's cross-pass quality
     // guard: the pass this one is judged against fits no surfaces of its own, so it is the uncorrected
@@ -7396,7 +7676,7 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
     decay_terms_kept = {};
     decay_terms_kept_valid = false;
     sub("relative B");
-    if (refine_surfaces && relative_b_deg > 0.0)
+    if (refine_surfaces && along_frames && relative_b_deg > 0.0)
         RefineRelativeB(n_groups);   // per-batch relative-B on top of the single decay slope
     // The three surfaces overlap, so the one fitted first gets the first claim and the others are tested
     // on what it left - the order decides what is adopted. Detector first: every frame of the sweep
@@ -7412,16 +7692,16 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
     // the grid is then tested on what the smooth surface left, and where that surface is refused the
     // grid sees exactly what it saw before.
     sub("modulation");
-    if (refine_surfaces && modulation_iter > 0)
+    if (refine_surfaces && along_frames && modulation_iter > 0)
         RefineModulation(modulation_iter, n_groups);
     sub("absorption time");
-    if (refine_surfaces && absorption_iter > 0)
+    if (refine_surfaces && along_frames && absorption_iter > 0)
         RefineAbsorptionTime(absorption_iter, n_groups);
     sub("absorption SH");
-    if (refine_surfaces && absorption_iter > 0)
+    if (refine_surfaces && along_frames && absorption_iter > 0)
         RefineAbsorptionSH(absorption_iter, n_groups);
     sub("absorption gonio");
-    if (refine_surfaces && absorption_iter > 0)
+    if (refine_surfaces && along_frames && absorption_iter > 0)
         RefineAbsorption(absorption_iter, n_groups);
     sub("disposition");
 
@@ -7462,6 +7742,17 @@ RotationScaleMerge::Result RotationScaleMerge::RunArm(bool scale_partials, bool
     // its CC1/2 is the answer without a second merge. A caller that did not ask keeps the NaN.
     if (measure_cc_before_corrections)
         r.cc_half_before_corrections = cc_half_before_corrections.value_or(r.statistics.overall.cc_half);
+    // Each segment's scale: the median of its frames' scales, at the run's gauge.
+    for (size_t k = 0; k < segments.size(); ++k) {
+        std::vector<double> g;
+        for (int32_t f = segments[k].first_frame; f < segments[k].first_frame + segments[k].n_frames; ++f)
+            if (fulls_frame_count[f] > 0 && std::isfinite(g_full[f]) && g_full[f] > 0.0)
+                g.push_back(g_full[f]);
+        if (!g.empty())
+            segment_results[k].scale = median_of(g);
+    }
+    r.segments = segment_results;
+    r.segment_decay_model = segment_decay_model;
     // The fulls and the merge's working arrays are rebuilt by every arm before they are read, so they go
     // with it: kept, they sat beside the next arm's combine and through everything the caller did between
     // merges - on a large cell in P1, the fulls and about 330 bytes per ASU group, gigabytes. Faulting them

image_analysis/scale_merge/RotationScaleMerge.h

@@ -11,6 +11,7 @@
 #include <memory>
 #include <new>
 #include <optional>
+#include <string>
 #include <type_traits>
 #include <utility>
 #include <vector>
@@ -139,6 +140,16 @@ public:
         // The fulls the merge kept, unmerged (see ScaledFull): what passed the outlier tests, inside the
         // resolution cut. Filled only when SetExportScaledFulls asked for it, and never on a search merge.
         std::vector<ScaledFull> scaled_fulls;
+        // Per sweep segment (SetSweepSegments), in segment order; empty without segments.
+        struct Segment {
+            double scale = NAN;         // median per-frame scale of the segment's fulls (run gauge)
+            double relative_b = NAN;    // A^2, against the joint merge, the least damaged segment at 0
+            int dose_rank = -1;         // the segment's place in the dose order those B factors decide
+        };
+        std::vector<Segment> segments;
+        // How the decay over the segments was modelled - "dose-linear" (along the dose order), "per-sweep B"
+        // (one B a segment) or "none" - as the held-out frames decided; empty without segments or decay.
+        std::string segment_decay_model;
     };
 
     // experiment: read live (its space group is changed by the caller between Run() calls).
@@ -288,6 +299,16 @@ public:
     // Merge only the first n frames of the outcomes (by default all of them). Set before Ingest().
     void SetFrameLimit(int n) { frame_limit = n; }
 
+    // The outcomes are several sweeps of pre-combined fulls laid end to end (SweepSegment; by default the
+    // experiment's ScalingSettings::GetSweepSegments). Set before Ingest(); the segments must cover the
+    // frames in order. Such an engine runs on the host: its fulls are a few per frame where a sweep's
+    // partials are millions, and the device pipeline is the partials'. Without segments nothing of this runs.
+    void SetSweepSegments(std::vector<SweepSegment> s) { segments = std::move(s); }
+
+    // Whether the fulls handed back (SetExportScaledFulls) are all the merge kept, past its resolution cut
+    // too: the fulls of one sweep of several, which the joint merge cuts again on its own.
+    void SetExportScaledFullsUncut(bool on) { export_scaled_fulls_uncut = on; }
+
     // The GPU the engine runs on (device 0 by default). Set before Ingest(). The device only says where
     // the merge runs: on cards of one model the result is the same on any of them.
     void SetGpuDevice(int device) { gpu_device = device; }
@@ -724,8 +745,18 @@ private:
     std::vector<int64_t> last_cc_n;
     std::vector<uint8_t> last_frame_scaled;
     int frame_limit = -1;                     // see SetFrameLimit
+    std::vector<SweepSegment> segments;       // see SetSweepSegments
+    std::vector<int32_t> frame_segment;       // each frame's segment (Ingest; empty without segments)
+    std::vector<Result::Segment> segment_results;   // the last arm's, copied into its Result
+    // Each frame's place on the dose axis the decay is fitted along: the frames of the segments laid end
+    // to end in the order the segments' relative B factors put them (MeasureSegmentDose).
+    std::vector<double> frame_dose;
+    // SmoothLogScale / ChooseLogScaleSmoothing within each segment on its own; lambda per segment.
+    std::vector<double> SmoothLogScaleBySegment(const std::vector<double> &y, const std::vector<double> &J,
+                                                std::vector<double> &lambda, bool choose) const;
     int gpu_device = 0;                       // see SetGpuDevice
     bool export_scaled_fulls = false;         // see SetExportScaledFulls
+    bool export_scaled_fulls_uncut = false;   // see SetExportScaledFullsUncut
     bool french_wilson = true;                // see SetFrenchWilson
 
     // --- helpers (each a flat pass; see the .cpp) ---
@@ -958,6 +989,12 @@ private:
         std::vector<int32_t> g_start;
     };
     void BuildDecayTerms(int n_groups, DecayTerms &t) const;
+    // Per segment: the relative B against the merged reference, and the dose order it puts the segments
+    // in (segment_results, frame_dose).
+    void MeasureSegmentDose(const DecayTerms &t, int n_groups);
+    // Each segment's B against a reference made of the same frames: those of `parity` (< 0: all).
+    std::vector<double> SegmentB(const DecayTerms &t, int n_groups, int parity) const;
+    std::string segment_decay_model;   // the last decay's choice (Result::segment_decay_model)
     // The terms the radiation-damage monitor built, handed on to RefineDecay: nothing between the two
     // touches corr, so they are the terms RefineDecay would build. Empty outside that stretch of RunArm.
     DecayTerms decay_terms_kept;

rugnux/CMakeLists.txt

@@ -17,6 +17,8 @@ ADD_LIBRARY(JFJochRugnux STATIC
         RugnuxFirstPassRescues.cpp
         RugnuxFirstPassAccept.cpp
         RugnuxScaleMerge.cpp
+        RugnuxJointSweeps.cpp
+        RugnuxJointSweeps.h
         RugnuxDefaults.cpp
         RugnuxDefaults.h
         RugnuxLive.cpp
@@ -40,7 +42,8 @@ TARGET_LINK_LIBRARIES(JFJochRugnux JFJochReader JFJochImageAnalysis JFJochPrevie
 
 # The files of a run, written by the programs that run one (rugnux, the viewer) from what JFJochRugnux
 # returns: the library itself writes none.
-ADD_LIBRARY(JFJochRugnuxOutput STATIC RugnuxOutputFiles.cpp RugnuxOutputFiles.h RugnuxInputKind.cpp RugnuxInputKind.h)
+ADD_LIBRARY(JFJochRugnuxOutput STATIC RugnuxOutputFiles.cpp RugnuxOutputFiles.h RugnuxInputKind.cpp RugnuxInputKind.h
+        RugnuxSweeps.cpp RugnuxSweeps.h)
 TARGET_LINK_LIBRARIES(JFJochRugnuxOutput JFJochRugnux JFJochWriter JFJochModelWriter ZLIB::ZLIB)
 
 # rugnux is the single offline analysis CLI; --mode picks what it does: mx (the full pipeline, the

rugnux/Rugnux.cpp

@@ -823,6 +823,21 @@ std::vector<bool> Rugnux::PrepassSubset(int prepass_images) const {
     return keep;
 }
 
+ProcessResult Rugnux::RunJoint(std::vector<IntegrationOutcome> outcomes, const UnitCell &cell,
+                               gemmi::CrystalSystem system, char centring, RugnuxObserver *observer) {
+    const auto run_start = std::chrono::steady_clock::now();
+    joint_outcomes_ = std::move(outcomes);
+    joint_cell_ = cell;
+    joint_system_ = system;
+    joint_centring_ = centring;
+    decide_late_ = true;
+    ProcessResult result = RunPipeline(observer, /*write_output=*/true, /*geometry_prepass=*/false);
+    joint_outcomes_.reset();
+    result.pass_count = 1;
+    result.total_time_s = std::chrono::duration<double>(std::chrono::steady_clock::now() - run_start).count();
+    return result;
+}
+
 ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, bool geometry_prepass) {
     Logger logger("Rugnux");
     ProcessResult result;
@@ -854,7 +869,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     if (config_.stride <= 0)
         throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image stride must be positive");
 
-    const auto total_images_in_file = static_cast<int>(reader_.GetNumberOfImages());
+    const bool joint = joint_outcomes_.has_value();
+    const auto total_images_in_file = joint ? static_cast<int>(joint_outcomes_->size())
+                                            : static_cast<int>(reader_.GetNumberOfImages());
     int end_image = config_.end_image;
     if (end_image < 0 || end_image > total_images_in_file)
         end_image = total_images_in_file;
@@ -943,7 +960,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
 
     // After any geometry refinement, so the shadow is found about the beam centre actually used,
     // and before the azimuthal mapping and the output mask below, which both read pixel_mask_.
-    if (config_.detect_beam_stop.has_value() || config_.adaptive_integration_radius)
+    if (!joint && (config_.detect_beam_stop.has_value() || config_.adaptive_integration_radius))
         PreScan(start_image, early_images,
                 std::min(config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES),
                          HotPixelFinder::MAX_SAMPLED_FRAMES), observer);
@@ -1156,6 +1173,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
         // (see the reflection mask in BraggIntegrationEngineCPU). A fixed group used to cost real
         // intensity accuracy for that reason.
         indexer->PredictAllCentringNodes(true);
+        if (joint) {
+            indexer->GetIntegrationOutcome() = std::move(*joint_outcomes_);
+            joint_outcomes_.reset();
+        }
     }
 
     const auto start_time = std::chrono::steady_clock::now();
@@ -1192,7 +1213,19 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
             .scan_result = scan_result, .unmerged_mtz = unmerged_mtz};
 
     // First pass of two-pass rotation indexing (full analysis only).
-    if (indexing && force_rotation_result_.has_value()) {
+    if (joint) {
+        // The sweeps were indexed, integrated and combined each on its own; their fulls are the frames.
+        result.rotation_lattice_found = true;
+        result.indexing_rate = 1.0f;
+        result.images_processed = static_cast<uint64_t>(images_to_process);
+        result.consensus_cell = joint_cell_;
+        result.consensus_centering = joint_centring_;
+        end_msg.rotation_lattice = CrystalLattice(*joint_cell_);
+        end_msg.rotation_lattice_type = LatticeMessage{.centering = joint_centring_, .niggli_class = 0,
+                                                       .crystal_system = joint_system_};
+        result.rotation_lattice_type = end_msg.rotation_lattice_type;
+        end_msg.unit_cell = result.consensus_cell;
+    } else if (indexing && force_rotation_result_.has_value()) {
         // Supercell-collapse fallback: force pass-1's WHOLE indexing result (lattice + refined orientation /
         // search metadata / axis), not just its lattice - a lattice-only force loses that metadata and
         // integrates markedly worse.
@@ -1230,7 +1263,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     //
     // Only an output pass probes: the answer the run returns is always one (RunAllPasses), and no other
     // pass's probe is read - the pre-pass, the beam-centre arm and the probe passes leave theirs unread.
-    if (indexing && indexer && experiment_.IsRotationIndexing() && write_output && !geometry_prepass) {
+    if (!joint && indexing && indexer && experiment_.IsRotationIndexing() && write_output && !geometry_prepass) {
         std::vector<int64_t> frames;
         for (const int o : select_equally_spaced_image_ordinals(images_to_process, 60))
             frames.push_back(o);
@@ -1241,7 +1274,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
     std::optional<GPUImagePass> gpu_image_pass;   // what a probe that does not fit beside it waits for
     if (!gpu_work_beside_)
         gpu_image_pass.emplace();
-    ImagePass(p);
+    if (!joint)
+        ImagePass(p);
     gpu_image_pass.reset();
     image_pass_mark.Stop();
 

rugnux/Rugnux.h

@@ -115,6 +115,11 @@ struct ProcessConfig {
     // scaling.
     bool write_merged = true;
 
+    // One sweep of several merged as one (--sweeps joint): the run is made as on its own, but its written
+    // merge hands back all of its fulls, past the resolution cut too (MergedOutput::scaled_fulls), for the
+    // joint merge to scale and cut again (Rugnux::RunJoint).
+    bool sweep_member = false;
+
     SpotFindingSettings spot_finding;   // FullAnalysis spot finding (and Calibration by spots)
 
     // Calibration: how the powder rings are measured, and the calibrant's rings as q = 2*pi/d (1/A).
@@ -1289,6 +1294,11 @@ class Rugnux {
     // Whether the first passes decide their near ties on all the frames: offline, on a re-check, and on
     // the run made again.
     bool decide_late_ = false;
+    // RunJoint's frames and cell, for the one pass it makes (RunPipeline); empty otherwise.
+    std::optional<std::vector<IntegrationOutcome>> joint_outcomes_;
+    std::optional<UnitCell> joint_cell_;
+    gemmi::CrystalSystem joint_system_ = gemmi::CrystalSystem::Triclinic;
+    char joint_centring_ = 'P';
     // Whether what the run otherwise decides on the pre-pass frames - the rotation-scale walk's refits,
     // the demoted metric-symmetry arm, the integration radius (RefinedPasses) - is decided on all the
     // frames: only on the run made again. Offline too these are decided on the pre-pass frames, and a near
@@ -1691,6 +1701,14 @@ public:
     // Runs the configured workflow to completion or until Cancel(). Throws on setup failure.
     ProcessResult Run(RugnuxObserver *observer = nullptr);
 
+    // The joint merge of several sweeps (--sweeps joint): `outcomes` are the frames of all of them laid
+    // end to end, each holding the fulls its own sweep combined (ScalingSettings::SweepSegments says which
+    // frames are whose), indexed on the axes of `cell`. Nothing is read from the images: the merge, the
+    // space-group search and everything after it run on these as on the outcomes of one sweep.
+    // `system` and `centring` are those of the lattice the cell is the conventional cell of.
+    ProcessResult RunJoint(std::vector<IntegrationOutcome> outcomes, const UnitCell &cell,
+                           gemmi::CrystalSystem system, char centring, RugnuxObserver *observer = nullptr);
+
     // The pinned host buffers the engines of a run kept for each other (CUDAMemHelpers.h), handed back to the
     // driver. Not part of Run: a process that outlives the run (the broker, the viewer) calls it once the
     // results are out, as it takes a tenth of a second the user would otherwise wait for; one that ends

rugnux/RugnuxJointSweeps.cpp

@@ -0,0 +1,89 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RugnuxJointSweeps.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "../common/CrystalLattice.h"
+#include "../image_analysis/scale_merge/ReindexAmbiguity.h"
+
+JointSweepRelation RelateJointSweeps(const std::vector<JointMemberSweep> &sweeps) {
+    JointSweepRelation out;
+    out.sweeps.resize(sweeps.size());
+    if (sweeps.empty())
+        return out;
+    out.reference = 0;
+    for (size_t k = 1; k < sweeps.size(); ++k)
+        if (sweeps[k].fulls.size() > sweeps[out.reference].fulls.size())
+            out.reference = static_cast<int>(k);
+    const JointMemberSweep &ref = sweeps[out.reference];
+    out.sweeps[out.reference] = {true, gemmi::Op::identity(), 1.0, 1.0, {}};
+    for (size_t k = 0; k < sweeps.size(); ++k) {
+        if (static_cast<int>(k) == out.reference)
+            continue;
+        JointSweepRelation::Sweep &s = out.sweeps[k];
+        if (sweeps[k].merged.empty()) {
+            s.why_left_out = "it merged nothing";
+            continue;
+        }
+        const ReferenceSettingChoice choice = ChooseReferenceSetting(sweeps[k].merged, sweeps[k].cell,
+                                                                     sweeps[k].space_group, ref.merged, ref.cell,
+                                                                     ref.space_group);
+        s.cob = choice.cob;
+        s.cc = choice.cc;
+        s.identity_cc = choice.identity_cc;
+        s.included = std::isfinite(choice.cc) && choice.cc >= MIN_REFERENCE_CC;
+        if (!s.included)
+            s.why_left_out = "its intensities correlate with the reference sweep's on no axes (best CC "
+                             + std::to_string(choice.cc) + ")";
+    }
+    return out;
+}
+
+JointMergeInput BuildJointMergeInput(const std::vector<JointMemberSweep> &sweeps,
+                                     const JointSweepRelation &relation) {
+    JointMergeInput out;
+    const CrystalLattice lattice(sweeps[relation.reference].cell);
+    int32_t first = 0;
+    for (size_t k = 0; k < sweeps.size(); ++k) {
+        const auto &sweep = sweeps[k];
+        const auto &rel = relation.sweeps[k];
+        if (!rel.included || sweep.n_frames <= 0)
+            continue;
+        out.segments.push_back({first, static_cast<int32_t>(sweep.n_frames)});
+        out.segment_sweep.push_back(static_cast<int>(k));
+        out.outcomes.resize(first + sweep.n_frames);
+        for (int64_t f = 0; f < sweep.n_frames; ++f)
+            out.outcomes[first + f].latt = lattice;
+        for (const ScaledFull &full : sweep.fulls) {
+            if (!std::isfinite(full.I) || !std::isfinite(full.var_bkg) || !(full.var_per_I > 0.0f)
+                || !(full.d > 0.0f))
+                continue;
+            const auto frame = std::clamp<int64_t>(std::llround(full.image_number), 0, sweep.n_frames - 1);
+            const gemmi::Op::Miller hkl = rel.cob.apply_to_hkl({{full.h, full.k, full.l}});
+            // A full as the merge takes a pre-combined observation, so that the combine rebuilds the variance
+            // its own sweep's merge gave it: a count c of intensity I/c at correction c = var_per_I is a
+            // variance c^2 var_bkg' + c I at I, which is the sweep's var_bkg + var_per_I I with
+            // var_bkg' = var_bkg / c^2 (the capture term folded in at the full's own intensity).
+            const float c = full.var_per_I;
+            const float capture_var = std::isfinite(full.capture) ? full.capture * full.I * full.capture * full.I : 0.0f;
+            Reflection r{};
+            r.h = hkl[0]; r.k = hkl[1]; r.l = hkl[2];
+            r.image_number = static_cast<float>(first + frame);
+            r.d = full.d;
+            r.I = full.I / c;
+            r.var_bkg = (full.var_bkg + capture_var) / (c * c);
+            r.sigma = std::sqrt(std::max(r.var_bkg + std::max(0.0f, r.I), 1e-12f));
+            r.prescaling_corr = c;
+            r.partiality = 1.0f;
+            r.zeta = 1.0f;
+            r.image_scale_corr = c;
+            r.observed = true;
+            out.outcomes[first + frame].reflections.push_back(r);
+        }
+        first += static_cast<int32_t>(sweep.n_frames);
+    }
+    return out;
+}

rugnux/RugnuxJointSweeps.h

@@ -0,0 +1,58 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "../common/Reflection.h"
+#include "../common/ScalingSettings.h"
+#include "../common/UnitCell.h"
+#include "../gemmi_gph/gemmi/symmetry.hpp"
+#include "../image_analysis/IntegrationOutcome.h"
+
+// Several sweeps of one crystal merged as one (--sweeps joint), in two levels: each sweep is processed on
+// its own - indexed, integrated, its partials combined into fulls and scaled - and the fulls of all of them
+// are then scaled and merged together (RotationScaleMerge with SweepSegments), as XSCALE merges the data
+// sets XDS made of each sweep. Between the two the sweeps are put on one set of axes: the indexing each
+// was given is its own, and an ambiguous lattice or a different choice of axes indexes the same crystal
+// differently from one sweep to the next.
+
+// What the joint merge takes from one sweep's own run.
+struct JointMemberSweep {
+    std::string label;
+    int64_t n_frames = 0;
+    UnitCell cell{};                                 // as its reflections are written
+    gemmi::SpaceGroup space_group;                   // the group its run determined
+    std::vector<MergedReflection> merged;
+    std::vector<ScaledFull> fulls;                   // all its fulls, unmerged, on the same axes
+};
+
+// How each sweep is put on the axes of the reference - the sweep with the most fulls: the change of basis
+// (an axis permutation or a twin law of the lattice, ChooseReferenceSetting) under which its merged
+// intensities correlate best with the reference's, and that correlation. A sweep whose best correlation is
+// below MIN_REFERENCE_CC is not the same crystal form on any axes, and is left out.
+struct JointSweepRelation {
+    struct Sweep {
+        bool included = false;
+        gemmi::Op cob = gemmi::Op::identity();
+        double cc = 0.0;              // with the reference, under cob
+        double identity_cc = 0.0;     // ... and on its own axes
+        std::string why_left_out;
+    };
+    int reference = -1;
+    std::vector<Sweep> sweeps;
+};
+JointSweepRelation RelateJointSweeps(const std::vector<JointMemberSweep> &sweeps);
+
+// The joint merge's input: the frames of the included sweeps laid end to end, in the order given, each
+// holding its sweep's fulls on the reference's axes as pre-combined observations, and the segments that
+// say whose frames are whose.
+struct JointMergeInput {
+    std::vector<IntegrationOutcome> outcomes;
+    std::vector<SweepSegment> segments;
+    std::vector<int> segment_sweep;   // the sweep each segment is
+};
+JointMergeInput BuildJointMergeInput(const std::vector<JointMemberSweep> &sweeps,
+                                     const JointSweepRelation &relation);

rugnux/RugnuxScaleMerge.cpp

@@ -1693,6 +1693,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
             // Only a merge that may be written hands its fulls back: they are a copy as long as the
             // fulls, and the search merges and the pre-pass never write a .hkl.
             rsm->SetExportScaledFulls(!for_search && !geometry_prepass && write_output && config_.write_merged);
+            rsm->SetExportScaledFullsUncut(config_.sweep_member);
             // Amplitudes are made once, on the merge that is written, where the anisotropy tensor
             // their prior uses is known (see below) - not on every merge the tail makes on its way.
             rsm->SetFrenchWilson(false);

rugnux/RugnuxSweeps.cpp

@@ -0,0 +1,140 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RugnuxSweeps.h"
+
+#include <algorithm>
+#include <cctype>
+#include <cmath>
+#include <filesystem>
+#include <map>
+#include <set>
+#include <stdexcept>
+
+#include <spdlog/fmt/fmt.h>
+
+namespace {
+    // How far two sweeps may state the wavelength apart and still be one dataset: JFJochMultiReader's
+    // tolerance for the inputs of one stills run.
+    constexpr double WAVELENGTH_TOLERANCE = 1e-3;   // relative
+
+    std::string Lower(std::string s) {
+        std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });
+        return s;
+    }
+}
+
+std::string RugnuxFrameTemplate(const std::string &path) {
+    const std::filesystem::path p(path);
+    std::string name = p.filename().string();
+    size_t end = name.size();
+    while (end > 0 && !std::isdigit(static_cast<unsigned char>(name[end - 1])))
+        end--;
+    if (end == 0)
+        return {};
+    size_t start = end;
+    while (start > 0 && std::isdigit(static_cast<unsigned char>(name[start - 1])))
+        start--;
+    name.replace(start, end - start, std::string(end - start, '#'));
+    return (p.parent_path() / name).string();
+}
+
+RugnuxSweepPlan PlanRugnuxSweeps(const std::vector<std::string> &inputs,
+                                 const std::function<RugnuxInputDescription(const std::string &)> &describe_file,
+                                 const std::function<RugnuxInputDescription(const std::string &)> &describe_sweep) {
+    RugnuxSweepPlan plan;
+    if (inputs.size() < 2)
+        return plan;
+
+    std::vector<RugnuxInputDescription> files;
+    for (const auto &f: inputs)
+        files.push_back(describe_file(f));
+    const bool any_foreign = std::any_of(files.begin(), files.end(), [](const auto &d) { return d.foreign; });
+    const size_t n_frames = std::count_if(files.begin(), files.end(),
+                                          [](const auto &d) { return !d.frame_template.empty(); });
+    // Inputs of the stills readers, or a mix of frames and masters: read as before, which takes the
+    // former and refuses the latter.
+    if (any_foreign || (n_frames != 0 && n_frames != inputs.size()))
+        return plan;
+
+    if (n_frames == inputs.size()) {
+        std::set<std::string> seen;
+        for (size_t i = 0; i < inputs.size(); i++)
+            if (seen.insert(files[i].frame_template).second)
+                plan.sweeps.push_back({inputs[i], files[i].frame_template, true});
+        if (plan.sweeps.size() == 1)
+            return plan;
+    } else {
+        // A first input of stills is a stills run, read together as before - and refused there if any
+        // other input turns. A rotation first input asks every other one.
+        if (!describe_sweep(inputs[0]).rotation)
+            return plan;
+        std::vector<std::string> rotation{inputs[0]}, stills;
+        for (size_t i = 1; i < inputs.size(); i++)
+            (describe_sweep(inputs[i]).rotation ? rotation : stills).push_back(inputs[i]);
+        if (!stills.empty())
+            throw std::runtime_error(fmt::format(
+                    "{} is a rotation sweep and {} a stills collection: rotation sweeps and stills are not "
+                    "processed together; give them in separate runs", rotation.front(), stills.front()));
+        for (const auto &f: inputs)
+            plan.sweeps.push_back({f, f, false});
+    }
+
+    std::vector<double> wavelength;
+    for (const auto &s: plan.sweeps)
+        wavelength.push_back(describe_sweep(s.input).wavelength_A);
+    for (size_t k = 1; k < plan.sweeps.size(); k++)
+        if (std::fabs(wavelength[k] - wavelength[0]) > WAVELENGTH_TOLERANCE * wavelength[0])
+            throw std::runtime_error(fmt::format(
+                    "{} is at {:.5f} A and {} at {:.5f} A: sweeps at different wavelengths are not processed "
+                    "together; give them in separate runs", plan.sweeps[0].label, wavelength[0],
+                    plan.sweeps[k].label, wavelength[k]));
+    plan.wavelength_A = wavelength[0];
+    return plan;
+}
+
+std::vector<std::string> RugnuxSiblingSweeps(const RugnuxSweepPlan &plan, const std::vector<std::string> &extensions) {
+    namespace fs = std::filesystem;
+    std::set<std::string> given, dirs_frames, dirs_masters;
+    for (const auto &s: plan.sweeps) {
+        const fs::path p(s.input);
+        const std::string dir = p.parent_path().empty() ? "." : p.parent_path().string();
+        if (s.frames) {
+            given.insert(fs::path(s.label).filename().string());
+            dirs_frames.insert(dir);
+        } else {
+            given.insert(p.filename().string());
+            dirs_masters.insert(dir);
+        }
+    }
+    auto has_extension = [&](std::string name) {
+        name = Lower(name);
+        if (name.size() > 3 && name.ends_with(".gz"))
+            name.resize(name.size() - 3);
+        return std::any_of(extensions.begin(), extensions.end(), [&](const auto &e) { return name.ends_with(e); });
+    };
+
+    // template (or master) name -> files, in name order
+    std::map<std::string, size_t> found;
+    std::error_code ec;
+    for (const auto &dir: dirs_frames)
+        for (const auto &e: fs::directory_iterator(dir, ec)) {
+            const std::string name = e.path().filename().string();
+            if (!e.is_regular_file() || !has_extension(name))
+                continue;
+            const std::string t = fs::path(RugnuxFrameTemplate(e.path().string())).filename().string();
+            if (!t.empty() && !given.contains(t))
+                found[(fs::path(dir) / t).string()]++;
+        }
+    for (const auto &dir: dirs_masters)
+        for (const auto &e: fs::directory_iterator(dir, ec)) {
+            const std::string name = e.path().filename().string();
+            if (e.is_regular_file() && name.ends_with("_master.h5") && !given.contains(name))
+                found[(fs::path(dir) / name).string()] = 0;
+        }
+
+    std::vector<std::string> out;
+    for (const auto &[name, n]: found)
+        out.push_back(n > 0 ? fmt::format("{} ({} frame{})", name, n, n == 1 ? "" : "s") : name);
+    return out;
+}

rugnux/RugnuxSweeps.h

@@ -0,0 +1,56 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <functional>
+#include <string>
+#include <vector>
+
+// The rotation sweeps a rugnux command line names, decided from the files in the order given:
+//   - frame files (CBF, marCCD, SMV): one sweep per file-name template - the name with its last run of
+//     digits, the frame number, taken out. Any frame of a sweep names the whole sweep, so several frames
+//     of one template are one sweep;
+//   - masters: one sweep per rotation master. Several stills masters stay one stills run, as before, and
+//     a mix of rotation and stills masters is refused.
+// Sweeps at different wavelengths are refused. Nothing is planned for a single input, or for inputs only
+// the stills readers take (XFEL HDF5, CXI, XTC): those are read as they always were.
+
+// What the planner needs to know about one input. `frame_template` is empty for anything but a frame
+// file; `rotation` and `wavelength_A` are asked only of the input that stands for a sweep.
+struct RugnuxInputDescription {
+    std::string frame_template;
+    bool foreign = false;        // read only by the stills readers (XFEL HDF5, CXI, XTC)
+    bool rotation = false;
+    double wavelength_A = 0.0;
+};
+
+struct RugnuxSweep {
+    std::string input;           // the file the sweep is read from: its first frame given, or its master
+    std::string label;           // the frame template or the master, for messages
+    bool frames = false;         // a sweep of frame files (label is its template)
+};
+
+struct RugnuxSweepPlan {
+    // Empty: the inputs are read together as before (one input, a stills run, or a mix the readers refuse).
+    std::vector<RugnuxSweep> sweeps;
+    double wavelength_A = 0.0;   // of every sweep, where there are several
+};
+
+// `describe_file` says what an input is from its name and contents without reading the sweep (its
+// frame template, or whether it is only for the stills readers); `describe_sweep` reads the input that
+// stands for a sweep for its rotation and wavelength. Throws std::runtime_error, with a message for the
+// user, on a mix of rotation sweeps and stills, or on sweeps at different wavelengths.
+RugnuxSweepPlan PlanRugnuxSweeps(const std::vector<std::string> &inputs,
+                                 const std::function<RugnuxInputDescription(const std::string &)> &describe_file,
+                                 const std::function<RugnuxInputDescription(const std::string &)> &describe_sweep);
+
+// The template of a frame file name, with its directory: the last run of digits in the name replaced by
+// '#' characters, one per digit ("/d/x_1_0001.cbf" -> "/d/x_1_####.cbf"). Empty for a name with no digits.
+std::string RugnuxFrameTemplate(const std::string &path);
+
+// The sweeps beside the given ones that were not given: frame templates of `extensions` (lower case, with
+// the dot; a trailing ".gz" is looked through) in the directories of the frame sweeps, and other
+// "*_master.h5" files beside the given masters - one line each, for the log. Only the inputs given are
+// processed; these are named so that a sweep left out by accident shows.
+std::vector<std::string> RugnuxSiblingSweeps(const RugnuxSweepPlan &plan, const std::vector<std::string> &extensions);

rugnux/rugnux_cli.cpp

@@ -86,6 +86,8 @@
 #include "RugnuxOutputFiles.h"
 #include "RugnuxConsole.h"
 #include "RugnuxInputKind.h"
+#include "RugnuxSweeps.h"
+#include "RugnuxJointSweeps.h"
 #ifdef __APPLE__
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineMetal.h"
 #endif
@@ -101,6 +103,9 @@ constexpr double RELATIVE_B_DEFAULT_DEG = 10.0;
 // mutually exclusive by construction.
 enum class RugnuxMode { MX, Azint, Scale, Calibration, Refine };
 
+// --sweeps: several rotation sweeps as one dataset (joint), or each alone into <prefix>_s<k> (separate).
+enum class RugnuxSweepsMode { Joint, Separate };
+
 // Azimuthal sectors used by --mode calibration --calibration rings when --azim-phi-bins was not given.
 // The default of 1 is a plain radial profile, in which a ring has been averaged over every direction so
 // nothing is left to say where its centre is (RingsFromAzimuthalProfile refuses below 4). Sweeping the
@@ -120,6 +125,7 @@ never a directory; a model or reflection file among them is recognised by its co
 Rotation data (--mode mx, the default; chosen when the file's goniometer axis turns)
   rugnux -o myrun dataset_master.h5              NXmx / EIGER / Jungfraujoch master file
   rugnux -o myrun sweep_00001.cbf                CBF, marCCD or SMV sweep: name one frame
+  rugnux -o myrun a_master.h5 b_master.h5      several sweeps of one crystal, merged as one
   rugnux -o myrun dataset_master.h5 model.cif    with a model (or --model model.cif)
   rugnux -o myrun x_master.h5 http://host:5232/  follow a broker's collection while it is written
   A model is placed and refined after merging; it does not change the processing.
@@ -167,6 +173,8 @@ Basic options
   --no-gpu                       CPU only, even where a GPU is present
   --no-merge                     No scaling or merging; write the per-image _process.h5
   --write-process-h5             Also write the (large) _process.h5 when merging
+  --sweeps separate              Several rotation sweeps (masters or frame templates): each alone,
+                                 into <prefix>_s<k>_*, instead of merged as one dataset
   -v, --verbose                  Whole log on the terminal (it is always in <prefix>_log.txt)
 
 All options (about 120, for unusual data and diagnosis), each with what it does:
@@ -280,7 +288,8 @@ enum {
     OPT_POLARIZATION,
     OPT_NO_REFINE_TILT,
     OPT_REFINE_SCORE_ONLY,
-    OPT_REFINE_D_MIN
+    OPT_REFINE_D_MIN,
+    OPT_SWEEPS
 };
 
 static option long_options[] = {
@@ -304,6 +313,7 @@ static option long_options[] = {
     {"no-refine-solvent", no_argument, nullptr, OPT_NO_REFINE_SOLVENT},
     {"refine-score-only", no_argument, nullptr, OPT_REFINE_SCORE_ONLY},
     {"refine-d-min", required_argument, nullptr, OPT_REFINE_D_MIN},
+    {"sweeps", required_argument, nullptr, OPT_SWEEPS},
     {"export-unmerged", no_argument, nullptr, OPT_EXPORT_UNMERGED},
     {"no-export-unmerged", no_argument, nullptr, OPT_NO_EXPORT_UNMERGED},
     {"export-unmerged-partials", no_argument, nullptr, OPT_EXPORT_UNMERGED_PARTIALS},
@@ -805,70 +815,27 @@ extern "C" void fftwf_assertion_failed(const char *s, int line, const char *file
     std::abort();
 }
 
-// The body of main. Settings setters and the pipeline itself throw JFJochException on input the parser
-// cannot reject on its own (a resolution limit of zero, a distance of zero, a polarisation above one,
-// an unreadable file), so main wraps this and reports rather than letting the exception terminate the
-// process with no diagnostic and exit code 134.
-static int RunRugnux(int argc, char **argv, bool no_gpu) {
-    // cuFFT plans are reused through the run and destroyed with it, ahead of the device reset in main.
-    FFTPlanReuse fft_plans;
-    // The whole invocation, reported at the end and written into the results report. Started before
-    // anything else so it covers reading the file and setting up as well as the processing.
-    const auto invocation_start = std::chrono::steady_clock::now();
-    RunProvenance provenance;
-    provenance.command_line = JoinCommandLine(argc, argv);
-    provenance.gpu_count = get_gpu_count();
-#ifdef __APPLE__
-    // Without CUDA, the Bragg integration's integer pixel work runs on the Apple GPU (MXAnalysisWithoutFPGA).
-    const std::string metal_device = BraggIntegrationEngineMetal::DeviceName();
-#else
-    const std::string metal_device;
-#endif
-    provenance.gpu_description = get_gpu_description();
+namespace {
+    // One input's run ended: what exit(EXIT_FAILURE) was where a command line ran one input only. A run of
+    // several sweeps goes on with the next sweep; RunRugnux returns the code.
+    struct RugnuxInputFailed { int code; };
+    // What Parse and PrepareRugnuxInput return when the run goes on.
+    constexpr int RUGNUX_CONTINUE = -1;
+}
 
-    // The terminal shows the run as steps (RugnuxConsole) and every Logger writes to the log file
-    // instead; with -v the whole log is on the terminal as well. -v is looked for here, before the
-    // options are parsed, because the parser already logs.
+// The command line as parsed, with what RunRugnux set up before it was parsed. The options are members so
+// that the code reading and running one input (RugnuxInputRun) reads them by name, as it did when it was
+// one function, and so that each sweep of a run of several is run from a copy of its own.
+struct RugnuxCli {
+    bool no_gpu;
     bool verbose = false;
-    for (int i = 1; i < argc; i++)
-        if (std::strcmp(argv[i], "-v") == 0 || std::strcmp(argv[i], "--verbose") == 0)
-            verbose = true;
-    console = std::make_shared<RugnuxConsole>(verbose);
-    std::atexit(CloseConsoleOnExit);
-    // Files with an older time than this are not this run's (the margin covers coarse file times).
-    const auto files_since = std::filesystem::file_time_type::clock::now() - std::chrono::seconds(2);
-    auto log_sinks = std::make_shared<spdlog::sinks::dist_sink_mt>();
-    // What is logged before the log file can be opened - its name comes from -o - waits here.
-    auto early_log = std::make_shared<spdlog::sinks::ringbuffer_sink_mt>(10000);
-    if (verbose)
-        log_sinks->add_sink(std::make_shared<spdlog::sinks::stdout_color_sink_mt>());
-    else
-        log_sinks->add_sink(MakeConsoleSink(console));
-    log_sinks->add_sink(early_log);
-    Logger::RouteAll(log_sinks);
-    // The HDF5 library prints its own error stack on stderr; every failure it reports reaches the
-    // reader as an error code anyway, which is turned into a message there.
-    H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr);
-
-    PrintBanner(*console);
-
-    Logger logger("rugnux");
-#ifdef JFJOCH_USE_CUDA
-    const char *build = "CUDA build";
-#else
-    const char *build = "CPU build";
-#endif
-    logger.Info("rugnux Version {}{} - {}", jfjoch_version(),
-                jfjoch_git_sha1().empty() ? std::string()
-                : fmt::format(" (git {} {})", jfjoch_git_sha1().substr(0, 6), jfjoch_git_date()), build);
-    logger.Info("Build flags: {}", jfjoch_build_cxx_flags().empty() ? "none" : jfjoch_build_cxx_flags());
-    logger.Info("Build: {}", jfjoch_build_info(true));
-    logger.Info("Command line: {}", provenance.command_line);
-    logger.Info("GPU: {}", no_gpu ? std::string("not used (--no-gpu) - running on the CPU")
-            : provenance.gpu_count > 0 ? fmt::format("{} ({} visible)", provenance.gpu_description, provenance.gpu_count)
-            : !metal_device.empty() ? fmt::format("{} through Metal, for the Bragg integration only - the rest runs on the CPU", metal_device)
-            : std::string("none visible - running on the CPU"));
-    RegisterHDF5Filter();
+    std::chrono::steady_clock::time_point invocation_start;
+    RunProvenance provenance;
+    std::string metal_device;
+    std::filesystem::file_time_type files_since;
+    std::shared_ptr<spdlog::sinks::dist_sink_mt> log_sinks;
+    std::shared_ptr<spdlog::sinks::ringbuffer_sink_mt> early_log;
+    Logger logger{"rugnux"};   // made after RunRugnux has routed every Logger to the log
 
     std::string input_file;
     std::string output_prefix = "output";
@@ -991,6 +958,79 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     bool simple_stills_flag = false;                   // --simple-stills: disable the default stills partiality post-refinement
     std::optional<double> outlier_reject_nsigma;       // merge per-observation outlier rejection
 
+    std::vector<std::string> input_files;
+    std::string broker_url;
+    std::shared_ptr<const RugnuxModel> model;
+    RugnuxSweepsMode sweeps_mode = RugnuxSweepsMode::Joint;   // --sweeps
+    // A log file for the whole run is open already (--sweeps joint): an input's own is not opened over it.
+    bool log_open = false;
+
+    RugnuxCli(bool no_gpu_, std::chrono::steady_clock::time_point invocation_start_, RunProvenance provenance_,
+              std::string metal_device_, std::filesystem::file_time_type files_since_,
+              std::shared_ptr<spdlog::sinks::dist_sink_mt> log_sinks_,
+              std::shared_ptr<spdlog::sinks::ringbuffer_sink_mt> early_log_)
+        : no_gpu(no_gpu_), invocation_start(invocation_start_), provenance(std::move(provenance_)),
+          metal_device(std::move(metal_device_)), files_since(files_since_), log_sinks(std::move(log_sinks_)),
+          early_log(std::move(early_log_)) {}
+
+    // The command line: RUGNUX_CONTINUE, or the exit code of a run that ends here.
+    int Parse(int argc, char **argv);
+    // --sweeps joint: the sweeps one dataset (RugnuxJointSweeps.h).
+    int RunJointSweeps(const RugnuxSweepPlan &plan);
+    // Each sweep of the command line, or its one input.
+    int Run();
+};
+
+// One input - one sweep, or one stills run - read, set up and processed. Copied from the command line it
+// belongs to, with its own input and output prefix.
+struct RugnuxInputRun : RugnuxCli {
+    explicit RugnuxInputRun(const RugnuxCli &cli) : RugnuxCli(cli) {}
+
+    // What PrepareRugnuxInput sets up and RunPreparedInput runs with.
+    std::future<void> cufft_warmup, context_warmup;
+    SimulatedLive<JFJochHDF5Reader> hdf5_reader;
+    SimulatedLive<JFJochCBFReader> cbf_reader;
+    SimulatedLive<JFJochMarCCDReader> marccd_reader;
+    SimulatedLive<JFJochSMVReader> smv_reader;
+    SimulatedLive<JFJochXFELReader> xfel_reader;
+    SimulatedLive<JFJochXTCReader> xtc_reader;
+    SimulatedLive<JFJochSFRMReader> sfrm_reader;
+    JFJochStartMessageReader live_reader;
+    std::unique_ptr<JFJochMultiReader> multi_reader;
+    std::unique_ptr<BrokerFeed> broker_feed;
+    JFJochReader *reader_ptr = nullptr;
+    std::shared_ptr<const JFJochReaderDataset> dataset;
+    OutputProvenance output_provenance;
+    DiffractionExperiment experiment;
+    ProcessConfig config;
+    std::optional<ProcessFileSink> sink;
+    // The log file of this input, taken out of the log again once it is done.
+    std::shared_ptr<spdlog::sinks::sink> log_file_sink;
+    // One sweep of a joint run (--sweeps joint): processed as on its own, its result kept here and nothing
+    // written.
+    bool sweep_member = false;
+    std::optional<ProcessResult> member_result;
+    // The joint merge of the sweeps' fulls, made in place of this input's own run (Rugnux::RunJoint).
+    struct JointRequest {
+        std::vector<IntegrationOutcome> outcomes;
+        UnitCell cell{};
+        gemmi::CrystalSystem system = gemmi::CrystalSystem::Triclinic;
+        char centring = 'P';
+    };
+    std::optional<JointRequest> joint_request;
+
+    // Reads the input and sets up the experiment and the run's configuration: RUGNUX_CONTINUE, or the
+    // exit code of a mode that ends here (refine, azint, calibration, scale).
+    int PrepareRugnuxInput();
+    // Runs the input PrepareRugnuxInput set up and writes its files.
+    int RunPreparedInput();
+    int Process() {
+        const int rc = PrepareRugnuxInput();
+        return rc != RUGNUX_CONTINUE ? rc : RunPreparedInput();
+    }
+};
+
+int RugnuxCli::Parse(int argc, char **argv) {
     if (argc == 1) {
         print_usage();
         exit(EXIT_FAILURE);
@@ -1198,6 +1238,16 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         case OPT_REFINE_D_MIN:
             refine_d_min = parse_number_arg<double>(optarg, "--refine-d-min", logger, 0.1, 100.0);
             break;
+        case OPT_SWEEPS:
+            if (std::strcmp(optarg, "joint") == 0)
+                sweeps_mode = RugnuxSweepsMode::Joint;
+            else if (std::strcmp(optarg, "separate") == 0)
+                sweeps_mode = RugnuxSweepsMode::Separate;
+            else {
+                logger.Error("--sweeps expects joint or separate, not {}", optarg);
+                exit(EXIT_FAILURE);
+            }
+            break;
         case OPT_EXPORT_UNMERGED:
             export_unmerged = true;
             break;
@@ -1625,8 +1675,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
 
     // The inputs, then a broker URL if the last word is one: told apart by the scheme, so that several
     // stills inputs can be given.
-    std::vector<std::string> input_files(argv + optind, argv + argc);
-    std::string broker_url;
+    input_files.assign(argv + optind, argv + argc);
+    broker_url.clear();
     if (!input_files.empty() && (input_files.back().starts_with("http://")
                                  || input_files.back().starts_with("https://"))) {
         broker_url = input_files.back();
@@ -1696,7 +1746,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         exit(EXIT_FAILURE);
     }
     // Read before anything runs: a file that holds no model ends the run here rather than after it.
-    std::shared_ptr<const RugnuxModel> model;
+    model.reset();
     if (const std::string path = mr_model.empty() ? model_pdb : mr_model; !path.empty()) {
         try {
             model = std::make_shared<const RugnuxModel>(ReadRugnuxModelFile(path));
@@ -1712,6 +1762,10 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         exit(EXIT_FAILURE);
     }
 
+    return RUGNUX_CONTINUE;
+}
+
+int RugnuxInputRun::PrepareRugnuxInput() {
     input_file = input_files[0];
 
     // -N defaults to 0 = "use every CPU the run may use" (AvailableCpus); resolve it to a concrete count here so every mode
@@ -1729,13 +1783,14 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
 
     // The log file, next to the outputs, and everything logged so far into it.
     std::string &log_file = g_log_file;
-    if (!output_prefix.empty()) {
+    if (!output_prefix.empty() && !log_open) {
         log_file = output_prefix + "_log.txt";
         try {
             auto file = std::make_shared<spdlog::sinks::basic_file_sink_mt>(log_file, /*truncate=*/true);
             for (const auto &m : early_log->last_raw())
                 file->log(m);
             log_sinks->add_sink(file);
+            log_file_sink = file;
         } catch (const std::exception &e) {
             logger.Warning("The log file {} cannot be written: {}", log_file, e.what());
             log_file.clear();
@@ -1768,7 +1823,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             data = ReadRugnuxReflectionFile(input_file);
         } catch (const std::exception &e) {
             logger.Error("Cannot read the reflections {}: {}", input_file, e.what());
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
         console->OnStep(refine_score_only ? "Scoring the model" : "Refining the model");
         RefinementSettings settings;
@@ -1815,27 +1870,16 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     // Off the critical path, so cuFFT's one-time initialisation overlaps the file open, the mask read
     // and the beam-stop pre-scan instead of landing inside the first pass. It also times a plan, which
     // shows whether the CUDA JIT cache works.
-    auto cufft_warmup = std::async(std::launch::async, [&logger] { LogCudaJitCache(logger, WarmUpCuFFT()); });
+    cufft_warmup = std::async(std::launch::async, [this] { LogCudaJitCache(logger, WarmUpCuFFT()); });
     // The same for the CUDA contexts of the other cards (the warm-up above creates card 0's), which
     // the pre-scan does not use: created one after another at start-up, they cost a fifth of a second
     // per card before the run began.
-    auto context_warmup = std::async(std::launch::async, create_gpu_contexts, 1);
+    context_warmup = std::async(std::launch::async, create_gpu_contexts, 1);
 
     // 1. Read Input File. An HDF5 master, one PILATUS miniCBF frame, or one marCCD frame - naming
     // any frame of a sweep reads the whole sweep, which is how these are archived and how a user
     // thinks of them.
     console->OnStep("Reading input");
-    SimulatedLive<JFJochHDF5Reader> hdf5_reader;
-    SimulatedLive<JFJochCBFReader> cbf_reader;
-    SimulatedLive<JFJochMarCCDReader> marccd_reader;
-    SimulatedLive<JFJochSMVReader> smv_reader;
-    SimulatedLive<JFJochXFELReader> xfel_reader;
-    SimulatedLive<JFJochXTCReader> xtc_reader;
-    SimulatedLive<JFJochSFRMReader> sfrm_reader;
-    JFJochStartMessageReader live_reader;
-    std::unique_ptr<JFJochMultiReader> multi_reader;
-    std::unique_ptr<BrokerFeed> broker_feed;
-    JFJochReader *reader_ptr = nullptr;
     // Which format this is. Asked INSIDE the try: these open and read the file, and an input the program
     // cannot read has to come out as the message below, not as an unhandled exception with nothing
     // printed at all.
@@ -1843,13 +1887,13 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     bool input_is_cxi = false, input_is_xtc = false, input_is_panel_table = false, input_is_sfrm = false;
     try {
         // Every input is the format the first one is.
-        const auto all_inputs = [&input_files](auto can_read) {
+        const auto all_inputs = [this](auto can_read) {
             return std::all_of(input_files.begin(), input_files.end(), can_read);
         };
         const auto refuse_mixed = [&]() {
             logger.Error("The inputs are not all of one format; several inputs must be files of the same "
                          "kind - Jungfraujoch/NXmx masters, SACLA HDF5, CXI or XTC");
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         };
         const auto is_foreign = [](const std::string &f) {
             return JFJochXFELReader::CanReadPanelTable(f) || JFJochXFELReader::CanReadCXI(f) || JFJochXTCReader::CanRead(f);
@@ -1861,7 +1905,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 || JFJochSMVReader::CanRead(input_file) || JFJochSFRMReader::CanRead(input_file)) {
                 logger.Error("{} is one frame of a sweep, which names the whole sweep: give just one of its "
                              "frames", input_file);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
             if (!all_inputs([&is_foreign](const std::string &f) {
                     return !is_foreign(f) && !JFJochCBFReader::CanRead(f) && !JFJochMarCCDReader::CanRead(f)
@@ -1874,7 +1918,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                     r->ReadFile(f);
                 } catch (const std::exception &e) {
                     logger.Error("{} is not a Jungfraujoch/NXmx master that can be read: {}", f, e.what());
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 inputs.push_back(std::move(r));
             }
@@ -1893,7 +1937,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             if (!extra_masks.empty() && !input_is_crystfel && !input_is_panel_table && !input_is_cxi) {
                 logger.Error("--mask-good and --mask-bad are for HDF5/CXI stills (SACLA, CXI or --crystfel-geometry "
                              "input); {} is not one", input_file);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
             xfel_reader.ExtraMasks(extra_masks);
             if (input_is_crystfel) {
@@ -1901,7 +1945,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 reader_ptr = &xfel_reader;
                 if (!xfel_reader.HasPhotonEnergy() && !wavelength_A) {
                     logger.Error("Nothing read states the photon energy of {}; give it with --wavelength", input_file);
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 logger.Info("Read {} frames from {} files with the CrystFEL geometry {}",
                             xfel_reader.GetNumberOfImages(), input_files.size(), crystfel_geometry);
@@ -1915,7 +1959,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 logger.Info("Read {} frames from {} CXI files", xfel_reader.GetNumberOfImages(), input_files.size());
             } else if (const auto refusal = JFJochXFELReader::CXIRefusal(input_file)) {
                 logger.Error("{}", *refusal);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
             input_is_cbf = !reader_ptr && JFJochCBFReader::CanRead(input_file);
             input_is_marccd = !reader_ptr && !input_is_cbf && JFJochMarCCDReader::CanRead(input_file);
@@ -1930,7 +1974,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 && !input_is_cbf) {
                 logger.Error("--gain is for input in detector units (SACLA, CXI, XTC, SMV, marCCD, SFRM, a CCD imgCIF or "
                              "--crystfel-geometry input); {} is read in the units it states", input_file);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
             if (reader_ptr) {
                 // read above
@@ -1972,7 +2016,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                     logger.Error("The XTC does not state the detector distance{}{}: give {} (the beam centre is "
                                  "measured; with -C the distance is too)", xtc_reader.HasPixelSize() ? "" : ", the pixel size",
                                  xtc_reader.HasWavelength() ? "" : ", the wavelength", options);
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 reader_ptr = &xtc_reader;
                 logger.Info("Read {} XTC still images from {} streams", xtc_reader.GetNumberOfImages(),
@@ -1981,13 +2025,13 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 // The collection is still being written: follow it from the broker.
                 if (mode != RugnuxMode::MX) {
                     logger.Error("Following a collection from the broker is for --mode mx only");
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 const char *token = std::getenv("JUNGFRAUJOCH_HTTP_TOKEN");
                 if (token == nullptr) {
                     logger.Error("{} does not exist yet, and following it from the broker needs the collection's token "
                                  "in the environment variable JUNGFRAUJOCH_HTTP_TOKEN", input_file);
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 broker_feed = std::make_unique<BrokerFeed>(broker_url, token);
                 // The broker's event stream starts with the start message: at once on a JUNGFRAU, when the
@@ -1999,11 +2043,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                 const std::string given = std::filesystem::path(input_file).lexically_normal().generic_string();
                 if (given != master && !given.ends_with("/" + master)) {
                     logger.Error("The broker at {} is writing {}, not {}", broker_feed->Host(), master, input_file);
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 if (!start.goniometer.has_value() || !start.goniometer->IsScanning()) {
                     logger.Error("Following a collection from the broker is for rotation data, and {} is not", master);
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
                 live_reader.Open(start, input_file);
                 broker_feed->Attach(live_reader);
@@ -2018,7 +2062,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         }
     } catch (const std::exception &e) {
         logger.Error("Error reading input file: {}", e.what());
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
     JFJochReader &reader = *reader_ptr;
 
@@ -2085,20 +2129,19 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                     ? fmt::format(", {:.3f} of them to a photon; intensities and sigmas in photons", adu_per_photon)
                     : std::string(", one to a photon"));
 
-    const auto dataset = reader.GetDataset();
+    dataset = reader.GetDataset();
     if (!dataset) {
         logger.Error("No experiment dataset found in the input file");
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
 
     if (rotation_indexing_image_count <= 0) {
         logger.Error("Invalid number of rotation indexing images: {}", rotation_indexing_image_count);
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
 
     logger.Info("Loaded dataset from {}", input_file);
     // Recorded in the reflection files.
-    OutputProvenance output_provenance;
     output_provenance.command_line = provenance.command_line;
     DescribeInput(input_file, output_provenance);
     output_provenance.collection_date = dataset->arm_date;
@@ -2192,7 +2235,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             }
         } catch (const std::exception &e) {
             logger.Error("Error reading reference {}: {}", ref_mtz, e.what());
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
     }
 
@@ -2203,7 +2246,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             const auto flags = LoadReferenceMtz(free_flags_file, std::nullopt);
             if (!flags.has_free_flags) {
                 logger.Error("{} carries no R-free flags (no FreeR_flag-like column, no _refln.status)", free_flags_file);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
             free_flags_data = flags.reflections;
             logger.Info("R-free flags from {} (column {}): {} of {} free, taken by the merge where these data match "
@@ -2211,7 +2254,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                         flags.reflections.size());
         } catch (const std::exception &e) {
             logger.Error("Error reading {}: {}", free_flags_file, e.what());
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
     }
 
@@ -2231,7 +2274,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             || input_is_xtc) {
             logger.Error("--mode scale needs the integrated reflections in a _process.h5; "
                          "a raw image sweep holds none. Run --mode mx on it first.");
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
         const auto total_images = static_cast<int>(reader.GetNumberOfImages());
         const int last_image = (end_image < 0 || end_image >= total_images) ? total_images - 1 : end_image;
@@ -2300,7 +2343,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
 
         if (!experiment.GetUnitCell()) {
             logger.Error("Experiment unit cell not found, cannot update reflection resolution");
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
         // The stored reflections have to be in the setting of the cell they are about to be merged in,
         // or the merge's per-frame cell gate rejects every observation and the run dies far downstream
@@ -2342,7 +2385,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                                         "written before that was recorded)",
                              file_cell->a, file_cell->b, file_cell->c, file_cell->alpha,
                              file_cell->beta, file_cell->gamma);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
         }
 
@@ -2366,7 +2409,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                              "has, or reprocess the images with --mode mx.",
                              merge_sg->xhm(), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma,
                              violation, MAX_METRIC_VIOLATION);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
         }
 
@@ -2804,12 +2847,12 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
 
     if (image_stride < 0) {
         logger.Error("Image stride cannot be negative");
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
 
     if (image_stride == 0) {
         logger.Error("Image stride cannot be zero");
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
 
     // Round UP: the worker loop walks ordinals until start + ordinal * stride reaches end_image, so a
@@ -2824,7 +2867,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     }
 
     // 2. Setup Experiment & Components
-    DiffractionExperiment experiment(dataset->experiment);
+    experiment = DiffractionExperiment(dataset->experiment);
 
     // Geometry overrides (default: keep the value stored in the input file). Applied before the
     // azimuthal-integration settings are derived, which depend on the geometry.
@@ -2860,20 +2903,20 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             }
             if (modules.empty()) {
                 logger.Error("--module-geometry {} describes fewer than two modules", module_geometry_file);
-                exit(EXIT_FAILURE);
+                throw RugnuxInputFailed{EXIT_FAILURE};
             }
             for (const auto &m: modules.Get())
                 if (m.x0 < 0 || m.y0 < 0 || m.x0 + m.width > experiment.GetXPixelsNum()
                     || m.y0 + m.height > experiment.GetYPixelsNum()) {
                     logger.Error("--module-geometry {}: a module lies outside the {} x {} image", module_geometry_file,
                                  experiment.GetXPixelsNum(), experiment.GetYPixelsNum());
-                    exit(EXIT_FAILURE);
+                    throw RugnuxInputFailed{EXIT_FAILURE};
                 }
             experiment.Modules(modules);
             logger.Info("--module-geometry: {} modules from {}", modules.size(), module_geometry_file);
         } catch (const std::exception &e) {
             logger.Error("Cannot read --module-geometry {}: {}", module_geometry_file, e.what());
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
     }
     if (!experiment.GetModules().empty())
@@ -2942,7 +2985,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             result = process.Run(console.get());
         } catch (const std::exception &e) {
             logger.Error("Processing failed: {}", FailureText(e));
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
         g_active_process = nullptr;
 
@@ -3228,7 +3271,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     if (force_still) {
         if (rotation_indexing) {
             logger.Error("--force-still conflicts with -R / --single-pass-rotation / --force-rotation-lattice");
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
         if (has_goniometer)
             logger.Info("--force-still: treating the rotation dataset as independent stills");
@@ -3250,7 +3293,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                      "measured over frames the sweep never recorded. Process the sweep whole, or select "
                      "a contiguous range with -s/-e. --mode azint and --force-still still take a stride.",
                      image_stride);
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
 
     // Scaling and merging are on by default (run_scaling initialised true); --no-merge turns them off
@@ -3510,7 +3553,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
 
     // Run the shared full-analysis workflow (rotation indexing + scaling/merging live in
     // Rugnux; the experiment above carries all algorithm settings).
-    ProcessConfig config;
     config.mode = ProcessMode::FullAnalysis;
     config.provenance = output_provenance;
     config.start_image = start_image;
@@ -3563,7 +3605,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     config.refine_solvent = refine_solvent;
     // When merging, the merged reflections (.mtz/.cif) are the wanted output; skip the large
     // _process.h5 unless explicitly requested. Without merging, the _process.h5 is the only output.
-    std::optional<ProcessFileSink> sink;
     if (!output_prefix.empty() && (run_scaling ? write_process_h5_flag : true))
         config.image_sink = &sink.emplace(output_prefix);
     config.finalist_ledger = finalist_ledger_flag;
@@ -3579,6 +3620,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         WarmUpModelValidation(model, static_cast<size_t>(config.nthreads), logger);
     simulated_start = std::chrono::steady_clock::now();
 
+    return RUGNUX_CONTINUE;
+}
+
+int RugnuxInputRun::RunPreparedInput() {
+    JFJochReader &reader = *reader_ptr;
     Rugnux process(reader, experiment, *dataset->pixel_mask, config);
 
     g_active_process = &process;
@@ -3586,7 +3632,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
 
     ProcessResult result;
     try {
-        result = process.Run(console.get());
+        result = joint_request ? process.RunJoint(std::move(joint_request->outcomes), joint_request->cell,
+                                                  joint_request->system, joint_request->centring, console.get())
+                               : process.Run(console.get());
     } catch (const NoLatticeFound &e) {
         // No crystal in these frames is a result, and the report says so - VERDICT= FAILED with
         // NO_LATTICE - before the run exits as failed.
@@ -3617,13 +3665,17 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                                                      provenance);
         console->Fail();
         PrintSummary(*console, doc, output_prefix, files_since, provenance.wall_time_s);
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     } catch (const std::exception &e) {
         logger.Error("Processing failed: {}", FailureText(e));
         sink.reset();
-        exit(EXIT_FAILURE);
+        throw RugnuxInputFailed{EXIT_FAILURE};
     }
     g_active_process = nullptr;
+    if (sweep_member) {
+        member_result = std::move(result);
+        return EXIT_SUCCESS;
+    }
 
     // Where the modules are, for a later run's --module-geometry (and any NXmx reader).
     if (!output_prefix.empty() && result.modules_placed) {
@@ -3641,7 +3693,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                                logger);
         } catch (const std::exception &e) {
             logger.Error("Processing failed: {}", FailureText(e));
-            exit(EXIT_FAILURE);
+            throw RugnuxInputFailed{EXIT_FAILURE};
         }
     }
 
@@ -3785,6 +3837,304 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     return EXIT_SUCCESS;
 }
 
+int RugnuxCli::Run() {
+    // The sweeps the inputs are (RugnuxSweeps.h). One input, a stills run and anything the readers refuse
+    // are read as they were.
+    const auto describe_file = [this](const std::string &f) {
+        RugnuxInputDescription d;
+        if (!crystfel_geometry.empty())
+            d.foreign = true;
+        else if (JFJochCBFReader::CanRead(f) || JFJochMarCCDReader::CanRead(f) || JFJochSMVReader::CanRead(f))
+            d.frame_template = RugnuxFrameTemplate(f);
+        else
+            d.foreign = JFJochXFELReader::CanReadPanelTable(f) || JFJochXFELReader::CanReadCXI(f)
+                        || JFJochXTCReader::CanRead(f);
+        return d;
+    };
+    // What a sweep is read as, for its rotation and wavelength: the reader the run itself would take.
+    const auto read_dataset = [](const std::string &f) -> std::shared_ptr<const JFJochReaderDataset> {
+        if (JFJochCBFReader::CanRead(f)) {
+            JFJochCBFReader r;
+            r.ReadFiles(f);
+            return r.GetDataset();
+        }
+        if (JFJochMarCCDReader::CanRead(f)) {
+            JFJochMarCCDReader r;
+            r.ReadFiles(f);
+            return r.GetDataset();
+        }
+        if (JFJochSMVReader::CanRead(f)) {
+            JFJochSMVReader r;
+            r.ReadFiles(f);
+            return r.GetDataset();
+        }
+        JFJochHDF5Reader r;
+        r.ReadFile(f);
+        return r.GetDataset();
+    };
+    const auto describe_sweep = [&read_dataset](const std::string &f) {
+        std::shared_ptr<const JFJochReaderDataset> dataset;
+        try {
+            dataset = read_dataset(f);
+        } catch (const std::exception &e) {
+            throw std::runtime_error(fmt::format("{} cannot be read: {}", f, e.what()));
+        }
+        RugnuxInputDescription d;
+        if (dataset) {
+            const auto g = dataset->experiment.GetGoniometer();
+            d.rotation = g.has_value() && g->GetIncrement_deg() != 0.0f;
+            d.wavelength_A = dataset->experiment.GetWavelength_A();
+        }
+        return d;
+    };
+    RugnuxSweepPlan plan;
+    try {
+        plan = PlanRugnuxSweeps(input_files, describe_file, describe_sweep);
+    } catch (const std::exception &e) {
+        logger.Error("{}", e.what());
+        return EXIT_FAILURE;
+    }
+
+    if (plan.sweeps.size() < 2) {
+        // Frames of one sweep: any of them names it.
+        if (plan.sweeps.size() == 1)
+            input_files = {plan.sweeps[0].input};
+        try {
+            return RugnuxInputRun(*this).Process();
+        } catch (const RugnuxInputFailed &f) {
+            return f.code;
+        }
+    }
+
+    for (size_t k = 0; k < plan.sweeps.size(); k++)
+        logger.Info("Sweep {}: {}", k + 1, plan.sweeps[k].label);
+    for (const auto &s: RugnuxSiblingSweeps(plan, {".cbf", ".img", ".mccd", ".osc"}))
+        logger.Info("Also beside these inputs, not given and so not processed: {}", s);
+    if (start_image != 0 || end_image >= 0) {
+        logger.Error("-s / -e select images of one sweep; with several sweeps give each sweep's frames instead");
+        return EXIT_FAILURE;
+    }
+    if (force_still) {
+        logger.Error("--force-still processes one rotation sweep as stills; give the sweeps one at a time");
+        return EXIT_FAILURE;
+    }
+    if (sweeps_mode == RugnuxSweepsMode::Joint)
+        return RunJointSweeps(plan);
+
+    // --sweeps separate: each sweep alone, as if it had been the one input, into <prefix>_s<k>.
+    int failed = 0;
+    for (size_t k = 0; k < plan.sweeps.size(); k++) {
+        RugnuxInputRun sweep(*this);
+        sweep.input_files = {plan.sweeps[k].input};
+        // Each sweep's report times that sweep; the first also the start-up before it.
+        if (k > 0)
+            sweep.invocation_start = std::chrono::steady_clock::now();
+        if (!output_prefix.empty())
+            sweep.output_prefix = output_prefix + "_s" + std::to_string(k + 1);
+        logger.Info("Sweep {} of {}: {} into {}_*", k + 1, plan.sweeps.size(), plan.sweeps[k].label,
+                    sweep.output_prefix);
+        int rc;
+        try {
+            rc = sweep.Process();
+        } catch (const RugnuxInputFailed &f) {
+            rc = f.code;
+        }
+        if (sweep.log_file_sink)
+            log_sinks->remove_sink(sweep.log_file_sink);
+        if (rc != EXIT_SUCCESS) {
+            logger.Error("Sweep {} ({}) failed", k + 1, plan.sweeps[k].label);
+            failed++;
+        }
+    }
+    if (failed > 0)
+        logger.Error("{} of {} sweeps failed", failed, plan.sweeps.size());
+    return failed > 0 ? EXIT_FAILURE : EXIT_SUCCESS;
+}
+
+
+int RugnuxCli::RunJointSweeps(const RugnuxSweepPlan &plan) {
+    // The whole run's log, from the start: each sweep writes its own beside it, and the joint merge goes on
+    // in this one.
+    if (!output_prefix.empty()) {
+        g_log_file = output_prefix + "_log.txt";
+        try {
+            auto file = std::make_shared<spdlog::sinks::basic_file_sink_mt>(g_log_file, /*truncate=*/true);
+            for (const auto &m : early_log->last_raw())
+                file->log(m);
+            log_sinks->add_sink(file);
+            log_open = true;
+        } catch (const std::exception &e) {
+            logger.Warning("The log file {} cannot be written: {}", g_log_file, e.what());
+            g_log_file.clear();
+        }
+    }
+    // Each sweep processed on its own first, writing nothing, its fulls kept: the first level.
+    std::vector<JointMemberSweep> members;
+    std::vector<std::string> inputs;
+    for (size_t k = 0; k < plan.sweeps.size(); k++) {
+        RugnuxInputRun sweep(*this);
+        sweep.input_files = {plan.sweeps[k].input};
+        sweep.log_open = false;
+        // A sweep of a joint run writes its log, and nothing else (sweep_member).
+        if (!output_prefix.empty())
+            sweep.output_prefix = output_prefix + "_s" + std::to_string(k + 1);
+        sweep.sweep_member = true;
+        logger.Info("Sweep {} of {}: {}, processed on its own", k + 1, plan.sweeps.size(), plan.sweeps[k].label);
+        int rc;
+        try {
+            rc = sweep.PrepareRugnuxInput();
+            if (rc == RUGNUX_CONTINUE) {
+                sweep.config.sweep_member = true;
+                sweep.config.model.reset();
+                sweep.config.force_mr = false;
+                sweep.config.image_sink = nullptr;
+                rc = sweep.RunPreparedInput();
+            }
+        } catch (const RugnuxInputFailed &f) {
+            rc = f.code;
+        }
+        if (sweep.log_file_sink)
+            log_sinks->remove_sink(sweep.log_file_sink);
+        JointMemberSweep m;
+        m.label = plan.sweeps[k].label;
+        if (rc == EXIT_SUCCESS && sweep.member_result && sweep.member_result->outputs
+            && sweep.member_result->outputs->merged) {
+            const MergedOutput &merged = *sweep.member_result->outputs->merged;
+            m.n_frames = static_cast<int64_t>(sweep.member_result->images_processed);
+            m.cell = merged.cell;
+            m.space_group = merged.experiment.GetSpaceGroupOrP1();
+            m.merged = merged.reflections;
+            m.fulls = merged.scaled_fulls;
+            logger.Info("Sweep {}: {} frames, {} in {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, {} fulls, {} merged",
+                        k + 1, m.n_frames, m.space_group.xhm(), m.cell.a, m.cell.b, m.cell.c, m.cell.alpha,
+                        m.cell.beta, m.cell.gamma, m.fulls.size(), m.merged.size());
+        } else
+            logger.Warning("Sweep {} ({}) gave no merge and is left out of the joint one", k + 1, plan.sweeps[k].label);
+        members.push_back(std::move(m));
+        inputs.push_back(plan.sweeps[k].input);
+    }
+
+    // On one set of axes, and laid end to end: the second level's input.
+    const JointSweepRelation relation = RelateJointSweeps(members);
+    if (relation.reference < 0 || members[relation.reference].fulls.empty()) {
+        logger.Error("No sweep merged on its own, so there is nothing to merge jointly");
+        return EXIT_FAILURE;
+    }
+    for (size_t k = 0; k < members.size(); k++) {
+        const auto &r = relation.sweeps[k];
+        if (static_cast<int>(k) == relation.reference)
+            logger.Info("Sweep {}: the reference the others are put on the axes of", k + 1);
+        else if (r.included)
+            logger.Info("Sweep {}: put on the reference's axes by {} (CC {:.3f}, {:.3f} as indexed)", k + 1,
+                        r.cob.as_hkl().triplet('h'), r.cc, r.identity_cc);
+        else
+            logger.Warning("Sweep {} is left out of the joint merge: {}", k + 1,
+                           r.why_left_out.empty() ? std::string("it gave no merge") : r.why_left_out);
+    }
+    JointMergeInput input = BuildJointMergeInput(members, relation);
+    if (input.segments.size() < 2)
+        logger.Warning("Only {} sweep goes into the joint merge", input.segments.size());
+
+    // The joint merge, in the reference sweep's setup, writing what a run writes.
+    const JointMemberSweep &ref = members[relation.reference];
+    RugnuxInputRun joint(*this);
+    joint.input_files = {inputs[relation.reference]};
+    if (log_open)
+        g_log_file = output_prefix + "_log.txt";
+    int rc;
+    try {
+        rc = joint.PrepareRugnuxInput();
+        if (rc != RUGNUX_CONTINUE)
+            return rc;
+        ScalingSettings s = joint.experiment.GetScalingSettings();
+        s.SweepSegments(input.segments);
+        joint.experiment.ImportScalingSettings(s);
+        joint.config.rotation_postrefine_geometry = false;
+        joint.config.export_unmerged = false;
+        joint.config.export_unmerged_partials = false;
+        joint.config.image_sink = nullptr;
+        joint.joint_request = RugnuxInputRun::JointRequest{std::move(input.outcomes), ref.cell,
+                                                           ref.space_group.crystal_system(),
+                                                           ref.space_group.centring_type()};
+        logger.Info("Joint merge of {} sweeps, {} frames", input.segments.size(), joint.joint_request->outcomes.size());
+        rc = joint.RunPreparedInput();
+    } catch (const RugnuxInputFailed &f) {
+        rc = f.code;
+    }
+    return rc;
+}
+
+// The body of main. Settings setters and the pipeline itself throw JFJochException on input the parser
+// cannot reject on its own (a resolution limit of zero, a distance of zero, a polarisation above one,
+// an unreadable file), so main wraps this and reports rather than letting the exception terminate the
+// process with no diagnostic and exit code 134.
+static int RunRugnux(int argc, char **argv, bool no_gpu) {
+    // cuFFT plans are reused through the run and destroyed with it, ahead of the device reset in main.
+    FFTPlanReuse fft_plans;
+    // The whole invocation, reported at the end and written into the results report. Started before
+    // anything else so it covers reading the file and setting up as well as the processing.
+    const auto invocation_start = std::chrono::steady_clock::now();
+    RunProvenance provenance;
+    provenance.command_line = JoinCommandLine(argc, argv);
+    provenance.gpu_count = get_gpu_count();
+#ifdef __APPLE__
+    // Without CUDA, the Bragg integration's integer pixel work runs on the Apple GPU (MXAnalysisWithoutFPGA).
+    const std::string metal_device = BraggIntegrationEngineMetal::DeviceName();
+#else
+    const std::string metal_device;
+#endif
+    provenance.gpu_description = get_gpu_description();
+
+    // The terminal shows the run as steps (RugnuxConsole) and every Logger writes to the log file
+    // instead; with -v the whole log is on the terminal as well. -v is looked for here, before the
+    // options are parsed, because the parser already logs.
+    bool verbose = false;
+    for (int i = 1; i < argc; i++)
+        if (std::strcmp(argv[i], "-v") == 0 || std::strcmp(argv[i], "--verbose") == 0)
+            verbose = true;
+    console = std::make_shared<RugnuxConsole>(verbose);
+    std::atexit(CloseConsoleOnExit);
+    // Files with an older time than this are not this run's (the margin covers coarse file times).
+    const auto files_since = std::filesystem::file_time_type::clock::now() - std::chrono::seconds(2);
+    auto log_sinks = std::make_shared<spdlog::sinks::dist_sink_mt>();
+    // What is logged before the log file can be opened - its name comes from -o - waits here.
+    auto early_log = std::make_shared<spdlog::sinks::ringbuffer_sink_mt>(10000);
+    if (verbose)
+        log_sinks->add_sink(std::make_shared<spdlog::sinks::stdout_color_sink_mt>());
+    else
+        log_sinks->add_sink(MakeConsoleSink(console));
+    log_sinks->add_sink(early_log);
+    Logger::RouteAll(log_sinks);
+    // The HDF5 library prints its own error stack on stderr; every failure it reports reaches the
+    // reader as an error code anyway, which is turned into a message there.
+    H5Eset_auto2(H5E_DEFAULT, nullptr, nullptr);
+
+    PrintBanner(*console);
+
+    Logger logger("rugnux");
+#ifdef JFJOCH_USE_CUDA
+    const char *build = "CUDA build";
+#else
+    const char *build = "CPU build";
+#endif
+    logger.Info("rugnux Version {}{} - {}", jfjoch_version(),
+                jfjoch_git_sha1().empty() ? std::string()
+                : fmt::format(" (git {} {})", jfjoch_git_sha1().substr(0, 6), jfjoch_git_date()), build);
+    logger.Info("Build flags: {}", jfjoch_build_cxx_flags().empty() ? "none" : jfjoch_build_cxx_flags());
+    logger.Info("Build: {}", jfjoch_build_info(true));
+    logger.Info("Command line: {}", provenance.command_line);
+    logger.Info("GPU: {}", no_gpu ? std::string("not used (--no-gpu) - running on the CPU")
+            : provenance.gpu_count > 0 ? fmt::format("{} ({} visible)", provenance.gpu_description, provenance.gpu_count)
+            : !metal_device.empty() ? fmt::format("{} through Metal, for the Bragg integration only - the rest runs on the CPU", metal_device)
+            : std::string("none visible - running on the CPU"));
+    RegisterHDF5Filter();
+
+    RugnuxCli cli(no_gpu, invocation_start, provenance, metal_device, files_since, log_sinks, early_log);
+    if (const int rc = cli.Parse(argc, argv); rc != RUGNUX_CONTINUE)
+        return rc;
+    return cli.Run();
+}
+
 int main(int argc, char **argv) {
     try {
         // Before anything creates a CUDA context, which fixes the flags for the device. A thread

tests/CMakeLists.txt

@@ -141,6 +141,7 @@ ADD_EXECUTABLE(jfjoch_test
         PostRefineTest.cpp
         CellEsdTest.cpp
         RotationScaleWalkTest.cpp
+        RotationScaleMergeSegmentsTest.cpp
         AndersonMixingTest.cpp
         RugnuxLargeTest.cpp
         TestData.h
@@ -214,6 +215,8 @@ ADD_EXECUTABLE(jfjoch_test
         WriteReflectionsTest.cpp
         RugnuxResultFilesTest.cpp
         RugnuxInputKindTest.cpp
+        RugnuxSweepsTest.cpp
+        RugnuxJointSweepsTest.cpp
         UnitCellTest.cpp
         CCTest.cpp
         MultiLatticeSearchTest.cpp

tests/RotationScaleMergeSegmentsTest.cpp

@@ -0,0 +1,231 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <cstring>
+#include <map>
+#include <random>
+#include <vector>
+
+#include "../common/CUDAWrapper.h"
+#include "../common/Logger.h"
+#include "../image_analysis/scale_merge/RotationScaleMerge.h"
+
+namespace {
+
+// A synthetic crystal measured as pre-combined fulls on several sweeps: a cubic 20 A cell, Wilson
+// intensities at B = 15 A^2 to 1.6 A, each sweep at its own scale and relative B, Poisson-like noise.
+struct SweepSpec {
+    int n_frames;
+    double scale;
+    double relative_b;                // at the sweep's first frame
+    double relative_b_end = NAN;      // at its last; NaN: the same B throughout
+};
+
+struct Synthetic {
+    std::vector<IntegrationOutcome> frames;
+    std::map<std::array<int, 3>, double> I_true;   // by the Friedel-merged index
+};
+
+constexpr double CELL = 20.0;
+
+std::array<int, 3> FriedelKey(int h, int k, int l) {
+    const std::array<int, 3> a{h, k, l}, b{-h, -k, -l};
+    return std::max(a, b);
+}
+
+Synthetic MakeSweeps(const std::vector<SweepSpec> &sweeps, uint32_t seed) {
+    std::mt19937 rng(seed);
+    std::vector<std::array<int, 3>> hkl;
+    for (int h = -12; h <= 12; ++h)
+        for (int k = -12; k <= 12; ++k)
+            for (int l = -12; l <= 12; ++l) {
+                const double d = CELL / std::sqrt(static_cast<double>(h * h + k * k + l * l));
+                if ((h || k || l) && d >= 1.6) hkl.push_back({h, k, l});
+            }
+    Synthetic out;
+    std::exponential_distribution<double> wilson(1.0);
+    for (const auto &m: hkl) {
+        const auto key = FriedelKey(m[0], m[1], m[2]);
+        if (!out.I_true.contains(key)) {
+            const double s2 = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]) / (4.0 * CELL * CELL);
+            out.I_true[key] = 2000.0 * std::exp(-2.0 * 15.0 * s2) * (0.05 + wilson(rng));
+        }
+    }
+    std::uniform_int_distribution<size_t> pick(0, hkl.size() - 1);
+    std::normal_distribution<double> noise(0.0, 1.0);
+    int frame = 0;
+    for (const auto &sweep: sweeps)
+        for (int f = 0; f < sweep.n_frames; ++f, ++frame) {
+            IntegrationOutcome o;
+            o.latt = CrystalLattice(UnitCell{static_cast<float>(CELL), static_cast<float>(CELL),
+                                             static_cast<float>(CELL), 90, 90, 90});
+            for (int j = 0; j < 400; ++j) {
+                const auto &m = hkl[pick(rng)];
+                const double s2 = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]) / (4.0 * CELL * CELL);
+                const double b = std::isfinite(sweep.relative_b_end)
+                                 ? sweep.relative_b + (sweep.relative_b_end - sweep.relative_b) * f / (sweep.n_frames - 1)
+                                 : sweep.relative_b;
+                const double expected = out.I_true.at(FriedelKey(m[0], m[1], m[2])) * sweep.scale
+                                        * std::exp(-2.0 * b * s2);
+                const double sigma = std::sqrt(expected + 20.0);
+                Reflection r{};
+                r.h = m[0]; r.k = m[1]; r.l = m[2];
+                r.image_number = static_cast<float>(frame);
+                r.d = static_cast<float>(1.0 / std::sqrt(4.0 * s2));
+                r.I = static_cast<float>(expected + sigma * noise(rng));
+                r.sigma = static_cast<float>(sigma);
+                // The combine rebuilds a full's variance as var_bkg + I: what is not the intensity's own
+                r.var_bkg = static_cast<float>(sigma * sigma - std::max(0.0, static_cast<double>(r.I)));
+                r.bkg = 20.0f;
+                r.prescaling_corr = 1.0f;
+                r.partiality = 1.0f;
+                r.zeta = 1.0f;
+                r.image_scale_corr = 1.0f;
+                r.predicted_x = 500.0f;
+                r.predicted_y = 500.0f;
+                o.reflections.push_back(r);
+            }
+            out.frames.push_back(std::move(o));
+        }
+    return out;
+}
+
+DiffractionExperiment SyntheticExperiment(bool correction_surfaces) {
+    DiffractionExperiment x;
+    x.IncidentEnergy_keV(12.4f);
+    x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(1, 0, 0), {}));
+    x.SetUnitCell(UnitCell{static_cast<float>(CELL), static_cast<float>(CELL), static_cast<float>(CELL), 90, 90, 90});
+    ScalingSettings s;
+    s.ScaleFulls(true).CorrectionSurfaces(correction_surfaces).ResolutionCutoff(ResolutionCutoffMethod::Off);
+    x.ImportScalingSettings(s);
+    return x;
+}
+
+std::vector<SweepSegment> Segments(const std::vector<SweepSpec> &sweeps) {
+    std::vector<SweepSegment> out;
+    int32_t first = 0;
+    for (const auto &s: sweeps) {
+        out.push_back({first, s.n_frames});
+        first += s.n_frames;
+    }
+    return out;
+}
+
+RotationScaleMerge::Result Merge(Synthetic &data, const DiffractionExperiment &x,
+                                 const std::vector<SweepSegment> &segments, bool export_fulls = false,
+                                 bool *on_gpu = nullptr) {
+    Logger logger("RotationScaleMergeSegmentsTest");
+    const UnitCell cell{static_cast<float>(CELL), static_cast<float>(CELL), static_cast<float>(CELL), 90, 90, 90};
+    RotationScaleMerge engine(x, data.frames, cell, 20, 4, logger);
+    if (!segments.empty())
+        engine.SetSweepSegments(segments);
+    engine.SetExportScaledFulls(export_fulls);
+    engine.Ingest();
+    if (on_gpu)
+        *on_gpu = engine.OnGPU();
+    return engine.Run(false, true, false);
+}
+
+// The median over a frame range of the scaled fulls' ratio to the true intensity (the strong ones).
+double MedianRatio(const RotationScaleMerge::Result &r, const Synthetic &data, int f0, int f1) {
+    std::vector<double> ratio;
+    for (const auto &o: r.scaled_fulls)
+        if (o.image_number >= static_cast<float>(f0) && o.image_number < static_cast<float>(f1)) {
+            const double t = data.I_true.at(FriedelKey(o.h, o.k, o.l));
+            if (t > 500.0)
+                ratio.push_back(o.I / t);
+        }
+    REQUIRE(ratio.size() > 20);
+    std::nth_element(ratio.begin(), ratio.begin() + ratio.size() / 2, ratio.end());
+    return ratio[ratio.size() / 2];
+}
+
+}  // namespace
+
+TEST_CASE("RotationScaleMergeSegments_PerSweepScaleRecovered", "[RotationScale]") {
+    const std::vector<SweepSpec> sweeps{{60, 1.0, 0.0}, {60, 0.5, 0.0}};
+    auto data = MakeSweeps(sweeps, 11);
+    const auto r = Merge(data, SyntheticExperiment(false), Segments(sweeps), true);
+    REQUIRE(r.segments.size() == 2);
+    CHECK(r.segments[1].scale / r.segments[0].scale == Catch::Approx(0.5).epsilon(0.03));
+    CHECK(r.statistics.overall.cc_half > 0.99);
+
+    // No coupling across the boundary: the frames on either side of it are put on the common scale as
+    // well as the middle of their sweeps are.
+    const double middle = MedianRatio(r, data, 20, 40);
+    CHECK(MedianRatio(r, data, 57, 60) / middle == Catch::Approx(1.0).epsilon(0.03));
+    CHECK(MedianRatio(r, data, 60, 63) / middle == Catch::Approx(1.0).epsilon(0.03));
+    CHECK(MedianRatio(r, data, 80, 100) / middle == Catch::Approx(1.0).epsilon(0.03));
+}
+
+TEST_CASE("RotationScaleMergeSegments_DoseOrderFromTheRelativeB", "[RotationScale]") {
+    // Given first, the sweep taken second: its B is 8 A^2 higher. The order comes from the data.
+    const std::vector<SweepSpec> sweeps{{60, 1.0, 8.0}, {60, 1.0, 0.0}, {60, 1.0, 4.0}};
+    auto data = MakeSweeps(sweeps, 12);
+    const auto r = Merge(data, SyntheticExperiment(true), Segments(sweeps));
+    REQUIRE(r.segments.size() == 3);
+    CHECK(r.segments[1].dose_rank == 0);
+    CHECK(r.segments[2].dose_rank == 1);
+    CHECK(r.segments[0].dose_rank == 2);
+    CHECK(r.segments[1].relative_b == Catch::Approx(0.0).margin(1e-9));
+    CHECK(r.segments[0].relative_b == Catch::Approx(8.0).margin(1.0));
+    CHECK(r.segments[2].relative_b == Catch::Approx(4.0).margin(1.0));
+    // Flat within each sweep: one B per sweep predicts the held-out frames better than a dose-linear ramp
+    CHECK(r.segment_decay_model == "per-sweep B");
+}
+
+TEST_CASE("RotationScaleMergeSegments_DoseLinearDecayAcrossSweeps", "[RotationScale]") {
+    // One crystal damaged steadily: the sweep given first was taken second, its B going on from where the
+    // other's ended. Along the dose order the B rises in a straight line, which the ramp describes.
+    const std::vector<SweepSpec> sweeps{{60, 1.0, 6.0, 12.0}, {60, 1.0, 0.0, 6.0}};
+    auto data = MakeSweeps(sweeps, 15);
+    const auto r = Merge(data, SyntheticExperiment(true), Segments(sweeps));
+    REQUIRE(r.segments.size() == 2);
+    CHECK(r.segments[1].dose_rank == 0);
+    CHECK(r.segments[0].dose_rank == 1);
+    CHECK(r.segments[0].relative_b == Catch::Approx(6.0).margin(1.5));
+    CHECK(r.segment_decay_model == "dose-linear");
+}
+
+TEST_CASE("RotationScaleMergeSegments_SegmentsMustCoverTheFrames", "[RotationScale]") {
+    const std::vector<SweepSpec> sweeps{{10, 1.0, 0.0}, {10, 1.0, 0.0}};
+    auto data = MakeSweeps(sweeps, 13);
+    const auto x = SyntheticExperiment(false);
+    CHECK_THROWS_AS(Merge(data, x, {{0, 10}, {11, 9}}), JFJochException);
+    CHECK_THROWS_AS(Merge(data, x, {{0, 10}}), JFJochException);
+}
+
+TEST_CASE("RotationScaleMergeSegments_NoSegmentsIsTheEngineAsBefore", "[RotationScale]") {
+    // The same synthetic sweep merged without segments gives, to the bit, what the engine gave before
+    // segments existed: FNV-1a over the merged I and sigma, recorded on the engine of rc181 e56b19abb -
+    // on the CPU, and on the GPU for the card model it was recorded on (the device path matches the host
+    // path to the bit only where every reduction does; this one carries its own, so it is checked against
+    // its own record).
+    const std::vector<SweepSpec> sweeps{{40, 1.0, 0.0}, {40, 0.7, 3.0}};
+    auto data = MakeSweeps(sweeps, 14);
+    bool on_gpu = false;
+    const auto r = Merge(data, SyntheticExperiment(true), {}, false, &on_gpu);
+    uint64_t hash = 1469598103934665603ULL;
+    auto mix = [&hash](float v) {
+        uint32_t u;
+        std::memcpy(&u, &v, sizeof u);
+        for (int b = 0; b < 4; ++b) {
+            hash ^= (u >> (8 * b)) & 0xffu;
+            hash *= 1099511628211ULL;
+        }
+    };
+    for (const auto &m: r.merged) {
+        mix(m.I);
+        mix(m.sigma);
+    }
+    CHECK(r.segments.empty());
+    CHECK(r.merged.size() == 4106);
+    if (!on_gpu)
+        CHECK(hash == 3376995616949734726ULL);
+    else if (get_gpu_description().find("RTX 5080") != std::string::npos)
+        CHECK(hash == 6594155396455850264ULL);
+}

tests/RugnuxJointSweepsTest.cpp

@@ -0,0 +1,108 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <cmath>
+#include <map>
+#include <random>
+
+#include "../rugnux/RugnuxJointSweeps.h"
+
+namespace {
+
+// A tetragonal P4 crystal: Laue class 4/m on a 4/mmm lattice, so (k, h, -l) indexes it a second way.
+const gemmi::SpaceGroup &P4() { return *gemmi::find_spacegroup_by_name("P 4"); }
+constexpr float A = 60.0f, C = 40.0f;
+
+std::array<int, 3> P4Key(int h, int k, int l) {
+    // the 4/m ASU representative, enough to give equivalents one intensity
+    std::array<int, 3> best{h, k, l};
+    const std::array<std::array<int, 3>, 8> eq{{{h, k, l}, {-k, h, l}, {-h, -k, l}, {k, -h, l},
+                                                {-h, -k, -l}, {k, -h, -l}, {h, k, -l}, {-k, h, -l}}};
+    for (const auto &e: eq) best = std::max(best, e);
+    return best;
+}
+
+JointMemberSweep MakeSweep(const std::map<std::array<int, 3>, double> &truth, const gemmi::Op &indexed_as,
+                           int n_frames, uint32_t seed) {
+    std::mt19937 rng(seed);
+    JointMemberSweep s;
+    s.n_frames = n_frames;
+    s.cell = UnitCell{A, A, C, 90, 90, 90};
+    s.space_group = P4();
+    for (const auto &[hkl, I]: truth) {
+        const auto m = indexed_as.apply_to_hkl({{hkl[0], hkl[1], hkl[2]}});
+        MergedReflection r;
+        r.h = m[0]; r.k = m[1]; r.l = m[2];
+        r.I = static_cast<float>(I);
+        r.sigma = 1.0f;
+        r.d = static_cast<float>(1.0 / std::sqrt((hkl[0] * hkl[0] + hkl[1] * hkl[1]) / (A * A) + hkl[2] * hkl[2] / (C * C)));
+        s.merged.push_back(r);
+        ScaledFull f;
+        f.h = m[0]; f.k = m[1]; f.l = m[2];
+        f.I = static_cast<float>(I);
+        f.sigma = 1.0f;
+        f.image_number = static_cast<float>(rng() % n_frames);
+        f.var_bkg = 4.0f;
+        f.var_per_I = 0.5f;
+        f.capture = 0.0f;
+        f.d = r.d;
+        s.fulls.push_back(f);
+    }
+    return s;
+}
+
+std::map<std::array<int, 3>, double> Truth(uint32_t seed) {
+    std::mt19937 rng(seed);
+    std::exponential_distribution<double> wilson(1.0);
+    std::map<std::array<int, 3>, double> t;
+    for (int h = -12; h <= 12; ++h)
+        for (int k = -12; k <= 12; ++k)
+            for (int l = -8; l <= 8; ++l) {
+                if (!(h || k || l)) continue;
+                const auto key = P4Key(h, k, l);
+                if (!t.contains(key)) t[key] = 1000.0 * wilson(rng);
+            }
+    return t;
+}
+
+}  // namespace
+
+TEST_CASE("RelateJointSweeps_AmbiguityRecoveredAndAStrangerLeftOut", "[Rugnux]") {
+    const auto truth = Truth(1);
+    const gemmi::Op twin = gemmi::parse_triplet("y,x,-z");   // h,k,l -> k,h,-l
+    std::vector<JointMemberSweep> sweeps;
+    sweeps.push_back(MakeSweep(truth, gemmi::Op::identity(), 100, 2));
+    sweeps.push_back(MakeSweep(truth, twin, 80, 3));
+    sweeps.push_back(MakeSweep(Truth(99), gemmi::Op::identity(), 50, 4));   // another crystal
+    sweeps[0].fulls.push_back(sweeps[0].fulls.front());                     // the reference: most fulls
+
+    const auto relation = RelateJointSweeps(sweeps);
+    REQUIRE(relation.reference == 0);
+    CHECK(relation.sweeps[0].included);
+    REQUIRE(relation.sweeps[1].included);
+    CHECK(relation.sweeps[1].cc > 0.99);
+    CHECK(relation.sweeps[1].identity_cc < 0.5);
+    // the operator found takes sweep 2's indices back onto the reference's (up to the Laue group)
+    const auto back = relation.sweeps[1].cob.apply_to_hkl(twin.apply_to_hkl({{3, 1, 2}}));
+    CHECK(P4Key(back[0], back[1], back[2]) == P4Key(3, 1, 2));
+    CHECK_FALSE(relation.sweeps[2].included);
+
+    const auto input = BuildJointMergeInput(sweeps, relation);
+    REQUIRE(input.segments.size() == 2);
+    CHECK(input.segments[0].first_frame == 0);
+    CHECK(input.segments[0].n_frames == 100);
+    CHECK(input.segments[1].first_frame == 100);
+    CHECK(input.segments[1].n_frames == 80);
+    CHECK(input.segment_sweep == std::vector<int>{0, 1});
+    REQUIRE(input.outcomes.size() == 180);
+    // a full as the combine takes it: count I/c at correction c = var_per_I, variance var_bkg / c^2
+    const Reflection *r = nullptr;
+    for (const auto &o: input.outcomes)
+        if (!o.reflections.empty()) { r = &o.reflections.front(); break; }
+    REQUIRE(r != nullptr);
+    CHECK(r->image_scale_corr == Catch::Approx(0.5));
+    CHECK(r->var_bkg == Catch::Approx(16.0));
+    CHECK(r->partiality == 1.0f);
+}

tests/RugnuxSweepsTest.cpp

@@ -0,0 +1,132 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <filesystem>
+#include <fstream>
+#include <map>
+#include <stdexcept>
+
+#include "../rugnux/RugnuxSweeps.h"
+
+namespace {
+
+// Inputs described by name alone: "<template>_<n>.cbf" frames, "<name>_rot.h5" / "<name>_still.h5"
+// masters, "<name>.cxi" for the stills readers; the wavelength from a table.
+struct FakeInputs {
+    std::map<std::string, double> wavelength;
+    int sweeps_described = 0;
+
+    RugnuxInputDescription File(const std::string &f) const {
+        RugnuxInputDescription d;
+        if (f.ends_with(".cbf"))
+            d.frame_template = RugnuxFrameTemplate(f);
+        d.foreign = f.ends_with(".cxi");
+        return d;
+    }
+    RugnuxInputDescription Sweep(const std::string &f) {
+        sweeps_described++;
+        RugnuxInputDescription d = File(f);
+        d.rotation = f.ends_with(".cbf") || f.ends_with("_rot.h5");
+        d.wavelength_A = wavelength.contains(f) ? wavelength.at(f) : 1.0;
+        return d;
+    }
+    RugnuxSweepPlan Plan(const std::vector<std::string> &inputs) {
+        return PlanRugnuxSweeps(inputs, [this](const std::string &f) { return File(f); },
+                                [this](const std::string &f) { return Sweep(f); });
+    }
+};
+
+}  // namespace
+
+TEST_CASE("RugnuxFrameTemplate", "[Rugnux][portable]") {
+    CHECK(RugnuxFrameTemplate("/d/x_1_0001.cbf") == "/d/x_1_####.cbf");
+    CHECK(RugnuxFrameTemplate("x_00017.cbf.gz") == "x_#####.cbf.gz");
+    CHECK(RugnuxFrameTemplate("formal.001") == "formal.###");
+    CHECK(RugnuxFrameTemplate("noframe.cbf").empty());
+}
+
+TEST_CASE("PlanRugnuxSweeps_single_input_is_not_planned", "[Rugnux][portable]") {
+    FakeInputs in;
+    CHECK(in.Plan({"a_0001.cbf"}).sweeps.empty());
+    CHECK(in.Plan({"run_rot.h5"}).sweeps.empty());
+    CHECK(in.sweeps_described == 0);   // nothing is opened for one input
+}
+
+TEST_CASE("PlanRugnuxSweeps_frames_by_template", "[Rugnux][portable]") {
+    FakeInputs in;
+    SECTION("frames of one template are its one sweep, named by the first given") {
+        const auto plan = in.Plan({"d/a_0003.cbf", "d/a_0001.cbf", "d/a_0002.cbf"});
+        REQUIRE(plan.sweeps.size() == 1);
+        CHECK(plan.sweeps[0].input == "d/a_0003.cbf");
+        CHECK(plan.sweeps[0].frames);
+    }
+    SECTION("two templates are two sweeps, in the order given") {
+        const auto plan = in.Plan({"d/b_0001.cbf", "d/b_0002.cbf", "d/a_0001.cbf"});
+        REQUIRE(plan.sweeps.size() == 2);
+        CHECK(plan.sweeps[0].input == "d/b_0001.cbf");
+        CHECK(plan.sweeps[0].label == std::filesystem::path("d/b_####.cbf").string());
+        CHECK(plan.sweeps[1].input == "d/a_0001.cbf");
+        CHECK(plan.wavelength_A == Catch::Approx(1.0));
+    }
+    SECTION("the same name in two directories is two sweeps") {
+        CHECK(in.Plan({"d1/a_0001.cbf", "d2/a_0001.cbf"}).sweeps.size() == 2);
+    }
+    SECTION("sweeps at different wavelengths are refused") {
+        in.wavelength["d/b_0001.cbf"] = 1.2;
+        CHECK_THROWS_AS(in.Plan({"d/a_0001.cbf", "d/b_0001.cbf"}), std::runtime_error);
+    }
+}
+
+TEST_CASE("PlanRugnuxSweeps_masters", "[Rugnux][portable]") {
+    FakeInputs in;
+    SECTION("stills masters stay one stills run") {
+        CHECK(in.Plan({"a_still.h5", "b_still.h5", "c_still.h5"}).sweeps.empty());
+        CHECK(in.sweeps_described == 1);   // the first decides; the stills reader refuses a rotation one
+    }
+    SECTION("rotation masters are a sweep each") {
+        const auto plan = in.Plan({"a_rot.h5", "b_rot.h5"});
+        REQUIRE(plan.sweeps.size() == 2);
+        CHECK(plan.sweeps[1].input == "b_rot.h5");
+        CHECK_FALSE(plan.sweeps[0].frames);
+    }
+    SECTION("rotation and stills are refused") {
+        CHECK_THROWS_AS(in.Plan({"a_rot.h5", "b_still.h5"}), std::runtime_error);
+    }
+    SECTION("rotation masters at different wavelengths are refused") {
+        in.wavelength["b_rot.h5"] = 0.98;
+        CHECK_THROWS_AS(in.Plan({"a_rot.h5", "b_rot.h5"}), std::runtime_error);
+    }
+    SECTION("a wavelength within the tolerance is the same") {
+        in.wavelength["b_rot.h5"] = 1.0 + 1e-4;
+        CHECK(in.Plan({"a_rot.h5", "b_rot.h5"}).sweeps.size() == 2);
+    }
+}
+
+TEST_CASE("PlanRugnuxSweeps_leaves_the_stills_readers_and_mixtures_alone", "[Rugnux][portable]") {
+    FakeInputs in;
+    CHECK(in.Plan({"a.cxi", "b.cxi"}).sweeps.empty());
+    CHECK(in.Plan({"a_0001.cbf", "b_rot.h5"}).sweeps.empty());   // the readers refuse the mix, as before
+    CHECK(in.sweeps_described == 0);
+}
+
+TEST_CASE("RugnuxSiblingSweeps", "[Rugnux][portable]") {
+    namespace fs = std::filesystem;
+    const fs::path dir = fs::path("rugnux_sweeps_siblings");
+    fs::remove_all(dir);
+    fs::create_directories(dir);
+    for (const char *name: {"a_0001.cbf", "a_0002.cbf", "b_0001.cbf", "b_0002.cbf", "b_0003.cbf", "c_0001.cbf.gz",
+                            "notes.txt", "x_master.h5", "y_master.h5"})
+        std::ofstream(dir / name) << "x";
+
+    RugnuxSweepPlan plan;
+    plan.sweeps.push_back({(dir / "a_0001.cbf").string(), RugnuxFrameTemplate((dir / "a_0001.cbf").string()), true});
+    plan.sweeps.push_back({(dir / "x_master.h5").string(), (dir / "x_master.h5").string(), false});
+    const auto siblings = RugnuxSiblingSweeps(plan, {".cbf"});
+    REQUIRE(siblings.size() == 3);
+    CHECK(siblings[0] == (dir / "b_####.cbf").string() + " (3 frames)");
+    CHECK(siblings[1] == (dir / "c_####.cbf.gz").string() + " (1 frame)");
+    CHECK(siblings[2] == (dir / "y_master.h5").string());
+    fs::remove_all(dir);
+}