leonarski_fandClaude Opus 5
6f4917dcee
Build Packages / build:viewer-tgz:cpu (push) Successful in 8m0s
Build Packages / build:viewer-tgz:cuda (push) Successful in 9m13s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 13m40s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 13m48s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 13m42s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 14m27s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 14m43s
Build Packages / build:rpm (rocky8) (push) Successful in 11m32s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 13m9s
Build Packages / XDS test (durin plugin) (push) Successful in 7m22s
Build Packages / Generate python client (push) Successful in 27s
Build Packages / Build documentation (push) Successful in 1m7s
Build Packages / Create release (push) Skipped
Build Packages / build:rpm (ubuntu2404) (push) Successful in 12m57s
Build Packages / build:rpm (rocky9) (push) Successful in 13m25s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 13m50s
Build Packages / DIALS test (push) Successful in 13m52s
Build Packages / XDS test (neggia plugin) (push) Successful in 8m0s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 8m55s
Build Packages / Unit tests (push) Successful in 1h2m59s
Build Packages / build:windows:nocuda (push) Canceled after 0s
Build Packages / build:windows:cuda (push) Canceled after 0s
It was held back because a 33-crystal rotation battery showed it breaking three
crystals deterministically - a lost space group, a halved indexing rate and a
collapsed merge. None of those causes turned out to be in detection.
The extra spots adaptive finds are real. Measured per spot against a
finder-neutral local background: 64% recur at the same position on the adjacent
frame (chance rate 0.5%) with 2-frame rocking curves, and 0.00% would fail a
conventional local SNR >= 4 test, median local SNR 34. What they include is
genuine peaks belonging to no lattice the indexer found, and the damage they did
scaled with their absolute COUNT (80.6 per frame against 36.8 for the fixed
finder), not with their quality - which is why nothing aimed at judging
individual spots ever worked.
The three failures fell to fixes elsewhere:
merge collapsed - a per-frame scale free to collapse toward zero amplified
two junk frames by 546x (704098712, ec7a82613). Not a
detection problem at all: the fixed-threshold finder trips
the same bug on a different frame range.
space group lost - the per-image geometry refinement was dragged 2.3-2.8 deg
off by the weak-spot tail in an unweighted fit; weighting
each spot by how strong it is FOR ITS RESOLUTION fixed it
(ae126c3d5), and gained a point group for the fixed finder
too.
indexing halved - gone with the same two; that crystal is now better under
adaptive (CC1/2 92.6 -> 95.9, high-shell 44.9 -> 56.4).
Battery, 33 crystals, adaptive vs the fixed finder:
exact space group matching XDS 26/33 vs 25/33
point group matching XDS 29/33 vs 30/33
ISa better on 6 crystals
recovers a screw axis the other misses (P321 -> P3121)
The one point group it loses is a tetragonal crystal where adaptive collects
2.4x the observations at better R_meas (21.4% vs 28.3%) and better ISa (4.58 vs
2.91), and forced to the right group gives CC1/2 99.2% at multiplicity 10.7 -
matching XDS. Only the automatic symmetry call fails there, and six candidate
causes have been measured and refuted (mixed indexing hands, off-crystal frames,
a badly integrated minority, radiation damage, pseudo-tetragonality, uncorrected
anisotropy). It is left as the subgroup, which is the recoverable direction: -S
gives XDS-quality data from the same run, whereas the failures this unblocks
were not recoverable.
--no-adaptive-spots reverts to the fixed-threshold finder.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
1597 lines
84 KiB
C++
1597 lines
84 KiB
C++
// 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 <optional>
|
|
#include <getopt.h>
|
|
#include <iostream>
|
|
#include <limits>
|
|
#include <sstream>
|
|
#include <string>
|
|
#include <thread>
|
|
#include <type_traits>
|
|
#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/RfreeFlags.h"
|
|
#include "../image_analysis/scale_merge/ScaleOnTheFly.h"
|
|
#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
|
|
#include "../image_analysis/scale_merge/RotationScaleMerge.h"
|
|
#include "../image_analysis/scale_merge/ResolutionCutoff.h"
|
|
#include "../image_analysis/scale_merge/TwinningAnalysis.h"
|
|
#include "../image_analysis/scale_merge/SearchSpaceGroup.h"
|
|
#include "Rugnux.h"
|
|
#include "ModelValidation.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-batch relative-B batch width (bare --relative-b), in degrees of rotation.
|
|
constexpr double RELATIVE_B_DEFAULT_DEG = 10.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 << " --min-pix-per-spot <num> Minimum connected strong pixels per spot. If omitted, min-pix is chosen PER IMAGE (stills indexing): the frame is indexed at min-pix 3/2/1 and the one maximising indexed count x indexed fraction is kept. Give an explicit value to force a fixed min-pix instead." << std::endl;
|
|
std::cout << " --adaptive-spots Self-calibrating detection (DEFAULT): the strong-pixel threshold comes from each image's own per-resolution-ring noise instead of the fixed --spot-threshold, so one setting adapts across datasets (no per-dataset --spot-threshold/--spot-sigma tuning)." << std::endl;
|
|
std::cout << " --no-adaptive-spots Turn adaptive detection off and use the fixed --spot-threshold / --spot-sigma finder instead" << std::endl;
|
|
std::cout << " --spot-false-pixels <num> Adaptive detection operating point: expected noise pixels tolerated per frame (default: 100; implies --adaptive-spots)" << std::endl;
|
|
std::cout << " --spot-high-resolution <num> High resolution limit for spot finding. If omitted, stills extend as far as the detector reaches (no resolution clipping) and rotation data keeps a 1.5 A limit." << std::endl;
|
|
std::cout << " --spot-low-resolution <num> Low resolution limit for spot finding, in A (default: 50; lower it, e.g. 24, to exclude the direct-beam halo on weakly-diffracting serial data)" << 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 << " --rotation-no-postrefine Disable the (default-on) two-pass rotation post-refine (post-refine detector distance/beam + cell/axis, then re-integrate; the refined pass is the canonical <prefix>_* output, the header-geometry pass is kept as <prefix>_01_*)" << std::endl;
|
|
std::cout << " -X, --indexing-algorithm <txt> Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl;
|
|
std::cout << " -S, --space-group <num|symbol> Space group number (92) or symbol (P43212) - for 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|flex); flex tries all three per image and keeps whichever indexes the most spots (alias: multi)" << std::endl;
|
|
std::cout << " --refine-geometry[=N|off] Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default: 200) then re-indexes (lifts weak-stills indexing). Default ON for stills when a reference cell is given (-C / reference MTZ); =off disables" << 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 << " --relative-b[=deg] rot3d: fit a per-batch relative-B (beyond the single decay slope) over deg-degree batches; cross-validated (default: 10 deg when bare; off otherwise)" << std::endl;
|
|
std::cout << " --no-scaling-corrections rot3d: disable the (default-on) decay + absorption correction surfaces fitted on the fulls after scale-fulls" << std::endl;
|
|
std::cout << " --no-expected-variance-merge stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting" << std::endl;
|
|
std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << 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 << " --min-captured-fraction <num> rot3d: drop a combined full whose rocking curve was captured below this fraction (edge-of-sweep truncated fulls) (default: 0.7 for rotation, 0 otherwise; 0 = off)" << 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 << " --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 << " -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, else FP/FOBS/F)" << std::endl;
|
|
std::cout << " --model <file.pdb> After merging, validate vs this model: R-free + 2Fo-Fc/Fo-Fc maps" << 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 << " --background-trim <f> Monochromatic (rotation + still): symmetric trimmed-mean fraction for the background ring (0<=f<0.5, default 0.10; 0 = plain mean). Removes the high-side bias that over-subtracts weak high-angle spots (broadband data keep the sigma-clip instead)" << std::endl;
|
|
std::cout << " --integrator <txt> Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
|
|
std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << 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). If omitted, integration extends to the highest Q the detector reaches." << 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_MIN_PIX_PER_SPOT,
|
|
OPT_ADAPTIVE_SPOTS,
|
|
OPT_NO_ADAPTIVE_SPOTS,
|
|
OPT_SPOT_FALSE_PIXELS,
|
|
OPT_SPOT_RESOLUTION,
|
|
OPT_SPOT_LOW_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_SINGLE_PASS_ROTATION,
|
|
OPT_REDO_ROTATION_SPOTS,
|
|
OPT_FORCE_ROTATION_LATTICE,
|
|
OPT_ROTATION_NO_POSTREFINE,
|
|
OPT_REFINE_GEOMETRY,
|
|
OPT_BANDWIDTH,
|
|
OPT_INTEGRATION_RADIUS,
|
|
OPT_BACKGROUND_TRIM,
|
|
OPT_REJECT_OUTLIERS,
|
|
OPT_REFERENCE_COLUMN,
|
|
OPT_MODEL,
|
|
OPT_DUMP_OBSERVATIONS,
|
|
OPT_INTEGRATOR,
|
|
OPT_SIMPLE_STILLS,
|
|
OPT_SCALE_FULLS,
|
|
OPT_CAPTURE_UNCERTAINTY,
|
|
OPT_MIN_CAPTURED_FRACTION,
|
|
OPT_MOSAICITY,
|
|
OPT_SMOOTH_G,
|
|
OPT_RELATIVE_B,
|
|
OPT_NO_SCALING_CORRECTIONS,
|
|
OPT_NO_EXPECTED_VARIANCE_MERGE,
|
|
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},
|
|
{"model", required_argument, nullptr, OPT_MODEL},
|
|
{"dump-observations", required_argument, nullptr, OPT_DUMP_OBSERVATIONS},
|
|
{"space-group", required_argument, nullptr, 'S'},
|
|
{"anomalous", no_argument, nullptr, 'A'},
|
|
{"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},
|
|
{"relative-b", optional_argument, nullptr, OPT_RELATIVE_B},
|
|
{"no-scaling-corrections", no_argument, nullptr, OPT_NO_SCALING_CORRECTIONS},
|
|
{"no-expected-variance-merge", no_argument, nullptr, OPT_NO_EXPECTED_VARIANCE_MERGE},
|
|
{"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},
|
|
{"rotation-no-postrefine", no_argument, nullptr, OPT_ROTATION_NO_POSTREFINE},
|
|
{"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY},
|
|
|
|
|
|
{"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA},
|
|
{"spot-threshold", required_argument, nullptr, OPT_SPOT_THRESHOLD},
|
|
{"min-pix-per-spot", required_argument, nullptr, OPT_MIN_PIX_PER_SPOT},
|
|
{"adaptive-spots", no_argument, nullptr, OPT_ADAPTIVE_SPOTS},
|
|
{"no-adaptive-spots", no_argument, nullptr, OPT_NO_ADAPTIVE_SPOTS},
|
|
{"spot-false-pixels", required_argument, nullptr, OPT_SPOT_FALSE_PIXELS},
|
|
{"spot-high-resolution", required_argument, nullptr, OPT_SPOT_RESOLUTION},
|
|
{"spot-low-resolution", required_argument, nullptr, OPT_SPOT_LOW_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},
|
|
{"min-captured-fraction", required_argument, nullptr, OPT_MIN_CAPTURED_FRACTION},
|
|
{"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},
|
|
{"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
|
|
{"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
|
|
{"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM},
|
|
{"integrator", required_argument, nullptr, OPT_INTEGRATOR},
|
|
{"simple-stills", no_argument, nullptr, OPT_SIMPLE_STILLS},
|
|
{"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS},
|
|
{"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
|
|
{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, optionally bounded to [min_value, max_value], or print a
|
|
// clear error and exit. getopt hands option arguments over as raw C strings; atoi()/atof() silently
|
|
// return 0 on non-numeric input (so a typo like "--min-pix-per-spot 2O" becomes 2 or 0) and std::sto*
|
|
// throws, which would terminate the program. This rejects non-numeric input, trailing garbage
|
|
// ("1.5foo"), integer overflow, and out-of-range values. T may be integral or floating-point; the
|
|
// bounds default to the full representable range (i.e. unbounded).
|
|
template <typename T>
|
|
T parse_number_arg(const char *arg, const char *option_name, Logger &logger,
|
|
T min_value = std::numeric_limits<T>::lowest(),
|
|
T max_value = std::numeric_limits<T>::max()) {
|
|
std::string s = arg ? arg : "";
|
|
trim_in_place(s);
|
|
T value{};
|
|
bool parsed = false;
|
|
if (!s.empty()) {
|
|
try {
|
|
size_t idx = 0;
|
|
if constexpr (std::is_integral_v<T>) {
|
|
const long long v = std::stoll(s, &idx);
|
|
value = static_cast<T>(v);
|
|
parsed = (idx == s.size()) && (static_cast<long long>(value) == v); // no overflow
|
|
} else {
|
|
value = static_cast<T>(std::stod(s, &idx));
|
|
parsed = (idx == s.size());
|
|
}
|
|
} catch (...) {}
|
|
}
|
|
if (!parsed) {
|
|
logger.Error("Invalid numeric value for {}: '{}'", option_name, arg ? arg : "<null>");
|
|
print_usage();
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
if (value < min_value || value > max_value) {
|
|
logger.Error("Value for {} out of range: {} (expected {} to {})",
|
|
option_name, value, min_value, max_value);
|
|
print_usage();
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
return value;
|
|
}
|
|
|
|
// Thin unbounded wrappers for the existing floating-point call sites.
|
|
double parse_double_arg(const char *arg, const char *option_name, Logger &logger) {
|
|
return parse_number_arg<double>(arg, option_name, logger);
|
|
}
|
|
|
|
float parse_float_arg(const char *arg, const char *option_name, Logger &logger) {
|
|
return parse_number_arg<float>(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);
|
|
// Offline CLI: the operator chose the output path, so allow an absolute -o (the multi-user guard
|
|
// that FilePrefix() applies is only for remotely-supplied prefixes in the broker/writer).
|
|
experiment.FilePrefixTrusted(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;
|
|
bool rotation_postrefine_geometry = true; // default on; --rotation-no-postrefine disables it
|
|
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
|
|
std::optional<double> relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given
|
|
bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption surfaces
|
|
bool no_expected_variance_merge = false; // --no-expected-variance-merge: restore observed-sigma stills merge weighting
|
|
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;
|
|
std::optional<int64_t> min_pix_per_spot; // unset -> adaptive per image; a value -> fixed min-pix
|
|
std::optional<bool> adaptive_spots; // unset -> per-workflow default (stills on, rotation off)
|
|
float false_pixels_per_frame = 100.0f;
|
|
std::string ref_mtz;
|
|
std::string ref_column;
|
|
std::string model_pdb; // --model: PDB to validate merged intensities against (R-free + maps)
|
|
std::string dump_observations; // diagnostic: dump unmerged -P rot3d fulls to this path
|
|
double min_partiality = 0.02;
|
|
std::optional<double> min_captured_fraction_arg; // explicit --min-captured-fraction; default depends on rotation
|
|
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<int> refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass
|
|
bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on
|
|
|
|
std::optional<float> bandwidth_fwhm; // relative FWHM of dlambda/lambda
|
|
|
|
IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto;
|
|
GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter;
|
|
|
|
std::optional<float> d_min_spot_finding; // unset -> as far as the detector reaches
|
|
float d_max_spot_finding = 0; // 0 = keep the SpotFindingSettings default (50 A)
|
|
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<double> background_trim_arg; // --background-trim: background-ring trimmed-mean fraction
|
|
std::optional<IntegratorMode> integrator_mode; // --integrator boxsum|gaussian|empirical
|
|
bool simple_stills_flag = false; // --simple-stills: disable the default stills partiality post-refinement
|
|
std::optional<double> outlier_reject_nsigma; // merge per-observation outlier rejection
|
|
|
|
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:FAS: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_ROTATION_NO_POSTREFINE:
|
|
rotation_postrefine_geometry = false;
|
|
break;
|
|
case OPT_REFINE_GEOMETRY: {
|
|
if (optarg && std::string(optarg) == "off") {
|
|
refine_geometry = std::nullopt;
|
|
refine_geometry_disabled = true; // opt out of the stills default-on
|
|
break;
|
|
}
|
|
// Positive frame count fed to the bundle adjust; upper-bounded so refine_frames * 50 in
|
|
// RefineStillsGeometry cannot overflow int (1e6 is already far past any real dataset).
|
|
refine_geometry = optarg
|
|
? parse_number_arg<int>(optarg, "--refine-geometry", logger, 1, 1000000)
|
|
: 200;
|
|
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 if (alg == "flex" || alg == "multi") // "multi" kept as a back-compat alias
|
|
refinement_algorithm = GeomRefinementAlgorithmEnum::Flex;
|
|
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_MODEL:
|
|
model_pdb = optarg;
|
|
break;
|
|
case OPT_DUMP_OBSERVATIONS:
|
|
dump_observations = optarg;
|
|
break;
|
|
case 'F':
|
|
indexing_algorithm = IndexingAlgorithmEnum::FFT;
|
|
break;
|
|
case 'A':
|
|
anomalous_mode = true;
|
|
break;
|
|
case 'S': {
|
|
// Accept a space-group number ("92") or a Hermann-Mauguin symbol ("P43212", "P 43 21 2").
|
|
char *end = nullptr;
|
|
const long as_number = strtol(optarg, &end, 10);
|
|
if (end != optarg && *end == '\0') {
|
|
space_group_number = as_number;
|
|
} else if (const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(optarg)) {
|
|
space_group_number = sg->number;
|
|
} else {
|
|
logger.Error("Unknown space group '{}' (use a number like 92 or a symbol like P43212)", optarg);
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
break;
|
|
}
|
|
case OPT_SPOT_SIGMA:
|
|
sigma_spot_finding = parse_number_arg<float>(optarg, "--spot-sigma", logger, 1.0f);
|
|
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 = parse_number_arg<int64_t>(optarg, "--spot-threshold", logger, 0);
|
|
logger.Info("Photon-count threshold level for spot finding set to {:d}",
|
|
photon_count_threshold_spot_finding);
|
|
break;
|
|
case OPT_MIN_PIX_PER_SPOT:
|
|
// Giving an explicit min-pix opts out of the per-image adaptive selection.
|
|
min_pix_per_spot = parse_number_arg<int64_t>(optarg, "--min-pix-per-spot", logger, 1);
|
|
logger.Info("Minimum pixels per spot fixed at {:d} (adaptive per-image min-pix off)", *min_pix_per_spot);
|
|
break;
|
|
case OPT_ADAPTIVE_SPOTS:
|
|
adaptive_spots = true;
|
|
logger.Info("Adaptive (self-calibrating) spot detection enabled");
|
|
break;
|
|
case OPT_NO_ADAPTIVE_SPOTS:
|
|
adaptive_spots = false;
|
|
logger.Info("Adaptive spot detection off: using the fixed --spot-threshold / --spot-sigma finder");
|
|
break;
|
|
case OPT_SPOT_FALSE_PIXELS:
|
|
false_pixels_per_frame = parse_number_arg<float>(optarg, "--spot-false-pixels", logger, 1.0f);
|
|
adaptive_spots = true;
|
|
logger.Info("Adaptive spot detection: expected false pixels/frame set to {:.0f}", false_pixels_per_frame);
|
|
break;
|
|
case OPT_SPOT_LOW_RESOLUTION:
|
|
d_max_spot_finding = parse_number_arg<float>(optarg, "--spot-low-resolution", logger, 0.0f);
|
|
logger.Info("Low resolution limit for spot finding set to {:.1f} A", d_max_spot_finding);
|
|
break;
|
|
case OPT_SPOT_RESOLUTION: {
|
|
// 0 has always meant "no limit" for this setting; keep that, but express it as the unset
|
|
// optional the rest of the code understands. Passing the 0 through instead reached
|
|
// ResolutionShells (via the spot plot), which rejects a zero d_min and threw away every image.
|
|
const auto d_min = parse_number_arg<float>(optarg, "--spot-high-resolution", logger, 0.0f);
|
|
if (d_min > 0.0f) {
|
|
d_min_spot_finding = d_min;
|
|
logger.Info("High resolution limit for spot finding set to {:.2f} A", d_min);
|
|
} else {
|
|
d_min_spot_finding.reset();
|
|
logger.Info("No high resolution limit for spot finding: as far as the detector reaches");
|
|
}
|
|
break;
|
|
}
|
|
case OPT_MAX_SPOTS:
|
|
max_spot_count_override = parse_number_arg<int64_t>(optarg, "--max-spots", logger, 1);
|
|
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_RELATIVE_B:
|
|
relative_b_deg_arg = optarg ? parse_double_arg(optarg, "--relative-b", logger) : RELATIVE_B_DEFAULT_DEG;
|
|
break;
|
|
case OPT_NO_EXPECTED_VARIANCE_MERGE:
|
|
no_expected_variance_merge = true;
|
|
break;
|
|
case OPT_NO_SCALING_CORRECTIONS:
|
|
no_scaling_corrections = true;
|
|
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_MIN_CAPTURED_FRACTION:
|
|
min_captured_fraction_arg = parse_double_arg(optarg, "--min-captured-fraction", 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_BACKGROUND_TRIM:
|
|
background_trim_arg = parse_double_arg(optarg, "--background-trim", logger);
|
|
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_SIMPLE_STILLS:
|
|
simple_stills_flag = true;
|
|
break;
|
|
case OPT_REJECT_OUTLIERS:
|
|
outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", 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_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;
|
|
bool reference_has_free_flags = false;
|
|
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;
|
|
reference_has_free_flags = reference.has_free_flags;
|
|
|
|
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.has_free_flags)
|
|
logger.Info("Reference carries R-free flags (column {}): {} of {} free; the merged "
|
|
"reflections will inherit this test set",
|
|
reference.free_column, reference.n_free, reference_data.size());
|
|
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);
|
|
|
|
// A reference MTZ fixes the space group and unit cell, unless -S / -C override them.
|
|
// (-S with the wrong enantiomorph, or -C with a different cell, is allowed - the explicit
|
|
// flag always wins.) The cell is a soft reference: indexing may drift within tolerance.
|
|
if (!space_group_number.has_value() && reference.space_group_number.has_value()) {
|
|
space_group_number = static_cast<int64_t>(*reference.space_group_number);
|
|
logger.Info("Fixing space group from reference MTZ: {} ({})",
|
|
reference.space_group_name, *space_group_number);
|
|
}
|
|
if (!fixed_reference_unit_cell.has_value() && reference.cell.has_value()) {
|
|
fixed_reference_unit_cell = reference.cell;
|
|
logger.Info("Fixing reference unit cell from reference MTZ (indexing may drift within tolerance)");
|
|
}
|
|
} 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.MinPartiality(min_partiality);
|
|
scaling_settings.MinCapturedFraction(min_captured_fraction_arg.value_or(
|
|
(experiment.GetGoniometer().has_value() && !force_still) ? 0.7 : 0.0));
|
|
scaling_settings.MinCCForImage(min_image_cc / 100.0); // --min-image-cc is percent; the setting is a fraction
|
|
scaling_settings.StillsPartialityRefine(!simple_stills_flag);
|
|
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
|
|
scaling_settings.OutlierRejectNsigma(
|
|
outlier_reject_nsigma.value_or(
|
|
(experiment.GetGoniometer().has_value() && !force_still) ? REJECT_OUTLIERS_DEFAULT_NSIGMA : 0.0));
|
|
// Mirror the de-novo rotation scaling defaults so --scale reproduces the full-pipeline merge on the
|
|
// same reflections (scale-fulls + smooth-G, and the opt-in per-batch relative-B). Without these the
|
|
// offline re-scale silently used a weaker model than the pipeline that wrote the _process.h5.
|
|
{
|
|
const bool rot = experiment.GetGoniometer().has_value() && !force_still;
|
|
scaling_settings.ScaleFulls(scale_fulls_arg.value_or(rot));
|
|
scaling_settings.SmoothGDegrees(smooth_g_deg_arg.value_or(rot ? SMOOTH_G_DEFAULT_DEG : 0.0));
|
|
scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0));
|
|
if (no_scaling_corrections)
|
|
scaling_settings.CorrectionSurfaces(false);
|
|
}
|
|
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, including the default-on decay + absorption correction surfaces. It does not support
|
|
// external-reference scaling 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.GetRefineRotationWedgeInScaling()
|
|
|| experiment.GetScalingSettings().GetRotationWedgeForScaling().has_value())
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Rotation scaling/merging (RotationScaleMerge) does not support reference "
|
|
"scaling 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 {
|
|
// Scaling self-references: the reference MTZ (if any) fixes the cell/space group, reports
|
|
// CCref and provides the R-free test set, but is NOT a scale anchor - scaling each image
|
|
// against a foreign dataset injects cross-dataset systematics and is a worse reference than
|
|
// the data's own merge. The per-image scale G is the exact one-pass solution, so one pass
|
|
// (iterating a self-rebuilt reference only re-fits the freshly-scaled noise on weak stills).
|
|
ScaleOnTheFly(experiment, MergeAll(experiment, reflections)).Scale(reflections, nthreads);
|
|
// Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
|
|
// orientation tilt against the merge and recompute each reflection's partiality + scale correction
|
|
// (no re-integration), then merge with the improved corrections.
|
|
if (experiment.GetScalingSettings().GetStillsPartialityRefine()) {
|
|
StillsPartialityRefine refiner(experiment);
|
|
const double mean_tilt = refiner.Run(reflections, nthreads);
|
|
logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
|
|
}
|
|
MergeOnTheFly merge_engine(experiment);
|
|
merge_engine.ReferenceCell(experiment.GetUnitCell());
|
|
// --min-image-cc has to hold for the merge itself, not only for the reported statistics.
|
|
merge_engine.FilterByImageCC(experiment.GetScalingSettings().GetMinCCForImage() > 0.0);
|
|
// Fit the (a, b) error model from symmetry-mate scatter before merging, exactly as the full
|
|
// pipeline does (Rugnux.cpp). Without this the offline --scale merge would use the identity
|
|
// model and produce much worse stills intensities (no (b*I)^2 systematic term, no sigma floor).
|
|
merge_engine.RefineErrorModel(reflections);
|
|
if (merge_engine.ErrorModelActive())
|
|
logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", merge_engine.ErrorModelA(),
|
|
merge_engine.ErrorModelB(),
|
|
merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0,
|
|
merge_engine.ErrorModelChi2());
|
|
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();
|
|
// The offline --scale path re-scales a stored _process.h5 and is never a P1 search merge.
|
|
const std::optional<double> effective_d_min = ApplyResolutionCutoff(
|
|
merged_reflections, cut_ss.GetHighResolutionLimit_A(), cut_ss.GetResolutionCutoff(),
|
|
cut_ss.GetResolutionCCTarget(), /*for_search=*/false, logger);
|
|
|
|
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());
|
|
|
|
// Inherit the campaign's shared R-free test set from the reference MTZ (overriding the
|
|
// per-hkl hash the merge assigned), so every dataset flags the same free reflections.
|
|
if (reference_has_free_flags && !reference_data.empty() && !merged_reflections.empty()) {
|
|
const auto sg = experiment.GetSpaceGroupNumber().value_or(1);
|
|
const size_t matched = ApplyReferenceFreeFlags(merged_reflections, static_cast<int32_t>(sg),
|
|
reference_data);
|
|
logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)",
|
|
matched, 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);
|
|
|
|
if (!output_prefix.empty() && !model_pdb.empty()) {
|
|
const auto data_sg = experiment.GetSpaceGroupNumber();
|
|
// With a reference MTZ the merohedral indexing was already resolved (stills per-image
|
|
// scaling); only probe indexing by R-free when model-only, with no reference.
|
|
ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb,
|
|
output_prefix, logger,
|
|
data_sg ? std::optional<int>(static_cast<int>(*data_sg)) : std::nullopt,
|
|
/*probe_indexing_ambiguity=*/reference_data.empty());
|
|
}
|
|
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;
|
|
}
|
|
|
|
// 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 ? max_q : azint_settings.GetRequestedHighQ_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 a cubic
|
|
// crystal). A user-supplied -C cell still takes effect (clears to nullopt when absent).
|
|
experiment.SetUnitCell(fixed_reference_unit_cell);
|
|
|
|
// --refine-geometry defaults ON for stills whenever a reference cell is available (-C or a
|
|
// reference MTZ): that is exactly when the geometry bundle-adjust can act (it anchors on a known
|
|
// cell) and where it helps weak/sparse stills. It is a no-op for rotation (which has its own
|
|
// two-pass) and for de-novo stills (no cell yet), so auto-enabling it only where it does something
|
|
// avoids spurious "skipping" warnings. Explicit --refine-geometry[=N] still forces it on;
|
|
// --refine-geometry=off opts out.
|
|
if (!refine_geometry.has_value() && !refine_geometry_disabled) {
|
|
const bool is_stills = !(experiment.GetGoniometer().has_value() && !force_still);
|
|
if (is_stills && experiment.GetUnitCell().has_value())
|
|
refine_geometry = 200;
|
|
}
|
|
|
|
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));
|
|
scaling_settings.RelativeBDegrees(relative_b_deg_arg.value_or(0.0)); // opt-in only; default off
|
|
if (no_scaling_corrections)
|
|
scaling_settings.CorrectionSurfaces(false);
|
|
scaling_settings.StillsPartialityRefine(!simple_stills_flag);
|
|
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
|
|
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.MinPartiality(min_partiality);
|
|
// Drop edge-of-sweep truncated fulls (rocking curve captured < this fraction) from the rot3d combine.
|
|
// Defaults ON (0.7) for rotation - removes the low-capture fulls that inflate low-res R-meas and
|
|
// slightly bias accuracy; off for non-rot3d (no combine). 0.7 (rather than 0.5) also strips the
|
|
// partiality-extrapolated fulls that dominate the intensity second moment on weakly-diffracting
|
|
// crystals, so the de-novo space-group search is no longer starved by the error-model I/sigma floor
|
|
// (e.g. a weakly-diffracting F-cubic or hexagonal crystal recovers its true space group instead of
|
|
// P1). An explicit --min-captured-fraction wins.
|
|
scaling_settings.MinCapturedFraction(min_captured_fraction_arg.value_or(rotation_indexing ? 0.7 : 0.0));
|
|
// 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));
|
|
|
|
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);
|
|
} else if (!rotation_indexing) {
|
|
// Stills spots span a range of crystal orientations captured in a single shot, so they land
|
|
// wider on the detector than the r1=4 monochromatic-rotation default assumes. A larger signal
|
|
// box lets the profile-fit integrator capture the whole spot while its profile weighting keeps
|
|
// the extra background from adding noise (a plain box-sum degrades with it). Measured R-free
|
|
// gains on serial stills. An explicit --integration-radius always wins.
|
|
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
|
|
bis.R1(6.0f).R2(8.0f).R3(12.0f);
|
|
experiment.ImportBraggIntegrationSettings(bis);
|
|
logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=12.0 (override with --integration-radius)");
|
|
}
|
|
|
|
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)");
|
|
}
|
|
|
|
if (background_trim_arg) {
|
|
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
|
|
bis.BackgroundTrimFraction(static_cast<float>(*background_trim_arg));
|
|
experiment.ImportBraggIntegrationSettings(bis);
|
|
logger.Info("Background-ring trimmed-mean fraction set to {:.2f} (rotation; stills keep the sigma-clip)",
|
|
*background_trim_arg);
|
|
}
|
|
|
|
SpotFindingSettings spot_settings;
|
|
spot_settings.enable = true;
|
|
spot_settings.indexing = true;
|
|
spot_settings.signal_to_noise_threshold = sigma_spot_finding;
|
|
spot_settings.photon_count_threshold = photon_count_threshold_spot_finding;
|
|
// Detection defaults differ by workflow; each is overridden by its flag, which always wins.
|
|
// - min-pix: choosing it per image (unset) only means something where each frame is indexed on its
|
|
// own. Rotation indexing builds ONE lattice from all frames, so it keeps the fixed value.
|
|
// - adaptive detection: a clear win across the stills battery, but on the 33-crystal rotation
|
|
// battery it deterministically breaks three (a lost space group, a halved indexing rate, a
|
|
// collapsed merge) while helping four, so rotation keeps the fixed-threshold finder.
|
|
// - high-resolution limit: unset means "as far as the detector reaches", which is what the stills
|
|
// path wants; on rotation the extra high-resolution spots cost indexing (measured 100.0 -> 96.8%
|
|
// on a strong rotation set), so rotation keeps the historical limit.
|
|
spot_settings.min_pix_per_spot = min_pix_per_spot;
|
|
if (rotation_indexing && !spot_settings.min_pix_per_spot.has_value())
|
|
spot_settings.min_pix_per_spot = 2;
|
|
if (rotation_indexing && !d_min_spot_finding.has_value())
|
|
d_min_spot_finding = 1.5f;
|
|
spot_settings.adaptive_threshold = adaptive_spots.value_or(true);
|
|
spot_settings.high_resolution_limit = d_min_spot_finding;
|
|
spot_settings.false_pixels_per_frame = false_pixels_per_frame;
|
|
if (d_max_spot_finding > 0.0f)
|
|
spot_settings.low_resolution_limit = d_max_spot_finding;
|
|
|
|
// Validate the assembled spot-finding settings the same way the online receivers do (broker and
|
|
// receiver call this same function). It enforces the cross-field constraints that per-argument
|
|
// bounds cannot express - in particular that the low-resolution limit is coarser than the
|
|
// high-resolution limit, so --spot-low-resolution below the high-res cut no longer silently
|
|
// rejects every pixel.
|
|
try {
|
|
DiffractionExperiment::CheckDataProcessingSettings(spot_settings);
|
|
} catch (const std::exception &e) {
|
|
logger.Error("Invalid spot-finding settings: {}", e.what());
|
|
return 1;
|
|
}
|
|
|
|
// 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_postrefine_geometry = rotation_postrefine_geometry;
|
|
config.rotation_indexing_image_count = rotation_indexing_image_count;
|
|
config.forced_rotation_lattice = forced_rotation_lattice;
|
|
config.refine_geometry = refine_geometry;
|
|
config.run_scaling = run_scaling;
|
|
config.scaling_iter = scaling_iter;
|
|
config.reference_data = reference_data;
|
|
config.reference_has_free_flags = reference_has_free_flags;
|
|
config.observation_dump_path = dump_observations;
|
|
config.model_path = model_pdb;
|
|
// 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;
|
|
|
|
// The space-group search is rendered here (not in the library) so the viewer does not emit it on
|
|
// stdout and the CLI owns the format.
|
|
if (result.space_group_search.has_value())
|
|
std::cout << std::endl << SearchSpaceGroupResultToText(*result.space_group_search) << std::endl;
|
|
|
|
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;
|
|
// Final one-line summary of the adopted crystal (whether de-novo determined or fixed with -S),
|
|
// so it is not buried in the space-group-search block (which is de-novo only) or only in the mmCIF.
|
|
if (result.space_group_number.has_value()) {
|
|
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_number(
|
|
static_cast<int>(result.space_group_number.value()));
|
|
std::cout << fmt::format("Space group: {} (No. {})", sg ? sg->short_name() : "?",
|
|
result.space_group_number.value()) << 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);
|
|
}
|