From 8a123e794de1806e148ea7308a79836bf990fc51 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 26 Oct 2026 11:34:56 +0100 Subject: [PATCH] Stills: the run's rocking curve goes to _process.h5 and --mode scale merges with it; MeasureRunWidth takes the window at each reflection's resolution --mode scale re-measured the rocking width on the stored reflections, so on any file whose frames were predicted through the run's curve it merged with a different curve than the mx run did. The curve the post-refinement used (formula, table, radial scale) is now carried in the EndMessage (StillsRunWidth) and written to the master: /entry/MX/stillsRockingWidthFormula [gamma0 1/A, gamma_e rad, bandwidth sigma], /entry/MX/stillsRadialScale, /entry/MX/stillsRockingWidthTableDstar and /entry/MX/stillsRockingWidthTable (1/A). The reader returns it, and --mode scale passes it to the post-refinement. A file without it is re-measured as before. MeasureRunWidth grouped each frame's reflections by the frame's window, taken as the largest excitation error it predicted. With a window that follows sigma(d*) that is the high-resolution window, applied to the frame's low-resolution reflections too. Given the curve the frames were predicted through, the window is now that curve at the reflection's own resolution; without one, the frame's largest excitation error as before. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit 992343d0206489991880bb927efc79131134a1d9) --- common/JFJochMessages.h | 3 ++ common/StillsRunWidth.h | 18 ++++++++++++ .../scale_merge/StillsPartialityRefine.cpp | 24 ++++++++++++--- .../scale_merge/StillsPartialityRefine.h | 29 ++++++++++--------- reader/HDF5MetadataSource.cpp | 12 ++++++++ reader/JFJochReaderDataset.h | 3 ++ rugnux/RugnuxScaleMerge.cpp | 2 +- rugnux/rugnux_cli.cpp | 5 +++- writer/HDF5NXmx.cpp | 13 +++++++++ 9 files changed, 89 insertions(+), 20 deletions(-) create mode 100644 common/StillsRunWidth.h diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index e0da5f334..ac9255090 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -14,6 +14,7 @@ #include #include "CompressedImage.h" +#include "StillsRunWidth.h" #include "SpotToSave.h" #include "UnitCell.h" #include "MtzTable.h" @@ -455,6 +456,8 @@ struct EndMessage { // setting the merge works in - is only chosen afterwards, so the two settings are not always the // same one. Absent means they are, i.e. the identity. std::optional> reindex_matrix; + // Stills: the rocking curve the partialities were computed from (and the frames predicted through). + std::optional stills_run_width; // Geometry refined by the offline analysis. When present these overwrite the master-file geometry // at finalization (rugnux writes the refined values; the broker leaves them empty and keeps the diff --git a/common/StillsRunWidth.h b/common/StillsRunWidth.h new file mode 100644 index 000000000..41a27570c --- /dev/null +++ b/common/StillsRunWidth.h @@ -0,0 +1,18 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +// A stills run's rocking curve as measured by StillsPartialityRefine::MeasureRunWidth (see there): the +// three-term width formula sigma^2 = gamma0^2 + (gamma_e d*)^2 + (bandwidth_sigma |q_z|)^2 (1/A, rad, +// relative), the width measured in shells of d* (table_s at the shells' mean d*, 1/A; empty when too few +// shells could be measured), and the radial scale of the reciprocal lattice the curve is centred on. +// Carried in the EndMessage so the _process.h5 holds the curve its frames were predicted through and +// --mode scale merges with it rather than measuring another. +struct StillsRunWidth { + double gamma0 = 0.0, gamma_e = 0.0, bandwidth_sigma = 0.0; + double radial_scale = 1.0; + std::vector table_dstar, table_s; +}; diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 35d55eb5f..6096da7f4 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -492,7 +492,7 @@ double StillsPartialityRefine::MeasureRadialScale(const std::vector StillsPartialityRefine::MeasureRunWidth( const std::vector &outcomes, const std::optional &reference_cell, - double radial_scale) const { + double radial_scale, const RunWidth *window_width) const { // The rocking width measured on the crystals rather than on the prediction window. A lattice point // sits anywhere across the window with equal likelihood, while what it records is I_full * p(de), so // the intensity-weighted second moment of the excitation error, sum I de^2 / sum I, is the width s^2 @@ -510,6 +510,11 @@ std::optional StillsPartialityRefine::MeasureR // WINDOW_GROUPS groups of similar window, each measured against its own, and the shell's width is the // median of the groups whose curve is narrower than their window; a shell where no group can see the // curve is left out. + // + // A frame's window is the largest excitation error it predicted where it predicted through one width. + // Predicted through the run's width (window_width), the window is its formula at the reflection's own + // resolution: the frame's largest excitation error is then the high-resolution window, and taken as the + // window of its low-resolution reflections it put them in the wrong group. const auto &indexing = experiment_.GetIndexingSettings(); struct Point { double x, de2, w, var_w, window; }; std::vector pts; @@ -519,10 +524,10 @@ std::optional StillsPartialityRefine::MeasureR if (reference_cell && !o.latt.GetUnitCell().is_close(*reference_cell, indexing.GetUnitCellDistTolerance(), indexing.GetUnitCellAngleTolerance_deg())) continue; - float window = 0.0f; // the largest excitation error the frame predicted + float frame_window = 0.0f; // the largest excitation error the frame predicted for (const auto &r: o.reflections) if (std::isfinite(r.dist_ewald)) - window = std::max(window, std::fabs(r.dist_ewald)); + frame_window = std::max(frame_window, std::fabs(r.dist_ewald)); for (const auto &r: o.reflections) { if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_)) continue; @@ -539,6 +544,14 @@ std::optional StillsPartialityRefine::MeasureR const Coord S0 = o.geom.GetScatteringVector(); de = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - 1.0 / o.geom.GetWavelength_A(); } + double window = frame_window; + if (window_width) { + // the formula at the reflection's d*, its bandwidth term at |q_z| = lambda d*^2 / 2 + const double dstar = 1.0 / r.d; + const double a = window_width->gamma_e * dstar; + const double b = window_width->bandwidth_sigma * o.geom.GetWavelength_A() * dstar * dstar / 2.0; + window = std::sqrt(window_width->gamma0 * window_width->gamma0 + a * a + b * b); + } pts.push_back({1.0 / (static_cast(r.d) * r.d), de * de, w, sigma_w * sigma_w, window}); } } @@ -707,7 +720,8 @@ std::optional StillsPartialityRefine::MeasureR double StillsPartialityRefine::Run(std::vector &outcomes, size_t nthreads, const std::optional &reference_cell, - const std::optional &given_width) const { + const std::optional &given_width, + std::optional *used_width) const { const auto &indexing = experiment_.GetIndexingSettings(); auto merged_crystal = [&](const IntegrationOutcome &o) { return !reference_cell || o.latt.GetUnitCell().is_close(*reference_cell, @@ -733,6 +747,8 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si run_width = MeasureRunWidth(outcomes, reference_cell, radial_scale); } const double radial_scale = run_width ? run_width->radial_scale : 1.0; + if (used_width) + *used_width = run_width; for (int outer = 0; outer < settings_.outer_iterations; ++outer) { // Reference full intensities from the current corrections. diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index 38c63f43c..9ba83f1bd 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -7,6 +7,7 @@ #include #include "../../common/DiffractionExperiment.h" +#include "../../common/StillsRunWidth.h" #include "../IntegrationOutcome.h" #include "HKLKey.h" @@ -45,28 +46,28 @@ public: explicit StillsPartialityRefine(const DiffractionExperiment &x); - // The run's rocking width, sigma^2 = gamma0^2 + (gamma_e d*)^2 + (bandwidth_sigma |q_z|)^2 (1/A, rad, - // relative), measured from the intensity-weighted excitation errors with the reciprocal lattice scaled - // by radial_scale (the curve's centre, MeasureRadialScale); empty when the run has too few reflections - // to measure it. - // table_s is the width measured in shells of equal width in d* (at the shells' mean d*, non-decreasing), - // which the partiality uses where it is not empty; the formula is what the prediction window uses. - struct RunWidth { - double gamma0, gamma_e, bandwidth_sigma; - double radial_scale = 1.0; - std::vector table_dstar, table_s; - }; + // The run's rocking width (common/StillsRunWidth.h): the three-term formula, which the prediction + // window uses, and the table in d* the partiality uses where it is not empty; measured from the + // intensity-weighted excitation errors with the reciprocal lattice scaled by radial_scale (the curve's + // centre, MeasureRadialScale); empty when the run has too few reflections to measure it. + using RunWidth = StillsRunWidth; // Refine all crystals in place. Returns the mean |dpsi| applied (degrees), for diagnostics. // Only the crystals whose cell is close to reference_cell (all without one) are refined and make the // reference: the others are not merged. The partiality width is run_width where one is given (the - // width the frames were predicted with), else measured here (MeasureRunWidth). + // width the frames were predicted with), else measured here (MeasureRunWidth); the one used is put in + // used_width where that is given. double Run(std::vector &outcomes, size_t nthreads, const std::optional &reference_cell, - const std::optional &run_width = std::nullopt) const; + const std::optional &run_width = std::nullopt, + std::optional *used_width = nullptr) const; + // window: the width the frames were predicted through (|dist_ewald| <= k times its formula at the + // reflection's d*), where they were; without one each frame's window is taken as the largest excitation + // error it predicted. std::optional MeasureRunWidth(const std::vector &outcomes, const std::optional &reference_cell, - double radial_scale = 1.0) const; + double radial_scale = 1.0, + const RunWidth *window = nullptr) const; // The run's radial scale of the reciprocal lattice (a wavelength or cell scale the spot positions cannot // see), measured from where the rocking curve is centred against d*^2: one step from radial_scale, so diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 8b2299cee..f6fea51fa 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -957,6 +957,18 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen if (const auto m = master_file->ReadOptVector("/entry/MX/reindexMatrix"); m.size() == 9) dataset->reindex_matrix = std::array{m[0], m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8]}; + if (const auto f = master_file->ReadOptVector("/entry/MX/stillsRockingWidthFormula"); f.size() == 3) { + StillsRunWidth w{.gamma0 = f[0], .gamma_e = f[1], .bandwidth_sigma = f[2]}; + if (const auto r = master_file->ReadOptVector("/entry/MX/stillsRadialScale"); r.size() == 1) + w.radial_scale = r[0]; + w.table_dstar = master_file->ReadOptVector("/entry/MX/stillsRockingWidthTableDstar"); + w.table_s = master_file->ReadOptVector("/entry/MX/stillsRockingWidthTable"); + if (w.table_dstar.size() != w.table_s.size()) { + w.table_dstar.clear(); + w.table_s.clear(); + } + dataset->stills_run_width = w; + } dataset->experiment.SampleName(master_file->GetString("/entry/sample/name")); diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 9d82b2b5a..5c01a5248 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -12,6 +12,7 @@ #include "../common/DiffractionGeometry.h" #include "../common/DiffractionExperiment.h" +#include "../common/StillsRunWidth.h" #include "../common/PixelMask.h" #include "../common/AzimuthalIntegrationMapping.h" @@ -41,6 +42,8 @@ struct JFJochReaderDataset { // /entry/MX/reindexMatrix; the reader applies it, so everything this class hands out is already in // the cell's setting. std::optional> reindex_matrix; + // Stills: the rocking curve the run's partialities were computed from (/entry/MX/stillsRockingWidth*). + std::optional stills_run_width; std::vector az_int_bin_to_q; std::vector az_int_bin_to_phi; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 2ef752eda..3e5c3ca89 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -1431,7 +1431,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { TimingMark refine_mark(logger, "stills partiality post-refine " + label); StillsPartialityRefine refiner(experiment_); const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads, result.consensus_cell, - stills_run_width_); + stills_run_width_, &end_msg.stills_run_width); logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); // The refine fits its OWN per-image scale, with no floor, and adopts it unconditionally // when the image has no prior CC - which is exactly the state the guard above leaves. So diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 7fd2e06be..6b1114868 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2586,7 +2586,10 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { // (no re-integration), then merge with the improved corrections. if (experiment.GetScalingSettings().GetStillsPartialityRefine()) { StillsPartialityRefine refiner(experiment); - const double mean_tilt = refiner.Run(reflections, nthreads, experiment.GetUnitCell()); + // The curve the mx run computed its partialities from (and predicted the frames through), + // where the file carries it: measured again on frames predicted through it, it is another. + const double mean_tilt = refiner.Run(reflections, nthreads, experiment.GetUnitCell(), + dataset->stills_run_width); logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt); // As in the full pipeline: the refine fits its own per-image scale, with no floor. ScaleOnTheFly::RejectCollapsedScales(reflections); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 205652468..550fae738 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1194,6 +1194,19 @@ void NXmx::Finalize(const EndMessage &end) { SaveVector(*hdf5_file, "/entry/MX/reindexMatrix", std::vector(end.reindex_matrix->begin(), end.reindex_matrix->end())); + // Stills: the run's rocking curve, so that --mode scale merges with the curve the frames were + // predicted through instead of measuring another on them. + if (end.stills_run_width) { + const auto &w = *end.stills_run_width; + SaveVector(*hdf5_file, "/entry/MX/stillsRockingWidthFormula", + std::vector{w.gamma0, w.gamma_e, w.bandwidth_sigma}); + SaveScalar(*hdf5_file, "/entry/MX/stillsRadialScale", w.radial_scale); + if (!w.table_s.empty()) { + SaveVector(*hdf5_file, "/entry/MX/stillsRockingWidthTableDstar", w.table_dstar)->Units("Angstrom^-1"); + SaveVector(*hdf5_file, "/entry/MX/stillsRockingWidthTable", w.table_s)->Units("Angstrom^-1"); + } + } + if (end.indexing_rate) { SaveScalar(*hdf5_file, "/entry/MX/imageIndexedMean", end.indexing_rate.value()); }