From 1239c497313019ded59854a3572c1cb31fecad97 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 9 Aug 2026 21:08:29 +0200 Subject: [PATCH] Bragg integration: separate the three things a bandwidth used to switch Setting a bandwidth flipped three unrelated switches at once: it changed the profile's radial capture term, it moved the width measurement from the signal disk to the whole fit grid, and it silently overrode the background clip and trim, so --background-clip under --bandwidth was ignored - the two runs were bit-identical. The width measurement was the damaging one. The fit grid is an azimuthally averaged stack, so its second moment is sigma_r^2 + sigma_t^2 and the radial smear of a bandwidth leaked into the tangential model - a tangential width of 3.04 px against a 1.06 px truth, inflating the effective background pixel count where the weak signal is. The result was a step rather than a slope: on genuinely monochromatic data, declaring a 0.2% bandwidth cost ISa 28.4 -> 22.2. Measure the two widths separately, accumulated in each spot's own radial/tangential frame over the signal disk, from the signed profile cells - away from the peak a learned cell is background noise centred on zero, so the signed sum is unbiased, while clamping it at zero turns that noise into a pedestal the r^2 weight reads as width. The radial term is then the measured excess or the analytic floor, whichever is larger. With the two widths separated there is nothing left for the broadband switch to select, so it is gone - which is the proof the three were independent. The background clip and trim now come from the settings in every case; the tuned 3-sigma broadband default moves to the rugnux front end, which is the only place that knows whether the user gave a value. Monochromatic data: declaring a 0.2% bandwidth now costs ISa 28.4 -> 27.9 rather than 22.2, and forcing the old 3-sigma clip in the new build reproduces the good result, so none of the step came from the clip. On large-bandwidth data CC1/2 improves in 8 of 10 shells. Across 12 monochromatic crystals the space groups are unchanged and CC1/2 moves by at most 0.2 points. Co-Authored-By: Claude Opus 5 (1M context) --- common/BraggIntegrationSettings.h | 5 +- docs/CHANGELOG.md | 2 + docs/CPU_DATA_ANALYSIS.md | 12 +- .../BraggIntegrationEngine.cpp | 23 ++-- .../BraggIntegrationEngine.h | 9 +- .../BraggIntegrationEngineCPU.cpp | 97 +++++++------- .../BraggIntegrationEngineGPU.cu | 126 +++++++++++------- .../BraggIntegrationEngineGPU.h | 3 +- rugnux/rugnux_cli.cpp | 21 ++- 9 files changed, 174 insertions(+), 124 deletions(-) diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 4886611a..1a43be2f 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -44,8 +44,9 @@ class BraggIntegrationSettings { // High-side-only sigma clip: reject ring pixels above mean + this many sqrt(mean). Rejects the same // contamination as the trim - measurably better, in fact - without cutting the low side, so it does // not carry the trim's skew bias. Measured empty-aperture pedestal, counts: plain mean -0.03..-0.20, - // 10% symmetric trim +5.05..+6.34, 4 sigma clip +0.02..+0.54. Applied to monochromatic data - // (rugnux --background-clip); broadband (non-zero bandwidth) data always clip, at their tuned 3 sigma. + // 10% symmetric trim +5.05..+6.34, 4 sigma clip +0.02..+0.54. Whatever is set here is what the + // engine applies (rugnux --background-clip); the front end picks the default, and rugnux lowers it + // to 3 sigma for broadband (non-zero bandwidth) data, where longer spots leak further into the ring. float bkg_clip_nsigma = 4.0f; // Radial background curvature correction. The signal disk and the background annulus are // concentric, so for ANY background linear in position their means are equal - a plane fit buys diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index d5d4ee96..7f3ce011 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -13,6 +13,8 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Bragg prediction: How far the predictor walks the lattice is a setting (`bragg_integration_settings.max_hkl`, `--max-hkl`) instead of a fixed 100, derived per crystal offline and held fixed online. * Bragg integration: The local background ring is made robust with a **high-side sigma clip** (`--background-clip `, default 4) instead of the symmetric trimmed mean; the trim stays reachable with `--background-trim `. * Bragg integration: The **uncertainty of the background estimate** is now propagated into `sigma`, which both engines previously omitted. +* Bragg integration: Profile-fit `sigma` is **no longer inflated on weak reflections** - the fit weights take the signal estimate as it is instead of clamping it at zero, and the per-pixel variance floor is 0.01 counts instead of 1/12. +* Bragg integration: The spot width is now measured **separately along the radial and tangential directions**, and `--bandwidth` no longer changes the background estimator: `--background-clip` and `--background-trim` apply whatever the bandwidth is (the clip still defaults to 3 sigma when a bandwidth is set). * Bragg integration: New **radial background correction** (`--background-radial[=on|off|auto]`, default off); `auto` applies it per image where that image's ice score shows a smooth powder ring. * rugnux: **Ice-ring handling now runs only when the crystal is measured to have ice**, detected on two channels - the spot finder's peak-excluded radial profile (`--ice-min-score`, default 1.5) and found spots on the rings against ice-free flanks (`--ice-min-spot-ratio`, default 2.0). * rugnux: The merge-time **ice-ring mask has been removed**; ice reflections are still kept out of the scale fit and the space-group search, and kept in the final merge. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index ce04b78d..212c653c 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -582,7 +582,7 @@ $ $ with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates $\sigma$ by $\sqrt{1+n_S/n_B}$ — 1.109 with the shipped stencil — uniformly, on every reflection of every dataset. The same term is carried into the profile fit (§9.3), where it adds $(\sum wP/\sum P^2/v)^2\,\mathrm{var}(\hat b)$; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum. This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below. -**High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ on monochromatic data (`--background-clip`; $n=0$ disables) and $n=3$ on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. +**High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ (`--background-clip`; $n=0$ disables), lowered by `rugnux` to $n=3$ on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. That is only a default — the flag sets $n$ whatever the bandwidth is. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. The clip cuts only the high tail, which matters: the **symmetric** trimmed mean it replaced (drop the lowest and highest fraction $f$ of ring pixels, $f=0.10$; still reachable with `--background-trim`, which switches the clip off) is *not* a consistent estimator of the mean of a right-skewed Poisson sample. It sits $\approx0.1$ ct/px **below** the true mean at every level, and with $n_S\approx50$ signal pixels in the $r_1$ disk that under-estimate adds $\approx5$ counts to **every** partial — 0.2 % of the mean partial at 4 Å but $\approx13$ % of the mean and $\approx50$ % of the median partial in the outermost shell. Measured empty-aperture pedestal, in counts: plain mean $-0.03\ldots-0.20$, $10\,\%$ symmetric trim $+5.05\ldots+6.34$, $4\sigma$ clip $+0.02\ldots+0.54$. The same contamination is rejected either way — better, in fact: the trim collapses once contamination exceeds $\approx10\,\%$ of the ring (the same neighbour core shifts it by $+10.1$ ct/px). Note that removing a positive background bias *lowers* $\langle I/\sigma\rangle$ and *raises* edge $R_\text{meas}$, because both are inflated by information-free counts — so neither may be read as evidence against the change. The accuracy gain shows up instead in per-shell agreement with independent processing of the same images, with $CC_{1/2}$ neutral to slightly positive. @@ -601,18 +601,18 @@ The model is a function of **radius alone**, so it is applied only where that is A fixed signal disk captures a *width-dependent* fraction of each spot, which puts a multiplicative floor on the per-observation precision of strong reflections and weights weak reflections poorly. Profile fitting removes this by extracting each intensity against a fitted spot shape, without needing reference intensities. Per frame: 1. **Seed.** Box-sum every reflection (§9.2) to get a rough intensity and observed centroid, and select strong spots (significance $\ge 5$). -2. **Build the profile.** For `gaussian` (the default) the width is taken **per resolution shell** from the measured second moment of the strong spots (shell-dependent because spot size grows with resolution); the intrinsic spot is treated as round in the detector plane, the crystal's own anisotropy lying in the rocking direction, which a 2D detector-plane profile does not sample. For `empirical` the profile is instead the averaged, background-subtracted pixel grid of the shell's strong spots, accumulated on their **rounded predicted** positions. For `gaussian` only, the profile is then **rebuilt for each reflection**, centred on its **sub-pixel predicted position** (the noise-free geometric centre, not the observed centroid) and, where needed, **elongated only along the radial direction** (away from the beam centre) — because two effects stretch a spot radially but not tangentially: +2. **Build the profile.** For `gaussian` (the default) the width is taken **per resolution shell** from the measured second moments of the strong spots (shell-dependent because spot size grows with resolution). The moments are **anisotropic**: each strong spot's pixels are rotated into its *own* radial/tangential frame before being accumulated, giving $\sigma^2_r$ and $\sigma^2_t$ separately. Stacking the spots in the detector frame instead — they sit at every azimuth — averages the two directions away, leaving only $\sigma_r^2+\sigma_t^2$, so radial smearing is read back as a wider *tangential* spot. For `empirical` the profile is instead the averaged, background-subtracted pixel grid of the shell's strong spots, accumulated in the detector frame on their **rounded predicted** positions. For `gaussian` only, the profile is then **rebuilt for each reflection**, centred on its **sub-pixel predicted position** (the noise-free geometric centre, not the observed centroid) and, where needed, **elongated only along the radial direction** (away from the beam centre) — because two effects stretch a spot radially but not tangentially: - a finite energy **bandwidth** smears each spot by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$ ($R_\mathrm{px}$ = distance from the beam centre, large at high resolution), and - - sensor **parallax** — the depth over which a photon converts in a thick Si/CdTe sensor — adds a term $\propto\tan^2(2\theta)$ (material- and energy-dependent), plus, on the monochromatic path, a small fixed weak-spot capture term. + - sensor **parallax** — the depth over which a photon converts in a thick Si/CdTe sensor — adds a term $\propto\tan^2(2\theta)$ (material- and energy-dependent), plus a small fixed weak-spot capture term. - These combine as $\sigma^2_\mathrm{radial}=\sigma^2_\mathrm{intrinsic}+\sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)$ (tangential unchanged), on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add. The `empirical` profile keeps the fixed per-shell grid and gets none of this. + The two enter as a floor on the measured radial excess: $\sigma^2_\mathrm{radial}=\sigma^2_t+\max\!\left(\sigma^2_r-\sigma^2_t,\ \sigma_\mathrm{bw}^2+c_\mathrm{par}\tan^2(2\theta)\right)$, tangential unchanged at $\sigma^2_t$. The measured excess is what the signal disk can resolve; the analytic term takes over for a streak too long to be measured there. The Gaussian is built on a grid grown to hold the streak — capturing it without the tangential background an isotropic widening would add. The `empirical` profile keeps the fixed per-shell grid and gets none of this. 3. **Fit (Kabsch).** With profile $P$, background $B$ and the shell variance model, the intensity and its uncertainty are $ I = \frac{\sum P\,(c-B)/v}{\sum P^2/v},\qquad \sigma = \sqrt{\frac{1}{\sum P^2/v}},\qquad -v = B + \max(I,0)\,P, +v = \max\!\left(B + I\,P,\ \tfrac{1}{2}B\right), $ -where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the variance floors the background at $1/12$ (the integer-binning pixel-variance floor). The rotation/excitation partiality is carried exactly as in the box-sum path. +where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). The plug-in $I$ enters **as it is**: half-wave rectifying it, $v=B+\max(I,0)P$, lets $v$ — and with it the reported $1/\sum P^2/v$ — respond only to *upward* fluctuations of a noisy estimate, which adds $\approx0.4\,\sigma\sum P^3/(\sum P^2)^2$ to every $\sigma$ whatever the count rate. That offset is invisible on strong reflections and a large fractional inflation on weak ones; the $\tfrac12 B$ clamp keeps $v$ positive without reintroducing it. As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the background term is floored at $0.01$ ct/px — enough to keep $P^2/v$ finite when the ring mean reads exactly zero, which a ring of $n_B$ pixels cannot distinguish from any background below $\approx1/n_B$. The rotation/excitation partiality is carried exactly as in the box-sum path. The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`). diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index c1e3b629..7d239599 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -57,26 +57,25 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe G = 2 * R + 1; GG = G * G; - // A set bandwidth (broadband / stills) vs monochromatic (rotation) splits the treatment: the - // background sigma-clip and radial-elongation terms are path-dependent (see ProfileIntegrate2D). + // The X-ray bandwidth enters ONE place: it smears a reflection radially by bw_sigma * Rpx, which + // the per-reflection Gaussian carries as part of its radial variance. It is not a mode switch - + // the background estimator and the parallax/capture term below are the same whatever the beam is. bw_sigma = experiment.GetBandwidthFWHM().value_or(0.0f) / 2.3548f; - broadband = bw_sigma > 0.0; const double c_par = parallax_var_px2(det.GetSensorMaterial(), det.GetSensorThickness_um(), geom.GetWavelength_A(), geom.GetPixelSize_mm() * 1000.0); - c_radial = c_par + (broadband ? 0.0 : bragg_engine::C_CAPTURE); + c_radial = c_par + bragg_engine::C_CAPTURE; F_px = geom.GetDetectorDistance_mm() / std::max(1e-6f, geom.GetPixelSize_mm()); beam_x = geom.GetBeamX_pxl(); beam_y = geom.GetBeamY_pxl(); - use_ellipse = !empirical && (bw_sigma > 0.0 || c_radial > 0.0); + use_ellipse = !empirical; - // Robust background ring, one estimator or the other (see BraggIntegrationSettings). Broadband - // (non-zero bandwidth: pink-beam / DMM) data keep their tuned 3 sigma high-side clip whatever the - // settings say; monochromatic data - rotation AND stills, the discriminator is the beam, not the - // acquisition mode - take the clip multiplier from settings, and fall back to the symmetric trim - // only when the clip is switched off (rugnux --background-trim). - bkg_clip_nsigma = broadband ? 3.0f : settings.GetBackgroundClipNSigma(); - bkg_trim = (broadband || bkg_clip_nsigma > 0.0f) ? 0.0f : settings.GetBackgroundTrimFraction(); + // Robust background ring, one estimator or the other (see BraggIntegrationSettings): a high-side + // sigma-clip (rugnux --background-clip, the default) or, when the clip is switched off, a + // symmetric trimmed mean (rugnux --background-trim). The caller owns the choice - the engine no + // longer overrides it for broadband data. + bkg_clip_nsigma = settings.GetBackgroundClipNSigma(); + bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction(); // Radial-offset kernels for the background curvature correction. A stencil pixel at (dx, dy) // sits at radial offset dx*cos(phi) + dy*sin(phi) from the reflection, where phi is the diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 38bf8339..bd5caebf 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -48,7 +48,7 @@ namespace bragg_engine { constexpr int N_SHELL = 6; // resolution shells for per-shell profile learning constexpr double STRONG_I_OVER_SIGMA = 5.0; // strong-spot threshold that seeds the profile constexpr int MIN_STRONG_PER_SHELL = 30; // below this a shell falls back to the global profile -constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term (monochromatic only) +constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term (coefficient of tan^2(2theta)) // Lower bound on the background term of the Kabsch fit weights (v = max(bkg, floor) + signal). It // guards the background ESTIMATE, not the detector: the r2..r3 ring mean of a high-angle reflection // can come out exactly zero, and v = 0 makes the weights P^2/v diverge. A ring of n pixels cannot @@ -101,7 +101,6 @@ protected: float r3, r3_sq; int R, G, GG; // profile-grid half-size, edge (2R+1) and area (G*G) - bool broadband; // a set bandwidth (stills) vs monochromatic (rotation) double bw_sigma; // bandwidth sigma [dimensionless, * Rpx -> px] float bkg_clip_nsigma; // high-outlier background sigma-clip multiplier (0 = no clip) bool use_ellipse; // radially elongate the per-reflection Gaussian @@ -109,9 +108,9 @@ protected: double c_radial; // radial variance coefficient of tan^2(2theta): parallax + capture double F_px; // detector distance expressed in pixels float beam_x, beam_y; - // Effective symmetric trimmed-mean background fraction (BraggIntegrationSettings): the configured - // fraction for monochromatic (rotation) data, forced to 0 for broadband (stills, which keep their - // high-side sigma-clip). 0 = plain ring mean. Read by both the CPU and GPU engines. + // Effective symmetric trimmed-mean background fraction (BraggIntegrationSettings), used only when + // the high-side clip is switched off - the two are alternatives. 0 = plain ring mean. Read by both + // the CPU and GPU engines. float bkg_trim = 0.0f; // --- radial background curvature correction (BraggIntegrationSettings) --- diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index d7ce3f0a..3ee24108 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -262,13 +262,25 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (rough[i].ok) rough[i].shell = shell_of(predicted[i].d); // --- Learn the profile per shell (+ global) from the strong spots. --- + // Two things are learned. The empirical profile is the average grid in the DETECTOR frame, which is + // where it is applied. The width is a pair of second moments taken in each spot's OWN radial / + // tangential frame: a grid stacked in the detector frame is azimuthally averaged, so its is + // sigma_r^2 + sigma_t^2 with no way back, and a radially smeared spot reads as a wide TANGENTIAL + // one. Rotating each contribution into the spot's frame keeps the two apart. + struct Moments { double rad = 0.0, tan = 0.0, w = 0.0; }; + struct Sigma2 { double rad = 1.0, tan = 1.0; }; std::vector> shell_grid(N_SHELL, std::vector(GG, 0.0)); + std::vector shell_mom(N_SHELL); std::vector shell_n(N_SHELL, 0); std::vector global_grid(GG, 0.0); + Moments global_mom; int global_n = 0; for (size_t i = 0; i < npredicted; ++i) { const auto &rh = rough[i]; if (!rh.ok || !rh.strong || rh.I <= 0.0) continue; + const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; + const double Rpx = std::hypot(rx, ry); + const double ux = Rpx > 1e-6 ? rx / Rpx : 1.0, uy = Rpx > 1e-6 ? ry / Rpx : 0.0; for (int dy = -R; dy <= R; ++dy) for (int dx = -R; dx <= R; ++dx) { const int x = rh.cx + dx, y = rh.cy + dy; @@ -278,65 +290,54 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const double v = (static_cast(px) - rh.bkg) / rh.I; shell_grid[rh.shell][grid_idx(dx, dy)] += v; global_grid[grid_idx(dx, dy)] += v; + if (dx * dx + dy * dy >= r1_sq) continue; + const double rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + shell_mom[rh.shell].rad += v * rad * rad; + shell_mom[rh.shell].tan += v * tn * tn; + shell_mom[rh.shell].w += v; + global_mom.rad += v * rad * rad; + global_mom.tan += v * tn * tn; + global_mom.w += v; } ++shell_n[rh.shell]; ++global_n; } - // Isotropic width (2nd moment) of a learned grid: over the r1 disk (monochromatic) or the full - // grid (broadband); = 2 sigma^2 in 2D. The cells are signed: away from the peak a learned - // cell is pure background noise centred on zero, and clamping it at zero turns that noise into a - // positive pedestal spread over the whole domain, which the r^2 weight then reads as extra width. - auto measure_sigma2 = [&](const std::vector &grid) { - double m2 = 0.0, m2w = 0.0; - for (int dy = -R; dy <= R; ++dy) - for (int dx = -R; dx <= R; ++dx) { - if (!broadband && dx * dx + dy * dy >= r1_sq) continue; - const double g = grid[grid_idx(dx, dy)]; - m2 += g * (dx * dx + dy * dy); - m2w += g; - } - return m2w > 0.0 ? std::max(0.25, (m2 / m2w) / 2.0) : 1.0; + // Radial and tangential variances from the moments. The domain is the r1 disk, which is + // azimuthally symmetric and so adds no anisotropy of its own. The cells are SIGNED: away from the + // peak a learned cell is background noise centred on zero, and clamping it at zero turns that + // noise into a positive pedestal that the rad^2 / tan^2 weights read as extra width. + auto widths = [](const Moments &m) { + Sigma2 s; + if (m.w > 0.0) { + s.rad = std::max(0.25, m.rad / m.w); + s.tan = std::max(0.25, m.tan / m.w); + } + return s; }; - // Normalised profile (sum = 1): empirical average grid, or an isotropic Gaussian of the measured - // 2nd moment (only used by ProfileEmpirical; ProfileGaussian rebuilds per reflection in Pass B). - auto build_profile = [&](const std::vector &grid, int n) { + // Normalised empirical profile (sum = 1), the average grid over the strong spots of a shell. + // ProfileGaussian does not use it - it rebuilds a per-reflection Gaussian in Pass B. + auto build_profile = [&](const std::vector &grid) { std::vector P(GG, 0.0); - if (n <= 0) return P; double sum = 0.0; for (int k = 0; k < GG; ++k) { - const double g = std::max(0.0, grid[k]); - sum += g; - if (empirical) P[k] = g; + P[k] = std::max(0.0, grid[k]); + sum += P[k]; } - if (sum <= 0.0) return P; - if (empirical) { + if (sum > 0.0) for (double &p : P) p /= sum; - } else { - const double sigma2 = measure_sigma2(grid); - double gsum = 0.0; - for (int dy = -R; dy <= R; ++dy) - for (int dx = -R; dx <= R; ++dx) { - const double g = std::exp(-(dx * dx + dy * dy) / (2.0 * sigma2)); - P[grid_idx(dx, dy)] = g; - gsum += g; - } - for (double &p : P) p /= gsum; - } return P; }; - const std::vector global_P = build_profile(global_grid, global_n); - const double global_sigma2 = global_n > 0 ? measure_sigma2(global_grid) : 1.0; - std::vector> shell_P(N_SHELL); - std::vector shell_sigma2(N_SHELL, global_sigma2); + const std::vector global_P = empirical && global_n > 0 ? build_profile(global_grid) + : std::vector(GG, 0.0); + const Sigma2 global_sigma2 = widths(global_mom); + std::vector> shell_P(N_SHELL, global_P); + std::vector shell_sigma2(N_SHELL, global_sigma2); for (int s = 0; s < N_SHELL; ++s) { - if (shell_n[s] >= MIN_STRONG_PER_SHELL) { - shell_P[s] = build_profile(shell_grid[s], shell_n[s]); - shell_sigma2[s] = measure_sigma2(shell_grid[s]); - } else { - shell_P[s] = global_P; - } + if (shell_n[s] < MIN_STRONG_PER_SHELL) continue; + if (empirical) shell_P[s] = build_profile(shell_grid[s]); + shell_sigma2[s] = widths(shell_mom[s]); } // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). --- @@ -352,12 +353,16 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; const double Rpx = std::hypot(rx, ry); const double tan2t = Rpx / F_px; - const double s2t = shell_sigma2[sh]; + const double s2t = shell_sigma2[sh].tan; double s2r = s2t, ux = 1.0, uy = 0.0; bool elong = false; if (use_ellipse) { + // Radial excess over the tangential width: measured where the peak is resolved inside + // the r1 disk, with the analytic bandwidth + parallax/capture term as the floor. The + // analytic term is what carries a streak the disk is too small to measure. const double sbw = bw_sigma * Rpx; - const double radial_extra = sbw * sbw + c_radial * tan2t * tan2t; + const double radial_extra = std::max(shell_sigma2[sh].rad - s2t, + sbw * sbw + c_radial * tan2t * tan2t); if (Rpx > 1e-6 && radial_extra > 0.25) { ux = rx / Rpx; uy = ry / Rpx; s2r = s2t + radial_extra; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index c132c6b3..63530fa1 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -19,7 +19,6 @@ struct BraggGpuParams { int R, G, GG; float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip) int empirical; // ProfileEmpirical vs ProfileGaussian - int broadband; int use_ellipse; float bw_sigma; float c_radial; @@ -33,6 +32,9 @@ __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MA // idea 1: max annulus pixels held in shared memory for the per-block trimmed-mean sort (else flat mean). constexpr int BKG_TRIM_MAX = 512; +// Learned second moments, (sum v*rad^2, sum v*tan^2, sum v) per shell plus one global slot at N_SHELL. +constexpr int MOM_STRIDE = 3; + // --- Mark the r2 signal disk of every predicted reflection (race-free: all writes are 1). --- __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, BraggGpuParams p, int n) { const int i = blockIdx.x; @@ -276,12 +278,13 @@ __device__ inline int compute_shell(float d, const unsigned long long *invd2mm) // --- Zero the profile accumulators and seed the inv-d^2 range, in one launch (replaces a handful // of small cudaMemsetAsync calls, which matter when kernel-launch latency is high). --- -__global__ void reset(float *shell_grid, float *global_grid, int *shell_n, int *global_n, +__global__ void reset(float *shell_grid, float *global_grid, float *mom, int *shell_n, int *global_n, unsigned long long *invd2mm, int GG) { for (int k = blockIdx.x * blockDim.x + threadIdx.x; k < N_SHELL * GG; k += blockDim.x * gridDim.x) shell_grid[k] = 0.0f; for (int k = blockIdx.x * blockDim.x + threadIdx.x; k < GG; k += blockDim.x * gridDim.x) global_grid[k] = 0.0f; + if (blockIdx.x == 0 && threadIdx.x < MOM_STRIDE * (N_SHELL + 1)) mom[threadIdx.x] = 0.0f; if (blockIdx.x == 0 && threadIdx.x < N_SHELL) shell_n[threadIdx.x] = 0; if (blockIdx.x == 0 && threadIdx.x == 0) { *global_n = 0; @@ -291,11 +294,16 @@ __global__ void reset(float *shell_grid, float *global_grid, int *shell_n, int * } // --- Learn the profile: each strong spot adds its bkg-subtracted, I-normalised grid to its shell -// (and the global grid). One block per reflection. --- -__global__ void learn_profile(const int32_t *img, const int *cx_a, const int *cy_a, const float *dd, +// (and the global grid), and its second moments to the same shell's moment slot. The grid stays +// in the DETECTOR frame (that is where the empirical profile is applied); the moments are taken in +// the spot's own radial/tangential frame, because a detector-frame stack is azimuthally averaged +// and cannot tell a radially smeared spot from a tangentially wide one (see the CPU engine). +// One block per reflection. --- +__global__ void learn_profile(const int32_t *img, const float *px_x, const float *px_y, + const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_a, const float *bkg_a, const uint8_t *ok_a, const uint8_t *strong_a, - float *shell_grid, float *global_grid, int *shell_n, int *global_n, + float *shell_grid, float *global_grid, float *mom, int *shell_n, int *global_n, BraggGpuParams p, int n) { const int i = blockIdx.x; if (i >= n || !ok_a[i] || !strong_a[i]) return; @@ -303,57 +311,75 @@ __global__ void learn_profile(const int32_t *img, const int *cx_a, const int *cy if (!(I > 0.0f)) return; const int cx = cx_a[i], cy = cy_a[i], sh = compute_shell(dd[i], invd2mm); const float bkg = bkg_a[i]; + const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; + const float Rpx = sqrtf(rx * rx + ry * ry); + const float ux = Rpx > 1e-6f ? rx / Rpx : 1.0f, uy = Rpx > 1e-6f ? ry / Rpx : 0.0f; float *sg = shell_grid + (size_t) sh * p.GG; + __shared__ float s_rad, s_tan, s_w; + if (threadIdx.x == 0) { s_rad = 0.0f; s_tan = 0.0f; s_w = 0.0f; } + __syncthreads(); + float l_rad = 0.0f, l_tan = 0.0f, l_w = 0.0f; for (int k = threadIdx.x; k < p.GG; k += blockDim.x) { - const int x = cx + (k % p.G - p.R), y = cy + (k / p.G - p.R); + const int dx = k % p.G - p.R, dy = k / p.G - p.R; + const int x = cx + dx, y = cy + dy; if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; const float v = ((float) px - bkg) / I; atomicAdd(&sg[k], v); atomicAdd(&global_grid[k], v); + if ((float) (dx * dx + dy * dy) >= p.r1_sq) continue; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + l_rad += v * rad * rad; l_tan += v * tn * tn; l_w += v; + } + // Reduce in shared memory first: every block adds to the same handful of global words, so one + // atomic per block instead of one per thread. + atomicAdd(&s_rad, l_rad); atomicAdd(&s_tan, l_tan); atomicAdd(&s_w, l_w); + __syncthreads(); + if (threadIdx.x == 0) { + atomicAdd(&mom[MOM_STRIDE * sh + 0], s_rad); + atomicAdd(&mom[MOM_STRIDE * sh + 1], s_tan); + atomicAdd(&mom[MOM_STRIDE * sh + 2], s_w); + atomicAdd(&mom[MOM_STRIDE * N_SHELL + 0], s_rad); + atomicAdd(&mom[MOM_STRIDE * N_SHELL + 1], s_tan); + atomicAdd(&mom[MOM_STRIDE * N_SHELL + 2], s_w); + atomicAdd(&shell_n[sh], 1); + atomicAdd(global_n, 1); } - if (threadIdx.x == 0) { atomicAdd(&shell_n[sh], 1); atomicAdd(global_n, 1); } } -// --- Reduce each learned grid to its 2nd-moment width (and, for empirical, a normalised profile). -// One block per grid: blocks [0,N_SHELL) are the shells, block N_SHELL is the global grid. --- -__global__ void build_profiles(const float *shell_grid, const float *global_grid, +// --- Turn the learned moments into radial/tangential variances and, for empirical, normalise the +// learned grid into a profile. One block per grid: blocks [0,N_SHELL) are the shells, block +// N_SHELL is the global one. --- +__global__ void build_profiles(const float *shell_grid, const float *global_grid, const float *mom, const int *shell_n, const int *global_n, - float *shell_P, float *global_P, float *shell_sigma2, float *global_sigma2, + float *shell_P, float *global_P, float *sigma2_r, float *sigma2_t, BraggGpuParams p) { const int b = blockIdx.x; - const float *grid; int nstrong; float *P; float *sig2; - if (b < N_SHELL) { grid = shell_grid + (size_t) b * p.GG; nstrong = shell_n[b]; P = shell_P + (size_t) b * p.GG; sig2 = &shell_sigma2[b]; } - else { grid = global_grid; nstrong = *global_n; P = global_P; sig2 = global_sigma2; } + const float *grid; int nstrong; float *P; + if (b < N_SHELL) { grid = shell_grid + (size_t) b * p.GG; nstrong = shell_n[b]; P = shell_P + (size_t) b * p.GG; } + else { grid = global_grid; nstrong = *global_n; P = global_P; } - __shared__ float s_m2, s_m2w, s_sum; - if (threadIdx.x == 0) { s_m2 = 0.0f; s_m2w = 0.0f; s_sum = 0.0f; } + if (threadIdx.x == 0) { + const float w = mom[MOM_STRIDE * b + 2]; + sigma2_r[b] = w > 0.0f ? fmaxf(0.25f, mom[MOM_STRIDE * b + 0] / w) : 1.0f; + sigma2_t[b] = w > 0.0f ? fmaxf(0.25f, mom[MOM_STRIDE * b + 1] / w) : 1.0f; + } + if (!p.empirical) return; + + __shared__ float s_sum; + if (threadIdx.x == 0) s_sum = 0.0f; __syncthreads(); - - float l_m2 = 0.0f, l_m2w = 0.0f, l_sum = 0.0f; + float l_sum = 0.0f; for (int k = threadIdx.x; k < p.GG; k += blockDim.x) { - const int dx = k % p.G - p.R, dy = k / p.G - p.R; - // The moment reads the SIGNED cell (see the CPU engine): clamping pure-background noise at - // zero leaves a positive pedestal that the r^2 weight reads as extra width. The empirical - // profile still takes the clamped grid - a profile has to be non-negative. - const float raw = grid[k], g = fmaxf(0.0f, raw); - const int r2i = dx * dx + dy * dy; - if (p.broadband || (float) r2i < p.r1_sq) { l_m2 += raw * (float) r2i; l_m2w += raw; } + const float g = fmaxf(0.0f, grid[k]); // a profile has to be non-negative + P[k] = g; l_sum += g; - if (p.empirical) P[k] = g; // pre-store clamped grid for in-place normalisation below } - atomicAdd(&s_m2, l_m2); atomicAdd(&s_m2w, l_m2w); atomicAdd(&s_sum, l_sum); + atomicAdd(&s_sum, l_sum); __syncthreads(); - - if (threadIdx.x == 0) - *sig2 = (nstrong > 0 && s_m2w > 0.0f) ? fmaxf(0.25f, (s_m2 / s_m2w) / 2.0f) : 1.0f; - - if (p.empirical) { - const float sum = s_sum; - const bool normalise = nstrong > 0 && sum > 0.0f; - for (int k = threadIdx.x; k < p.GG; k += blockDim.x) P[k] = normalise ? P[k] / sum : 0.0f; - } + const bool normalise = nstrong > 0 && s_sum > 0.0f; + for (int k = threadIdx.x; k < p.GG; k += blockDim.x) P[k] = normalise ? P[k] / s_sum : 0.0f; } // --- Radial background curvature correction: one thread per reflection, no pixel reads. @@ -388,7 +414,7 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const float *I_seed, const float *sigma_seed, const float *bkg_a, const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, const float *shell_P, const float *global_P, - const float *shell_sigma2, const float *global_sigma2, const int *shell_n, + const float *sigma2_r, const float *sigma2_t, const int *shell_n, float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n) { const int i = blockIdx.x; if (i >= n) return; @@ -410,15 +436,18 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, for (int k = threadIdx.x; k < p.GG; k += blockDim.x) Pbuf[k] = Psrc[k]; __syncthreads(); } else { - const float s2t = use_shell ? shell_sigma2[sh] : *global_sigma2; + const int si = use_shell ? sh : N_SHELL; // N_SHELL is the global slot + const float s2t = sigma2_t[si]; const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; const float Rpx = sqrtf(rx * rx + ry * ry); const float tan2t = Rpx / p.F_px; float s2r = s2t, ux = 1.0f, uy = 0.0f; bool elong = false; if (p.use_ellipse) { + // Measured radial excess over the tangential width, floored by the analytic bandwidth + + // parallax/capture term (see the CPU engine). const float sbw = p.bw_sigma * Rpx; - const float radial_extra = sbw * sbw + p.c_radial * tan2t * tan2t; + const float radial_extra = fmaxf(sigma2_r[si] - s2t, sbw * sbw + p.c_radial * tan2t * tan2t); if (Rpx > 1e-6f && radial_extra > 0.25f) { ux = rx / Rpx; uy = ry / Rpx; s2r = s2t + radial_extra; elong = true; } } const int Rf = elong ? min(3 * p.R, (int) ceilf(p.r2 + 2.0f * sqrtf(s2r))) : p.R; @@ -498,8 +527,9 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment d_global_grid(GG), d_shell_P(static_cast(bragg_engine::N_SHELL) * GG), d_global_P(GG), - d_shell_sigma2(bragg_engine::N_SHELL), - d_global_sigma2(1), + d_mom(MOM_STRIDE * (bragg_engine::N_SHELL + 1)), + d_sigma2_r(bragg_engine::N_SHELL + 1), + d_sigma2_t(bragg_engine::N_SHELL + 1), d_shell_n(bragg_engine::N_SHELL), d_global_n(1), d_invd2(2) { @@ -598,7 +628,6 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .R = R, .G = G, .GG = GG, .bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f, .empirical = empirical ? 1 : 0, - .broadband = broadband ? 1 : 0, .use_ellipse = use_ellipse ? 1 : 0, .bw_sigma = static_cast(bw_sigma), .c_radial = static_cast(c_radial), .F_px = static_cast(F_px), @@ -616,7 +645,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(float) * n_rad, *stream)); cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream)); } - reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_invd2, GG); + reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG); mark_mask<<>>(d_px_x, d_px_y, d_mask, p, n); boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, p, n, d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y, @@ -633,14 +662,15 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (mode != IntegratorMode::BoxSum) { // Pass B: learn (shell computed inline) -> build -> fit. - learn_profile<<>>(img, d_cx, d_cy, d_d, d_invd2, d_I, d_bkg, d_ok, d_strong, - d_shell_grid, d_global_grid, d_shell_n, d_global_n, p, n); + learn_profile<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, + d_I, d_bkg, d_ok, d_strong, + d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, p, n); build_profiles<<>>( - d_shell_grid, d_global_grid, d_shell_n, d_global_n, - d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, p); + d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, + d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, p); fit<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_shell_P, d_global_P, - d_shell_sigma2, d_global_sigma2, d_shell_n, + d_sigma2_r, d_sigma2_t, d_shell_n, d_I, d_sigma, d_var_bkg, d_ok, p, n); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index a9d0e4eb..e698c941 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -45,7 +45,8 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaDevicePtr d_mask; // per-pixel r2-disk reflection mask CudaDevicePtr d_shell_grid, d_global_grid; // learned profile accumulators (N_SHELL*GG, GG) CudaDevicePtr d_shell_P, d_global_P; // normalised profiles (empirical mode) - CudaDevicePtr d_shell_sigma2, d_global_sigma2; + CudaDevicePtr d_mom; // learned 2nd moments, 3 per shell + global + CudaDevicePtr d_sigma2_r, d_sigma2_t; // radial/tangential widths, N_SHELL + global CudaDevicePtr d_shell_n, d_global_n; CudaDevicePtr d_invd2; // [min,max] inv-d^2 as monotonic bit patterns diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e7a39f82..71aee9ab 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -47,6 +47,10 @@ // in its favour and no measurable cost in wall clock. It is also what jfjoch_viewer already sends. constexpr int64_t RUGNUX_MAX_SPOT_COUNT = 1000; +// Default background-ring high-side clip for broadband (non-zero bandwidth) data, tighter than the +// monochromatic default because a pink-beam reflection is long and its wings reach into the ring. +constexpr float BROADBAND_BACKGROUND_CLIP_NSIGMA = 3.0f; + // Default rot3d per-frame scale-G smoothing range (XDS DELPHI-like), in degrees of rotation. constexpr double SMOOTH_G_DEFAULT_DEG = 5.0; @@ -156,9 +160,9 @@ void print_usage() { std::cout << " --integration-radius Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl; std::cout << " --integration-high-resolution High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl; std::cout << " --max-hkl Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl; - std::cout << " --background-clip Monochromatic (rotation + still): high-side clip of the background ring at mean + n*sqrt(mean) (default 4; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Broadband data always clip, at 3 sigma; ignored by --integrator boxsum" << std::endl; + std::cout << " --background-clip High-side clip of the background ring at mean + n*sqrt(mean) (default 4, or 3 when --bandwidth is set; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl; std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default off). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl; - std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean" << std::endl; + std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean. Applies whatever --bandwidth is set to" << std::endl; std::cout << " --integrator 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 Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl; @@ -1847,8 +1851,8 @@ static int RunRugnux(int argc, char **argv) { BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); bis.BackgroundTrimFraction(static_cast(*background_trim_arg)); experiment.ImportBraggIntegrationSettings(bis); - logger.Info("Background ring: symmetric trimmed mean at {:.2f} instead of the default high-side clip " - "(monochromatic data; broadband always clips)", *background_trim_arg); + logger.Info("Background ring: symmetric trimmed mean at {:.2f} instead of the default high-side clip", + *background_trim_arg); } if (background_clip_arg) { @@ -1856,6 +1860,15 @@ static int RunRugnux(int argc, char **argv) { bis.BackgroundClipNSigma(static_cast(*background_clip_arg)); experiment.ImportBraggIntegrationSettings(bis); logger.Info("Background ring: high-side clip at {:.1f} sigma", *background_clip_arg); + } else if (!background_trim_arg && experiment.GetBandwidthFWHM().value_or(0.0f) > 0.0f) { + // Broadband (pink-beam / DMM) reflections are long and their wings reach into the background + // ring, so the default clip is tighter there. It is only a default - --background-clip and + // --background-trim override it as they do for any other data. + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.BackgroundClipNSigma(BROADBAND_BACKGROUND_CLIP_NSIGMA); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Background ring: high-side clip at {:.1f} sigma (broadband default)", + BROADBAND_BACKGROUND_CLIP_NSIGMA); } if (background_radial_given) {