Skip to content

Commit 93eff4b

Browse files
jl-wynenjl-wynen
committed
Start using monitor masks
1 parent 489c65e commit 93eff4b

File tree

2 files changed

+240
-3
lines changed

2 files changed

+240
-3
lines changed

‎src/ess/reduce/normalization.py‎

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

5+
import functools
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
*,
1415
monitor: sc.DataArray,
1516
uncertainty_broadcast_mode: UncertaintyBroadcastMode,
17+
combine_monitor_masks: str | None = None,
1618
) -> sc.DataArray:
1719
"""Normalize detector data by a histogrammed monitor.
1820
@@ -71,6 +73,11 @@ def normalize_by_monitor_histogram(
7173
Must be one-dimensional and have a dimension coordinate, typically "wavelength".
7274
uncertainty_broadcast_mode:
7375
Choose how uncertainties of the monitor are broadcast to the sample data.
76+
combine_monitor_masks:
77+
If specified, all masks of the monitor are combined and stored in the
78+
output under this name.
79+
Otherwise, the masks are stored separately, or-ed with
80+
detector masks of the same name.
7481
7582
Returns
7683
-------
@@ -84,14 +91,17 @@ def normalize_by_monitor_histogram(
8491
"""
8592
dim = monitor.dim
8693

94+
detector = _mask_detector_for_norm(
95+
detector=detector, monitor=monitor, combine_monitor_masks=combine_monitor_masks
96+
)
8797
clipped = _clip_monitor_to_detector_range(monitor=monitor, detector=detector)
8898
coord = clipped.coords[dim]
89-
delta_w = coord[1:] - coord[:-1]
99+
delta_w = sc.DataArray(coord[1:] - coord[:-1], masks=clipped.masks)
90100
total_monitor_weight = broadcast_uncertainties(
91101
clipped.sum() / delta_w.sum(),
92102
prototype=clipped,
93103
mode=uncertainty_broadcast_mode,
94-
)
104+
).data
95105
delta_w *= total_monitor_weight
96106
norm = broadcast_uncertainties(
97107
clipped / delta_w, prototype=detector, mode=uncertainty_broadcast_mode
@@ -107,6 +117,7 @@ def normalize_by_monitor_integrated(
107117
*,
108118
monitor: sc.DataArray,
109119
uncertainty_broadcast_mode: UncertaintyBroadcastMode,
120+
combine_monitor_masks: str | None = None,
110121
) -> sc.DataArray:
111122
"""Normalize detector data by an integrated monitor.
112123
@@ -137,6 +148,11 @@ def normalize_by_monitor_integrated(
137148
Must be one-dimensional and have a dimension coordinate, typically "wavelength".
138149
uncertainty_broadcast_mode:
139150
Choose how uncertainties of the monitor are broadcast to the sample data.
151+
combine_monitor_masks:
152+
If specified, all masks of the monitor are combined and stored in the
153+
output under this name.
154+
Otherwise, the masks are stored separately, or-ed with
155+
detector masks of the same name.
140156
141157
Returns
142158
-------
@@ -148,6 +164,9 @@ def normalize_by_monitor_integrated(
148164
normalize_by_monitor_histogram:
149165
Normalize by a monitor histogram without integration.
150166
"""
167+
detector = _mask_detector_for_norm(
168+
detector=detector, monitor=monitor, combine_monitor_masks=combine_monitor_masks
169+
)
151170
clipped = _clip_monitor_to_detector_range(monitor=monitor, detector=detector)
152171
coord = clipped.coords[clipped.dim]
153172
norm = (clipped * (coord[1:] - coord[:-1])).data.sum()
@@ -201,3 +220,32 @@ def _clip_monitor_to_detector_range(
201220
# But integration would pick up all monitor bins.
202221
return monitor.rebin({dim: det_coord})
203222
return monitor[dim, lo:hi]
223+
224+
225+
def _mask_detector_for_norm(
226+
*, detector: sc.DataArray, monitor: sc.DataArray, combine_monitor_masks: str | None
227+
) -> sc.DataArray:
228+
masks = {}
229+
for name, mask in _monitor_masks(monitor, combine_monitor_masks).items():
230+
if name in detector.masks:
231+
masks[name] = detector.masks[name] | mask
232+
else:
233+
masks[name] = mask
234+
return detector.assign_masks(masks)
235+
236+
237+
def _monitor_masks(
238+
monitor: sc.DataArray, combine_masks: str | None
239+
) -> dict[str, sc.Variable]:
240+
if combine_masks is None:
241+
masks = dict(monitor.masks)
242+
else:
243+
masks = {
244+
combine_masks: functools.reduce(lambda a, b: a | b, monitor.masks.values())
245+
}
246+
247+
finite = sc.isfinite(monitor.data)
248+
if not finite.all():
249+
masks["_scipp_nonfinite_monitor"] = ~finite
250+
251+
return masks

‎tests/normalization_test.py‎

Lines changed: 189 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,191 @@ 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(m=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(
544+
m1=sc.array(dims=['w'], values=[False, True, False]),
545+
m2=sc.array(dims=['w'], values=[False, True, True]),
546+
)
547+
)
548+
549+
sc.testing.assert_identical(normalized, expected)
550+
551+
552+
def test_normalize_by_monitor_histogram_monitor_mask_multiple_combine() -> None:
553+
detector = sc.DataArray(
554+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
555+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
556+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
557+
monitor = sc.DataArray(
558+
sc.array(dims=['w'], values=[5.0, 6.0, 7.0], unit='counts'),
559+
coords={'w': sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å')},
560+
masks={
561+
'm1': sc.array(dims=['w'], values=[False, True, False]),
562+
'm2': sc.array(dims=['w'], values=[False, True, True]),
563+
},
564+
)
565+
normalized = normalize_by_monitor_histogram(
566+
detector,
567+
monitor=monitor,
568+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
569+
combine_monitor_masks="MM",
570+
)
571+
572+
expected = (
573+
sc.DataArray(
574+
sc.array(dims=['w'], values=[0.0, 10, 0, 0], unit='counts'),
575+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
576+
)
577+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
578+
.assign_masks(MM=sc.array(dims=['w'], values=[False, True, True]))
579+
)
580+
581+
sc.testing.assert_identical(normalized, expected)
582+
583+
584+
def test_normalize_by_monitor_histogram_monitor_mask_unaligned_bins() -> None:
585+
detector = sc.DataArray(
586+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
587+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
588+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
589+
monitor = sc.DataArray(
590+
sc.array(dims=['w'], values=[5.0, 6.0, 7.0, 8.0], unit='counts'),
591+
coords={'w': sc.array(dims=['w'], values=[0.0, 2, 3.5, 4, 7], unit='Å')},
592+
masks={'m': sc.array(dims=['w'], values=[False, True, False, False])},
593+
)
594+
normalized = normalize_by_monitor_histogram(
595+
detector,
596+
monitor=monitor,
597+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
598+
)
599+
600+
expected = (
601+
sc.DataArray(
602+
sc.array(
603+
dims=['w'], values=[0.0, 100 / 11, 200 / 11, 450 / 11], unit='counts'
604+
),
605+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
606+
)
607+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
608+
.assign_masks(m=sc.array(dims=['w'], values=[False, True, True, False]))
609+
)
610+
611+
sc.testing.assert_identical(normalized, expected)
612+
613+
614+
def test_normalize_by_monitor_histogram_monitor_mask_at_edge() -> None:
615+
detector = sc.DataArray(
616+
sc.array(dims=['w'], values=[0, 10, 30], unit='counts'),
617+
coords={'w': sc.arange('w', 3.0, unit='Å')},
618+
).bin(w=sc.array(dims=['w'], values=[0.0, 2, 3], unit='Å'))
619+
monitor = sc.DataArray(
620+
sc.array(dims=['w'], values=[5.0, 6.0], unit='counts'),
621+
coords={'w': sc.array(dims=['w'], values=[0.0, 2, 3], unit='Å')},
622+
masks={'m': sc.array(dims=['w'], values=[False, True])},
623+
)
624+
normalized = normalize_by_monitor_histogram(
625+
detector,
626+
monitor=monitor,
627+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
628+
)
629+
630+
expected = (
631+
sc.DataArray(
632+
sc.array(dims=['w'], values=[0, 10, 0], unit='counts'),
633+
coords={'w': sc.arange('w', 3.0, unit='Å')},
634+
)
635+
.bin(w=sc.array(dims=['w'], values=[0.0, 2, 3], unit='Å'))
636+
.assign_masks(m=sc.array(dims=['w'], values=[False, True]))
637+
)
638+
639+
sc.testing.assert_identical(normalized, expected)
640+
641+
642+
@pytest.mark.parametrize("nonfinite_value", [np.nan, np.inf])
643+
def test_normalize_by_monitor_histogram_nonfinite_in_monitor_is_masked(
644+
nonfinite_value: float,
645+
) -> None:
646+
detector = sc.DataArray(
647+
sc.array(dims=['w'], values=[0, 10, 20, 30], unit='counts'),
648+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
649+
).bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
650+
monitor = sc.DataArray(
651+
sc.array(dims=['w'], values=[nonfinite_value, 6.0, 7.0], unit='counts'),
652+
coords={'w': sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å')},
653+
)
654+
normalized = normalize_by_monitor_histogram(
655+
detector,
656+
monitor=monitor,
657+
uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
658+
)
659+
660+
expected = (
661+
sc.DataArray(
662+
sc.array(
663+
dims=['w'], values=[np.nan, np.nan, 65 / 6, 585 / 14], unit='counts'
664+
),
665+
coords={'w': sc.arange('w', 1.0, 5.0, unit='Å')},
666+
)
667+
.bin(w=sc.array(dims=['w'], values=[1.0, 3, 4, 7], unit='Å'))
668+
.assign_masks(
669+
_scipp_nonfinite_monitor=sc.array(
670+
dims=['w'], values=[True, True, False, False]
671+
)
672+
)
673+
)
674+
675+
sc.testing.assert_identical(normalized, expected)

0 commit comments

Comments
 (0)