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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit 992343d0206489991880bb927efc79131134a1d9)

common/JFJochMessages.h

@@ -14,6 +14,7 @@
 #include <nlohmann/json.hpp>
 
 #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<std::array<int32_t, 9>> reindex_matrix;
+    // Stills: the rocking curve the partialities were computed from (and the frames predicted through).
+    std::optional<StillsRunWidth> 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

common/StillsRunWidth.h

@@ -0,0 +1,18 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <vector>
+
+// 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<double> table_dstar, table_s;
+};

image_analysis/scale_merge/StillsPartialityRefine.cpp

@@ -492,7 +492,7 @@ double StillsPartialityRefine::MeasureRadialScale(const std::vector<IntegrationO
 
 std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureRunWidth(
         const std::vector<IntegrationOutcome> &outcomes, const std::optional<UnitCell> &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::RunWidth> 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<Point> pts;
@@ -519,10 +524,10 @@ std::optional<StillsPartialityRefine::RunWidth> 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::RunWidth> 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<double>(r.d) * r.d), de * de, w, sigma_w * sigma_w, window});
         }
     }
@@ -707,7 +720,8 @@ std::optional<StillsPartialityRefine::RunWidth> StillsPartialityRefine::MeasureR
 
 double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, size_t nthreads,
                                    const std::optional<UnitCell> &reference_cell,
-                                   const std::optional<RunWidth> &given_width) const {
+                                   const std::optional<RunWidth> &given_width,
+                                   std::optional<RunWidth> *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<IntegrationOutcome> &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.

image_analysis/scale_merge/StillsPartialityRefine.h

@@ -7,6 +7,7 @@
 #include <vector>
 
 #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<double> 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<IntegrationOutcome> &outcomes, size_t nthreads,
                const std::optional<UnitCell> &reference_cell,
-               const std::optional<RunWidth> &run_width = std::nullopt) const;
+               const std::optional<RunWidth> &run_width = std::nullopt,
+               std::optional<RunWidth> *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<RunWidth> MeasureRunWidth(const std::vector<IntegrationOutcome> &outcomes,
                                             const std::optional<UnitCell> &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

reader/HDF5MetadataSource.cpp

@@ -957,6 +957,18 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen
         if (const auto m = master_file->ReadOptVector<int32_t>("/entry/MX/reindexMatrix"); m.size() == 9)
             dataset->reindex_matrix = std::array<int32_t, 9>{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<double>("/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<double>("/entry/MX/stillsRadialScale"); r.size() == 1)
+                w.radial_scale = r[0];
+            w.table_dstar = master_file->ReadOptVector<double>("/entry/MX/stillsRockingWidthTableDstar");
+            w.table_s = master_file->ReadOptVector<double>("/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"));
 
 

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<std::array<int32_t, 9>> reindex_matrix;
+    // Stills: the rocking curve the run's partialities were computed from (/entry/MX/stillsRockingWidth*).
+    std::optional<StillsRunWidth> stills_run_width;
 
     std::vector<float> az_int_bin_to_q;
     std::vector<float> az_int_bin_to_phi;

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

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

writer/HDF5NXmx.cpp

@@ -1194,6 +1194,19 @@ void NXmx::Finalize(const EndMessage &end) {
             SaveVector(*hdf5_file, "/entry/MX/reindexMatrix",
                        std::vector<int32_t>(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<double>{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());
         }