Jungfraujoch/tools/rugnux_cli.cpp

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <csignal>
#include <getopt.h>
#include <iostream>
#include <sstream>
#include <string>
#include <thread>
#include <vector>

#include "../reader/JFJochHDF5Reader.h"
#include "../common/Logger.h"
#include "../common/Definitions.h"
#include "../common/DiffractionExperiment.h"
#include "../common/PixelMask.h"
#include "../common/print_license.h"
#include "../image_analysis/LoadFCalcFromMtz.h"
#include "../image_analysis/UpdateReflectionResolution.h"
#include "../image_analysis/WriteReflections.h"
#include "../image_analysis/scale_merge/Merge.h"
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
#include "../rugnux/Rugnux.h"

// Default rot3d per-frame scale-G smoothing range (XDS DELPHI-like), in degrees of rotation.
constexpr double SMOOTH_G_DEFAULT_DEG = 5.0;

// Default rot3d per-observation merge outlier rejection (sigma from the per-reflection median, XDS-like).
// Neutral-or-better across the test crystals; pass --reject-outliers 0 to disable.
constexpr double REJECT_OUTLIERS_DEFAULT_NSIGMA = 6.0;

void print_usage() {
    std::cout << "Usage rugnux {<options>} <input.h5>" << std::endl;
    std::cout << "Options:" << std::endl;
    std::cout << "   -o, --output-prefix <txt>        Output file prefix (default: output)" << std::endl;
    std::cout << "   -N, --threads <num>              Number of threads (default: all hardware threads)" << std::endl;
    std::cout << "   -s, --start-image <num>          Start image number (default: 0)" << std::endl;
    std::cout << "   -e, --end-image <num>            End image number (default: all)" << std::endl;
    std::cout << "   -t, --stride <num>               Image stride (default: 1)" << std::endl;
    std::cout << "   -v, --verbose                    Verbose output" << std::endl;
    std::cout << std::endl;

    std::cout << "  Modes (default: full analysis - spot finding, indexing, integration and merging)" << std::endl;
    std::cout << "   --azint-only                     Only run azimuthal integration (no spot finding/indexing); writes <prefix>_process.h5" << std::endl;
    std::cout << "   --scale                          Only re-scale/merge the already-integrated reflections in <input> (no re-integration)" << std::endl;
    std::cout << std::endl;

    std::cout << "  Spot finding" << std::endl;
    std::cout << "   --spot-sigma <num>               Noise sigma level for spot finding (default: 3.0)" << std::endl;
    std::cout << "   --spot-threshold <num>           Photon count threshold for spot finding (default: 10)" << std::endl;
    std::cout << "   --spot-high-resolution <num>     High resolution limit for spot finding (default: 1.5)" << std::endl;
    std::cout << "   --max-spots <num>                Max spot count (default: 250)" << std::endl;
    std::cout << "   --detect-ice-rings[=on|off]      Flag ice-ring spots (de-prioritised in indexing) and exclude ice-ring reflections from scaling/merging; overrides the dataset/master-file setting (default: use dataset value)" << std::endl;
    std::cout << std::endl;

    std::cout << "  Indexing" << std::endl;
    std::cout << "   (A dataset with a rotation goniometer axis is processed as rotation data by default; use --force-still to override)" << std::endl;
    std::cout << "   --force-still                    Process a rotation (goniometer) dataset as independent stills (still indexing + per-image ScaleOnTheFly) instead of rotation" << std::endl;
    std::cout << "   -R, --two-pass-rotation[=num]    Two-pass offline rotation indexing (default for goniometer data; optional first-pass image count, default: 100)" << std::endl;
    std::cout << "   --single-pass-rotation[=num]     Use online-like single-pass rotation indexing (optional: min angular range deg)" << std::endl;
    std::cout << "   --redo-rotation-spots            Redo spot finding for two-pass rotation indexing" << std::endl;
    std::cout << "   --force-rotation-lattice <vec>   Force rotation indexer with external lattice (in Angstrom) : \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, skips first pass)" << std::endl;
    std::cout << "   -X, --indexing-algorithm <txt>   Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl;
    std::cout << "   -S, --space-group <num>          Space group number - used for both indexing and scaling" << std::endl;
    std::cout << "   -C, --unit-cell <cell>           Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl;
    std::cout << "   -r, --refine <txt>               Geometry refinement algorithm (none|orientation|beam_and_lattice)" << std::endl;
    std::cout << std::endl;

    std::cout << "  Scaling and merging (on by default)" << std::endl;
    std::cout << "   --no-merge                      Skip scaling and merging; write only the per-image _process.h5" << std::endl;
    std::cout << "   --scale-fulls                   rot3d: after the 3D combine, refit a per-frame scale on the fulls (XDS order, Unity model). Default ON for rot3d" << std::endl;
    std::cout << "   --no-scale-fulls                Disable the rot3d scale-fulls refit (it is on by default for rot3d)" << std::endl;
    std::cout << "   --write-process-h5              Also write the (large) _process.h5 when merging (default: only .mtz/.cif when merging)" << std::endl;
    std::cout << "   --smooth-g[=deg]                rot3d: smooth per-frame scale G over a deg-degree rotation range (XDS DELPHI-like) before the combine (default: 5 for rot3d; 0 = off)" << std::endl;
    std::cout << "   -A, --anomalous                 Anomalous mode (don't merge Friedel pairs)" << std::endl;
    std::cout << "   -B, --refine-bfactor            Refine per image B-factor" << std::endl;
    std::cout << "   --scaling-high-resolution <num> High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl;
    std::cout << "   --resolution-cutoff <txt>       Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl;
    std::cout << "   --resolution-cc-target <num>    CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl;
    std::cout << "   --resolution-shells <num>       Number of resolution shells in the reported statistics table (default: 10)" << std::endl;
    std::cout << "   --min-partiality <num>          Minimum partiality to accept reflection (default: 0.02)" << std::endl;
    std::cout << "   --capture-uncertainty <num>     rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl;
    std::cout << "   --mosaicity <num>               Diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed" << std::endl;
    std::cout << "   --reject-outliers <num>         Per-observation merge outlier rejection, N sigma from the per-reflection median (default: 6 for rot3d, XDS/DIALS-style; 0 = off)" << std::endl;
    std::cout << "   --reject-delta-cchalf <num>     Per-crystal CC1/2-delta rejection: drop images with deltaCChalf below mean - N*stddev (default: off; e.g. 2.5)" << std::endl;
    std::cout << "   --min-image-cc <num>            Per-image CC limit in percent (default: no limit)" << std::endl;
    std::cout << "   --scaling-iterations <num>      Number of scaling iterations with no reference data (default: 3)" << std::endl;
    std::cout << "   --scaling-output <txt>              Output format for scaling results mtz|cif|txt (default: cif)" << std::endl;
    std::cout << "   -z, --reference-mtz <file>      Reference MTZ file" << std::endl;
    std::cout << "   --reference-column <label>      Reference MTZ column to use (default: auto - F-model, else IMEAN/I/...)" << std::endl;
    std::cout << std::endl;

    std::cout << "  Integration" << std::endl;
    std::cout << "   --bandwidth <num>               Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
    std::cout << "   --integration-radius <r>        Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
    std::cout << "   --integrator <txt>              Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
    std::cout << "   -q, --azim-q-spacing <num>      Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;
    std::cout << "   --azim-min-q <num>              Azimuthal-integration minimum Q (1/A)" << std::endl;
    std::cout << "   --azim-max-q <num>              Azimuthal-integration maximum Q (1/A)" << std::endl;
    std::cout << "   --azim-phi-bins <num>           Number of azimuthal (phi) bins (default: 1)" << std::endl;
    std::cout << "   --polarization-correction <on|off>  Enable/disable azimuthal polarization correction" << std::endl;
    std::cout << "   --solid-angle-correction <on|off>   Enable/disable azimuthal solid angle correction" << std::endl;
    std::cout << std::endl;

    std::cout << "  Geometry overrides (defaults taken from the input file)" << std::endl;
    std::cout << "   --beam-x <num>                   Beam center X (pixel)" << std::endl;
    std::cout << "   --beam-y <num>                   Beam center Y (pixel)" << std::endl;
    std::cout << "   --detector-distance <num>        Detector distance (mm)" << std::endl;
    std::cout << "   --wavelength <num>               Wavelength (A)" << std::endl;
    std::cout << "   --rot1 <num>                     PONI rotation 1 (rad)" << std::endl;
    std::cout << "   --rot2 <num>                     PONI rotation 2 (rad)" << std::endl;
    std::cout << "   --polarization <num>             Polarization factor" << std::endl;
    }

enum {
    OPT_SPOT_SIGMA = 1000,
    OPT_SPOT_THRESHOLD,
    OPT_SPOT_RESOLUTION,
    OPT_MAX_SPOTS,
    OPT_MIN_PARTIALITY,
    OPT_MIN_IMAGE_CC,
    OPT_SCALING_ITERATIONS,
    OPT_SCALING_HIGH_RESOLUTION,
    OPT_RESOLUTION_CUTOFF,
    OPT_RESOLUTION_CC_TARGET,
    OPT_RESOLUTION_SHELLS,
    OPT_SCALING_OUTPUT,
    OPT_SINGLE_PASS_ROTATION,
    OPT_REDO_ROTATION_SPOTS,
    OPT_FORCE_ROTATION_LATTICE,
    OPT_BANDWIDTH,
    OPT_INTEGRATION_RADIUS,
    OPT_REJECT_OUTLIERS,
    OPT_REJECT_DELTA_CCHALF,
    OPT_REFERENCE_COLUMN,
    OPT_DUMP_OBSERVATIONS,
    OPT_INTEGRATOR,
    OPT_SCALE_FULLS,
    OPT_CAPTURE_UNCERTAINTY,
    OPT_MOSAICITY,
    OPT_SMOOTH_G,
    OPT_DETECT_ICE_RINGS,
    OPT_NO_SCALE_FULLS,
    OPT_WRITE_PROCESS_H5,
    OPT_FORCE_STILL,
    OPT_AZIM_MIN_Q,
    OPT_AZIM_MAX_Q,
    OPT_AZIM_PHI_BINS,
    OPT_AZINT_ONLY,
    OPT_SCALE,
    OPT_NO_MERGE,
    OPT_POLARIZATION_CORRECTION,
    OPT_SOLID_ANGLE_CORRECTION,
    OPT_BEAM_X,
    OPT_BEAM_Y,
    OPT_DETECTOR_DISTANCE,
    OPT_WAVELENGTH,
    OPT_ROT1,
    OPT_ROT2,
    OPT_POLARIZATION
};

static option long_options[] = {
    {"verbose", no_argument, nullptr, 'v'},
    {"output-prefix", required_argument, nullptr, 'o'},
    {"threads", required_argument, nullptr, 'N'},
    {"start-image", required_argument, nullptr, 's'},
    {"end-image", required_argument, nullptr, 'e'},
    {"stride", required_argument, nullptr, 't'},
    {"indexing-algorithm", required_argument, nullptr, 'X'},
    {"unit-cell", required_argument, nullptr, 'C'},
    {"reference-mtz", required_argument, nullptr, 'z'},
    {"reference-column", required_argument, nullptr, OPT_REFERENCE_COLUMN},
    {"dump-observations", required_argument, nullptr, OPT_DUMP_OBSERVATIONS},
    {"space-group", required_argument, nullptr, 'S'},
    {"anomalous", no_argument, nullptr, 'A'},
    {"refine-bfactor", no_argument, nullptr, 'B'},
    {"azint-only", no_argument, nullptr, OPT_AZINT_ONLY},
    {"scale", no_argument, nullptr, OPT_SCALE},
    {"no-merge", no_argument, nullptr, OPT_NO_MERGE},
    {"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS},
    {"no-scale-fulls", no_argument, nullptr, OPT_NO_SCALE_FULLS},
    {"write-process-h5", no_argument, nullptr, OPT_WRITE_PROCESS_H5},
    {"smooth-g", optional_argument, nullptr, OPT_SMOOTH_G},
    {"refine", required_argument, nullptr, 'r'},

    {"two-pass-rotation", optional_argument, nullptr, 'R'},
    {"single-pass-rotation", optional_argument, nullptr, OPT_SINGLE_PASS_ROTATION},
    {"force-still", no_argument, nullptr, OPT_FORCE_STILL},
    {"azim-q-spacing", required_argument, nullptr, 'q'},
    {"azim-min-q", required_argument, nullptr, OPT_AZIM_MIN_Q},
    {"azim-max-q", required_argument, nullptr, OPT_AZIM_MAX_Q},
    {"azim-phi-bins", required_argument, nullptr, OPT_AZIM_PHI_BINS},
    {"polarization-correction", required_argument, nullptr, OPT_POLARIZATION_CORRECTION},
    {"solid-angle-correction", required_argument, nullptr, OPT_SOLID_ANGLE_CORRECTION},
    {"beam-x", required_argument, nullptr, OPT_BEAM_X},
    {"beam-y", required_argument, nullptr, OPT_BEAM_Y},
    {"detector-distance", required_argument, nullptr, OPT_DETECTOR_DISTANCE},
    {"wavelength", required_argument, nullptr, OPT_WAVELENGTH},
    {"rot1", required_argument, nullptr, OPT_ROT1},
    {"rot2", required_argument, nullptr, OPT_ROT2},
    {"polarization", required_argument, nullptr, OPT_POLARIZATION},
    {"redo-rotation-spots", no_argument, nullptr, OPT_REDO_ROTATION_SPOTS},
    {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE},


    {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA},
    {"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD},
    {"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION},
    {"max-spots", required_argument, nullptr, OPT_MAX_SPOTS},
    {"min-partiality", required_argument, nullptr, OPT_MIN_PARTIALITY},
    {"capture-uncertainty", required_argument, nullptr, OPT_CAPTURE_UNCERTAINTY},
    {"mosaicity", required_argument, nullptr, OPT_MOSAICITY},
    {"min-image-cc", required_argument, nullptr, OPT_MIN_IMAGE_CC},
    {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS},
    {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION},
    {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF},
    {"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET},
    {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS},
    {"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT},
    {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
    {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
    {"integrator", required_argument, nullptr, OPT_INTEGRATOR},
    {"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS},
    {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
    {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF},
    {nullptr, 0, nullptr, 0}
};

void trim_in_place(std::string &t) {
    size_t b = 0;
    while (b < t.size() && std::isspace(static_cast<unsigned char>(t[b]))) b++;
    size_t e = t.size();
    while (e > b && std::isspace(static_cast<unsigned char>(t[e - 1]))) e--;
    t = t.substr(b, e - b);
};

bool parse_float_strict(const std::string &t, float &out) {
    try {
        size_t idx = 0;
        out = std::stof(t, &idx);
        return idx == t.size();
    } catch (...) {
        return false;
    }
};

// Parse a required numeric option argument, or print a clear error and exit. getopt hands option
// arguments over as raw C strings; atof() silently returns 0 on non-numeric input (so a typo like
// "--beam-x 21OO" would move the beam to pixel 0) and std::stod() throws std::invalid_argument,
// which terminates the program. This rejects both, plus trailing garbage like "1.5foo".
double parse_double_arg(const char *arg, const char *option_name, Logger &logger) {
    std::string s = arg ? arg : "";
    trim_in_place(s);
    if (!s.empty()) {
        try {
            size_t idx = 0;
            double out = std::stod(s, &idx);
            if (idx == s.size())
                return out;
        } catch (...) {}
    }
    logger.Error("Invalid numeric value for {}: '{}'", option_name, arg ? arg : "<null>");
    print_usage();
    exit(EXIT_FAILURE);
}

float parse_float_arg(const char *arg, const char *option_name, Logger &logger) {
    return static_cast<float>(parse_double_arg(arg, option_name, logger));
}

bool parse_on_off(const char *arg, bool &out) {
    std::string s = arg ? arg : "";
    std::transform(s.begin(), s.end(), s.begin(),
                   [](unsigned char c) { return static_cast<char>(std::tolower(c)); });
    if (s == "on" || s == "1" || s == "true" || s == "yes") {
        out = true;
        return true;
    }
    if (s == "off" || s == "0" || s == "false" || s == "no") {
        out = false;
        return true;
    }
    return false;
}

std::optional<UnitCell> parse_unit_cell_arg(const char *arg) {
    if (!arg)
        return std::nullopt;

    std::string s(arg);


    trim_in_place(s);

    if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
        s = s.substr(1, s.size() - 2);
        trim_in_place(s);
    }

    std::vector<std::string> parts;
    parts.reserve(6);
    size_t start = 0;
    while (true) {
        size_t pos = s.find(',', start);
        if (pos == std::string::npos) {
            parts.push_back(s.substr(start));
            break;
        }
        parts.push_back(s.substr(start, pos - start));
        start = pos + 1;
    }

    if (parts.size() != 6)
        return std::nullopt;



    UnitCell uc{};
    if (!parse_float_strict(parts[0], uc.a)) return std::nullopt;
    if (!parse_float_strict(parts[1], uc.b)) return std::nullopt;
    if (!parse_float_strict(parts[2], uc.c)) return std::nullopt;
    if (!parse_float_strict(parts[3], uc.alpha)) return std::nullopt;
    if (!parse_float_strict(parts[4], uc.beta)) return std::nullopt;
    if (!parse_float_strict(parts[5], uc.gamma)) return std::nullopt;

    return uc;
}

std::optional<CrystalLattice> parse_lattice_arg(const char *arg) {
    if (!arg)
        return std::nullopt;

    std::string s(arg);
    trim_in_place(s);

    if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
        s = s.substr(1, s.size() - 2);
        trim_in_place(s);
    }

    std::vector<std::string> parts;
    parts.reserve(9);
    size_t start = 0;
    while (true) {
        size_t pos = s.find(',', start);
        if (pos == std::string::npos) {
            parts.push_back(s.substr(start));
            break;
        }
        parts.push_back(s.substr(start, pos - start));
        start = pos + 1;
    }

    if (parts.size() != 9)
        return std::nullopt;

    std::vector<float> vals(9);
    for (int i = 0; i < 9; i++) {
        if (!parse_float_strict(parts[i], vals[i]))
            return std::nullopt;
    }

    return CrystalLattice(vals);
}

// Shared offline-output settings for the _process.h5 / reflection writer, used by both the --scale
// path and the full-analysis path (each then sets its own space group and images-per-trigger).
void configure_offline_output(DiffractionExperiment &experiment, const std::string &output_prefix) {
    experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4);
    experiment.FilePrefix(output_prefix);
    experiment.Mode(DetectorMode::Standard);   // full image analysis
    experiment.PixelSigned(true);
    experiment.OverwriteExistingFiles(true);
    experiment.PolarizationFactor(0.99);
    experiment.SetFileWriterFormat(FileWriterFormat::NXmxLegacy);
    experiment.NumTriggers(1);
}

namespace {
    std::atomic<Rugnux *> g_active_process{nullptr};
    void handle_sigint(int) {
        if (auto *p = g_active_process.load())
            p->Cancel();
    }
}

int main(int argc, char **argv) {
    for (int i = 0; i < argc; i++) {
        std::cout << argv[i] << " ";
    }
    std::cout << std::endl << std::endl;


    RegisterHDF5Filter();

    print_license("rugnux");

    Logger logger("rugnux");

    std::string input_file;
    std::string output_prefix = "output";
    int nthreads = 0;   // 0 = auto: resolved to all hardware threads after parsing (see below)
    int start_image = 0;
    int end_image = -1; // -1 indicates process until end
    int image_stride = 1;

    bool verbose = false;
    bool azint_only = false;      // --azint-only: azimuthal integration only (no spot finding/indexing)
    bool scale_only = false;      // --scale: re-scale/merge stored reflections only (no re-integration)
    bool rotation_indexing = false;
    bool force_still = false;     // --force-still: process a rotation dataset as stills (indexing + scaling)
    bool two_pass_rotation = true;
    bool reuse_rotation_spots = true;
    int rotation_indexing_image_count = 100;
    std::optional<float> rotation_indexing_range;
    bool run_scaling = true;      // merge is on by default; --no-merge turns it off
    std::optional<bool> scale_fulls_arg;   // --scale-fulls / --no-scale-fulls; default on for rot3d
    bool write_process_h5_flag = false;    // --write-process-h5; also write _process.h5 when merging
    std::optional<bool> detect_ice_rings;  // --detect-ice-rings[=on|off]; unset => use the dataset (file) value
    std::optional<float> min_q, max_q, q_spacing;   // azimuthal integration range / -q spacing (1/A)
    std::optional<int32_t> azimuthal_bins;          // --azimuthal-bins
    std::optional<bool> polarization_correction;    // --polarization-correction (azimuthal integration)
    std::optional<bool> solid_angle_correction;     // --solid-angle-correction (azimuthal integration)

    // Geometry overrides (default: keep the value stored in the input file)
    std::optional<float> beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, polarization_factor;
    std::optional<double> smooth_g_deg_arg;   // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise
    bool anomalous_mode = false;
    std::optional<int64_t> space_group_number;
    std::optional<UnitCell> fixed_reference_unit_cell;
    std::optional<int64_t> max_spot_count_override;
    float sigma_spot_finding = 3.0;
    int64_t photon_count_threshold_spot_finding = 10;
    bool refine_bfactor = false;
    std::string ref_mtz;
    std::string ref_column;
    std::string dump_observations;   // diagnostic: dump unmerged -P rot3d fulls to this path
    double min_partiality = 0.02;
    std::optional<double> capture_uncertainty_arg;   // explicit --capture-uncertainty; default depends on rot3d
    std::optional<double> forced_mosaicity_arg;      // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed
    double min_image_cc = 0.0;
    int64_t scaling_iter = 3;
    std::optional<CrystalLattice> forced_rotation_lattice;

    std::optional<float> bandwidth_fwhm; // relative FWHM of dlambda/lambda

    IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
    GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;

    std::optional<IntensityFormat> intensity_format;   // --scaling-output override; default lives in ScalingSettings (mmCIF)

    float d_min_spot_finding = 1.5;
    std::optional<float> d_min_scale_merge;
    std::optional<ResolutionCutoffMethod> resolution_cutoff_method;   // --resolution-cutoff cc-logistic|off
    std::optional<double> resolution_cc_target;                      // --resolution-cc-target
    std::optional<int> report_shell_count;                           // --resolution-shells
    std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
    std::optional<IntegratorMode> integrator_mode;     // --integrator boxsum|gaussian|empirical
    std::optional<double> outlier_reject_nsigma;       // merge per-observation outlier rejection
    std::optional<double> delta_cchalf_nsigma;         // per-crystal CC1/2-delta rejection

    if (argc == 1) {
        print_usage();
        exit(EXIT_FAILURE);
    }

    int opt;
    int option_index = 0;
    const char *short_opts = "vo:N:s:e:t:R::X:C:z:FABS:r:q:";

  while ((opt = getopt_long(argc, argv, short_opts, long_options, &option_index)) != -1) {
    switch (opt) {
        case 'o':
            output_prefix = optarg;
            break;
        case 'v':
            verbose = true;
            break;
        case 'N':
            nthreads = atoi(optarg);
            break;
        case 's':
            start_image = atoi(optarg);
            break;
        case 'e':
            end_image = atoi(optarg);
            break;
        case 't':
            image_stride = atoi(optarg);
            break;
        case 'R':
            if (rotation_indexing) {
                logger.Error("Rotation indexing already enabled");
                exit(EXIT_FAILURE);
            }
            rotation_indexing = true;
            two_pass_rotation = true;
            if (optarg)
                rotation_indexing_image_count = atoi(optarg);

            break;
        case OPT_SINGLE_PASS_ROTATION:
            if (rotation_indexing) {
                logger.Error("Rotation indexing already enabled");
                exit(EXIT_FAILURE);
            }
            rotation_indexing = true;
            two_pass_rotation = false;

            if (optarg)
                rotation_indexing_range = atof(optarg);
            break;
        case OPT_REDO_ROTATION_SPOTS:
            reuse_rotation_spots = false;
            break;
        case OPT_FORCE_ROTATION_LATTICE: {
            if (rotation_indexing) {
                logger.Error("Rotation indexing already enabled");
                exit(EXIT_FAILURE);
            }
            rotation_indexing = true;

            auto latt = parse_lattice_arg(optarg);
            if (!latt.has_value()) {
                logger.Error(
                    "Invalid rotation lattice. Expected: \"a0x,a0y,a0z,a1x,a1y,a1z,a2x,a2y,a2z\" (9 floats, comma-separated). Got: {}",
                    optarg ? optarg : "<null>");
                print_usage();
                exit(EXIT_FAILURE);
            }
            forced_rotation_lattice = latt;
            auto uc = latt->GetUnitCell();
            logger.Info(
                "Forced rotation lattice set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}",
                uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
            break;
        }
        case 'X': {
            std::string alg = optarg ? optarg : "";
            std::transform(alg.begin(), alg.end(), alg.begin(),
                           [](unsigned char c) { return static_cast<char>(std::tolower(c)); });

            if (alg == "ffbidx")
                indexing_algorithm = IndexingAlgorithmEnum::FFBIDX;
            else if (alg == "fft")
                indexing_algorithm = IndexingAlgorithmEnum::FFT;
            else if (alg == "fftw")
                indexing_algorithm = IndexingAlgorithmEnum::FFTW;
            else if (alg == "auto")
                indexing_algorithm = IndexingAlgorithmEnum::Auto;
            else if (alg == "none")
                indexing_algorithm = IndexingAlgorithmEnum::None;
            else {
                logger.Error("Invalid indexing algorithm: {}", alg);
                print_usage();
                exit(EXIT_FAILURE);
            }
            break;
        }
        case 'r': {
            std::string alg = optarg ? optarg : "";
            std::transform(alg.begin(), alg.end(), alg.begin(),
                           [](unsigned char c) { return static_cast<char>(std::tolower(c)); });
            if (alg == "none")
                refinement_algorithm = GeomRefinementAlgorithmEnum::None;
            else if (alg == "beam_and_lattice")
                refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
            else if (alg == "orientation")
                refinement_algorithm = GeomRefinementAlgorithmEnum::OrientationOnly;
            else {
                logger.Error("Invalid geom refinement algorithm: {}", alg);
                print_usage();
                exit(EXIT_FAILURE);
            }
            break;
        }
        case 'C': {
            auto uc = parse_unit_cell_arg(optarg);
            if (!uc.has_value()) {
                logger.Error(
                    "Invalid unit cell. Expected: \"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}",
                    optarg ? optarg : "<null>");
                print_usage();
                exit(EXIT_FAILURE);
            }
            fixed_reference_unit_cell = uc;
            logger.Info(
                "Fixed reference unit cell set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}",
                uc->a, uc->b, uc->c, uc->alpha, uc->beta, uc->gamma);
            break;
        }
        case 'z':
            ref_mtz = optarg;
            break;
        case OPT_REFERENCE_COLUMN:
            ref_column = optarg;
            break;
        case OPT_DUMP_OBSERVATIONS:
            dump_observations = optarg;
            break;
        case 'F':
            indexing_algorithm = IndexingAlgorithmEnum::FFT;
            break;
        case 'A':
            anomalous_mode = true;
            break;
        case 'B':
            refine_bfactor = true;
            break;
        case 'S':
            space_group_number = atoi(optarg);
            break;
        case OPT_SPOT_SIGMA:
            sigma_spot_finding = atof(optarg);
            logger.Info("Noise threshold level for spot finding set to {:.2f} sigma", sigma_spot_finding);
            break;
        case OPT_SPOT_THRESHOLD:
            photon_count_threshold_spot_finding = atoi(optarg);
            logger.Info("Photon-count threshold level for spot finding set to {:d}",
                        photon_count_threshold_spot_finding);
            break;
        case OPT_SPOT_RESOLUTION:
            d_min_spot_finding = atof(optarg);
            logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min_spot_finding);
            break;
        case OPT_MAX_SPOTS:
            max_spot_count_override = atoll(optarg);
            logger.Info("Max spot count overridden to {}", max_spot_count_override.value());
            break;
        case OPT_AZINT_ONLY:
            azint_only = true;
            break;
        case OPT_SCALE:
            scale_only = true;
            break;
        case OPT_NO_MERGE:
            run_scaling = false;
            break;
        case OPT_SCALE_FULLS:
            scale_fulls_arg = true;
            break;
        case OPT_NO_SCALE_FULLS:
            scale_fulls_arg = false;
            break;
        case OPT_DETECT_ICE_RINGS:
            if (optarg == nullptr || strcmp(optarg, "on") == 0)
                detect_ice_rings = true;
            else if (strcmp(optarg, "off") == 0)
                detect_ice_rings = false;
            else {
                logger.Error("Invalid --detect-ice-rings value: {} (expected on|off)", optarg);
                print_usage();
                exit(EXIT_FAILURE);
            }
            break;
        case OPT_WRITE_PROCESS_H5:
            write_process_h5_flag = true;
            break;
        case OPT_SMOOTH_G:
            smooth_g_deg_arg = optarg ? parse_double_arg(optarg, "--smooth-g", logger) : SMOOTH_G_DEFAULT_DEG;
            break;
        case OPT_MIN_PARTIALITY:
            min_partiality = parse_double_arg(optarg, "--min-partiality", logger);
            break;
        case OPT_CAPTURE_UNCERTAINTY:
            capture_uncertainty_arg = parse_double_arg(optarg, "--capture-uncertainty", logger);
            break;
        case OPT_MOSAICITY:
            forced_mosaicity_arg = parse_double_arg(optarg, "--mosaicity", logger);
            break;
        case OPT_INTEGRATION_RADIUS:
            integration_radius_arg = optarg;
            break;
        case OPT_INTEGRATOR:
            if (strcmp(optarg, "boxsum") == 0) integrator_mode = IntegratorMode::BoxSum;
            else if (strcmp(optarg, "gaussian") == 0) integrator_mode = IntegratorMode::ProfileGaussian;
            else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical;
            else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; }
            break;
        case OPT_REJECT_OUTLIERS:
            outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger);
            break;
        case OPT_REJECT_DELTA_CCHALF:
            delta_cchalf_nsigma = parse_double_arg(optarg, "--reject-delta-cchalf", logger);
            break;
        case OPT_MIN_IMAGE_CC:
            min_image_cc = parse_double_arg(optarg, "--min-image-cc", logger);
            break;
        case OPT_SCALING_HIGH_RESOLUTION:
            d_min_scale_merge = atof(optarg);
            break;
        case OPT_RESOLUTION_CUTOFF:
            if (strcmp(optarg, "cc-logistic") == 0)
                resolution_cutoff_method = ResolutionCutoffMethod::CCHalfLogistic;
            else if (strcmp(optarg, "off") == 0)
                resolution_cutoff_method = ResolutionCutoffMethod::Off;
            else {
                logger.Error("Invalid --resolution-cutoff value: {} (expected cc-logistic|off)", optarg);
                print_usage();
                exit(EXIT_FAILURE);
            }
            break;
        case OPT_RESOLUTION_CC_TARGET:
            resolution_cc_target = parse_double_arg(optarg, "--resolution-cc-target", logger);
            break;
        case OPT_RESOLUTION_SHELLS:
            report_shell_count = atoi(optarg);
            if (report_shell_count.value() < 1) {
                logger.Error("Invalid --resolution-shells value: {} (must be >= 1)", report_shell_count.value());
                exit(EXIT_FAILURE);
            }
            break;
        case OPT_SCALING_OUTPUT:
            if (strcmp(optarg, "mtz") == 0) {
                intensity_format = IntensityFormat::MTZ;
            } else if (strcmp(optarg, "cif") == 0) {
                intensity_format = IntensityFormat::mmCIF;
            } else if (strcmp(optarg, "txt") == 0) {
                intensity_format = IntensityFormat::Text;
            } else {
                logger.Error("Invalid intensity format: {}", optarg);
                exit(EXIT_FAILURE);
            }
            break;
        case OPT_FORCE_STILL:
            force_still = true;
            break;
        case 'q':
            q_spacing = atof(optarg);
            break;
        case OPT_AZIM_MIN_Q:
            min_q = atof(optarg);
            break;
        case OPT_AZIM_MAX_Q:
            max_q = atof(optarg);
            break;
        case OPT_AZIM_PHI_BINS:
            azimuthal_bins = atoi(optarg);
            break;
        case OPT_POLARIZATION_CORRECTION: {
            bool value;
            if (!parse_on_off(optarg, value)) {
                logger.Error("Invalid polarization correction value (expected on|off): {}", optarg);
                exit(EXIT_FAILURE);
            }
            polarization_correction = value;
            break;
        }
        case OPT_SOLID_ANGLE_CORRECTION: {
            bool value;
            if (!parse_on_off(optarg, value)) {
                logger.Error("Invalid solid angle correction value (expected on|off): {}", optarg);
                exit(EXIT_FAILURE);
            }
            solid_angle_correction = value;
            break;
        }
        case OPT_BEAM_X: beam_x = parse_float_arg(optarg, "--beam-x", logger); break;
        case OPT_BEAM_Y: beam_y = parse_float_arg(optarg, "--beam-y", logger); break;
        case OPT_DETECTOR_DISTANCE: detector_distance_mm = parse_float_arg(optarg, "--detector-distance", logger); break;
        case OPT_WAVELENGTH: {
            // Guard > 0: wavelength is used as a divisor (WVL_1A_IN_KEV / wavelength) below, and a 0
            // would produce a non-finite incident energy that throws unguarded and aborts the process.
            float w = parse_float_arg(optarg, "--wavelength", logger);
            if (!(w > 0.0f)) {
                logger.Error("Invalid wavelength (must be > 0 A): {}", optarg);
                exit(EXIT_FAILURE);
            }
            wavelength_A = w;
            break;
        }
        case OPT_ROT1: rot1_rad = parse_float_arg(optarg, "--rot1", logger); break;
        case OPT_ROT2: rot2_rad = parse_float_arg(optarg, "--rot2", logger); break;
        case OPT_POLARIZATION: polarization_factor = parse_float_arg(optarg, "--polarization", logger); break;
        case OPT_SCALING_ITERATIONS:
            scaling_iter = atoi(optarg);
            if (scaling_iter <= 0) {
                logger.Error("Invalid scaling iteration count: {}", scaling_iter);
                exit(EXIT_FAILURE);
            }
            break;
        case OPT_BANDWIDTH:
            bandwidth_fwhm = atof(optarg);
            if (!(bandwidth_fwhm.value() >= 0.0f)) {
                logger.Error("Invalid bandwidth: {}", optarg);
                exit(EXIT_FAILURE);
            }
            break;

        default:
            print_usage();
            exit(EXIT_FAILURE);
    }
}

    if (optind != argc - 1) {
        logger.Error("Input file not specified");
        print_usage();
        exit(EXIT_FAILURE);
    }

    input_file = argv[optind];
    logger.Verbose(verbose);

    // -N defaults to 0 = "use all hardware threads"; resolve it to a concrete count here so every mode
    // behaves the same. The scale/merge engines expand 0 on their own, but the per-image processing
    // loop (Rugnux) spawns exactly nthreads workers, so passing 0 there would spawn none and process
    // nothing - hence resolving it centrally rather than relying on each consumer.
    if (nthreads <= 0) {
        unsigned int hw = std::thread::hardware_concurrency();
        nthreads = hw > 0 ? static_cast<int>(hw) : 1;
    }

    if (azint_only && scale_only) {
        logger.Error("--azint-only and --scale are mutually exclusive");
        exit(EXIT_FAILURE);
    }

    // Validate space group number early
    const gemmi::SpaceGroup *space_group = nullptr;
    if (space_group_number.has_value()) {
        space_group = gemmi::find_spacegroup_by_number(space_group_number.value());
        if (!space_group) {
            logger.Error("Unknown space group number {}", space_group_number.value());
            exit(EXIT_FAILURE);
        }
        logger.Info("Using space group {} (number {})", space_group->hm, space_group_number.value());
    }

    // 1. Read Input File
    JFJochHDF5Reader reader;
    try {
        reader.ReadFile(input_file);
    } catch (const std::exception &e) {
        logger.Error("Error reading input file: {}", e.what());
        exit(EXIT_FAILURE);
    }

    const auto dataset = reader.GetDataset();
    if (!dataset) {
        logger.Error("No experiment dataset found in the input file");
        exit(EXIT_FAILURE);
    }

    if (rotation_indexing_image_count <= 0) {
        logger.Error("Invalid number of rotation indexing images: {}", rotation_indexing_image_count);
        exit(EXIT_FAILURE);
    }

    logger.Info("Loaded dataset from {}", input_file);

    std::vector<MergedReflection> reference_data;
    if (!ref_mtz.empty()) {
        try {
            const auto reference = LoadReferenceMtz(
                ref_mtz, ref_column.empty() ? std::nullopt : std::optional<std::string>(ref_column));
            reference_data = reference.reflections;

            logger.Info("Loaded {} reference reflections from {} (column {}{}{})",
                        reference_data.size(), ref_mtz, reference.used_column,
                        reference.squared ? ", squared to intensity" : "",
                        reference.default_column ? ", auto-selected" : ", user-specified");
            if (reference.d_max > 0.0)
                logger.Info("Reference resolution range {:.2f} - {:.2f} A", reference.d_max, reference.d_min);
            if (reference.cell.has_value())
                logger.Info("Reference unit cell: a={:.3f} b={:.3f} c={:.3f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
                            reference.cell->a, reference.cell->b, reference.cell->c,
                            reference.cell->alpha, reference.cell->beta, reference.cell->gamma);
            if (!reference.space_group_name.empty())
                logger.Info("Reference space group: {} (number {})",
                            reference.space_group_name, reference.space_group_number.value_or(0));

            // Check the reference against the cell that will actually drive the merge. --scale merges
            // in the cell stored in the input file (as the former jfjoch_scale did); the -C override
            // only takes effect on the full-analysis path, which otherwise determines its cell later by
            // indexing (unknown here, so nothing can be checked yet).
            const std::optional<UnitCell> data_cell =
                scale_only ? dataset->experiment.GetUnitCell() : fixed_reference_unit_cell;
            const auto warning = ReferenceConsistencyWarning(
                reference, data_cell,
                space_group_number.has_value() ? std::optional<int>(static_cast<int>(*space_group_number))
                                               : std::nullopt);
            if (!warning.empty())
                logger.Warning("{}", warning);
        } catch (const std::exception &e) {
            logger.Error("Error reading reference MTZ {}: {}", ref_mtz, e.what());
            exit(EXIT_FAILURE);
        }
    }

    // --scale: re-scale and merge the already-integrated reflections stored in the input file,
    // without re-running spot finding or integration (folded in from the former rugnux_scale tool).
    if (scale_only) {
        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;
        auto reflections = reader.ReadReflections(start_image, last_image);

        DiffractionExperiment experiment(dataset->experiment);
        configure_offline_output(experiment, output_prefix);
        // Keep the space group stored in the input file (written by the full pipeline) unless -S overrides.
        if (space_group_number.has_value())
            experiment.SpaceGroupNumber(space_group_number);

        // A rotation (goniometer) dataset uses RotationScaleMerge unless --force-still asks for stills scaling.
        IndexingSettings indexing_settings;
        indexing_settings.RotationIndexing(experiment.GetGoniometer().has_value() && !force_still);
        experiment.ImportIndexingSettings(indexing_settings);

        ScalingSettings scaling_settings;
        if (d_min_scale_merge)
            scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
        if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
        if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target);
        if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count);
        scaling_settings.MergeFriedel(!anomalous_mode);
        scaling_settings.RefineB(refine_bfactor);
        scaling_settings.MinPartiality(min_partiality);
        scaling_settings.MinCCForImage(min_image_cc / 100.0);   // --min-image-cc is percent; the setting is a fraction
        if (intensity_format)
            scaling_settings.FileFormat(intensity_format.value());
        experiment.ImportScalingSettings(scaling_settings);

        if (!experiment.GetUnitCell()) {
            logger.Error("Experiment unit cell not found, cannot update reflection resolution");
            exit(EXIT_FAILURE);
        }
        auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections);
        logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images);
        experiment.ImagesPerTrigger(refl_stats.n_images);

        const auto scale_start = std::chrono::steady_clock::now();
        std::vector<MergedReflection> merged_reflections;
        MergeStatistics merged_statistics;
        double error_model_isa = 0.0;

        // Rotation (rot3d): the dedicated RotationScaleMerge does the whole self-scale -> 3D combine ->
        // merge. It does not support external-reference scaling, B-factor or wedge refinement. Everything
        // else (stills, reference scaling) uses ScaleOnTheFly + MergeOnTheFly.
        const bool is_rotation = experiment.IsRotationIndexing();
        if (is_rotation) {
            if (!reference_data.empty() || experiment.GetScalingSettings().GetRefineB()
                || experiment.GetRefineRotationWedgeInScaling()
                || experiment.GetScalingSettings().GetRotationWedgeForScaling().has_value())
                throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                      "Rotation scaling/merging (RotationScaleMerge) does not support reference "
                                      "scaling, B-factor refinement or wedge refinement");
            RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(),
                                   scaling_iter, 0.0f, nthreads, logger);
            rsm.Ingest();
            auto r = rsm.Run(false);
            merged_reflections = std::move(r.merged);
            merged_statistics = std::move(r.statistics);
            error_model_isa = r.isa;
        } else {
            for (int i = 0; i < scaling_iter; i++) {
                if (reference_data.empty())
                    ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads);
                else
                    ScaleOnTheFly(experiment, reference_data).Scale(reflections, nthreads);
            }
            MergeOnTheFly merge_engine(experiment);
            merge_engine.ReferenceCell(experiment.GetUnitCell());
            for (size_t i = 0; i < reflections.size(); ++i)
                merge_engine.AddImage(reflections[i], static_cast<int64_t>(i));
            merged_reflections = merge_engine.ExportReflections();

            // Automatic high-resolution cutoff (post-merge), matching the full-analysis path: a manual
            // --scaling-high-resolution wins, otherwise trim the written reflections + reported shells
            // to the CC1/2 fall-off. (Rotation is cut inside RotationScaleMerge above.)
            const auto &cut_ss = experiment.GetScalingSettings();
            std::optional<double> effective_d_min = cut_ss.GetHighResolutionLimit_A();
            if (!effective_d_min
                && cut_ss.GetResolutionCutoff() == ResolutionCutoffMethod::CCHalfLogistic) {
                const auto rc = ComputeCCHalfLogisticCutoff(merged_reflections, cut_ss.GetResolutionCCTarget(),
                                                            logger);
                if (rc.d_cut) {
                    effective_d_min = rc.d_cut;
                    logger.Info("Auto resolution cutoff: {:.2f} A ({}; override with --scaling-high-resolution)",
                                *rc.d_cut, rc.note);
                }
            }
            if (effective_d_min)
                merged_reflections.erase(std::remove_if(merged_reflections.begin(), merged_reflections.end(),
                    [&](const MergedReflection &m) { return std::isfinite(m.d) && m.d < *effective_d_min; }),
                    merged_reflections.end());

            merged_statistics = merge_engine.MergeStats(merged_reflections, reflections, reference_data,
                                                        effective_d_min);
            error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0;
        }

        logger.Info("Scale + merge completed in {:.2f} s ({} unique reflections)",
                    std::chrono::duration<double>(std::chrono::steady_clock::now() - scale_start).count(),
                    merged_reflections.size());

        std::cout << merged_statistics;

        // Space-group determination lives in the full rugnux pipeline; --scale only consumes a space
        // group (from the file or -S) and merges in it.
        const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
        if (!fixed_space_group)
            logger.Warning("No space group in the input file or on the command line - merged in P1. "
                           "Re-run rugnux (which determines and stores the space group) or pass "
                           "-S to scale and merge in the correct symmetry.");

        const auto twin_sg_number = experiment.GetSpaceGroupNumber();
        const gemmi::SpaceGroup *twin_sg = twin_sg_number
            ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr;
        const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg);
        std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl;

        if (!output_prefix.empty())
            WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics,
                             error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?",
                             twinning, output_prefix);
        return 0;
    }

    uint64_t total_images_in_file = reader.GetNumberOfImages();
    if (end_image < 0 || end_image > total_images_in_file)
        end_image = total_images_in_file;

    if (image_stride < 0) {
        logger.Error("Image stride cannot be negative");
        exit(EXIT_FAILURE);
    }

    if (image_stride == 0) {
        logger.Error("Image stride cannot be zero");
        exit(EXIT_FAILURE);
    }

    int images_to_process = (end_image - start_image) / image_stride;

    if (images_to_process <= 0) {
        logger.Warning("No images to process (Start: {}, End: {} Stride: {}, Total: {})", start_image, end_image,
            image_stride, total_images_in_file);
        return 0;
    }

    logger.Info("Starting analysis of {} images (range {}-{}) using {} threads",
                images_to_process, start_image, end_image, nthreads);

    // 2. Setup Experiment & Components
    DiffractionExperiment experiment(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.
    if (beam_x) experiment.BeamX_pxl(beam_x.value());
    if (beam_y) experiment.BeamY_pxl(beam_y.value());
    if (detector_distance_mm) experiment.DetectorDistance_mm(detector_distance_mm.value());
    if (wavelength_A) experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / wavelength_A.value());
    if (rot1_rad) experiment.PoniRot1_rad(rot1_rad.value());
    if (rot2_rad) experiment.PoniRot2_rad(rot2_rad.value());
    if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value());

    // Azimuthal integration (default q-spacing 0.01 1/A, from AzimuthalIntegrationSettings): the profile
    // resolves the narrow ice rings for the ice-ring score. Shared by --azint-only and full analysis.
    // -q / --azim-* / correction flags override; defaults come from the input file.
    {
        AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings();
        if (min_q || max_q)
            azint_settings.QRange_recipA(min_q.value_or(azint_settings.GetLowQ_recipA()),
                                         max_q.value_or(azint_settings.GetHighQ_recipA()));
        if (q_spacing)
            azint_settings.QSpacing_recipA(q_spacing.value());
        if (azimuthal_bins)
            azint_settings.AzimuthalBinCount(azimuthal_bins.value());
        if (polarization_correction)
            azint_settings.PolarizationCorrection(polarization_correction.value());
        if (solid_angle_correction)
            azint_settings.SolidAngleCorrection(solid_angle_correction.value());
        experiment.ImportAzimuthalIntegrationSettings(azint_settings);
        logger.Info("Azimuthal integration: Q [{:.4f}, {:.4f}] 1/A, spacing {:.4f}, {} Q x {} azimuthal bins",
                    azint_settings.GetLowQ_recipA(), azint_settings.GetHighQ_recipA(),
                    azint_settings.GetQSpacing_recipA(), azint_settings.GetQBinCount(),
                    azint_settings.GetAzimuthalBinCount());
    }

    // --azint-only: azimuthal integration only (no spot finding / indexing / scaling). Rugnux reads
    // the geometry and azimuthal-integration settings configured above off the experiment.
    if (azint_only) {
        ProcessConfig config;
        config.mode = ProcessMode::AzimuthalIntegration;
        config.start_image = start_image;
        config.end_image = end_image;
        config.stride = image_stride;
        config.nthreads = nthreads;
        config.output_prefix = output_prefix;

        Rugnux process(reader, experiment, *dataset->pixel_mask, config);
        g_active_process = &process;
        std::signal(SIGINT, handle_sigint);

        ProcessResult result;
        try {
            result = process.Run();
        } catch (const std::exception &e) {
            logger.Error("Processing failed: {}", e.what());
            exit(EXIT_FAILURE);
        }
        g_active_process = nullptr;

        std::cout << fmt::format("Processing time:  {:.2f} s", result.processing_time_s) << std::endl;
        std::cout << fmt::format("Frame rate:       {:.2f} Hz", result.frame_rate_hz) << std::endl;
        std::cout << fmt::format("Total throughput: {:.2f} MB/s", result.throughput_MBs) << std::endl;
        if (result.cancelled)
            logger.Warning("Processing was cancelled after {} images", result.images_processed);
        return 0;
    }

    configure_offline_output(experiment, output_prefix);
    experiment.SpaceGroupNumber(space_group_number);
    experiment.ImagesPerTrigger(images_to_process);

    // Re-determine the unit cell from scratch: discard any cell stored in the input file so
    // indexing is not biased by it. A stale or wrong stored cell otherwise resolves the indexing
    // algorithm to FFBIDX and drives it to the wrong lattice (e.g. a non-cubic cell for cubic
    // insulin). A user-supplied -C cell still takes effect (clears to nullopt when absent).
    experiment.SetUnitCell(fixed_reference_unit_cell);

    if (max_spot_count_override.has_value()) {
        experiment.MaxSpotCount(max_spot_count_override.value());
        logger.Info("Max spot count overridden to {}", max_spot_count_override.value());
    }

    // X-ray bandwidth: CLI overrides the value carried in the dataset; otherwise
    // keep whatever the dataset provided (0 / none -> monochromatic).
    if (bandwidth_fwhm)
        experiment.BandwidthFWHM(bandwidth_fwhm);
    if (experiment.GetBandwidthFWHM())
        logger.Info("X-ray bandwidth FWHM set to {:.4f}", experiment.GetBandwidthFWHM().value());

    // Rotation vs stills. A dataset collected on a rotation goniometer is processed as rotation data
    // (two-pass indexing) by default; --force-still forces per-frame stills. The rotation flags
    // (-R / --single-pass-rotation / --force-rotation-lattice) still request rotation explicitly and
    // choose the pass/lattice; at this point they show up as rotation_indexing already being set.
    const bool has_goniometer = experiment.GetGoniometer().has_value();
    if (force_still) {
        if (rotation_indexing) {
            logger.Error("--force-still conflicts with -R / --single-pass-rotation / --force-rotation-lattice");
            exit(EXIT_FAILURE);
        }
        if (has_goniometer)
            logger.Info("--force-still: treating the rotation dataset as independent stills");
    } else if (!rotation_indexing && has_goniometer) {
        rotation_indexing = true;
        two_pass_rotation = true;
        logger.Info("Dataset has a rotation goniometer axis: processing as rotation data (two-pass "
                    "indexing). Use --force-still to treat it as stills.");
    }

    // Scaling and merging are on by default (run_scaling initialised true); --no-merge turns them off
    // for both rotation and stills, in which case only the per-image _process.h5 is written.

    // Configure Indexing
    IndexingSettings indexing_settings;
    indexing_settings.Algorithm(indexing_algorithm);
    indexing_settings.RotationIndexing(rotation_indexing);
    if (rotation_indexing_range.has_value())
        indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value());
    indexing_settings.GeomRefinementAlgorithm(refinement_algorithm);
    experiment.ImportIndexingSettings(indexing_settings);

    // --detect-ice-rings[=on|off] overrides the value carried in from the dataset (HDF5MetadataSource
    // sets DetectIceRings from the master file's detect_ice_rings key); with no flag the dataset stands.
    if (detect_ice_rings.has_value())
        experiment.DetectIceRings(detect_ice_rings.value());

    // Scale-fulls refits the per-frame scale on the rotation combined fulls; on by default for rotation
    // data (where it lifts ISa substantially) and off for stills. --no-scale-fulls overrides.
    const bool scale_fulls = scale_fulls_arg.value_or(rotation_indexing);

    ScalingSettings scaling_settings;
    scaling_settings.ScaleFulls(scale_fulls);
    scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rotation_indexing ? SMOOTH_G_DEFAULT_DEG : 0.0));
    if (d_min_scale_merge)
        scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
    if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
    if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target);
    if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count);
    scaling_settings.MergeFriedel(!anomalous_mode);
    scaling_settings.RefineB(refine_bfactor);
    scaling_settings.MinPartiality(min_partiality);
    // Capture-aware systematic sigma defaults ON (1.0) for the rot3d combine - it down-weights the
    // over-extrapolated under-captured fulls and, with the mosaicity fix, lifts rotation ISa/anomalous
    // substantially. Off for non-rot3d (no combine). An explicit --capture-uncertainty always wins.
    scaling_settings.CaptureUncertaintyCoeff(capture_uncertainty_arg.value_or(rotation_indexing ? 1.0 : 0.0));
    scaling_settings.ForcedMosaicity(forced_mosaicity_arg);
    scaling_settings.MinCCForImage(min_image_cc / 100.0);   // --min-image-cc is in percent; the setting is a fraction
    scaling_settings.OutlierRejectNsigma(
        outlier_reject_nsigma.value_or(rotation_indexing ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));
    if (intensity_format)
        scaling_settings.FileFormat(intensity_format.value());

    experiment.ImportScalingSettings(scaling_settings);

    // Integration radii: r1 (signal box), r2/r3 (background annulus).
    if (integration_radius_arg) {
        std::vector<float> rr;
        std::stringstream ss(*integration_radius_arg);
        std::string tok;
        while (std::getline(ss, tok, ',')) {
            trim_in_place(tok);
            if (!tok.empty()) rr.push_back(std::stof(tok));
        }
        float r1, r2, r3;
        if (rr.size() == 1) { r1 = rr[0]; r2 = r1 + 2.0f; r3 = r1 + 4.0f; }
        else if (rr.size() == 3) { r1 = rr[0]; r2 = rr[1]; r3 = rr[2]; }
        else { logger.Error("--integration-radius expects r1 or r1,r2,r3"); return 1; }
        BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
        bis.R1(r1).R2(r2).R3(r3);
        experiment.ImportBraggIntegrationSettings(bis);
        logger.Info("Integration radii set to r1={:.1f} r2={:.1f} r3={:.1f}", r1, r2, r3);
    }

    if (integrator_mode) {
        BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
        bis.Integrator(*integrator_mode);
        experiment.ImportBraggIntegrationSettings(bis);
        logger.Info("Integrator set to {}", *integrator_mode == IntegratorMode::BoxSum ? "box-sum"
                    : *integrator_mode == IntegratorMode::ProfileGaussian ? "profile (gaussian)"
                    : "profile (empirical)");
    }

    SpotFindingSettings spot_settings;
    spot_settings.enable = true;
    spot_settings.indexing = true;
    spot_settings.high_resolution_limit = d_min_spot_finding;
    spot_settings.signal_to_noise_threshold = sigma_spot_finding;
    spot_settings.photon_count_threshold = photon_count_threshold_spot_finding;
    if (d_min_spot_finding > 0.0f)
        spot_settings.high_resolution_limit = d_min_spot_finding;

    // 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.start_image = start_image;
    config.end_image = end_image;
    config.stride = image_stride;
    config.nthreads = nthreads;
    config.output_prefix = output_prefix;
    config.spot_finding = spot_settings;
    config.rotation_indexing = rotation_indexing;
    config.two_pass_rotation = two_pass_rotation;
    config.reuse_rotation_spots = reuse_rotation_spots;
    config.rotation_indexing_image_count = rotation_indexing_image_count;
    config.forced_rotation_lattice = forced_rotation_lattice;
    config.run_scaling = run_scaling;
    config.scaling_iter = scaling_iter;
    config.reference_data = reference_data;
    config.observation_dump_path = dump_observations;
    // 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.
    config.write_process_h5 = run_scaling ? write_process_h5_flag : true;

    Rugnux process(reader, experiment, *dataset->pixel_mask, config);

    g_active_process = &process;
    std::signal(SIGINT, handle_sigint);

    ProcessResult result;
    try {
        result = process.Run();
    } catch (const std::exception &e) {
        logger.Error("Processing failed: {}", e.what());
        exit(EXIT_FAILURE);
    }
    g_active_process = nullptr;

    if (!result.merge_statistics_text.empty())
        std::cout << std::endl << result.merge_statistics_text << std::endl;

    // Report statistics
    std::cout << fmt::format("Processing time: {:.2f} s", result.processing_time_s) << std::endl;
    std::cout << fmt::format("Frame rate:      {:.2f} Hz", result.frame_rate_hz) << std::endl;
    std::cout << fmt::format("Total throughput:{:.2f} MB/s", result.throughput_MBs) << std::endl;
    if (result.indexing_rate.has_value())
        std::cout << fmt::format("Indexing rate:   {:.2f}%", result.indexing_rate.value() * 100.0) << std::endl;
    if (result.consensus_cell.has_value()) {
        const auto &c = result.consensus_cell.value();
        std::cout << fmt::format("Unit cell:       a={:.2f} b={:.2f} c={:.2f} alpha={:.2f} beta={:.2f} gamma={:.2f}",
                                 c.a, c.b, c.c, c.alpha, c.beta, c.gamma) << std::endl;
    }

    const auto &t = result.mean_processing_time;
    std::cout << fmt::format(
        "Per-image time (mean; ms): decompress {:.2f} preprocess {:.2f} azint {:.2f} spot finding {:.2f} "
        "indexing {:.2f} refinement {:.2f} indexing analysis {:.2f} prediction {:.2f} integration {:.2f} "
        "scaling {:.2f} total {:.2f}",
        t.compression * 1e3, t.preprocessing * 1e3, t.azint * 1e3, t.spot_finding * 1e3,
        t.indexing * 1e3, t.refinement * 1e3, t.indexing_analysis * 1e3, t.bragg_prediction * 1e3,
        t.integration * 1e3, t.image_scale * 1e3, t.processing * 1e3) << std::endl;

    if (result.cancelled)
        logger.Warning("Processing was cancelled after {} images", result.images_processed);
}