Add physical partiality post-refinement for stills (default on)

Replace the frozen scalar-sigma stills partiality with a physical, refined model.
Per crystal, refine an orientation tilt (dpsi_x, dpsi_y) against the running merge
and recompute each reflection's partiality analytically from the refined geometry
(angular Ewald-proximity model, sigma(d*) = gamma_e*d*), with the per-crystal scale
G profiled out by the existing robust IRLS - no re-integration. A soft Gaussian
prior on dpsi tames weak-data overfit while staying inert on strong data. The
merge <-> refine loop iterates a few times.

This is now the stills default via ScalingSettings::stills_partiality_refine (on).
A single opt-out flag `--simple-stills` reverts to treating every reflection as a
full (p=1, single pass). Retires the experimental `--still-partiality` flag. The
viewer gains a "Partiality post-refinement (stills)" checkbox in Scaling settings.

Validated (integrate-once / --scale): CC1/2 and R_meas both improve on lyso8
(+2.8 / -10), KR2 (+5.6 / -3.4) and LOV (+2.1 / -4); neutral on a pink-beam DMM
set (already-full reflections); R-free/R-work down vs a fixed model; competitive
with CrystFEL partialator on matched frames.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

common/ScalingSettings.cpp

@@ -160,6 +160,15 @@ bool ScalingSettings::GetStillsModulation() const {
     return stills_modulation;
 }
 
+ScalingSettings &ScalingSettings::StillsPartialityRefine(bool input) {
+    stills_partiality_refine = input;
+    return *this;
+}
+
+bool ScalingSettings::GetStillsPartialityRefine() const {
+    return stills_partiality_refine;
+}
+
 ScalingSettings &ScalingSettings::ExpectedVarianceMerge(bool input) {
     expected_variance_merge = input;
     return *this;

common/ScalingSettings.h

@@ -51,6 +51,13 @@ class ScalingSettings {
     // rotation path fits its own modulation via RotationScaleMerge). Enabled by rugnux --stills-modulation.
     bool stills_modulation = false;
 
+    // Physical partiality post-refinement for the STILLS merge (StillsPartialityRefine): refine a per-crystal
+    // orientation tilt against the running merge, recompute each reflection's partiality from the refined
+    // geometry (angular Ewald-proximity model), and re-scale/merge - the "full model" for stills. ON by
+    // default (helps mono stills, neutral on pink beam, tames weak data via a soft prior). rugnux
+    // --simple-stills turns it OFF, reverting to treating every reflection as a full (p = 1, single pass).
+    bool stills_partiality_refine = true;
+
     // Expected-variance merge weighting for the STILLS merge (MergeOnTheFly). When combining a reflection's
     // redundant observations by inverse variance, rebuild the Poisson signal part of each observation's
     // variance at the reflection's EXPECTED <I> instead of the observation's own intensity. Weighting by an
@@ -101,6 +108,7 @@ public:
     ScalingSettings& AbsorptionIter(int input);
     ScalingSettings& CorrectionSurfaces(bool input);
     ScalingSettings& StillsModulation(bool input);
+    ScalingSettings& StillsPartialityRefine(bool input);
     ScalingSettings& ExpectedVarianceMerge(bool input);
     ScalingSettings& SmoothGDegrees(double input);
     ScalingSettings& RelativeBDegrees(double input);
@@ -141,6 +149,7 @@ public:
     [[nodiscard]] int GetAbsorptionIter() const;
     [[nodiscard]] bool GetCorrectionSurfaces() const;
     [[nodiscard]] bool GetStillsModulation() const;
+    [[nodiscard]] bool GetStillsPartialityRefine() const;
     [[nodiscard]] bool GetExpectedVarianceMerge() const;
     [[nodiscard]] double GetSmoothGDegrees() const;
     [[nodiscard]] double GetRelativeBDegrees() const;

image_analysis/scale_merge/CMakeLists.txt

@@ -7,6 +7,8 @@ ADD_LIBRARY(JFJochScaleMerge
         Merge.h
         ScaleOnTheFly.cpp
         ScaleOnTheFly.h
+        StillsPartialityRefine.cpp
+        StillsPartialityRefine.h
         RotationScaleMerge.cpp
         RotationScaleMerge.h
         ResolutionCutoff.cpp

image_analysis/scale_merge/StillsPartialityRefine.cpp

@@ -0,0 +1,329 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "StillsPartialityRefine.h"
+
+#include <atomic>
+#include <cmath>
+#include <future>
+#include <thread>
+#include <vector>
+
+#include <ceres/ceres.h>
+#include <ceres/rotation.h>
+
+#include "Merge.h"
+
+namespace {
+    constexpr size_t MIN_FIT_REFLECTIONS = 20;
+    constexpr double kRadToDeg = 180.0 / 3.14159265358979323846;
+
+    double SafeInv(double x, double fallback) {
+        if (!std::isfinite(x) || x == 0.0)
+            return fallback;
+        return 1.0 / x;
+    }
+
+    // One accepted reflection reduced to the physical partiality fit: the base reciprocal vector q (crystal
+    // frame at the stored per-image orientation), the reference full intensity, the measured intensity, the
+    // Lorentz factor and the weight. dist_ewald / partiality are recomputed from q under the refined tilt.
+    struct FitObs {
+        double qx, qy, qz;
+        double Iref;
+        double Iobs;
+        double lp;           // 1 / rlp
+        double weight;       // 1 / sigma
+    };
+
+    double VecLen(double x, double y, double z) { return std::sqrt(x * x + y * y + z * z); }
+
+    // Analytic partiality for a reflection whose base reciprocal vector is q, tilted by (psi_x, psi_y).
+    // Mirrors BraggPrediction: dist_ewald = |S| - 1/lambda with S = q_rot + S0, and
+    // p = exp(-dist_ewald^2 / 2 sigma^2), sigma^2 = (gamma0 + gamma_e*d*)^2 + (bw*|q_z|)^2.
+    double ComputeP(double qx, double qy, double qz,
+                    double psi_x, double psi_y,
+                    double s0x, double s0y, double s0z, double inv_lambda,
+                    double gamma0, double gamma_e, double bw) {
+        double aa[3] = {psi_x, psi_y, 0.0};
+        double q[3] = {qx, qy, qz};
+        double qr[3];
+        ceres::AngleAxisRotatePoint(aa, q, qr);
+        const double Sx = qr[0] + s0x, Sy = qr[1] + s0y, Sz = qr[2] + s0z;
+        const double de = std::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - inv_lambda;
+        const double dstar = std::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]);
+        const double sig = gamma0 + gamma_e * dstar;
+        const double sbw = bw * std::fabs(qr[2]);
+        const double sig2 = sig * sig + sbw * sbw;
+        if (!(sig2 > 0.0))
+            return 1.0;
+        return std::exp(-0.5 * de * de / sig2);
+    }
+
+    // Robust per-crystal scale G (linear in G given the model coefficients), identical objective to
+    // ScaleOnTheFly::SolveScaleIRLS: minimise sum Cauchy_k( w (G*coeff - Iobs) ) over G >= 0.
+    double SolveScaleIRLS(const std::vector<double> &coeff, const std::vector<double> &Iobs,
+                          const std::vector<double> &weight, double robust_k) {
+        auto weighted_scale = [&](auto robust_weight) {
+            double num = 0.0, den = 0.0;
+            for (size_t i = 0; i < coeff.size(); ++i) {
+                const double rw = robust_weight(i);
+                const double w2 = weight[i] * weight[i];
+                num += rw * w2 * coeff[i] * Iobs[i];
+                den += rw * w2 * coeff[i] * coeff[i];
+            }
+            return den > 0.0 ? num / den : NAN;
+        };
+
+        double G = weighted_scale([](size_t) { return 1.0; });
+        if (!std::isfinite(G))
+            return 1.0;
+        G = std::max(0.0, G);
+
+        const double k2 = robust_k * robust_k;
+        for (int iter = 0; iter < 30; ++iter) {
+            const double G_prev = G;
+            const double G_next = weighted_scale([&](size_t i) {
+                const double res = weight[i] * (G * coeff[i] - Iobs[i]);
+                return 1.0 / (1.0 + res * res / k2);
+            });
+            if (!std::isfinite(G_next))
+                break;
+            G = std::max(0.0, G_next);
+            if (std::abs(G - G_prev) <= 1e-7 * std::max(G, 1.0))
+                break;
+        }
+        return G;
+    }
+
+    // Ceres residual: refine the orientation tilt (psi_x, psi_y) holding the scale G fixed. The tilt
+    // rotates the base reciprocal vector q; partiality follows analytically. Residual is the intensity
+    // mismatch weighted by 1/sigma, exactly matching ScaleOnTheFly's intensity-space objective.
+    struct PsiResidual {
+        double qx, qy, qz;
+        double s0x, s0y, s0z, inv_lambda;
+        double gamma0, gamma_e, bw;
+        double G, lp, Iref, Iobs, weight;
+
+        template<typename T>
+        bool operator()(const T *const psi, T *residual) const {
+            T q[3] = {T(qx), T(qy), T(qz)};
+            T aa[3] = {psi[0], psi[1], T(0.0)};
+            T qr[3];
+            ceres::AngleAxisRotatePoint(aa, q, qr);
+            const T Sx = qr[0] + T(s0x), Sy = qr[1] + T(s0y), Sz = qr[2] + T(s0z);
+            const T de = ceres::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - T(inv_lambda);
+            const T dstar = ceres::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]);
+            const T sig = T(gamma0) + T(gamma_e) * dstar;
+            const T sbw = T(bw) * ceres::abs(qr[2]);
+            const T sig2 = sig * sig + sbw * sbw;
+            const T p = ceres::exp(T(-0.5) * de * de / sig2);
+            residual[0] = T(weight) * (T(G) * p * T(lp) * T(Iref) - T(Iobs));
+            return true;
+        }
+    };
+
+    // Gaussian prior N(0, sigma_prior^2) on the tilt. The data residuals above are (model-obs)/sigma, a
+    // proper chi^2, so the MAP prior residual is simply dpsi/sigma_prior - no scale calibration needed.
+    struct PsiPrior {
+        double inv_sigma;
+        template<typename T>
+        bool operator()(const T *const psi, T *residual) const {
+            residual[0] = T(inv_sigma) * psi[0];
+            residual[1] = T(inv_sigma) * psi[1];
+            return true;
+        }
+    };
+}
+
+StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x)
+    : StillsPartialityRefine(x, Settings{}) {}
+
+StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x, Settings settings)
+    : experiment_(x),
+      settings_(settings),
+      hkl_key_generator_(x.GetScalingSettings().GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)),
+      d_min_limit_(x.GetScalingSettings().GetHighResolutionLimit_A()),
+      bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {}
+
+double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome,
+                                         const std::map<HKLKey, double> &reference) const {
+    if (outcome.reflections.empty())
+        return 0.0;
+
+    const Coord Astar = outcome.latt.Astar();
+    const Coord Bstar = outcome.latt.Bstar();
+    const Coord Cstar = outcome.latt.Cstar();
+    const Coord S0 = outcome.geom.GetScatteringVector();
+    const double inv_lambda = 1.0 / outcome.geom.GetWavelength_A();
+    const double bw = bandwidth_sigma_;
+    const double gamma0 = 0.0;   // width is purely angular: sigma(d*) = gamma_e * d* (set per crystal below)
+
+    auto base_q = [&](const Reflection &r) {
+        return Astar * static_cast<float>(r.h) + Bstar * static_cast<float>(r.k)
+               + Cstar * static_cast<float>(r.l);
+    };
+
+    // Collect the reflections that constrain the fit (accepted, non-ice, finite, present in the reference).
+    std::vector<FitObs> obs;
+    obs.reserve(outcome.reflections.size());
+    double sum_ang2 = 0.0;   // RMS angular excitation error (dist_ewald / d*) -> per-crystal mosaic width
+    size_t n_de = 0;
+    for (const Reflection &r: outcome.reflections) {
+        if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_))
+            continue;
+        if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f)
+            continue;
+        const auto it = reference.find(hkl_key_generator_(r));
+        if (it == reference.end() || !std::isfinite(it->second))
+            continue;
+
+        const Coord q = base_q(r);
+        obs.push_back(FitObs{
+            .qx = q.x, .qy = q.y, .qz = q.z,
+            .Iref = it->second,
+            .Iobs = static_cast<double>(r.I),
+            .lp = SafeInv(r.rlp, 1.0),
+            .weight = SafeInv(r.sigma, 1.0),
+        });
+
+        // Angular excitation error delta_psi = dist_ewald / d* at the stored orientation (psi = 0). Using
+        // the ANGULAR distance (not the linear reciprocal-space distance) makes the partiality width
+        // resolution-clean: a fixed mosaic angle smears high-resolution rlps more in reciprocal space, so a
+        // constant linear width computes p too small at high resolution and over-divides those shells.
+        const double dstar = VecLen(q.x, q.y, q.z);
+        const double de0 = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - inv_lambda;
+        if (dstar > 1e-9) {
+            const double dpsi = de0 / dstar;
+            sum_ang2 += dpsi * dpsi;
+            ++n_de;
+        }
+    }
+
+    if (obs.size() < MIN_FIT_REFLECTIONS || n_de == 0)
+        return 0.0;
+
+    // Per-crystal angular mosaic width from the RMS angular excitation error. sigma(d*) = gamma_e * d*
+    // (gamma0 = 0), i.e. p = exp(-0.5 (delta_psi / gamma_e)^2) is a Gaussian in the angular distance from
+    // the Ewald sphere - the physical mosaic/divergence model, independent of resolution. A positive
+    // settings_.gamma_e overrides the per-crystal estimate with a shared (pooled) width.
+    const double gamma_e_ang = std::max(std::sqrt(sum_ang2 / static_cast<double>(n_de)), 1e-9);
+    const double gamma_e = settings_.gamma_e > 0.0 ? settings_.gamma_e : gamma_e_ang;
+
+    double psi[2] = {0.0, 0.0};
+    double G = 1.0;
+
+    const bool refine_tilt = obs.size() >= settings_.min_reflections;
+    const int inner = refine_tilt ? settings_.inner_iterations : 1;
+
+    for (int it = 0; it < inner; ++it) {
+        // (1) Solve G given the current partialities.
+        std::vector<double> coeff(obs.size()), Iobs(obs.size()), weight(obs.size());
+        for (size_t j = 0; j < obs.size(); ++j) {
+            const double p = ComputeP(obs[j].qx, obs[j].qy, obs[j].qz, psi[0], psi[1],
+                                      S0.x, S0.y, S0.z, inv_lambda, gamma0, gamma_e, bw);
+            coeff[j] = p * obs[j].lp * obs[j].Iref;
+            Iobs[j] = obs[j].Iobs;
+            weight[j] = obs[j].weight;
+        }
+        G = SolveScaleIRLS(coeff, Iobs, weight, settings_.robust_k);
+        if (!(G > 0.0) || !std::isfinite(G))
+            return 0.0;
+
+        if (!refine_tilt)
+            break;
+
+        // (2) Refine the tilt holding G fixed.
+        ceres::Problem problem;
+        for (const auto &o: obs) {
+            auto *cost = new ceres::AutoDiffCostFunction<PsiResidual, 1, 2>(new PsiResidual{
+                .qx = o.qx, .qy = o.qy, .qz = o.qz,
+                .s0x = S0.x, .s0y = S0.y, .s0z = S0.z, .inv_lambda = inv_lambda,
+                .gamma0 = gamma0, .gamma_e = gamma_e, .bw = bw,
+                .G = G, .lp = o.lp, .Iref = o.Iref, .Iobs = o.Iobs, .weight = o.weight});
+            problem.AddResidualBlock(cost, new ceres::CauchyLoss(settings_.robust_k), psi);
+        }
+        if (settings_.prior_sigma_deg > 0.0) {
+            const double inv_sigma = kRadToDeg / settings_.prior_sigma_deg;  // 1 / sigma_prior (rad)
+            problem.AddResidualBlock(new ceres::AutoDiffCostFunction<PsiPrior, 2, 2>(
+                                         new PsiPrior{inv_sigma}), nullptr, psi);
+        }
+        problem.SetParameterLowerBound(psi, 0, -settings_.max_tilt_rad);
+        problem.SetParameterUpperBound(psi, 0, settings_.max_tilt_rad);
+        problem.SetParameterLowerBound(psi, 1, -settings_.max_tilt_rad);
+        problem.SetParameterUpperBound(psi, 1, settings_.max_tilt_rad);
+
+        ceres::Solver::Options options;
+        options.linear_solver_type = ceres::DENSE_QR;
+        options.minimizer_progress_to_stdout = false;
+        options.num_threads = 1;
+        options.max_num_iterations = 25;
+        ceres::Solver::Summary summary;
+        ceres::Solve(options, &problem, &summary);
+    }
+
+    // Write the refined partiality + scale correction onto every reflection of the crystal (not only the
+    // fit subset), so the merge sees a consistent model. image_scale_corr = rlp / (partiality * G).
+    for (auto &r: outcome.reflections) {
+        const Coord q = base_q(r);
+        const double p = ComputeP(q.x, q.y, q.z, psi[0], psi[1], S0.x, S0.y, S0.z, inv_lambda,
+                                  gamma0, gamma_e, bw);
+        r.partiality = static_cast<float>(p);
+        const double denom = p * G;
+        r.image_scale_corr = (std::isfinite(r.rlp) && denom > 0.0)
+                                 ? static_cast<float>(r.rlp / denom)
+                                 : NAN;
+    }
+    outcome.image_scale_g = static_cast<float>(G);
+
+    const double tilt_deg = std::sqrt(psi[0] * psi[0] + psi[1] * psi[1]) * kRadToDeg;
+    return tilt_deg;
+}
+
+double StillsPartialityRefine::Run(std::vector<IntegrationOutcome> &outcomes, size_t nthreads) const {
+    if (nthreads == 0)
+        nthreads = std::thread::hardware_concurrency();
+    nthreads = std::max<size_t>(1, nthreads);
+
+    double last_mean_tilt = 0.0;
+
+    for (int outer = 0; outer < settings_.outer_iterations; ++outer) {
+        // Reference full intensities from the current corrections.
+        const std::vector<MergedReflection> merged = MergeAll(experiment_, outcomes, false);
+        std::map<HKLKey, double> reference;
+        for (const auto &m: merged)
+            reference[hkl_key_generator_(m)] = m.I;
+
+        std::atomic<double> tilt_sum{0.0};
+        std::atomic<size_t> tilt_n{0};
+        std::atomic<size_t> next{0};
+
+        auto worker = [&]() {
+            size_t i = next.fetch_add(1);
+            while (i < outcomes.size()) {
+                const double t = RefineOne(outcomes[i], reference);
+                if (t > 0.0) {
+                    double prev = tilt_sum.load();
+                    while (!tilt_sum.compare_exchange_weak(prev, prev + t)) {}
+                    tilt_n.fetch_add(1);
+                }
+                i = next.fetch_add(1);
+            }
+        };
+
+        const size_t nt = std::min(nthreads, std::max<size_t>(1, outcomes.size()));
+        if (nt <= 1) {
+            worker();
+        } else {
+            std::vector<std::future<void>> futures;
+            futures.reserve(nt);
+            for (size_t t = 0; t < nt; ++t)
+                futures.emplace_back(std::async(std::launch::async, worker));
+            for (auto &f: futures)
+                f.get();
+        }
+
+        last_mean_tilt = tilt_n > 0 ? tilt_sum.load() / static_cast<double>(tilt_n.load()) : 0.0;
+    }
+
+    return last_mean_tilt;
+}

image_analysis/scale_merge/StillsPartialityRefine.h

@@ -0,0 +1,59 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <vector>
+
+#include "../../common/DiffractionExperiment.h"
+#include "../IntegrationOutcome.h"
+#include "HKLKey.h"
+
+// Experimental physical partiality post-refinement for STILLS (env JFJOCH_STILL_POSTREFINE).
+//
+// The default stills partiality is a frozen scalar-sigma Gaussian (or p == 1): p is set once at
+// prediction and never optimised, and any attempt to free a per-image sigma jointly with the per-image
+// scale G collapses (within one still the excitation-error spread is narrow, so a scalar sigma is
+// degenerate with G). This class instead parametrises partiality by the crystal ORIENTATION - a small
+// tilt (dpsi_x, dpsi_y) about the two axes perpendicular to the beam - shared by all of a crystal's
+// reflections. A tilt moves each reflection's excitation error by an amount that depends on where the
+// reflection sits on the pattern (one side of the Ewald sphere approaches, the opposite recedes), so it
+// reshapes the SPATIAL pattern of partialities in a way a single G cannot mimic. That breaks the
+// degeneracy that killed the scalar-sigma fit.
+//
+// Because the integrated intensity I_obs is fixed, refining the tilt only recomputes the partiality
+// analytically from the stored per-image lattice/geometry (q = A*.h + B*.k + C*.l, then
+// dist_ewald = | |q + S0| - 1/lambda |) - NO pixel re-integration. The loop is: merge -> per-crystal
+// refine tilt (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and
+// image_scale_corr -> re-merge, iterated a few times. Mutates each reflection's `partiality` and
+// `image_scale_corr` in place; the existing MergeOnTheFly then consumes the improved corrections.
+class StillsPartialityRefine {
+public:
+    struct Settings {
+        int outer_iterations = 2;      // merge <-> refine cycles
+        int inner_iterations = 3;      // (solve G) <-> (refine tilt) alternations per crystal
+        size_t min_reflections = 40;   // skip tilt refinement below this (anti-overfit on sparse crystals)
+        double max_tilt_rad = 0.0175;  // hard bound on |dpsi| (~1 deg); indexing already refined orientation
+        double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0 (0 = off); tames weak-data overfit,
+                                       // inert on strong data (well-supported tilts overcome it)
+        double robust_k = 3.0;         // Cauchy loss scale (sigma units)
+        double gamma_e = 0.0;          // angular width sigma(d*) = gamma_e*d* (0 = estimate per crystal from data)
+    };
+
+    explicit StillsPartialityRefine(const DiffractionExperiment &x);
+    StillsPartialityRefine(const DiffractionExperiment &x, Settings settings);
+
+    // Refine all crystals in place. Returns the mean |dpsi| applied (degrees), for diagnostics.
+    double Run(std::vector<IntegrationOutcome> &outcomes, size_t nthreads = 0) const;
+
+private:
+    const DiffractionExperiment experiment_;
+    const Settings settings_;
+    const HKLKeyGenerator hkl_key_generator_;
+    const std::optional<double> d_min_limit_;
+    const float bandwidth_sigma_;
+
+    // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped).
+    double RefineOne(IntegrationOutcome &outcome,
+                     const std::map<HKLKey, double> &reference) const;
+};

rugnux/Rugnux.cpp

@@ -40,6 +40,7 @@
 #include "../image_analysis/lattice_search/LatticeSearch.h"
 #include "../image_analysis/scale_merge/TwinningAnalysis.h"
 #include "../image_analysis/scale_merge/HKLKey.h"
+#include "../image_analysis/scale_merge/StillsPartialityRefine.h"
 #include "../image_analysis/WriteReflections.h"
 #include "../image_analysis/bragg_integration/CalcISigma.h"
 #include "../common/Definitions.h"
@@ -970,6 +971,15 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
                 auto merge_result = MergeAll(experiment_, indexer->GetIntegrationOutcome(), false);
                 indexer->ScaleAllImages(merge_result);
             }
+            // Physical partiality post-refinement (default on; --simple-stills disables): refine a per-crystal
+            // orientation tilt against the running merge and recompute each reflection's partiality + scale
+            // correction (no re-integration). Never for the P1 search pass.
+            if (!for_search && experiment_.GetScalingSettings().GetStillsPartialityRefine()) {
+                phase("Partiality post-refine (" + label + ")");
+                StillsPartialityRefine refiner(experiment_);
+                const double mean_tilt = refiner.Run(indexer->GetIntegrationOutcome(), config_.nthreads);
+                logger.Info("Stills partiality post-refine: mean |dpsi| = {:.3f} deg", mean_tilt);
+            }
             const std::vector<IntegrationOutcome> &merge_input = indexer->GetIntegrationOutcome();
 
             phase("Merging");

rugnux/RugnuxCommandLine.cpp

@@ -106,8 +106,6 @@ std::string RugnuxCommandLine(const ProcessConfig &config,
         std::ostringstream radii;
         radii << bragg.GetR1() << "," << bragg.GetR2() << "," << bragg.GetR3();
         add("--integration-radius", radii.str());
-        if (bragg.GetStillPartiality())
-            args.emplace_back("--still-partiality");
         // Background trim defaults to 0.10 in the CLI, so emit it whenever the GUI value differs (a
         // custom fraction, or 0 when the box is unchecked) to reproduce the GUI's choice faithfully.
         if (bragg.GetBackgroundTrimFraction() != 0.10f)
@@ -153,6 +151,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config,
                 add("--partiality-uncertainty", num(sc.GetPartialityUncertaintyCoeff()));
             if (sc.GetStillsModulation())
                 args.emplace_back("--stills-modulation");
+            if (!sc.GetStillsPartialityRefine())
+                args.emplace_back("--simple-stills");
             if (!sc.GetExpectedVarianceMerge())
                 args.emplace_back("--no-expected-variance-merge");
             // When merging, the CLI skips the large _process.h5 unless asked; emit the flag when it is

rugnux/rugnux_cli.cpp

@@ -27,6 +27,7 @@
 #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"
@@ -103,7 +104,7 @@ void print_usage() {
     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 << "   --partiality-uncertainty <num>  stills: extra merge sigma ~num*(1-partiality)*<I> on partials (use with --still-partiality; auto-gated to error-model b>1 / ISa<1; default 0, ~2.5 recommended)" << std::endl;
+    std::cout << "   --partiality-uncertainty <num>  stills: extra merge sigma ~num*(1-partiality)*<I> on partials (auto-gated to error-model b>1 / ISa<1; default 0, ~2.5 recommended)" << 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;
@@ -120,7 +121,7 @@ void print_usage() {
     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 << "   --still-partiality              Experimental: weight stills reflections by a Gaussian excitation-error partiality exp(-dist_ewald^2/2sigma^2) instead of treating each as a full" << 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)" << std::endl;
@@ -167,7 +168,7 @@ enum {
     OPT_MODEL,
     OPT_DUMP_OBSERVATIONS,
     OPT_INTEGRATOR,
-    OPT_STILL_PARTIALITY,
+    OPT_SIMPLE_STILLS,
     OPT_SCALE_FULLS,
     OPT_CAPTURE_UNCERTAINTY,
     OPT_PARTIALITY_UNCERTAINTY,
@@ -271,7 +272,7 @@ static option long_options[] = {
     {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
     {"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM},
     {"integrator", required_argument, nullptr, OPT_INTEGRATOR},
-    {"still-partiality", no_argument, nullptr, OPT_STILL_PARTIALITY},
+    {"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},
     {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF},
@@ -547,7 +548,7 @@ int main(int argc, char **argv) {
     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 still_partiality_flag = false;                // --still-partiality (experimental stills partiality)
+    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
     std::optional<double> delta_cchalf_nsigma;         // per-crystal CC1/2-delta rejection
 
@@ -829,8 +830,8 @@ int main(int argc, char **argv) {
             else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical;
             else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; }
             break;
-        case OPT_STILL_PARTIALITY:
-            still_partiality_flag = true;
+        case OPT_SIMPLE_STILLS:
+            simple_stills_flag = true;
             break;
         case OPT_REJECT_OUTLIERS:
             outlier_reject_nsigma = parse_double_arg(optarg, "--reject-outliers", logger);
@@ -1081,6 +1082,7 @@ int main(int argc, char **argv) {
             (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.StillsModulation(stills_modulation_flag);
+        scaling_settings.StillsPartialityRefine(!simple_stills_flag);
         scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
         scaling_settings.OutlierRejectNsigma(
             outlier_reject_nsigma.value_or(
@@ -1137,6 +1139,14 @@ int main(int argc, char **argv) {
             // 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());
             // Optional detector-plane modulation (flat-field) correction, folded into each reflection's
@@ -1394,6 +1404,7 @@ int main(int argc, char **argv) {
     if (no_scaling_corrections)
         scaling_settings.CorrectionSurfaces(false);
     scaling_settings.StillsModulation(stills_modulation_flag);
+    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());
@@ -1469,13 +1480,6 @@ int main(int argc, char **argv) {
                     *background_trim_arg);
     }
 
-    if (still_partiality_flag) {
-        BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
-        bis.StillPartiality(true);
-        experiment.ImportBraggIntegrationSettings(bis);
-        logger.Info("Stills partiality enabled (experimental Gaussian excitation-error weighting)");
-    }
-
     SpotFindingSettings spot_settings;
     spot_settings.enable = true;
     spot_settings.indexing = true;

viewer/widgets/JFJochViewerSettingsDock.cpp

@@ -446,11 +446,6 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() {
     auto *r3 = new NumberLineEdit(1.0f, 30.0f, bragg_.GetR3(), 1, "px", this);
     auto *radii = new QHBoxLayout();
     radii->addWidget(r1); radii->addWidget(r2); radii->addWidget(r3);
-    auto *stillPartiality = new QCheckBox("Stills partiality (experimental)", this);
-    stillPartiality->setChecked(bragg_.GetStillPartiality());
-    stillPartiality->setToolTip("Experimental, stills only: weight each reflection by a Gaussian "
-                                "excitation-error partiality exp(-d_ewald²/2σ²) instead of treating it as a "
-                                "full. Pairs with \"Partiality uncertainty\" in Scaling.");
     // Background trim: replace the r2..r3 ring mean with a symmetric trimmed mean (drop the lowest and
     // highest fraction of ring pixels), which removes the high-side bias that over-subtracts weak
     // high-angle reflections. Checkbox + fraction; 0 = plain mean.
@@ -468,7 +463,6 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() {
     trimRow->addWidget(bkgTrimFrac, 1);
     form->addRow("", gaussian);
     form->addRow("Radii r1/r2/r3", radii);
-    form->addRow("", stillPartiality);
     form->addRow("", trimRow);
     section->setContentLayout(form);
     section->setExpanded(false);   // folded on start (only geometry + unit cell start open)
@@ -477,12 +471,10 @@ QWidget *JFJochViewerSettingsDock::BuildBraggSection() {
         bragg_.Integrator(gaussian->isChecked() ? IntegratorMode::ProfileGaussian : IntegratorMode::BoxSum);
         bragg_.R1(static_cast<float>(r1->value())).R2(static_cast<float>(r2->value()))
               .R3(static_cast<float>(r3->value()));
-        bragg_.StillPartiality(stillPartiality->isChecked());
         bragg_.BackgroundTrimFraction(bkgTrim->isChecked() ? static_cast<float>(bkgTrimFrac->value()) : 0.0f);
         emit braggChanged(bragg_);
     };
     connect(gaussian, &QCheckBox::toggled, this, [emitBragg] { emitBragg(); });
-    connect(stillPartiality, &QCheckBox::toggled, this, [emitBragg] { emitBragg(); });
     connect(bkgTrim, &QCheckBox::toggled, this,
             [bkgTrim, bkgTrimFrac, emitBragg] { bkgTrimFrac->setEnabled(bkgTrim->isChecked()); emitBragg(); });
     connect(bkgTrimFrac, &NumberLineEdit::newValue, this, [emitBragg] { emitBragg(); });
@@ -512,20 +504,27 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() {
     modulation->setToolTip("Stills: fit a detector-plane modulation (flat-field) surface over where each "
                            "reflection lands, cross-validated so it no-ops when the systematic is absent. "
                            "For rotation, modulation is part of \"Correction surfaces\" above.");
+    auto *partRefine = new QCheckBox("Partiality post-refinement (stills)", this);
+    partRefine->setChecked(scaling_.GetStillsPartialityRefine());
+    partRefine->setToolTip("Stills: refine a per-crystal orientation tilt against the running merge and "
+                           "recompute each reflection's partiality (physical Ewald-proximity model), then "
+                           "re-scale/merge. On by default; uncheck for the simple model (each reflection a "
+                           "full, single-pass). No effect on rotation data.");
     auto *limitRes = new QCheckBox("High-resolution limit", this);
     limitRes->setChecked(scaling_.GetHighResolutionLimit_A().has_value());
     auto *highRes = new NumberLineEdit(0.3f, 5.0f, scaling_.GetHighResolutionLimit_A().value_or(2.0), 1, "Å", this);
     highRes->setEnabled(limitRes->isChecked());
 
     // Stills partiality-uncertainty merge term: adds a systematic sigma ~c*(1-partiality)*<I> on partials,
-    // so strong low-partiality partials are not over-trusted. Pairs with "Stills partiality" (Bragg); the
-    // library auto-gates it to strong/medium data. 0 = off; ~2.5 recommended.
+    // so strong low-partiality partials are not over-trusted. Relevant once partials exist (partiality
+    // post-refinement on, the default); the library auto-gates it to strong/medium data. 0 = off; ~2.5 rec.
     const double part_unc = scaling_.GetPartialityUncertaintyCoeff();
     auto *partUncertain = new QCheckBox("Partiality uncertainty", this);
     partUncertain->setChecked(part_unc > 0.0);
     partUncertain->setToolTip("Stills: add a systematic merge σ ~c·(1−partiality)·⟨I⟩ to partials so strong "
-                              "low-partiality partials are not over-trusted. Use with \"Stills partiality\"; "
-                              "auto-gated to strong/medium data. ~2.5 recommended; unchecked = off.");
+                              "low-partiality partials are not over-trusted. Relevant with partiality "
+                              "post-refinement (the default); auto-gated to strong/medium data. ~2.5 "
+                              "recommended; unchecked = off.");
     auto *partUncertainCoeff = new NumberLineEdit(0.1f, 10.0f, part_unc > 0.0 ? part_unc : 2.5, 1, "", this);
     partUncertainCoeff->setEnabled(partUncertain->isChecked());
 
@@ -533,6 +532,7 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() {
     form->addRow("", refineB);
     form->addRow("", corrections);
     form->addRow("", modulation);
+    form->addRow("", partRefine);
     // Compact, and aligned with the checkboxes above: the limit checkbox + value sit together in the
     // field column (not as a row label, which would indent it differently).
     auto *resRow = new QHBoxLayout();
@@ -551,6 +551,7 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() {
         scaling_.RefineB(refineB->isChecked());
         scaling_.CorrectionSurfaces(corrections->isChecked());
         scaling_.StillsModulation(modulation->isChecked());
+        scaling_.StillsPartialityRefine(partRefine->isChecked());
         scaling_.HighResolutionLimit_A(limitRes->isChecked()
             ? std::optional<double>(highRes->value()) : std::nullopt);
         scaling_.PartialityUncertaintyCoeff(partUncertain->isChecked() ? partUncertainCoeff->value() : 0.0);
@@ -560,6 +561,7 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() {
     connect(refineB, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); });
     connect(corrections, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); });
     connect(modulation, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); });
+    connect(partRefine, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); });
     connect(limitRes, &QCheckBox::toggled, this, [emitScaling, highRes](bool on) {
         highRes->setEnabled(on); emitScaling(); });
     connect(highRes, &NumberLineEdit::newValue, this, [emitScaling] { emitScaling(); });