Skip to content

Commit 8d45af3

Browse files
jl-wynenjl-wynen
committed
Apply monitor masks
1 parent 489c65e commit 8d45af3

File tree

2 files changed

+191
-3
lines changed

2 files changed

+191
-3
lines changed

‎src/ess/reduce/normalization.py‎

Lines changed: 44 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -2,12 +2,13 @@
22
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
33
"""Normalization routines for neutron data reduction."""
44

5+
import itertools
6+
57
import scipp as sc
68

79
from .uncertainty import UncertaintyBroadcastMode, broadcast_uncertainties
810

911

10-
# TODO explain impact of masking, esp multi-dim
1112
def normalize_by_monitor_histogram(
1213
detector: sc.DataArray,
1314
*,
@@ -84,14 +85,15 @@ def normalize_by_monitor_histogram(
8485
"""
8586
dim = monitor.dim
8687

88+
detector = _mask_detector_for_norm(detector=detector, monitor=monitor)
8789
clipped = _clip_monitor_to_detector_range(monitor=monitor, detector=detector)
8890
coord = clipped.coords[dim]
89-
delta_w = coord[1:] - coord[:-1]
91+
delta_w = sc.DataArray(coord[1:] - coord[:-1], masks=clipped.masks)
9092
total_monitor_weight = broadcast_uncertainties(
9193
clipped.sum() / delta_w.sum(),
9294
prototype=clipped,
9395
mode=uncertainty_broadcast_mode,
94-
)
96+
).data
9597
delta_w *= total_monitor_weight
9698
norm = broadcast_uncertainties(
9799
clipped / delta_w, prototype=detector, mode=uncertainty_broadcast_mode
@@ -148,6 +150,7 @@ def normalize_by_monitor_integrated(
148150
normalize_by_monitor_histogram:
149151
Normalize by a monitor histogram without integration.
150152
"""
153+
detector = _mask_detector_for_norm(detector=detector, monitor=monitor)
151154
clipped = _clip_monitor_to_detector_range(monitor=monitor, detector=detector)
152155
coord = clipped.coords[clipped.dim]
153156
norm = (clipped * (coord[1:] - coord[:-1])).data.sum()
@@ -201,3 +204,41 @@ def _clip_monitor_to_detector_range(
201204
# But integration would pick up all monitor bins.
202205
return monitor.rebin({dim: det_coord})
203206
return monitor[dim, lo:hi]
207+
208+
209+
def _mask_detector_for_norm(
210+
*, detector: sc.DataArray, monitor: sc.DataArray
211+
) -> sc.DataArray:
212+
"""Mask the detector where the monitor is masked.
213+
214+
For performance, this applies the monitor mask to the detector bins.
215+
This can lead to masking more vents than strictly necessary if we
216+
used an event mask.
217+
"""
218+
if (monitor_mask := _monitor_mask(monitor)) is None:
219+
return detector
220+
221+
# Use rebin to reshape the mask to the detector:
222+
dim = monitor.dim
223+
mask = sc.DataArray(monitor_mask, coords={dim: monitor.coords[dim]}).rebin(
224+
{dim: detector.coords[dim]}
225+
).data != sc.scalar(0, unit=None)
226+
return detector.assign_masks({"_monitor_mask": mask})
227+
228+
229+
def _monitor_mask(monitor: sc.DataArray) -> sc.Variable | None:
230+
"""Mask nonfinite monitor values and combine all masks."""
231+
masks = monitor.masks.values()
232+
233+
finite = sc.isfinite(monitor.data)
234+
if not finite.all():
235+
masks = itertools.chain(masks, (~finite,))
236+
237+
mask = None
238+
for m in masks:
239+
if mask is None:
240+
mask = m
241+
else:
242+
mask |= m
243+
244+
return mask

‎tests/normalization_test.py‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
# SPDX-License-Identifier: BSD-3-Clause
22
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
33

4+
import numpy as np
45
import pytest
56
import scipp as sc
67
import scipp.testing
@@ -484,3 +485,149 @@ def test_normalize_by_monitor_integrated_raises_if_monitor_range_too_narrow() ->
484485
monitor=monitor,
485486
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
486487
)
488+
489+
490+
def test_normalize_by_monitor_histogram_monitor_mask_aligned_bins() -> None:
491+
detector = sc.DataArray(
492+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
493+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
494+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
495+
monitor = sc.DataArray(
496+
sc.array(dims=['w'], values=[5.0, 6.0, 7.0], unit='counts'),
497+
coords={'w': sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å')},
498+
masks={'m': sc.array(dims=['w'], values=[False, True, False])},
499+
)
500+
normalized = normalize_by_monitor_histogram(
501+
detector,
502+
monitor=monitor,
503+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
504+
)
505+
506+
expected = (
507+
sc.DataArray(
508+
sc.array(dims=['w'], values=[0.0, 9.6, 0, 216 / 7], unit='counts'),
509+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
510+
)
511+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
512+
.assign_masks(_monitor_mask=sc.array(dims=['w'], values=[False, True, False]))
513+
)
514+
515+
sc.testing.assert_allclose(normalized, expected)
516+
517+
518+
def test_normalize_by_monitor_histogram_monitor_mask_multiple() -> None:
519+
detector = sc.DataArray(
520+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
521+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
522+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
523+
monitor = sc.DataArray(
524+
sc.array(dims=['w'], values=[5.0, 6.0, 7.0], unit='counts'),
525+
coords={'w': sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å')},
526+
masks={
527+
'm1': sc.array(dims=['w'], values=[False, True, False]),
528+
'm2': sc.array(dims=['w'], values=[False, True, True]),
529+
},
530+
)
531+
normalized = normalize_by_monitor_histogram(
532+
detector,
533+
monitor=monitor,
534+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
535+
)
536+
537+
expected = (
538+
sc.DataArray(
539+
sc.array(dims=['w'], values=[0.0, 10, 0, 0], unit='counts'),
540+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
541+
)
542+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
543+
.assign_masks(_monitor_mask=sc.array(dims=['w'], values=[False, True, True]))
544+
)
545+
546+
sc.testing.assert_identical(normalized, expected)
547+
548+
549+
def test_normalize_by_monitor_histogram_monitor_mask_unaligned_bins() -> None:
550+
detector = sc.DataArray(
551+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
552+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
553+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
554+
monitor = sc.DataArray(
555+
sc.array(dims=['w'], values=[5.0, 6.0, 7.0, 8.0], unit='counts'),
556+
coords={'w': sc.array(dims=['w'], values=[0.0, 2, 3.5, 4, 7], unit='Å')},
557+
masks={'m': sc.array(dims=['w'], values=[False, True, False, False])},
558+
)
559+
560+
normalized = normalize_by_monitor_histogram(
561+
detector,
562+
monitor=monitor,
563+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
564+
)
565+
566+
expected = (
567+
sc.DataArray(
568+
sc.array(dims=['w'], values=[0.0, 0, 0, 30], unit='counts'),
569+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
570+
)
571+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
572+
.assign_masks(_monitor_mask=sc.array(dims=['w'], values=[True, True, False]))
573+
)
574+
575+
sc.testing.assert_identical(normalized, expected)
576+
577+
578+
def test_normalize_by_monitor_histogram_monitor_mask_at_edge() -> None:
579+
detector = sc.DataArray(
580+
sc.array(dims=['w'], values=[0, 10, 30], unit='counts'),
581+
coords={'w': sc.arange('w', 3.0, unit='Å')},
582+
).bin(w=sc.array(dims=['w'], values=[0.0, 2, 3], unit='Å'))
583+
monitor = sc.DataArray(
584+
sc.array(dims=['w'], values=[5.0, 6.0], unit='counts'),
585+
coords={'w': sc.array(dims=['w'], values=[0.0, 2, 3], unit='Å')},
586+
masks={'m': sc.array(dims=['w'], values=[False, True])},
587+
)
588+
normalized = normalize_by_monitor_histogram(
589+
detector,
590+
monitor=monitor,
591+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
592+
)
593+
594+
expected = (
595+
sc.DataArray(
596+
sc.array(dims=['w'], values=[0, 10, 0], unit='counts'),
597+
coords={'w': sc.arange('w', 3.0, unit='Å')},
598+
)
599+
.bin(w=sc.array(dims=['w'], values=[0.0, 2, 3], unit='Å'))
600+
.assign_masks(_monitor_mask=sc.array(dims=['w'], values=[False, True]))
601+
)
602+
603+
sc.testing.assert_identical(normalized, expected)
604+
605+
606+
@pytest.mark.parametrize("nonfinite_value", [np.nan, np.inf])
607+
def test_normalize_by_monitor_histogram_nonfinite_in_monitor_is_masked(
608+
nonfinite_value: float,
609+
) -> None:
610+
detector = sc.DataArray(
611+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
612+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
613+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
614+
monitor = sc.DataArray(
615+
sc.array(dims=['w'], values=[nonfinite_value, 6.0, 7.0], unit='counts'),
616+
coords={'w': sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å')},
617+
)
618+
normalized = normalize_by_monitor_histogram(
619+
detector,
620+
monitor=monitor,
621+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
622+
)
623+
624+
expected = (
625+
sc.DataArray(
626+
sc.array(dims=['w'], values=[0.0, 0.0, 65 / 6, 585 / 14], unit='counts'),
627+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
628+
)
629+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
630+
.assign_masks(_monitor_mask=sc.array(dims=['w'], values=[True, False, False]))
631+
)
632+
633+
sc.testing.assert_allclose(normalized, expected)

0 commit comments

Comments
 (0)