v1.0.0-rc.164 (#74)

* rugnux now tells you whether a crystal diffracts anisotropically and how far it reaches in each direction, without a second program: a new `9. DIFFRACTION ANISOTROPY` section in `<prefix>_report.txt` and matching `_reflns.pdbx_aniso_B_tensor_*` / `_reflns.jfjoch_aniso_*` items in the merged mmCIF report the anisotropic deltaB, the diffraction limit along each principal direction, and a `NOT DETECTED` / `DETECTED` / `CANNOT DETERMINE` verdict measured against the data set's own systematic error. It is a description only - no intensity is corrected, no reflection is removed, and the merged data do not depend on direction.
* rugnux can hand its integrated observations to another scaling program: `--export-unmerged` writes `<prefix>_unmerged.mtz`, an unmerged MTZ readable by aimless, pointless, careless and `iotbx.merging_statistics`, in `--mode mx` and `--mode scale` alike. Each rotation reflection's partials are summed into one full; `--export-unmerged-partials` writes one row per image instead. Intensities carry the Lorentz-polarization factor and nothing else, since those programs scale the data themselves. Lattice-centring absences are not written; screw and glide absences are.
* rugnux integrates crystals with broad spots better - where it changes anything, per-shell mean I/sigma improves by up to 31% and R_meas by up to 24% - because on rotation data the integration signal radius is now taken from the crystal's own measured spot width instead of a fixed 4 px. `--adaptive-integration-radius=off` restores the fixed radius and an explicit `--integration-radius` still overrides both. The widened radius applies to the final integration pass only, and a pattern too dense for it is re-integrated at 4 px with a note in the log.
* rugnux discards fewer stills reflections for want of a background ring, improving per-shell R_meas over most of the signal-bearing range: the stills background ring now runs to 14 px instead of 12. The gain reverses in shells below a mean I/sigma of about 4.
* rugnux determines the space group with thresholds that mean the same thing on a weak crystal as on a strong one: symmetry operators are scored on resolution-normalised intensities (E squared) instead of raw merged intensities, and a reflection counts as genuinely present on its counting significance instead of on the merged I/sigma, which saturates at the merge's own ISa. The search resolution cut is no longer able to move the answer, and the twin-law H bound moves from 1.70 to 1.85, which stops one class of correct high-symmetry assignment being refused as twinning.
* rugnux says what the space-group search tested and what it could not: the twin-law disagreement H is printed for every operator together with the adopted point group's H ratio and its bound; alternatives that are not on the reported lattice are named with how their cell differs; and a lattice centring the data could not test - the crystal having been integrated on the primitive sub-cell, so the reflections it extinguishes were never measured - is marked `UNTESTED` and warned about where it is adopted, as coming from the lattice metric rather than from the intensities.
* rugnux `--mode scale` re-merges a `_process.h5` in the right symmetry without being told it: the file now records the space group on every run - a two-pass rotation run wrote none before, so re-merging defaulted to P1 - together with the change of basis under `/entry/MX/reindexMatrix` where the lattice was re-seated, and `--mode scale` also reports the Wilson B-factor estimate instead of `WILSON_B= nan`. A file written before this stops with a message naming the two cells and the override to use, instead of failing inside the merge. A third-party reader of a `_process.h5` must apply `reindexMatrix` where it is present.
* rugnux installs on its own, as a package called `rugnux` - `dnf install rugnux` or `apt install rugnux` - instead of arriving inside `jfjoch-viewer`. It pulls in none of the acquisition stack, so a machine that only processes data no longer has to carry the broker, the detector libraries or Qt to get it. Installing it over a `jfjoch-viewer` from rc.163 or earlier, which still owns `/usr/bin/rugnux`, upgrades cleanly rather than failing on the duplicate file.
* rugnux is also a standalone download, built for arm64 as well as x86_64: `rugnux-<version>-linux-{x86_64|aarch64}-cuda<major>.tgz` and `rugnux-<version>-win64-cuda<major>.zip` on the release page, for machines that are not managed by a package manager. The aarch64 build targets GH200 and DGX Spark, and is untested on hardware.
* Every portable Linux binary is now a single self-contained file: cuFFT is linked statically instead of being shipped beside the executable and found through an rpath, so `rugnux` and `jfjoch_viewer` need nothing but an NVIDIA driver, and only to use the GPU. The `.rpm`/`.deb` continue to take cuFFT from the distribution. The developer utilities `jfjoch_extract_hkl` and `jfjoch_recompress` are no longer packaged anywhere.
* Jungfraujoch needs six fewer shared libraries on the machine - libopenblas and libmetis, and libgfortran, libquadmath, libgomp and libz behind them - because the Ceres LAPACK, METIS and SuiteSparse back-ends are no longer built. Nothing in the code ever selected them, and results are unchanged.
* The PCIe driver DKMS package builds for the kernel it is being installed for instead of the running one, so a module built while a kernel update is being applied loads after the reboot.
* The PCIe driver builds on RHEL 9.5 and later, and on their CentOS Stream, Rocky and AlmaLinux equivalents, where the `vm_flags` kernel interface was backported into the 5.14 kernel.
* A data collection started with `async_start` that fails to start - a writer refusing to overwrite an existing file, for instance - is reported as an error by `/wait_until_running` and `/wait_till_done` instead of as a timeout and a successful collection respectively. The error message is the one the writer gave.
* A calibration that is cancelled or that fails to collect its pedestals is no longer reported as a successful one. The broker goes to `Inactive` with an error message and has to be initialized again, instead of sitting in `Idle` looking ready to measure while holding partial pedestals - data collected in that state was silently mis-converted.
* A failed `/initialize` is reported to `/wait_until_running` and `/wait_till_done` as soon as it happens, instead of when their timeout expires.
* `space_group_number` accepts space groups up to 230 in the API schema, so cubic space groups can be recorded. The broker always accepted them; the generated clients rejected them before the request was sent.
* The results report's `REPORT_VERSION` is 3, two sections having been added. Existing key names and table columns are unchanged.
* The merged statistics table has **9** resolution shells instead of 10, which is what XDS reports. The bins were already XDS's - equal steps in 1/d^2 between the lowest- and the highest-resolution reflection the merge kept - so at the same resolution limits the two tables now have the same shell boundaries and can be read row for row. `--resolution-shells` sets a different count.
* `rugnux --model` now settles the frame the merged reflections are written in, not only the frame the R-factors and the maps are computed in: the `.mtz`/`.cif`/`.hkl` come out in the model's indexing, and where the data were merged in the model's enantiomorph they take the model's hand and space group - which on anomalous data puts I(+) and I(-) the right way round. The indexing choice is logged with the winning R-free and the runner-up, so a decision made within noise is visible.
* `rugnux --model` can resolve the indexing ambiguity of a **serial stills** run, which a model could not do before: structure factors computed from the model become the per-image reference, the same role a reference MTZ plays. It needs the cell and space group up front (`-C` / `-S`). Without one or the other, a merohedral serial run still merges both hands together and says so.
* The rugnux documentation opens with a quick start - the default run, and runs with a reference MTZ, with a model, or with the space group and cell pinned - and explains the indexing ambiguity: what it costs on rotation and on serial data, and which of `-z` / `--model` resolves it in each case. The long reference pages now carry a table of contents.

Reviewed-on: #74
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>

image_analysis/WriteReflections.cpp

@@ -5,12 +5,14 @@
 #include "scale_merge/Merge.h"
 #include "scale_merge/HKLKey.h"
 #include "scale_merge/TwinningAnalysis.h"
+#include "bragg_integration/SystematicAbsence.h"
 #include "../common/ParallelFor.h"
 
 #include <algorithm>
 #include <cmath>
 #include <cstdio>
 #include <map>
+#include <set>
 #include <tuple>
 #include <fstream>
 #include <iomanip>
@@ -272,6 +274,51 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
         if (std::isfinite(statistics.radiation_damage_delta_b))
             out << "_reflns.jfjoch_radiation_damage_relative_B  " << Fmt(statistics.radiation_damage_delta_b, 2)
                 << "   # relative-B first->last over the run (A^2); + = high-res fades with dose\n";
+        // Diffraction anisotropy. The eigen-decomposition of the anisotropy tensor has standard PDBx
+        // items; the eigenvalues there must be non-negative, and only the deviatoric part of the tensor
+        // is determined at all (its isotropic part is degenerate with the overall scale), so they are
+        // written relative to the weakest direction - eigenvalue_3 is 0 by construction and
+        // eigenvalue_1 is the anisotropic deltaB. The eigenvectors are in the PDB orthogonalisation
+        // convention, which is the one gemmi (and hence rugnux) uses throughout.
+        //
+        // The directional diffraction LIMITS are deliberately NOT written as
+        // _reflns.pdbx_aniso_diffraction_limit_*: the dictionary defines those as the semi-axes of an
+        // ellipsoid fitted to a diffraction cut-off surface, which is a different construction from the
+        // one below and one rugnux does not perform - it cuts nothing on a directional criterion. They
+        // go under the jfjoch local prefix with their own definition instead.
+        const auto &an = statistics.anisotropy;
+        if (an.n_reflections > 0 && an.n_cells > 0 && std::isfinite(an.delta_b)) {
+            out << "_reflns.pdbx_orthogonalization_convention  pdb\n";
+            for (int i = 0; i < 3; ++i) {
+                out << "_reflns.pdbx_aniso_B_tensor_eigenvalue_" << (i + 1) << "  "
+                    << Fmt(an.eigenvalue[i] - an.eigenvalue[2], 2) << "\n";
+                for (int j = 0; j < 3; ++j)
+                    out << "_reflns.pdbx_aniso_B_tensor_eigenvector_" << (i + 1) << "_ortho[" << (j + 1)
+                        << "]  " << Fmt(an.eigenvector[i][j], 4) << "\n";
+            }
+            out << "_reflns.jfjoch_aniso_delta_B          " << Fmt(an.delta_b, 2)
+                << "   # range of the principal components (A^2), fitted on intensities\n";
+            if (std::isfinite(an.delta_b_linear))
+                out << "_reflns.jfjoch_aniso_delta_B_linear   " << Fmt(an.delta_b_linear, 2)
+                    << "   # deltaB implied by the s^2 slope alone: what a Debye-Waller B accounts for\n";
+            out << "_reflns.jfjoch_aniso_shape            " << AnisotropyShapeCode(an.shape)
+                << "   # resolution signature of the directional deficit\n";
+            if (std::isfinite(an.floor))
+                out << "_reflns.jfjoch_aniso_floor            " << Fmt(an.floor, 3)
+                    << "   # deltaB this data set's own systematic error could manufacture (A^2)\n";
+            if (std::isfinite(an.significance))
+                out << "_reflns.jfjoch_aniso_significance     " << Fmt(an.significance, 2)
+                    << "   # deltaB(linear) / floor\n";
+            out << "_reflns.jfjoch_aniso_verdict          " << AnisotropyVerdictCode(an.verdict)
+                << "   # DETECTED / NOT_DETECTED / CANNOT_DETERMINE, at deltaB(linear)/floor > 3.5\n";
+            for (int i = 0; i < 3; ++i)
+                if (std::isfinite(an.d_min_axis[i]))
+                    out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << "          "
+                        << Fmt(an.d_min_axis[i], 2)
+                        << "   # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1)
+                        << (an.d_min_censored[i] ? "; at the edge of the measured data, so a bound" : "")
+                        << "\n";
+        }
         out << "#\n";
 
         // Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B
@@ -537,6 +584,297 @@ void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
     out.close();
 }
 
+namespace {
+
+// The detector position an observation is written at: a reflection too weak to have a measured
+// centroid still has a predicted one.
+float DetectorX(const Reflection &r) { return std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x; }
+float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y; }
+
+// Sum each rocking event into one full observation. A rotation reflection is integrated image by
+// image, so it arrives here as a run of partials over consecutive frames; the run is cut where the
+// 3D combine cuts it - same raw hkl, frames no further apart than MAX_FRAME_GAP - so the exported
+// file and rugnux's own merge see exactly the same events.
+// The parts are added, plainly, with their variances in quadrature, which is what every other
+// rotation program writes as a full. Nothing is divided by the partiality: FRACTIONCALC carries the
+// summed rocking-curve fraction, and an event the sweep cut short says so there rather than being
+// scaled up or dropped here.
+// The metadata a full carries is the partiality-weighted mean of its parts - the centroid of the
+// rocking curve, which is what XD/YD/ZD mean in an XDS file and what a scale model downstream reads
+// these columns for. An intensity weighting is the other candidate and it does not survive weak
+// data: the parts of a weak reflection scatter about zero, so on the crystals this was measured on
+// 8-41% of events have parts summing to zero or less and 10-21% of the intensity-weighted centroids
+// fall outside the event's own frame range. Partialities are the rocking curve's own weights and are
+// positive by construction.
+// min_partiality is the combine's own floor on the assembled full (--min-partiality, default 0.02,
+// 0 = off). An event that caught a thousandth of its rocking curve is not a measurement of that
+// reflection - rugnux does not merge one either - and writing it as a full hands the reading program
+// a whole observation whose intensity is noise and whose sigma is small, so it is weighted heavily.
+// On a crystal whose rocking curves span twenty frames those events are 7% of the file and cost
+// POINTLESS the point group. Everything above the floor is written with its honest FRACTIONCALC,
+// truncated edge-of-sweep events included, for the reader to cut where it wants.
+std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &outcomes,
+                                         double min_partiality) {
+    constexpr float MAX_FRAME_GAP = 2.0f;   // == RotationScaleMerge's: what makes one rocking event
+
+    // The sort key travels with the part instead of being read back through the pointer, the way the
+    // merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
+    // and an indirect compare is a cache miss on every one of them. The keys are the same values in
+    // the same order, so introsort makes the same comparisons and the same swaps and leaves the same
+    // order - which matters, because two parts can genuinely share (h,k,l) and image_number and the
+    // event sums below are floating point.
+    struct Part {
+        int32_t h, k, l;
+        float image_number;
+        const Reflection *r;
+    };
+    std::vector<Part> parts;
+    size_t n_parts = 0;
+    for (const auto &outcome : outcomes)
+        n_parts += outcome.reflections.size();
+    parts.reserve(n_parts);
+    for (const auto &outcome : outcomes)
+        for (const auto &r : outcome.reflections)
+            parts.push_back({r.h, r.k, r.l, r.image_number, &r});
+    std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &b) {
+        return std::tie(a.h, a.k, a.l, a.image_number) < std::tie(b.h, b.k, b.l, b.image_number);
+    });
+
+    std::vector<Reflection> fulls;
+    for (size_t i = 0; i < parts.size(); ) {
+        size_t j = i + 1;
+        while (j < parts.size() && parts[j].h == parts[i].h && parts[j].k == parts[i].k
+               && parts[j].l == parts[i].l
+               && parts[j].image_number - parts[j - 1].image_number <= MAX_FRAME_GAP)
+            ++j;
+
+        double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0;
+        double p_rlp = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0, p_delta_phi = 0.0, p_zeta = 0.0,
+               p_bkg = 0.0;
+        for (size_t m = i; m < j; ++m) {
+            const Reflection &r = *parts[m].r;
+            const double p = r.partiality;
+            sum_p += p;
+            sum_I += static_cast<double>(r.I) * r.rlp;
+            sum_var += static_cast<double>(r.sigma) * r.sigma * r.rlp * r.rlp;
+            sum_var_bkg += static_cast<double>(r.var_bkg) * r.rlp * r.rlp;
+            p_rlp += p * r.rlp;
+            p_frame += p * r.image_number;
+            p_x += p * DetectorX(r);
+            p_y += p * DetectorY(r);
+            p_delta_phi += p * r.delta_phi_deg;
+            p_zeta += p * r.zeta;
+            p_bkg += p * r.bkg;
+        }
+        Reflection full = *parts[i].r;
+        i = j;
+        if (sum_p < min_partiality)
+            continue;
+
+        // The Lorentz-polarization factor is applied by the writer, which multiplies I by rlp, so
+        // divide the event's own factor back out of the sums here. LP is the same geometry for every
+        // part of one event to a median 2e-4, so the file's I/LP is still the raw count sum.
+        full.rlp = static_cast<float>(p_rlp / sum_p);
+        full.I = static_cast<float>(sum_I / full.rlp);
+        full.sigma = static_cast<float>(std::sqrt(sum_var) / full.rlp);
+        full.var_bkg = static_cast<float>(sum_var_bkg / (static_cast<double>(full.rlp) * full.rlp));
+        full.partiality = static_cast<float>(sum_p);
+        full.image_number = static_cast<float>(p_frame / sum_p);
+        full.observed_x = static_cast<float>(p_x / sum_p);
+        full.observed_y = static_cast<float>(p_y / sum_p);
+        full.delta_phi_deg = static_cast<float>(p_delta_phi / sum_p);
+        full.zeta = static_cast<float>(p_zeta / sum_p);
+        full.bkg = static_cast<float>(p_bkg / sum_p);
+        fulls.push_back(full);
+    }
+    return fulls;
+}
+
+} // namespace
+
+void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes,
+                                 const UnitCell &unitCell,
+                                 const DiffractionExperiment &experiment,
+                                 bool sum_partials,
+                                 const std::string &filename) {
+    gemmi::Mtz mtz;
+    mtz.spacegroup = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1));
+    mtz.set_cell_for_all(unitCell);
+    mtz.title = "Rugnux unmerged observations";
+    mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
+    mtz.add_base();                 // the HKL_base dataset and the H K L columns
+    gemmi::Mtz::Dataset &ds = mtz.add_dataset("unmerged");
+    ds.crystal_name = experiment.GetSampleName();
+    ds.wavelength = experiment.GetWavelength_A();
+
+    // Every data column belongs to HKL_base and only the batches carry a dataset id, as in a
+    // POINTLESS file. Down to FLAG this is POINTLESS's own column set; the four after it are what
+    // rugnux measures beyond it - the offset of the reflection from the centre of its rocking curve,
+    // the Lorentz geometry of that curve, and the background that was subtracted.
+    mtz.add_column("M/ISYM", 'Y', 0, -1, false);
+    mtz.add_column("BATCH", 'B', 0, -1, false);
+    mtz.add_column("I", 'J', 0, -1, false);
+    mtz.add_column("SIGI", 'Q', 0, -1, false);
+    mtz.add_column("FRACTIONCALC", 'R', 0, -1, false);
+    mtz.add_column("XDET", 'R', 0, -1, false);
+    mtz.add_column("YDET", 'R', 0, -1, false);
+    mtz.add_column("ROT", 'R', 0, -1, false);
+    mtz.add_column("LP", 'R', 0, -1, false);
+    mtz.add_column("FLAG", 'I', 0, -1, false);
+    mtz.add_column("DELPHI", 'R', 0, -1, false);
+    mtz.add_column("ZETA", 'R', 0, -1, false);
+    mtz.add_column("BGMEAN", 'R', 0, -1, false);
+    mtz.add_column("BGVAR", 'R', 0, -1, false);
+
+    const auto gon = experiment.GetGoniometer();
+    // A scanning goniometer means rugnux integrated a rocking curve image by image, so each
+    // observation is one PART of a reflection unless the parts are summed here; on a still each is a
+    // whole measurement either way, and consecutive stills are different crystals, so there is
+    // nothing to sum there. Which of the two the file holds has to be declared in two places, because
+    // that is where POINTLESS and AIMLESS look for it: the batch header's data type, and the M flag
+    // that M/ISYM packs as 256*M + ISYM. Told nothing, they take each partial for a whole reflection
+    // and neither the symmetry determination nor the scaling survives it.
+    const bool scanning = gon && gon->IsScanning();
+    const bool partials = scanning && !sum_partials;
+    const float wedge_deg = gon ? gon->GetWedge_deg() : 0.0f;
+    const auto phi_start_deg = [&](float image_number) {
+        return gon ? gon->GetAngle_deg(image_number) : 0.0f;
+    };
+
+    // H K L are the ASU indices and M/ISYM says which symmetry operation (and which Friedel hand) got
+    // them there, so the index the reflection was actually measured at is recoverable - that is the
+    // crystal-frame information careless scales on, and what makes the file unmerged rather than a
+    // merge waiting to happen.
+    // I and SIGI are the integrated intensity with the Lorentz-polarization factor applied and
+    // nothing else, which is what IOBS means in every unmerged format (LP records the factor, so the
+    // raw counts are I/LP). LP is geometry, not a scale, and a program that reads this file has no
+    // way to recover it. The partiality is NOT divided out - that is a scale, FRACTIONCALC carries
+    // it, and every program this file is for wants to handle it its own way.
+    gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
+    std::set<int> batch_numbers;
+
+    // Lattice-centring absences are integrated on purpose - prediction runs in P so the space-group
+    // search can confirm or disprove the centring - but they are not written here. POINTLESS reads
+    // the interstitial nodes, takes the lattice for primitive and demotes the space group, and they
+    // roughly halve AIMLESS's <I/sigma>. Screw and glide absences ARE written: they are the evidence
+    // the space group was chosen on, so a program reading this file can check that choice, and
+    // removing them turns its test into an assumption. XDS and DIALS draw the line in the same place.
+    const char centering = mtz.spacegroup ? mtz.spacegroup->hm[0] : 'P';
+
+    const auto add_row = [&](const Reflection &r) {
+        if (systematic_absence(r.h, r.k, r.l, centering))
+            return;
+        std::array<int, 3> hkl{r.h, r.k, r.l};
+        const int isym = hkl_mover.move_to_asu(hkl);
+        // A summed full's image_number is its rocking-curve centroid, so this is the batch the
+        // reflection is centred on - which is what a batch means for a full everywhere else.
+        const int batch = 1 + static_cast<int>(std::lround(r.image_number));
+        batch_numbers.insert(batch);
+        mtz.data.push_back(static_cast<float>(hkl[0]));
+        mtz.data.push_back(static_cast<float>(hkl[1]));
+        mtz.data.push_back(static_cast<float>(hkl[2]));
+        mtz.data.push_back(static_cast<float>((partials ? 256 : 0) + isym));
+        mtz.data.push_back(static_cast<float>(batch));
+        mtz.data.push_back(r.I * r.rlp);
+        mtz.data.push_back(r.sigma * r.rlp);
+        mtz.data.push_back(r.partiality);
+        mtz.data.push_back(DetectorX(r));
+        mtz.data.push_back(DetectorY(r));
+        mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f);
+        mtz.data.push_back(r.rlp);
+        mtz.data.push_back(0.0f);            // FLAG: nothing here is a rejected observation
+        mtz.data.push_back(r.delta_phi_deg);
+        mtz.data.push_back(r.zeta);
+        mtz.data.push_back(r.bkg);
+        mtz.data.push_back(r.var_bkg);
+    };
+    if (scanning && sum_partials) {
+        for (const auto &r : SumRockingEvents(outcomes,
+                                              experiment.GetScalingSettings().GetMinPartiality()))
+            add_row(r);
+    } else {
+        for (const auto &outcome : outcomes)
+            for (const auto &r : outcome.reflections)
+                add_row(r);
+    }
+    mtz.nreflections = static_cast<int>(mtz.data.size() / mtz.columns.size());
+
+    // The batch header's orientation matrix is the crystal at rotation angle zero - each batch's own
+    // PHISTT is applied on top of it - but the lattice stored with an outcome is the crystal as it
+    // stood on that image. Turn the first indexed one back by its own angle to get the orientation of
+    // the sweep, which is the one matrix POINTLESS also writes into every batch.
+    std::optional<CrystalLattice> lattice_at_zero;
+    std::optional<float> mosaicity_deg;
+    for (const auto &outcome : outcomes) {
+        if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f)
+            continue;
+        const float mid_deg = phi_start_deg(outcome.reflections.front().image_number) + wedge_deg / 2.0f;
+        lattice_at_zero = gon ? outcome.latt.Multiply(gon->GetTransformationAngle(mid_deg)) : outcome.latt;
+        mosaicity_deg = outcome.mosaicity_deg;
+        break;
+    }
+
+    // The batch header is written in the "Cambridge" frame - z along the principal rotation axis, x
+    // along the beam - while the jfjoch lab frame has the beam along +z, so the two are related by a
+    // rotation. These three lab-frame vectors are the Cambridge axes; a lab vector's components in
+    // that frame are its dot products with them. A still has no rotation axis, and any axis across
+    // the beam then defines the frame just as consistently.
+    // The axis is NEGATED: rugnux turns an observation made at angle phi back to phi = 0 by rotating
+    // it by +phi about the goniometer axis, so the crystal itself turns by -phi about it, while the
+    // MTZ batch header's scan axis is the one a batch's own increasing PHI turns the crystal about.
+    const Coord beam = experiment.GetDiffractionGeometry().GetScatteringVector().Normalize();
+    const Coord z_cam = gon ? -gon->GetAxis().Normalize() : Coord(0, 1, 0);
+    const Coord y_cam = (z_cam % beam).Normalize();
+    const Coord x_cam = (y_cam % z_cam).Normalize();
+
+    gemmi::Mtz::Batch batch;
+    batch.title = "Rugnux";
+    batch.axes.emplace_back("PHI");
+    batch.set_dataset_id(ds.id);
+    batch.ints[12] = 1;   // ncryst
+    batch.ints[14] = partials ? 1 : 2;   // ldtype: oscillation data (2D spots) / area detector (3D)
+    batch.ints[15] = 1;   // jsaxs: the goniostat scan axis
+    batch.ints[17] = 1;   // ngonax
+    batch.ints[19] = 1;   // ndet
+    batch.set_cell(unitCell);
+    if (lattice_at_zero) {
+        // Orientation matrix U, built from the reciprocal axes and stored column by column in
+        // Cambridge components, as gemmi's XDS_ASCII converter builds it (gemmi/xds2mtz.hpp).
+        const Coord ar = lattice_at_zero->Astar().Normalize();
+        const Coord cr = (ar % lattice_at_zero->Bstar()).Normalize();
+        const Coord u[3] = {ar, cr % ar, cr};
+        for (int i = 0; i < 3; i++) {
+            batch.floats[6 + 3 * i] = u[i] * x_cam;
+            batch.floats[7 + 3 * i] = u[i] * y_cam;
+            batch.floats[8 + 3 * i] = u[i] * z_cam;
+        }
+    }
+    batch.floats[21] = mosaicity_deg.value_or(0.0f);   // crydat(0), the reflecting range
+    batch.floats[40] = 1.0f;    // scanax = [0, 0, 1]: the rotation axis IS z in the Cambridge frame
+    batch.floats[47] = wedge_deg;
+    batch.floats[61] = 1.0f;    // e1 = scanax, the only goniostat axis
+    batch.floats[80] = -1.0f;   // idealised source vector, antiparallel to the beam
+    batch.floats[83] = -(beam * x_cam);   // s0, the source vector of the geometry as it really stands
+    batch.floats[84] = -(beam * y_cam);
+    batch.floats[85] = -(beam * z_cam);
+    batch.set_wavelength(experiment.GetWavelength_A());
+    batch.floats[111] = experiment.GetDetectorDistance_mm();
+    batch.floats[113] = 1.0f;                                             // detector limits, in pixels
+    batch.floats[114] = static_cast<float>(experiment.GetXPixelsNum());
+    batch.floats[115] = 1.0f;
+    batch.floats[116] = static_cast<float>(experiment.GetYPixelsNum());
+
+    for (const int number : batch_numbers) {
+        batch.number = number;
+        batch.floats[36] = phi_start_deg(static_cast<float>(number - 1));   // phistt
+        batch.floats[37] = batch.floats[36] + wedge_deg;                    // phiend
+        mtz.batches.push_back(batch);
+    }
+
+    mtz.sort(5);   // by H K L M/ISYM BATCH, the order POINTLESS leaves an unmerged file in
+    mtz.write_to_file(filename);
+}
+
 void WriteReflections(const std::vector<MergedReflection> &reflections,
                       const UnitCell &unitCell,
                       const DiffractionExperiment &experiment,