diff --git a/broker/gen/model/Spot_finding_settings.cpp b/broker/gen/model/Spot_finding_settings.cpp
index 1d2001f6..d85be081 100644
--- a/broker/gen/model/Spot_finding_settings.cpp
+++ b/broker/gen/model/Spot_finding_settings.cpp
@@ -32,7 +32,7 @@ Spot_finding_settings::Spot_finding_settings()
m_Low_resolution_limit = 0.0f;
m_High_resolution_limit_for_spot_count_low_res = 0.0f;
m_Quick_integration = false;
- m_Ice_ring_width_q_recipA = 0.02f;
+ m_Ice_ring_width_q_recipA = 0.03f;
m_High_res_gap_Q_recipA = 1.5f;
m_High_res_gap_Q_recipAIsSet = false;
m_Adaptive_threshold = false;
diff --git a/broker/gen/model/Spot_finding_settings.h b/broker/gen/model/Spot_finding_settings.h
index 7de307f9..63838b1b 100644
--- a/broker/gen/model/Spot_finding_settings.h
+++ b/broker/gen/model/Spot_finding_settings.h
@@ -110,7 +110,7 @@ public:
bool isQuickIntegration() const;
void setQuickIntegration(bool const value);
///
- /// Width of ice ring in q-space in reciprocal space
+ /// Half-width of the ice ring band in q (1/A). Matches the offline default in image_analysis/spot_finding/SpotFindingSettings.h, which was set from a measured ring FWHM of ~0.06; the two must agree or the same data gets a narrower ice band online.
///
float getIceRingWidthQRecipA() const;
void setIceRingWidthQRecipA(float const value);
diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml
index 011a3eac..eb75e3ea 100644
--- a/broker/jfjoch_api.yaml
+++ b/broker/jfjoch_api.yaml
@@ -1098,8 +1098,10 @@ components:
format: float
minimum: 0.0
maximum: 1.0
- default: 0.02
- description: Width of ice ring in q-space in reciprocal space
+ default: 0.03
+ description: Half-width of the ice ring band in q (1/A). Matches the offline default in
+ image_analysis/spot_finding/SpotFindingSettings.h, which was set from a measured ring
+ FWHM of ~0.06; the two must agree or the same data gets a narrower ice band online.
high_res_gap_Q_recipA:
type: number
format: float
diff --git a/common/AzimuthalIntegrationProfile.cpp b/common/AzimuthalIntegrationProfile.cpp
index 3b271293..412eea01 100644
--- a/common/AzimuthalIntegrationProfile.cpp
+++ b/common/AzimuthalIntegrationProfile.cpp
@@ -180,13 +180,30 @@ float AzimuthalIntegrationProfile::GetBkgEstimate(const AzimuthalIntegrationSett
float AzimuthalIntegrationProfile::GetIceRingScore(const AzimuthalIntegrationSettings &settings,
float half_width_q) const {
+ return IceRingScore(GetResult1D(), q_bins, settings, half_width_q);
+}
+
+float AzimuthalIntegrationProfile::IceRingScore(const std::vector &profile, int32_t q_bins,
+ const AzimuthalIntegrationSettings &settings,
+ float half_width_q) {
// Strongest hexagonal-ice ring's intensity relative to the background *under* it (1 = no ice). The
// background is a smooth whole-profile estimate: a running median of the NON-ice bins, interpolated to
// each ring position - not a couple of adjacent shoulder bins (the azint binning is coarser than the
// ring width, so a local shoulder is only ~1 bin and a narrow ratio is noisy and can double-count the
// ring's own edge). Clean profiles then sit at ~1 at every ring; ice makes the ring bin stand out.
constexpr float two_pi = 6.283185307f;
- const std::vector prof = GetResult1D();
+ // Average over the azimuthal bins, if there are any; a profile that is already 1-D passes through.
+ std::vector prof(std::max(q_bins, 0), 0.0f);
+ std::vector nbin(prof.size(), 0);
+ for (size_t i = 0; i < profile.size() && !prof.empty(); i++) {
+ const size_t q_bin = i % prof.size();
+ if (std::isfinite(profile[i])) {
+ prof[q_bin] += profile[i];
+ nbin[q_bin]++;
+ }
+ }
+ for (size_t i = 0; i < prof.size(); i++)
+ prof[i] = nbin[i] ? prof[i] / static_cast(nbin[i]) : NAN;
const int nq = static_cast(prof.size());
const float low_q = settings.GetLowQ_recipA();
const float dq = settings.GetQSpacing_recipA();
diff --git a/common/AzimuthalIntegrationProfile.h b/common/AzimuthalIntegrationProfile.h
index f7d998f3..7eb6fcd2 100644
--- a/common/AzimuthalIntegrationProfile.h
+++ b/common/AzimuthalIntegrationProfile.h
@@ -47,6 +47,11 @@ public:
// Single per-image ice indicator: the strongest hexagonal-ice ring's intensity relative to the smooth
// radial background interpolated under it (1 = no ice, >1 = ice). See the .cpp for the background fit.
float GetIceRingScore(const AzimuthalIntegrationSettings& settings, float half_width_q) const;
+ // The same score for a radial profile that is not this object's: the adaptive spot finder computes a
+ // peak-excluded per-ring background in these same bins, which is a much cleaner input (see the .cpp).
+ // profile is q_bins long, or q_bins x azimuthal bins, in which case it is averaged over azimuth first.
+ static float IceRingScore(const std::vector& profile, int32_t q_bins,
+ const AzimuthalIntegrationSettings& settings, float half_width_q);
MultiLinePlot GetPlot(bool force_1d = false, PlotAzintUnit plot_unit = PlotAzintUnit::Q_recipA) const;
AzimuthalIntegrationProfile& operator+=(const AzimuthalIntegrationProfile& profile); // Not thread safe
};
diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h
index 940bd37c..1c1787ac 100644
--- a/common/JFJochMessages.h
+++ b/common/JFJochMessages.h
@@ -108,6 +108,12 @@ struct DataMessage {
std::vector spots;
std::optional spot_count;
std::optional spot_count_ice_rings;
+ // Spots in the ICE-FREE control flanks either side of the hexagonal rings, rescaled to the ring
+ // bands' own q width - so spot_count_ice_rings / spot_count_ice_control is 1 when spots are spread
+ // evenly and > 1 when they pile up on the rings. Pooled over the run (a ratio of two per-image
+ // counts this small is meaningless on its own image), it is the second ice channel: TEXTURED ice
+ // arrives as discrete spots and leaves the radial profile - and so ice_ring_score - flat.
+ std::optional spot_count_ice_control;
std::optional spot_count_low_res;
std::vector spot_plot_count;
@@ -366,6 +372,7 @@ struct EndMessage {
std::vector data_collection_efficiency;
std::vector spot_count;
std::vector spot_count_ice_ring;
+ std::vector spot_count_ice_control;
std::vector spot_count_low_res;
std::vector spot_count_indexed;
std::vector image_indexed;
@@ -386,7 +393,12 @@ struct EndMessage {
std::vector image_scale_factor;
std::vector image_scale_cc;
std::vector image_scale_mosaicity;
+ // Per-image ice strength. Note the name carries no v_ prefix, unlike v_bkg_estimate above - it
+ // was shipped that way and the CBOR key is part of the stream format.
std::vector ice_ring_score;
+ // Run mean of the above, the single "how icy was this dataset" number (1 = no ice). The
+ // bkg_estimate scalar's counterpart; written to /entry/MX/iceRingScoreMean.
+ std::optional ice_ring_score_mean;
};
struct MetadataMessage {
diff --git a/common/JFJochReceiverPlots.cpp b/common/JFJochReceiverPlots.cpp
index 7802dadc..17116054 100644
--- a/common/JFJochReceiverPlots.cpp
+++ b/common/JFJochReceiverPlots.cpp
@@ -66,6 +66,7 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A
spot_count_low_res.Clear(r);
spot_count_indexed.Clear(r);
spot_count_ice.Clear(r);
+ spot_count_ice_control.Clear(r);
indexing_solution.Clear(r);
indexing_uc_a.Clear(r);
@@ -133,6 +134,7 @@ void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegration
spot_count_low_res.AddElement(msg.number, msg.spot_count_low_res);
spot_count_indexed.AddElement(msg.number, msg.spot_count_indexed);
spot_count_ice.AddElement(msg.number, msg.spot_count_ice_rings);
+ spot_count_ice_control.AddElement(msg.number, msg.spot_count_ice_control);
error_pixels.AddElement(msg.number, msg.error_pixel_count);
saturated_pixels.AddElement(msg.number, msg.saturated_pixel_count);
pixel_sum.AddElement(msg.number, msg.pixel_sum);
@@ -470,6 +472,33 @@ std::optional JFJochReceiverPlots::GetBkgEstimate() const {
return {};
}
+std::optional JFJochReceiverPlots::GetIceRingSpotRatio() const {
+ // A ratio of MEANS, not a mean of ratios: one image holds a handful of control spots, so a
+ // per-image ratio is dominated by its own denominator. Pooling over the run is the measurement.
+ const float ring = spot_count_ice.Mean();
+ const float control = spot_count_ice_control.Mean();
+ if (!std::isfinite(ring) || !std::isfinite(control))
+ return std::nullopt;
+ // An empty control with spots on the rings is the STRONGEST evidence of ice there is, not the
+ // absence of it - a crystal whose found spots are all ice leaves nothing in the flanks. Report a
+ // large finite ratio rather than dividing by zero.
+ if (!(control > 0.0f))
+ return ring > 0.0f ? 1.0e3f : std::optional{};
+ return ring / control;
+}
+
+std::optional JFJochReceiverPlots::GetIceRingScore() const {
+ auto tmp = ice_ring_score.Mean();
+ if (std::isfinite(tmp))
+ return tmp;
+ else
+ return {};
+}
+
+std::vector JFJochReceiverPlots::GetIceRingScoreArray() const {
+ return ice_ring_score.ExportArray();
+}
+
MeanProcessingTime JFJochReceiverPlots::GetMeanProcessingTime() const {
MeanProcessingTime ret{};
ret.compression = compression_time.Mean();
diff --git a/common/JFJochReceiverPlots.h b/common/JFJochReceiverPlots.h
index 4af96d24..23d1fc23 100644
--- a/common/JFJochReceiverPlots.h
+++ b/common/JFJochReceiverPlots.h
@@ -47,6 +47,7 @@ class JFJochReceiverPlots {
StatusVector spot_count_low_res;
StatusVector spot_count_indexed;
StatusVector spot_count_ice;
+ StatusVector spot_count_ice_control;
StatusVector indexing_solution;
StatusVector indexing_lattice_count;
@@ -121,6 +122,11 @@ public:
std::optional GetIndexingRate() const;
std::optional GetBkgEstimate() const;
+ std::optional GetIceRingScore() const;
+ // Pooled over the run: spots on the hexagonal rings over the same q width of ice-free control
+ // flanks. 1 = spots spread evenly, > 1 = they pile up on the rings (textured ice).
+ [[nodiscard]] std::optional GetIceRingSpotRatio() const;
+ std::vector GetIceRingScoreArray() const;
std::vector GetAzIntProfile() const;
MultiLinePlot GetAzIntProfilePlot(bool force_1d = false, PlotAzintUnit azint_unit = PlotAzintUnit::Q_recipA) const;
diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp
index c1562f61..9717ad84 100644
--- a/common/ScalingSettings.cpp
+++ b/common/ScalingSettings.cpp
@@ -173,6 +173,28 @@ bool ScalingSettings::GetIceRingMergeMask() const {
return ice_ring_merge_mask;
}
+float ScalingSettings::GetIceMinScore() const {
+ return ice_min_score;
+}
+
+float ScalingSettings::GetIceMinSpotRatio() const {
+ return ice_min_spot_ratio;
+}
+
+ScalingSettings &ScalingSettings::IceMinSpotRatio(float input) {
+ if (input < 0)
+ throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Ice spot-ratio gate must be non-negative");
+ ice_min_spot_ratio = input;
+ return *this;
+}
+
+ScalingSettings &ScalingSettings::IceMinScore(float input) {
+ if (input < 0)
+ throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Ice score gate must be non-negative");
+ ice_min_score = input;
+ return *this;
+}
+
ScalingSettings &ScalingSettings::SmoothGDegrees(double input) {
if (input < 0)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Smooth-G range must be non-negative");
diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h
index a7a0df32..15ab255f 100644
--- a/common/ScalingSettings.h
+++ b/common/ScalingSettings.h
@@ -68,6 +68,23 @@ class ScalingSettings {
// SCALE fit - so this can be turned off on its own to de-confound a merge experiment without changing
// how the data were indexed and scaled. On by default; only consulted when detect_ice_rings is on.
bool ice_ring_merge_mask = true;
+ // Minimum measured ice strength (iceRingScore, 1 = no ice) before any ice-ring handling is applied
+ // at all. The eleven fixed hexagonal bands cover 16-26% of the unique reflections at typical
+ // resolutions REGARDLESS of whether the crystal has ice, so flagging unconditionally taxes clean
+ // data for nothing - and the merge-time ring mask has been observed to fire on crystals with no
+ // measurable ice. 0 disables the gate (always handle ice, the previous behaviour).
+ // 1.5 is measured, not guessed: over 37 rotation crystals the score lands at 1.00-1.22 on the
+ // thirty with no ice, 1.28-1.44 on four borderline ones whose ice-ring positions show no
+ // azimuthally smooth elevation, and 2.08-2.37 on the three with confirmed ice - and a decoy null
+ // (the identical statistic at ring positions where hexagonal ice cannot be) never exceeded 1.29.
+ float ice_min_score = 1.5f;
+ // The same gate on the SECOND ice channel: spots found on the hexagonal rings over the same q width
+ // of ice-free flanks beside them (1 = spots spread evenly). This is what catches ice in large
+ // crystallites, which diffracts as discrete spots and leaves the radial profile - and so
+ // ice_min_score - flat. Also measured, not guessed: over 36 rotation crystals thirty read
+ // 0.65-1.37 and a clean control 1.04, then 1.63/1.78 and a gap to 2.18-14.6 on the five whose
+ // spots really do pile up on the rings. 0 disables this channel.
+ float ice_min_spot_ratio = 2.0f;
// Smooth the per-frame scale G across frames (centered moving average of log G) before the rot3d
// combine, so a rocking event's partials share a consistent scale. Given as a ROTATION RANGE in
@@ -110,6 +127,8 @@ public:
ScalingSettings& StillsPartialityRefine(bool input);
ScalingSettings& ExpectedVarianceMerge(bool input);
ScalingSettings& IceRingMergeMask(bool input);
+ ScalingSettings& IceMinScore(float input);
+ ScalingSettings& IceMinSpotRatio(float input);
ScalingSettings& SmoothGDegrees(double input);
ScalingSettings& RelativeBDegrees(double input);
@@ -146,6 +165,8 @@ public:
[[nodiscard]] bool GetStillsPartialityRefine() const;
[[nodiscard]] bool GetExpectedVarianceMerge() const;
[[nodiscard]] bool GetIceRingMergeMask() const;
+ [[nodiscard]] float GetIceMinScore() const;
+ [[nodiscard]] float GetIceMinSpotRatio() const;
[[nodiscard]] double GetSmoothGDegrees() const;
[[nodiscard]] double GetRelativeBDegrees() const;
diff --git a/common/ScanResult.h b/common/ScanResult.h
index bd705165..92a409bd 100644
--- a/common/ScanResult.h
+++ b/common/ScanResult.h
@@ -31,6 +31,7 @@ struct ScanResultElem {
std::optional spot_count_low_res;
std::optional spot_count_indexed;
std::optional spot_count_ice;
+ std::optional spot_count_ice_control;
std::optional indexing_solution;
std::optional profile_radius;
std::optional b_factor;
diff --git a/common/ScanResultGenerator.cpp b/common/ScanResultGenerator.cpp
index f9256ffa..1510da68 100644
--- a/common/ScanResultGenerator.cpp
+++ b/common/ScanResultGenerator.cpp
@@ -56,6 +56,7 @@ void ScanResultGenerator::Add(const DataMessage &message) {
v[image_number].max_viable_pixel = message.max_viable_pixel_value;
v[image_number].sat_pixels = message.saturated_pixel_count;
v[image_number].spot_count_ice = message.spot_count_ice_rings;
+ v[image_number].spot_count_ice_control = message.spot_count_ice_control;
v[image_number].spot_count_low_res = message.spot_count_low_res;
v[image_number].spot_count_indexed = message.spot_count_indexed;
v[image_number].res = message.resolution_estimate;
@@ -93,6 +94,7 @@ void ScanResultGenerator::FillEndMessage(EndMessage &message) const {
message.data_collection_efficiency.resize(n);
message.spot_count.resize(n);
message.spot_count_ice_ring.resize(n);
+ message.spot_count_ice_control.resize(n);
message.spot_count_low_res.resize(n);
message.spot_count_indexed.resize(n);
message.image_indexed.resize(n);
@@ -124,6 +126,7 @@ void ScanResultGenerator::FillEndMessage(EndMessage &message) const {
message.data_collection_efficiency[number] = e.collection_efficiency;
message.spot_count[number] = static_cast(value_or_zero(e.spot_count));
message.spot_count_ice_ring[number] = static_cast(value_or_zero(e.spot_count_ice));
+ message.spot_count_ice_control[number] = e.spot_count_ice_control.value_or(NAN);
message.spot_count_low_res[number] = static_cast(value_or_zero(e.spot_count_low_res));
message.spot_count_indexed[number] = static_cast(value_or_zero(e.spot_count_indexed));
message.image_indexed[number] = static_cast(e.indexing_solution.value_or(0));
diff --git a/docs/CBOR.md b/docs/CBOR.md
index 0dbcd207..3536f9d6 100644
--- a/docs/CBOR.md
+++ b/docs/CBOR.md
@@ -216,6 +216,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s
| packets_received | uint64 | Number of packets received per image (in units of 2 kB) | | |
| bkg_estimate | float | Mean value for pixels in resolution range from 3.0 to 5.0 A \[photons\] | | |
| ice_ring_score | float | Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) | | |
+| spot_count_ice_control | float | Spots in the ice-free flanks beside the hexagonal rings, rescaled to the ring bands' own q width (control for spot_count_ice_rings) | | |
| beam_corr_x | float | Beam center correction X applied during processing \[pixel\] | | X |
| beam_corr_y | float | Beam center correction Y applied during processing \[pixel\] | | X |
| image_scale_factor | float | Scaling result: Image scale factor (g) | | X |
@@ -306,6 +307,8 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s
| image_indexed | Array(uint8) | Per-image indexing result; 0 = not indexed, nonzero = indexed | |
| v_bkg_estimate | Array(float) | Per-image background estimate | |
| ice_ring_score | Array(float) | Per-image strongest ice-ring intensity over the smooth radial background (1 = no ice) | |
+| spot_count_ice_control | Array(float) | Per-image spot count in the ice-free flanks beside the hexagonal rings, rescaled to the ring bands' q width | |
+| ice_ring_score_mean | float | Mean ice-ring score for the whole run (1 = no ice) | |
| profile_radius | Array(float) | Per-image profile radius \[Angstrom^-1\] | |
| mosaicity | Array(float) | Per-image mosaicity \[degree\] | |
| bFactor | Array(float) | Per-image estimated B-factor \[Angstrom^2\] | |
diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md
index e564c9df..4651331a 100644
--- a/docs/CPU_DATA_ANALYSIS.md
+++ b/docs/CPU_DATA_ANALYSIS.md
@@ -258,7 +258,11 @@ Because detection reads the pixel's ring, a pixel that falls outside the azimuth
Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, spots in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6).
-A single per-image **ice-ring score** is derived from the azimuthally-integrated radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins so the estimate is robust to the radial binning: at a coarse Q-spacing a local shoulder can be only ~1 bin and would double-count the ring's own edge (offline processing defaults to a fine 0.01 1/Å spacing, `--azim-q-spacing`, so the rings are well resolved). The reported quantity is the ice *magnitude* rather than a significance: with many photons any real ice ring is statistically significant, so significance does not discriminate. It is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity, distinct from the merge-time ice masking, which is data-driven from the per-ring merged CC1/2.
+A single per-image **ice-ring score** is derived from a radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins so the estimate is robust to the radial binning: at a coarse Q-spacing a local shoulder can be only ~1 bin and would double-count the ring's own edge (offline processing defaults to a fine 0.01 1/Å spacing, `--azim-q-spacing`, so the rings are well resolved). The reported quantity is the ice *magnitude* rather than a significance: with many photons any real ice ring is statistically significant, so significance does not discriminate.
+
+The profile the score is read off is the **peak-excluded** one, not the plain azimuthal integration: where adaptive spot finding runs (§3.2 — the offline and viewer default), the score uses the sigma-clipped per-resolution-ring background that finder already computes for its threshold. This matters more than it sounds. A plain azimuthal profile is a per-ring *mean*, so a few strong low-resolution reflections landing in a ring's bin raise it exactly as ice would; measured over 37 rotation crystals that alone put ice-free crystals at scores of 1.5–4.2, above crystals that really are iced, and the strongest apparent "ice" in the set was a crystal with none. An ice ring is azimuthally smooth and survives the sigma clip, while Bragg peaks do not, so on the clipped profile the same 30 ice-free crystals sit at 1.00–1.22 and the three with confirmed ice at 2.08–2.37. Only where no adaptive finder ran (the FPGA workflow) does the score fall back to the plain profile.
+
+The score is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity, and offline it also **gates** ice handling: `--ice-min-score` (default 1.5) is the score a run must reach before ice-ring flagging, the exclusion from the scale fit and the merge-time ice mask are applied at all. The eleven fixed bands cover 16–26 % of the unique reflections at typical resolutions whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. Which rings are then dropped from the merge remains data-driven, from the per-ring merged CC1/2.
A further optional safeguard removes isolated high-resolution “spur” spots by detecting large gaps in $1/d$ (or $q$) space and discarding spots beyond the gap. This is intended for macromolecular diffraction where edge-of-detector backgrounds can be extremely low.
diff --git a/docs/HDF5.md b/docs/HDF5.md
index e75db39c..9aaba923 100644
--- a/docs/HDF5.md
+++ b/docs/HDF5.md
@@ -245,6 +245,7 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i
| `peakCountUnfiltered` | | spots found before filtering |
| `peakCountLowRes` | | low-resolution spots |
| `peakCountIceRingRes` | | spots inside ice-ring bands |
+| `peakCountIceRingControl` | | spots in the ice-free flanks beside those bands, rescaled to their q width - the control for the count above (their ratio, pooled over the run, is the spot-based ice indicator) |
| `peakCountIndexed` | | spots fitting the indexing solution |
| `imageIndexed` | | image was indexed (0/1) |
| `indexingLatticeCount` | | number of lattices found for the image |
@@ -277,6 +278,7 @@ variants.
| `rotationLatticeNiggliClass` | | Niggli class of the run lattice |
| `imageIndexedMean` | | mean indexing rate over the run |
| `bkgEstimateMean` | photons | mean background over the run |
+| `iceRingScoreMean` | ratio | mean `iceRingScore` over the run — the single "how icy was this dataset" number (1 = no ice) |
| `indexedLatticeCount` | | per-image lattice count summary (master). *Note: data files use `indexingLatticeCount`; readers accept either.* |
CrystFEL can read the spots directly with:
diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md
index 9c878f40..7bec03f2 100644
--- a/docs/RUGNUX.md
+++ b/docs/RUGNUX.md
@@ -269,6 +269,8 @@ Scaling and merging:
| `--resolution-shells ` | Number of resolution shells in the reported statistics table (default: 10) |
| `--min-partiality ` | Minimum partiality to accept a reflection (default: 0.02) |
| `--ice-ring-mask[=on\|off]` | Merge-time ice-ring mask: after a first merge, drop a hexagonal-ice ring whose merged half-set CC1/2 has collapsed below its resolution shoulders and merge again (default: on; only acts when `--detect-ice-rings` is on). `off` disables **only** this mask — ice-spot flagging and the ice exclusion from scaling stay as `--detect-ice-rings` set them |
+| `--ice-min-score ` | Ice-presence gate: the measured per-run ice score (1 = no ice) a dataset must reach before **any** ice handling is applied — the flagging, the exclusion from scaling and the merge-time mask (default: 1.5; 0 = no gate). The eleven fixed hexagonal bands cover 16–26 % of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal only costs completeness |
+| `--ice-min-spot-ratio ` | The second ice-presence channel: found **spots** on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so `--ice-min-score` alone is blind to it (default: 2.0; 0 disables this channel) |
| `--reject-outliers ` | Per-observation outlier rejection, N σ from the per-reflection median (default: 6 for `rot3d`, off otherwise) |
| `--min-image-cc ` | Per-image CC limit, percent (default: no limit) |
| `--search-min-zeta ` | De-novo space-group search only: also search a merge of just the observations whose Lorentz geometry \|ζ\| reaches this, and keep whichever search found more symmetry (default: 0.85 for rotation, 0 = single search). Reflections crossing the Ewald sphere near-tangentially are measured worst and can make a real symmetry operator look like a twin law. The point group only — the systematic absences always come from the merge of all the observations |
diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp
index 2593d9ee..35aaa7cf 100644
--- a/frame_serialize/CBORStream2Deserializer.cpp
+++ b/frame_serialize/CBORStream2Deserializer.cpp
@@ -717,6 +717,8 @@ namespace {
message.spot_count_indexed = GetCBORUInt(value);
else if (key == "spot_count_ice_rings")
message.spot_count_ice_rings = GetCBORUInt(value);
+ else if (key == "spot_count_ice_control")
+ message.spot_count_ice_control = GetCBORFloat(value);
else if (key == "az_int_profile")
GetCBORFloatArray(value, message.az_int_profile);
else if (key == "az_int_profile_std")
@@ -1427,6 +1429,10 @@ namespace {
GetCBORFloatArray(value, message.v_bkg_estimate);
else if (key == "ice_ring_score")
GetCBORFloatArray(value, message.ice_ring_score);
+ else if (key == "spot_count_ice_control")
+ GetCBORFloatArray(value, message.spot_count_ice_control);
+ else if (key == "ice_ring_score_mean")
+ message.ice_ring_score_mean = GetCBORFloat(value);
else if (key == "profile_radius")
GetCBORFloatArray(value, message.profile_radius);
else if (key == "mosaicity")
diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp
index 229b8510..70e29d9f 100644
--- a/frame_serialize/CBORStream2Serializer.cpp
+++ b/frame_serialize/CBORStream2Serializer.cpp
@@ -774,6 +774,8 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) {
CBOR_ENC(mapEncoder, "image_indexed", message.image_indexed);
CBOR_ENC(mapEncoder, "v_bkg_estimate", message.v_bkg_estimate);
CBOR_ENC(mapEncoder, "ice_ring_score", message.ice_ring_score);
+ CBOR_ENC(mapEncoder, "ice_ring_score_mean", message.ice_ring_score_mean);
+ CBOR_ENC(mapEncoder, "spot_count_ice_control", message.spot_count_ice_control);
CBOR_ENC(mapEncoder, "profile_radius", message.profile_radius);
CBOR_ENC(mapEncoder, "mosaicity", message.mosaicity);
CBOR_ENC(mapEncoder, "bFactor", message.bFactor);
@@ -807,6 +809,7 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons
CBOR_ENC(mapEncoder, "spot_count", message.spot_count);
CBOR_ENC(mapEncoder, "spot_count_ice_rings", message.spot_count_ice_rings);
+ CBOR_ENC(mapEncoder, "spot_count_ice_control", message.spot_count_ice_control);
CBOR_ENC(mapEncoder, "spot_count_low_res", message.spot_count_low_res);
CBOR_ENC(mapEncoder, "spot_count_indexed", message.spot_count_indexed);
CBOR_ENC(mapEncoder, "az_int_profile", message.az_int_profile);
diff --git a/frontend/src/client/types.gen.ts b/frontend/src/client/types.gen.ts
index d09a64a6..460d6004 100644
--- a/frontend/src/client/types.gen.ts
+++ b/frontend/src/client/types.gen.ts
@@ -733,7 +733,7 @@ export type spot_finding_settings = {
*/
quick_integration: boolean;
/**
- * Width of ice ring in q-space in reciprocal space
+ * Half-width of the ice ring band in q (1/A). Matches the offline default in image_analysis/spot_finding/SpotFindingSettings.h, which was set from a measured ring FWHM of ~0.06; the two must agree or the same data gets a narrower ice band online.
*/
ice_ring_width_q_recipA: number;
/**
diff --git a/frontend/src/client/zod.gen.ts b/frontend/src/client/zod.gen.ts
index 0ba3c5f0..13b1a2e5 100644
--- a/frontend/src/client/zod.gen.ts
+++ b/frontend/src/client/zod.gen.ts
@@ -304,7 +304,7 @@ export const zSpotFindingSettings = z.object({
low_resolution_limit: z.number(),
high_resolution_limit_for_spot_count_low_res: z.number().gte(2).lte(8),
quick_integration: z.boolean().default(false),
- ice_ring_width_q_recipA: z.number().gte(0).lte(1).default(0.02),
+ ice_ring_width_q_recipA: z.number().gte(0).lte(1).default(0.03),
high_res_gap_Q_recipA: z.number().gte(0.1).lte(5).optional().default(1.5),
adaptive_threshold: z.boolean().optional().default(false),
false_pixels_per_frame: z.number().gte(1).lte(100000).optional().default(100)
diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp
index 385e7da8..40b81d1b 100644
--- a/image_analysis/MXAnalysisWithoutFPGA.cpp
+++ b/image_analysis/MXAnalysisWithoutFPGA.cpp
@@ -246,8 +246,21 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
output.az_int_profile_std = profile.GetStd();
output.bkg_estimate = profile.GetBkgEstimate(integration.Settings());
- output.ice_ring_score = profile.GetIceRingScore(integration.Settings(),
- spot_finding_settings.ice_ring_width_Q_recipA);
+
+ // The ice score wants a radial profile with the Bragg peaks taken OUT of it. The azimuthal profile
+ // is a plain per-ring mean, so a strong low-resolution reflection landing in a ring's bin is
+ // indistinguishable from ice sitting there - measured, that alone lifts clean crystals to a score
+ // of 1.5-4.2, right into the range real ice occupies. The adaptive spot finder already computes
+ // exactly what is wanted: a sigma-clipped per-ring background, in the same bins, from which the
+ // peaks have been removed (an ice ring is azimuthally smooth, so it survives the clip). It is in
+ // raw counts rather than corrected ones, which the score does not care about - it is a ratio to the
+ // background interpolated under the ring, and the corrections are smooth in radius.
+ const std::vector &ring_bkg = adaptiveSpotFinder->GetRingBackground();
+ const bool have_ring_bkg = spot_finding_settings.enable && spot_finding_settings.adaptive_threshold
+ && !ring_bkg.empty();
+ output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore(
+ have_ring_bkg ? ring_bkg : profile.GetResult1D(), integration.GetQBinCount(),
+ integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA);
}
void MXAnalysisWithoutFPGA::RebuildROI() {
diff --git a/image_analysis/scale_merge/IceRingMask.cpp b/image_analysis/scale_merge/IceRingMask.cpp
index eb511f98..6c46173b 100644
--- a/image_analysis/scale_merge/IceRingMask.cpp
+++ b/image_analysis/scale_merge/IceRingMask.cpp
@@ -26,7 +26,13 @@ std::vector FindDecorrelatedIceRings(const std::vector &
continue;
const float dq = std::fabs(two_pi / m.d - q_ring);
if (dq < w) { ring.Add(m.I_half[0], m.I_half[1]); ++n_ring; }
- else if (dq < 3.0f * w) { shoulder.Add(m.I_half[0], m.I_half[1]); ++n_shoulder; }
+ // The shoulder is the control, so it has to be free of ice itself. The hexagonal rings are
+ // not evenly spaced - 1.947/1.916/1.882 A sit 0.05-0.06 apart in q - so for those three the
+ // [w, 3w) band around one ring lands squarely on its neighbours, and the test ends up
+ // comparing ice against ice. Measured, that is the ONLY thing this exclusion changes: over
+ // the battery it removes firings on those three rings and leaves every other firing's CC
+ // pair identical to three decimals.
+ else if (dq < 3.0f * w && !IsOnIceRing(m.d, w)) { shoulder.Add(m.I_half[0], m.I_half[1]); ++n_shoulder; }
}
// The 0.10 margin is the 99th percentile of this statistic measured on DECOY bands - the same
// ring/shoulder geometry evaluated at q positions carrying no ice ring - over the 37-crystal
diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
index ca63d69a..a60e7832 100644
--- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
+++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp
@@ -19,6 +19,7 @@ AdaptiveSpotFinderCPU::AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &
ring_mean.assign(nbins, 0.0f);
ring_sigma.assign(nbins, 0.0f);
ring_thr.assign(nbins, 0.0f);
+ ring_bkg.assign(nbins, NAN);
}
// Accumulate per-ring mean/variance from the raw (photon) image. clip_k <= 0 -> use every valid
@@ -80,9 +81,13 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image,
if (n_total == 0) {
// Nothing valid to threshold against: leave no strong pixels for ExtractSpots to build on.
std::fill(output_buffer.begin(), output_buffer.end(), 0);
+ std::fill(ring_bkg.begin(), ring_bkg.end(), NAN);
return;
}
+ for (size_t b = 0; b < nbins; ++b)
+ ring_bkg[b] = (ring_cnt[b] < adaptive_threshold::MIN_RING_PIXELS) ? NAN : ring_mean[b];
+
const double E = std::max(1.0f, settings.false_pixels_per_frame);
double p = E / static_cast(n_total);
p = std::min(std::max(p, 1e-9), 0.1);
diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h
index fa67ebfc..c78d732b 100644
--- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h
+++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.h
@@ -37,10 +37,14 @@ class AdaptiveSpotFinderCPU : public ImageSpotFinder {
std::vector ring_mean;
std::vector ring_sigma;
std::vector ring_thr;
+ // ring_mean of the last Detect(), NaN where the ring holds too few pixels to be its own background.
+ // Kept separately because ring_mean carries the previous frame's value for an empty ring.
+ std::vector ring_bkg;
void AccumulateRings(const ImagePreprocessorBuffer &image, float clip_k);
public:
explicit AdaptiveSpotFinderCPU(const AzimuthalIntegrationMapping &mapping);
void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override;
+ [[nodiscard]] const std::vector &GetRingBackground() const override { return ring_bkg; }
};
diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu
index aa403eb0..a20caa01 100644
--- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu
+++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu
@@ -187,6 +187,7 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
host_sum(nbins),
host_sum2(nbins),
host_count(nbins),
+ host_bkg(nbins, NAN),
prof_sum(nbins),
prof_sum2(nbins),
prof_count(nbins),
@@ -256,6 +257,7 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting
}
if (n_total == 0) {
host_thr.clear();
+ std::fill(host_bkg.begin(), host_bkg.end(), NAN);
return;
}
@@ -273,11 +275,13 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting
for (int b = 0; b < nbins; ++b) {
if (host_count[b] < adaptive_threshold::MIN_RING_PIXELS) {
host_thr[b] = g_thr;
+ host_bkg[b] = NAN;
} else {
const double m = static_cast(static_cast(host_sum[b])) / host_count[b];
const double var = std::max(0.0, static_cast(host_sum2[b]) / host_count[b] - m * m);
host_thr[b] = adaptive_threshold::RingThreshold(static_cast(m),
static_cast(std::sqrt(var)), p, z);
+ host_bkg[b] = static_cast(m);
}
}
}
diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h
index fb7dc749..5eecba20 100644
--- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h
+++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h
@@ -75,6 +75,7 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder {
std::vector host_sum2; // clipped raw sum^2 } (exact integers - see the kernel)
std::vector host_count; // clipped raw count }
std::vector host_thr; // per-ring threshold (empty -> frame had no valid pixels)
+ std::vector host_bkg; // clipped per-ring mean, NaN where the ring is too sparse to trust
std::vector prof_sum; // plain corrected sum } azimuthal-integration profile
std::vector prof_sum2; // plain corrected sum^2 }
std::vector prof_count; // plain pixel count }
@@ -106,4 +107,5 @@ public:
// The azimuthal profile computed as a byproduct of the last Detect() - lets this engine replace the
// separate azint pass in the analysis pipeline.
[[nodiscard]] const AzimuthalIntegrationProfile &GetProfile() const { return last_profile; }
+ [[nodiscard]] const std::vector &GetRingBackground() const override { return host_bkg; }
};
diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp
index 74f57d5e..d50063f6 100644
--- a/image_analysis/spot_finding/ImageSpotFinder.cpp
+++ b/image_analysis/spot_finding/ImageSpotFinder.cpp
@@ -41,6 +41,11 @@ void ImageSpotFinder::SetResolutionMask(const std::vector &mask) {
res_mask_bits.back() |= ~((1u << (npixel % 32)) - 1u);
}
+const std::vector &ImageSpotFinder::GetRingBackground() const {
+ static const std::vector none;
+ return none;
+}
+
void ImageSpotFinder::ExtractComponentsHost(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings) {
// Collect the strong pixels first and read their values afterwards, instead of reading the image
diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h
index d03dc9db..781939a2 100644
--- a/image_analysis/spot_finding/ImageSpotFinder.h
+++ b/image_analysis/spot_finding/ImageSpotFinder.h
@@ -45,6 +45,11 @@ public:
// components from those pixels.
virtual void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) = 0;
+ // Peak-excluded per-ring background of the last Detect(), in the bins of the azimuthal-integration
+ // mapping and in raw photon counts. Only the adaptive finders build one (it is what sets their
+ // threshold); empty for everyone else, and for a frame with nothing valid to reduce.
+ [[nodiscard]] virtual const std::vector &GetRingBackground() const;
+
// Pixels to ignore, one bool per pixel (true = ignore). Set when the resolution limits change,
// not per image: the GPU finders keep a bit-packed device copy of it, and re-uploading that for
// every image would cost more than the extraction it feeds.
diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp
index 81dc4f4f..304653cf 100644
--- a/image_analysis/spot_finding/SpotUtils.cpp
+++ b/image_analysis/spot_finding/SpotUtils.cpp
@@ -22,6 +22,39 @@ void CountSpots(DataMessage &msg,
msg.spot_count_ice_rings = ice_ring;
}
+// Spots in the ice-free control flanks either side of the hexagonal rings, rescaled to the ring bands'
+// own q width. The control for one ring is the two intervals [w, 2w) beside it - same total width as
+// the ring band, and symmetric, so the fall-off of spot density with resolution cancels to first
+// order. A flank that lands on another ring is not a control and is dropped, its width with it; the
+// three rings at 1.947/1.916/1.882 A are 0.05-0.06 apart in q and usually lose both.
+float CountIceRingControlSpots(const std::vector &spots, float w) {
+ if (!(w > 0.0f))
+ return 0.0f;
+ float control = 0.0f;
+ for (const float d : ICE_RING_RES_A) {
+ const float q_ring = 2 * PI / d;
+ bool lo_free = true, hi_free = true;
+ for (const float other : ICE_RING_RES_A) {
+ const float q_other = 2 * PI / other;
+ if (q_other > q_ring && q_other < q_ring + 3 * w) hi_free = false;
+ if (q_other < q_ring && q_other > q_ring - 3 * w) lo_free = false;
+ }
+ const int free_flanks = (lo_free ? 1 : 0) + (hi_free ? 1 : 0);
+ if (free_flanks == 0)
+ continue;
+ int64_t n = 0;
+ for (const auto &s: spots) {
+ if (!(s.d_A > 0.0f)) continue;
+ const float dq = 2 * PI / s.d_A - q_ring;
+ if (hi_free && dq >= w && dq < 2 * w) n++;
+ if (lo_free && dq <= -w && dq > -2 * w) n++;
+ }
+ // One free flank covers half the ring band's width, so it counts double.
+ control += static_cast(n) * 2.0f / static_cast(free_flanks);
+ }
+ return control;
+}
+
void MarkIceRings(std::vector &spots, float tolerance_q_recipA) {
std::vector ice_rings_q;
@@ -146,8 +179,12 @@ void SpotAnalyze(const DiffractionExperiment &experiment,
if (spot_finding_settings.high_res_gap_Q_recipA.has_value())
FilterSpuriousHighResolutionSpots(spots_out, spot_finding_settings.high_res_gap_Q_recipA.value());
- if (experiment.GetDatasetSettings().IsDetectIceRings() && spot_finding_settings.ice_ring_width_Q_recipA > 0.0f)
+ if (experiment.GetDatasetSettings().IsDetectIceRings() && spot_finding_settings.ice_ring_width_Q_recipA > 0.0f) {
MarkIceRings(spots_out, spot_finding_settings.ice_ring_width_Q_recipA);
+ // Before FilterSpotsByCount below, which orders ice spots LAST and would throw them away first.
+ output.spot_count_ice_control =
+ CountIceRingControlSpots(spots_out, spot_finding_settings.ice_ring_width_Q_recipA);
+ }
CountSpots(output, spots_out, spot_finding_settings.cutoff_spot_count_low_res);
diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h
index 9ba458be..23ca156d 100644
--- a/image_analysis/spot_finding/SpotUtils.h
+++ b/image_analysis/spot_finding/SpotUtils.h
@@ -16,6 +16,8 @@ void CountSpots(DataMessage &msg,
const std::vector &spots,
float d_min_A);
+float CountIceRingControlSpots(const std::vector &spots, float half_width_q_recipA);
+
void MarkIceRings(std::vector &spots, float tolerance_q_recipA);
void FilterSpotsByCount(std::vector &input, int64_t count);
diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp
index 27eb96d3..c4205dd4 100644
--- a/reader/HDF5MetadataSource.cpp
+++ b/reader/HDF5MetadataSource.cpp
@@ -360,6 +360,7 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen
dataset->spot_count_low_res = master_file->ReadOptVector("/entry/MX/peakCountLowRes");
dataset->spot_count_indexed = master_file->ReadOptVector("/entry/MX/peakCountIndexed");
dataset->spot_count_ice_rings = master_file->ReadOptVector("/entry/MX/peakCountIceRingRes");
+ dataset->spot_count_ice_control = master_file->ReadOptVector("/entry/MX/peakCountIceRingControl");
dataset->indexing_result = master_file->ReadOptVector("/entry/MX/imageIndexed");
dataset->bkg_estimate = master_file->ReadOptVector("/entry/MX/bkgEstimate");
@@ -468,6 +469,9 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen
data_file, "/entry/MX/nPeaks",
number_of_images, fimages);
+ ReadVector(dataset->spot_count_ice_control,
+ data_file, "/entry/MX/peakCountIceRingControl",
+ number_of_images, fimages);
ReadVector(dataset->spot_count_ice_rings,
data_file, "/entry/MX/peakCountIceRingRes",
number_of_images, fimages);
diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h
index c38ea2c1..3f2f0427 100644
--- a/reader/JFJochReaderDataset.h
+++ b/reader/JFJochReaderDataset.h
@@ -36,6 +36,7 @@ struct JFJochReaderDataset {
std::vector spot_count_indexed;
std::vector spot_count_low_res;
std::vector spot_count_ice_rings;
+ std::vector spot_count_ice_control;
std::vector indexing_result;
std::vector indexing_lattice_count;
diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp
index 17ef151d..7c1a2179 100644
--- a/receiver/JFJochReceiver.cpp
+++ b/receiver/JFJochReceiver.cpp
@@ -162,6 +162,7 @@ void JFJochReceiver::SendEndMessage() {
message.run_name = experiment.GetRunName();
message.bkg_estimate = plots.GetBkgEstimate();
+ message.ice_ring_score_mean = plots.GetIceRingScore();
message.indexing_rate = plots.GetIndexingRate();
message.az_int_result["dataset"] = plots.GetAzIntProfile();
diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp
index 9b7fecb7..94d3abf9 100644
--- a/rugnux/Rugnux.cpp
+++ b/rugnux/Rugnux.cpp
@@ -1054,6 +1054,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
end_msg.run_number = experiment_.GetRunNumber();
end_msg.run_name = experiment_.GetRunName();
end_msg.bkg_estimate = plots.GetBkgEstimate();
+ end_msg.ice_ring_score = plots.GetIceRingScoreArray();
+ end_msg.ice_ring_score_mean = plots.GetIceRingScore();
end_msg.az_int_result["dataset"] = plots.GetAzIntProfile();
end_msg.indexing_rate = result.indexing_rate;
@@ -1094,7 +1096,32 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// they drag the per-image scale fit. Flag them here so scaling (the per-image G fit and the
// fulls refit) skips them, while the combine, the merge and the statistics keep them - dropping
// them outright guts low/mid-resolution completeness on crystals that merge them fine.
- if (experiment_.IsDetectIceRings()) {
+ // The eleven hexagonal-ice bands cover 16-26% of the unique reflections at typical resolutions
+ // whether or not the crystal has any ice, so handling ice unconditionally taxes clean data for
+ // nothing. Gate the whole mechanism - this flagging and the merge-time CC1/2 ring mask below -
+ // on the run's own measured ice strength (1 = no ice; the strongest ring over the smooth radial
+ // background, measured on the peak-excluded profile so Bragg peaks cannot pose as ice). The
+ // score is logged either way, so a run that gates off still reports why.
+ // Ice is looked for on TWO channels, because it arrives in two forms. Fine-grained ice is a
+ // smooth powder ring and shows up in the radial profile (ice_ring_score). Ice in large
+ // crystallites arrives as discrete reflections instead, leaves the profile flat, and is only
+ // visible as a pile-up of found SPOTS on the ring positions - measured against the same q width
+ // of ice-free flanks beside them, so the fall-off of spot density with resolution cancels. On
+ // the rotation battery the two channels barely overlap: the crystals with smooth ice read
+ // 2.1-2.4 / ~1.0, the ones whose spots are mostly ice read ~1.1 / 3.8-14.6, and a clean crystal
+ // reads 1.04 on both. Either channel alone is blind to half the cases.
+ const auto ice_score = plots.GetIceRingScore();
+ const auto ice_spot_ratio = plots.GetIceRingSpotRatio();
+ const float ice_min_score = experiment_.GetScalingSettings().GetIceMinScore();
+ const float ice_min_spot_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio();
+ const bool ice_present = (!ice_score.has_value() || *ice_score >= ice_min_score)
+ || (ice_min_spot_ratio > 0.0f && ice_spot_ratio.has_value()
+ && *ice_spot_ratio >= ice_min_spot_ratio);
+ if (experiment_.IsDetectIceRings() && !ice_present) {
+ logger.Info("Ice-ring handling: measured ice score {:.2f} < {:.2f} and spot ratio {:.2f} < "
+ "{:.2f}, no ice detected - ice-ring handling skipped entirely",
+ *ice_score, ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio);
+ } else if (experiment_.IsDetectIceRings()) {
const float ice_width = config_.spot_finding.ice_ring_width_Q_recipA;
size_t total = 0, flagged = 0;
for (auto &outcome : indexer->GetIntegrationOutcome()) {
@@ -1105,8 +1132,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
++flagged;
}
}
- logger.Info("Ice-ring handling: flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); "
- "excluded from scaling, kept for merging", flagged, total, ice_width);
+ logger.Info("Ice-ring handling: ice score {:.2f} (gate {:.2f}), spot ratio {:.2f} (gate "
+ "{:.2f}); flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); "
+ "excluded from scaling, kept for merging",
+ ice_score.value_or(NAN), ice_min_score, ice_spot_ratio.value_or(NAN),
+ ice_min_spot_ratio, flagged, total, ice_width);
}
// Scale the images and merge. Factored so it can run twice: first in P1 to give the
@@ -1696,7 +1726,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// re-merge. Weak/absent rings track their neighbours and stay, so completeness on clean crystals
// is untouched. Uses the ring band (+/- ice width in q=2pi/d) vs the shoulders either side.
if (experiment_.IsDetectIceRings() && experiment_.GetScalingSettings().GetIceRingMergeMask()
- && !sm.merged.empty()) {
+ && ice_present && !sm.merged.empty()) {
auto mask = FindDecorrelatedIceRings(sm.merged, config_.spot_finding.ice_ring_width_Q_recipA,
logger);
if (!mask.empty()) {
diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp
index 9672cd8b..aa115b6e 100644
--- a/rugnux/rugnux_cli.cpp
+++ b/rugnux/rugnux_cli.cpp
@@ -110,6 +110,8 @@ void print_usage() {
std::cout << " --resolution-cutoff Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl;
std::cout << " --resolution-cc-target CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl;
std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 10)" << std::endl;
+ std::cout << " --ice-min-score Ice-presence gate: measured ice score (1 = no ice) a run must reach before ANY ice handling is applied - the flagging, the exclusion from scaling and the merge-time mask (default: 1.5). The eleven hexagonal bands cover 16-26% of the unique reflections whether or not the crystal has ice, so handling ice on a clean crystal is a pure loss. 0 = no gate (always handle ice)" << std::endl;
+ std::cout << " --ice-min-spot-ratio Second ice-presence channel: found spots on the hexagonal rings over the same q width of ice-free flanks beside them (1 = spots spread evenly). Ice in large crystallites diffracts as discrete spots and leaves the radial profile flat, so --ice-min-score alone is blind to it. Default 2.0; 0 disables this channel" << std::endl;
std::cout << " --ice-ring-mask[=on|off] Drop a hexagonal-ice ring from the merge when its merged half-set CC1/2 has collapsed below its resolution shoulders, then re-merge (default: on; needs --detect-ice-rings on). Off disables ONLY this merge-time mask - ice-spot flagging and the ice exclusion from scaling stay as --detect-ice-rings sets them" << std::endl;
std::cout << " --min-partiality Minimum partiality to accept reflection (default: 0.02)" << std::endl;
std::cout << " --capture-uncertainty rot3d: systematic sigma ~num*(1-captured_fraction)*I on under-captured fulls (default: 1.0 for rot3d, 0 otherwise)" << std::endl;
@@ -198,6 +200,8 @@ enum {
OPT_NO_EXPECTED_VARIANCE_MERGE,
OPT_DETECT_ICE_RINGS,
OPT_ICE_RING_MASK,
+ OPT_ICE_MIN_SCORE,
+ OPT_ICE_MIN_SPOT_RATIO,
OPT_NO_SCALE_FULLS,
OPT_WRITE_PROCESS_H5,
OPT_FORCE_STILL,
@@ -298,6 +302,8 @@ static option long_options[] = {
{"simple-stills", no_argument, nullptr, OPT_SIMPLE_STILLS},
{"detect-ice-rings", optional_argument, nullptr, OPT_DETECT_ICE_RINGS},
{"ice-ring-mask", optional_argument, nullptr, OPT_ICE_RING_MASK},
+ {"ice-min-score", required_argument, nullptr, OPT_ICE_MIN_SCORE},
+ {"ice-min-spot-ratio", required_argument, nullptr, OPT_ICE_MIN_SPOT_RATIO},
{"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
{nullptr, 0, nullptr, 0}
};
@@ -531,6 +537,8 @@ static int RunRugnux(int argc, char **argv) {
bool write_process_h5_flag = false; // --write-process-h5; also write _process.h5 when merging
std::optional detect_ice_rings; // --detect-ice-rings[=on|off]; unset => use the dataset (file) value
bool ice_ring_mask = true; // --ice-ring-mask[=on|off]; merge-time CC1/2 ice-ring mask
+ std::optional ice_min_score_arg; // --ice-min-score: ice-presence gate on the measured score
+ std::optional ice_min_spot_ratio_arg; // --ice-min-spot-ratio: the same gate on the spot channel
std::optional min_q, max_q, q_spacing; // azimuthal integration range / -q spacing (1/A)
std::optional azimuthal_bins; // --azimuthal-bins
std::optional polarization_correction; // --polarization-correction (azimuthal integration)
@@ -862,6 +870,12 @@ static int RunRugnux(int argc, char **argv) {
exit(EXIT_FAILURE);
}
break;
+ case OPT_ICE_MIN_SCORE:
+ ice_min_score_arg = parse_double_arg(optarg, "--ice-min-score", logger);
+ break;
+ case OPT_ICE_MIN_SPOT_RATIO:
+ ice_min_spot_ratio_arg = parse_double_arg(optarg, "--ice-min-spot-ratio", logger);
+ break;
case OPT_WRITE_PROCESS_H5:
write_process_h5_flag = true;
break;
@@ -1188,6 +1202,10 @@ static int RunRugnux(int argc, char **argv) {
scaling_settings.StillsPartialityRefine(!simple_stills_flag);
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
scaling_settings.IceRingMergeMask(ice_ring_mask);
+ if (ice_min_score_arg)
+ scaling_settings.IceMinScore(static_cast(*ice_min_score_arg));
+ if (ice_min_spot_ratio_arg)
+ scaling_settings.IceMinSpotRatio(static_cast(*ice_min_spot_ratio_arg));
scaling_settings.OutlierRejectNsigma(
outlier_reject_nsigma.value_or(scaling_settings.GetOutlierRejectNsigma()));
scaling_settings.ScaleFulls(scale_fulls_arg.value_or(scaling_settings.GetScaleFulls()));
@@ -1211,7 +1229,33 @@ static int RunRugnux(int argc, char **argv) {
// otherwise --scale re-scales a dataset the writing run had scaled without those reflections,
// and the per-image scales come out of a different fit than the ones in the file.
const float ice_width = SpotFindingSettings().ice_ring_width_Q_recipA;
- if (experiment.IsDetectIceRings()) {
+ // ...and gated the same way, on the per-image ice score the writing run stored in the file, so
+ // --scale reaches the same verdict on the same data as the pipeline that produced it.
+ double ice_sum = 0.0;
+ size_t ice_n = 0;
+ for (const float s : dataset->ice_ring_score)
+ if (std::isfinite(s)) {
+ ice_sum += s;
+ ++ice_n;
+ }
+ const float ice_min_score = experiment.GetScalingSettings().GetIceMinScore();
+ // ...and the spot channel, pooled over the run exactly as the full pipeline pools it.
+ double ring_sum = 0.0, ctrl_sum = 0.0;
+ for (const float v : dataset->spot_count_ice_rings)
+ if (std::isfinite(v)) ring_sum += v;
+ for (const float v : dataset->spot_count_ice_control)
+ if (std::isfinite(v)) ctrl_sum += v;
+ // Empty control + spots on the rings = the strongest ice evidence, not its absence.
+ const double ice_spot_ratio = ctrl_sum > 0.0 ? ring_sum / ctrl_sum : (ring_sum > 0.0 ? 1.0e3 : 0.0);
+ const float ice_min_spot_ratio = experiment.GetScalingSettings().GetIceMinSpotRatio();
+ const bool ice_present = (ice_n == 0 || ice_sum / static_cast(ice_n) >= ice_min_score)
+ || (ice_min_spot_ratio > 0.0f && ice_spot_ratio >= ice_min_spot_ratio);
+ if (experiment.IsDetectIceRings() && !ice_present) {
+ logger.Info("Ice-ring handling: measured ice score {:.2f} and spot ratio {:.2f} below the "
+ "gates ({:.2f} / {:.2f}), no ice detected - ice-ring handling skipped entirely",
+ ice_sum / static_cast(ice_n), ice_spot_ratio, ice_min_score,
+ ice_min_spot_ratio);
+ } else if (experiment.IsDetectIceRings()) {
size_t total = 0, flagged = 0;
for (auto &outcome : reflections) {
for (auto &r : outcome.reflections) {
@@ -1260,7 +1304,7 @@ static int RunRugnux(int argc, char **argv) {
error_model_isa = r.isa;
};
run({});
- if (experiment.IsDetectIceRings() && experiment.GetScalingSettings().GetIceRingMergeMask()) {
+ if (experiment.IsDetectIceRings() && experiment.GetScalingSettings().GetIceRingMergeMask() && ice_present) {
masked_ice_rings = FindDecorrelatedIceRings(merged_reflections, ice_width, logger);
if (!masked_ice_rings.empty())
run(masked_ice_rings);
@@ -1317,7 +1361,7 @@ static int RunRugnux(int argc, char **argv) {
error_model_isa = merge_engine.ErrorModelB() > 0 ? 1.0 / merge_engine.ErrorModelB() : 0.0;
};
merge({});
- if (experiment.IsDetectIceRings() && experiment.GetScalingSettings().GetIceRingMergeMask()) {
+ if (experiment.IsDetectIceRings() && experiment.GetScalingSettings().GetIceRingMergeMask() && ice_present) {
masked_ice_rings = FindDecorrelatedIceRings(merged_reflections, ice_width, logger);
if (!masked_ice_rings.empty())
merge(masked_ice_rings);
@@ -1547,6 +1591,10 @@ static int RunRugnux(int argc, char **argv) {
scaling_settings.StillsPartialityRefine(!simple_stills_flag);
scaling_settings.ExpectedVarianceMerge(!no_expected_variance_merge);
scaling_settings.IceRingMergeMask(ice_ring_mask);
+ if (ice_min_score_arg)
+ scaling_settings.IceMinScore(static_cast(*ice_min_score_arg));
+ if (ice_min_spot_ratio_arg)
+ scaling_settings.IceMinSpotRatio(static_cast(*ice_min_spot_ratio_arg));
if (d_min_scale_merge)
scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value());
if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method);
diff --git a/tests/AdaptiveSpotFinderCPUTest.cpp b/tests/AdaptiveSpotFinderCPUTest.cpp
index 321a3776..819ccbf1 100644
--- a/tests/AdaptiveSpotFinderCPUTest.cpp
+++ b/tests/AdaptiveSpotFinderCPUTest.cpp
@@ -1,6 +1,8 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-only
+#include
+
#include
#include "../common/AzimuthalIntegrationMapping.h"
@@ -130,3 +132,61 @@ TEST_CASE("AdaptiveSpotFinderCPU_ThresholdTracksBackground", "[AdaptiveSpotFinde
CHECK(std::lround(scaled[0].RawCoord().x) == std::lround(plain[0].RawCoord().x));
CHECK(std::lround(scaled[0].RawCoord().y) == std::lround(plain[0].RawCoord().y));
}
+
+// The per-ring background the finder hands out is the peak-EXCLUDED one: a handful of very bright
+// pixels planted on a ring must not move it, which is the property that lets the ice score be read
+// off it instead of off the plain azimuthal profile. Rings with too few pixels to be their own
+// background come back as NaN rather than as a stale value from the previous frame.
+TEST_CASE("AdaptiveSpotFinderCPU_RingBackgroundExcludesPeaks", "[AdaptiveSpotFinder]") {
+ DiffractionExperiment x(DetJF4M());
+ x.DetectorDistance_mm(80).BeamX_pxl(1030).BeamY_pxl(1080);
+ x.QSpacingForAzimInt_recipA(0.05).QRangeForAzimInt_recipA(0.05, 5.0);
+ x.GeometryTransformation(false);
+
+ PixelMask pixel_mask(x);
+ AzimuthalIntegrationMapping mapping(x, pixel_mask);
+ const auto &pixel_to_bin = mapping.GetPixelToBin();
+
+ const size_t w = x.GetXPixelsNum();
+ const size_t h = x.GetYPixelsNum();
+
+ ImagePreprocessorBuffer buffer(x.GetPixelsNum());
+ for (size_t i = 0; i < w * h; i++)
+ buffer[i] = 8 + static_cast(i % 5); // background 8..12
+
+ std::vector res_mask(x.GetPixelsNum(), false);
+ AdaptiveSpotFinderCPU finder(mapping);
+ finder.SetResolutionMask(res_mask);
+
+ finder.Detect(buffer, AdaptiveSettings());
+ const std::vector clean = finder.GetRingBackground();
+ REQUIRE(clean.size() == mapping.GetBinNumber());
+
+ std::vector pixels_per_bin(mapping.GetBinNumber(), 0);
+ for (size_t i = 0; i < w * h; i++)
+ if (pixel_to_bin[i] < mapping.GetBinNumber())
+ pixels_per_bin[pixel_to_bin[i]]++;
+
+ // Plant 200 bright pixels spread over one well-populated ring - far more than a real spot, so a
+ // plain mean would move visibly.
+ const uint16_t ring = static_cast(
+ std::max_element(pixels_per_bin.begin(), pixels_per_bin.end()) - pixels_per_bin.begin());
+ REQUIRE(pixels_per_bin[ring] > 10000);
+ int planted = 0;
+ for (size_t i = 0; i < w * h && planted < 200; i++)
+ if (pixel_to_bin[i] == ring) {
+ buffer[i] = 100000;
+ planted++;
+ }
+ REQUIRE(planted == 200);
+
+ finder.Detect(buffer, AdaptiveSettings());
+ const std::vector spiked = finder.GetRingBackground();
+
+ REQUIRE(std::isfinite(clean[ring]));
+ REQUIRE(std::isfinite(spiked[ring]));
+ CHECK(spiked[ring] == Catch::Approx(clean[ring]).epsilon(0.01));
+
+ for (size_t b = 0; b < clean.size(); b++)
+ CHECK(std::isfinite(clean[b]) == (pixels_per_bin[b] >= 40));
+}
diff --git a/tests/AzimuthalIntegrationTest.cpp b/tests/AzimuthalIntegrationTest.cpp
index 11235f1d..98dc3e8d 100644
--- a/tests/AzimuthalIntegrationTest.cpp
+++ b/tests/AzimuthalIntegrationTest.cpp
@@ -545,4 +545,52 @@ TEST_CASE("AzimuthalIntegrationMapping_Threading_ConvertedGeometry_18Modules", "
REQUIRE(x.GetModulesNum() == 18);
CheckAzimuthalIntegrationMappingThreadingExact(x);
-}
\ No newline at end of file
+}
+
+// The ice score on an explicit profile: flat means no ice (1.0), and a bump planted on a hexagonal
+// ring is reported at its own height over the background. Uses the static entry point, which is what
+// feeds the score the peak-excluded per-ring background instead of the plain profile.
+TEST_CASE("AzimuthalIntegrationProfile_IceRingScore","[AzimuthalIntegration]") {
+ AzimuthalIntegrationSettings settings;
+ settings.QSpacing_recipA(0.01f).QRange_recipA(0.1f, 4.5f);
+
+ const int q_bins = settings.GetQBinCount();
+ REQUIRE(q_bins > 400);
+
+ auto bin_of = [&](float d_A) {
+ const float q = 6.283185307f / d_A;
+ return static_cast(std::lround((q - settings.GetLowQ_recipA()) / settings.GetQSpacing_recipA() - 0.5f));
+ };
+
+ std::vector profile(q_bins, 100.0f);
+ CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(1.0f));
+
+ // A bump at the 2.249 A hexagonal ring, on a background the running median still reads as 100.
+ profile[bin_of(2.249f)] = 250.0f;
+ CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(2.5f));
+
+ // A bump of the same size well away from every ring is not ice and must not be reported.
+ std::vector off_ring(q_bins, 100.0f);
+ off_ring[bin_of(2.500f)] = 250.0f;
+ CHECK(AzimuthalIntegrationProfile::IceRingScore(off_ring, q_bins, settings, 0.03f) == Catch::Approx(1.0f));
+}
+
+// A profile given per (q, azimuth) bin is averaged over azimuth first, so a ring seen in every sector
+// scores the same as the equivalent 1-D profile.
+TEST_CASE("AzimuthalIntegrationProfile_IceRingScore_Azimuthal","[AzimuthalIntegration]") {
+ AzimuthalIntegrationSettings settings;
+ settings.QSpacing_recipA(0.01f).QRange_recipA(0.1f, 4.5f).AzimuthalBinCount(4);
+
+ const int q_bins = settings.GetQBinCount();
+ const float q = 6.283185307f / 2.249f;
+ const int ring = static_cast(std::lround((q - settings.GetLowQ_recipA()) / settings.GetQSpacing_recipA() - 0.5f));
+
+ std::vector flat(q_bins, 100.0f);
+ std::vector sectors(static_cast(q_bins) * 4, 100.0f);
+ flat[ring] = 250.0f;
+ for (int az = 0; az < 4; az++)
+ sectors[static_cast(az) * q_bins + ring] = 250.0f;
+
+ CHECK(AzimuthalIntegrationProfile::IceRingScore(sectors, q_bins, settings, 0.03f)
+ == Catch::Approx(AzimuthalIntegrationProfile::IceRingScore(flat, q_bins, settings, 0.03f)));
+}
diff --git a/writer/HDF5DataFilePluginMX.cpp b/writer/HDF5DataFilePluginMX.cpp
index 752e45e9..808a2006 100644
--- a/writer/HDF5DataFilePluginMX.cpp
+++ b/writer/HDF5DataFilePluginMX.cpp
@@ -86,6 +86,7 @@ void HDF5DataFilePluginMX::OpenFile(HDF5File &data_file, const DataMessage &msg,
spot_count_total.reserve(images_per_file);
spot_count_ice.reserve(images_per_file);
+ spot_count_ice_control.reserve(images_per_file);
spot_count_indexed.reserve(images_per_file);
spot_count_low_res.reserve(images_per_file);
integrated_reflections.reserve(images_per_file);
@@ -145,6 +146,7 @@ void HDF5DataFilePluginMX::Write(const DataMessage &msg, uint64_t image_number)
strong_pixel_count[image_number] = msg.strong_pixel_count.value_or(0);
spot_count_total[image_number] = msg.spot_count.value_or(0);
spot_count_ice[image_number] = msg.spot_count_ice_rings.value_or(0);
+ spot_count_ice_control[image_number] = msg.spot_count_ice_control.value_or(NAN);
spot_count_low_res[image_number] = msg.spot_count_low_res.value_or(0);
if (indexing) {
@@ -225,6 +227,8 @@ void HDF5DataFilePluginMX::WriteFinal(HDF5File &data_file) {
if (!spot_count_ice.empty())
data_file.SaveVector("/entry/MX/peakCountIceRingRes", spot_count_ice.vec());
+ if (!spot_count_ice_control.empty())
+ data_file.SaveVector("/entry/MX/peakCountIceRingControl", spot_count_ice_control.vec());
if (!spot_count_indexed.empty())
data_file.SaveVector("/entry/MX/peakCountIndexed", spot_count_indexed.vec());
if (!spot_count_low_res.empty())
diff --git a/writer/HDF5DataFilePluginMX.h b/writer/HDF5DataFilePluginMX.h
index 3f0265b1..68a568ed 100644
--- a/writer/HDF5DataFilePluginMX.h
+++ b/writer/HDF5DataFilePluginMX.h
@@ -30,6 +30,7 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin {
AutoIncrVector spot_count_total;
AutoIncrVector spot_count_ice;
+ AutoIncrVector spot_count_ice_control{NAN};
AutoIncrVector spot_count_indexed;
AutoIncrVector spot_count_low_res;
diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp
index 8cdffec1..fd779752 100644
--- a/writer/HDF5NXmx.cpp
+++ b/writer/HDF5NXmx.cpp
@@ -960,6 +960,9 @@ void NXmx::Finalize(const EndMessage &end) {
if (end.bkg_estimate) {
SaveScalar(*hdf5_file, "/entry/MX/bkgEstimateMean", end.bkg_estimate.value());
}
+ if (end.ice_ring_score_mean) {
+ SaveScalar(*hdf5_file, "/entry/MX/iceRingScoreMean", end.ice_ring_score_mean.value());
+ }
hdf5_file->Close();
hdf5_file.reset();
@@ -1029,6 +1032,7 @@ void NXmx::EndResultVectors(const EndMessage &end) {
mx_group.NXClass("NXcollection");
SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountIceRingRes", end.spot_count_ice_ring);
+ SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountIceRingControl", end.spot_count_ice_control);
SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountLowRes", end.spot_count_low_res);
SaveVectorIfMissing(*hdf5_file, "/entry/MX/peakCountIndexed", end.spot_count_indexed);
SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageIndexed", end.image_indexed);