rugnux: --sweeps joint merges several sweeps of one crystal as one dataset

Two levels, as XSCALE merges the data sets XDS made of each sweep. Each sweep is processed on its own
first (a member run: everything a run does, writing only <prefix>_s<k>_log.txt, without the model),
and its written merge hands back all of its fulls - past its resolution cut, with the reflections its
group extinguishes, and with the counting variance the merge rebuilds them from (ScaledFull var_bkg,
var_per_I, capture, d). RelateJointSweeps puts every sweep on the axes of the one with the most fulls
by ChooseReferenceSetting (axis permutations and the lattice's twin laws, scored by the merged-intensity
CC); a sweep below MIN_REFERENCE_CC on every axes is left out and named. BuildJointMergeInput lays the
sweeps' frames end to end, each full a pre-combined observation whose count, correction and variance
reproduce what its own merge gave it, and Rugnux::RunJoint runs the one pass of the pipeline on them:
no image is read, the outcomes are the indexer's, and the unchanged ScaleMergeAndSymmetry makes the
merge (RotationScaleMerge with the sweep segments, now taken from ScalingSettings::SweepSegments), the
space-group search, the cut and everything after it, into <prefix>.*. The run's log is opened first and
holds the sweeps' relation and the joint merge. Not yet: the unmerged MTZ of a joint merge, the per-sweep
report keys, and joint geometry (P1d / phase 2).

Single-sweep runs are unchanged (RunJoint's branches in RunPipeline are only taken by it; the exports
are only asked for by member runs).

Test pairs (rugnux joint against the best sweep alone, and against XSCALE, at the XDS range):
- in-house lysozyme, two 450 x 0.2 deg sweeps: P41212 (both sweeps alone too); to 1.85 A fitted against
  1.91-1.94; at 50-1.97 A R_meas 14.5% against 14.9%, I/sigma 14.1 against 9.4, ISa 24.4 against 20.8;
  XSCALE of the two XDS sweeps: R_meas 17.4%, I/sigma 9.7, ISa 16.0 / 15.2.
- 8C3E, four short CCD sweeps of a twinned P3121 crystal: the joint merge determines P3121 (the
  deposition), where the sweeps alone give P321, P6222, P321 and P622; completeness 98% against 96%.
- 5JVN, two 3600-frame sweeps: P622; to 2.32 A against 2.35; at 50-2.62 A R_meas 16.6% against 13.7%,
  I/sigma 19.7 against 17.0, ISa 11.2 against 15.3; the depositors' XSCALE: R_meas 15.5%, ISa 15.3 / 13.5.

Tests: RugnuxJointSweepsTest (a sweep indexed by the twin law put back, another crystal left out, the
segments and the pre-combined observation's count and variance).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

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/RUGNUX_CLI.md

@@ -30,11 +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 --sweeps separate a_master.h5 b_master.h5   # several sweeps, each alone
+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. Several sweeps can be given (below, Inputs); for now
-each is processed alone with `--sweeps separate`. 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)).
 
@@ -105,9 +104,19 @@ template (the name with its frame number taken out - any frame of a sweep names
 `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.
-`--sweeps separate` processes each sweep alone, as if it were the only input, into
-`<prefix>_s<k>_*` (k = 1, 2, ... in the order given); processing several sweeps as one dataset
-(`--sweeps joint`, the default) is not implemented yet and refuses.
+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
@@ -142,7 +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): `separate` processes each alone into `<prefix>_s<k>_*`; `joint` (the default), one dataset of all of them, is not implemented yet |
+| `--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)
@@ -6113,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};
             }
         });
     }

image_analysis/scale_merge/RotationScaleMerge.h

@@ -81,17 +81,6 @@ std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes
                                        const std::vector<SweepEndTally> &tally);
 std::pair<int, int> TrimBlankSweepEnds(std::vector<IntegrationOutcome> &outcomes, size_t nthreads);
 
-// One sweep of several whose frames are laid end to end in one engine (the joint merge of several
-// sweeps): frames [first_frame, first_frame + n_frames) of the outcomes. 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 nothing is combined again, nothing is recomputed from the rocking geometry, and nothing that runs
-// along the frames - the per-frame scale's smoothness above all - reaches across from one sweep into the
-// next.
-struct SweepSegment {
-    int32_t first_frame = 0;
-    int32_t n_frames = 0;
-};
-
 class RotationScaleMerge {
 public:
     // The log scale x minimising sum_f J_f (x_f - y_f)^2 + lambda * sum_f (x_f - 2 x_{f+1} + x_{f+2})^2:
@@ -310,12 +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). 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 (the default) nothing of this runs.
+    // 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; }
@@ -763,6 +756,7 @@ private:
                                                 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) ---

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

rugnux/Rugnux.cpp

@@ -787,6 +787,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;
@@ -818,7 +833,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;
@@ -907,7 +924,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);
@@ -1120,6 +1137,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();
@@ -1156,7 +1177,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.
@@ -1194,7 +1227,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);
@@ -1205,7 +1238,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

@@ -114,6 +114,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).
@@ -1261,6 +1266,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
@@ -1656,6 +1666,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

@@ -1688,6 +1688,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/rugnux_cli.cpp

@@ -86,6 +86,7 @@
 #include "RugnuxConsole.h"
 #include "RugnuxInputKind.h"
 #include "RugnuxSweeps.h"
+#include "RugnuxJointSweeps.h"
 #ifdef __APPLE__
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineMetal.h"
 #endif
@@ -123,7 +124,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 --sweeps separate a_master.h5 b_master.h5    several sweeps, each alone
+  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.
@@ -172,7 +173,7 @@ Basic options
   --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>_*; processing them as one dataset is to come
+                                 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:
@@ -960,6 +961,8 @@ struct RugnuxCli {
     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_,
@@ -971,6 +974,8 @@ struct RugnuxCli {
 
     // 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();
 };
@@ -999,6 +1004,18 @@ struct RugnuxInputRun : RugnuxCli {
     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).
@@ -1764,7 +1781,7 @@ int RugnuxInputRun::PrepareRugnuxInput() {
 
     // 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);
@@ -3600,7 +3617,9 @@ int RugnuxInputRun::RunPreparedInput() {
 
     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.
@@ -3638,6 +3657,10 @@ int RugnuxInputRun::RunPreparedInput() {
         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) {
@@ -3880,11 +3903,8 @@ int RugnuxCli::Run() {
         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) {
-        logger.Error("{} sweeps: processing them as one dataset is not implemented yet. --sweeps separate "
-                     "processes each sweep alone, into <prefix>_s<k>_*", plan.sweeps.size());
-        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;
@@ -3917,6 +3937,118 @@ int RugnuxCli::Run() {
 }
 
 
+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

tests/CMakeLists.txt

@@ -214,6 +214,7 @@ ADD_EXECUTABLE(jfjoch_test
         RugnuxResultFilesTest.cpp
         RugnuxInputKindTest.cpp
         RugnuxSweepsTest.cpp
+        RugnuxJointSweepsTest.cpp
         UnitCellTest.cpp
         CCTest.cpp
         MultiLatticeSearchTest.cpp

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);
+}