Reconstruct wavelength as bin edges

Author: SimonHeybrock

src/ess/powder/correction.py

@@ -4,6 +4,7 @@
 
 import enum
 from typing import TypeVar
+from uuid import uuid4
 
 import sciline
 import scipp as sc
@@ -58,7 +59,9 @@ def normalize_by_monitor_histogram(
     )
     lut = sc.lookup(norm, dim="wavelength")
     if detector.bins is None:
-        result = detector / lut[detector.coords['wavelength']]
+        result = (
+            detector / lut[sc.midpoints(detector.coords['wavelength'], dim='dspacing')]
+        )
     else:
         result = detector.bins / lut
     return ScaledCountsDspacing[RunType](result)
@@ -104,19 +107,30 @@ def normalize_by_monitor_integrated(
 
     # Clip `monitor` to the range of `detector`, where the bins at the boundary
     # may extend past the detector range (how label-based indexing works).
-    if dim not in detector.coords:
+    if detector.bins is not None and dim in detector.bins.coords:
         det_coord = detector.bins.coords.get(dim)
+        lo = det_coord.nanmin()
+        hi = det_coord.nanmax()
     else:
-        # In this case, we have the wavelength at the bin centers. There is thus some
-        # imprecision in the definition of the wavelength limits, since the detector
-        # bins strictly speaking extend beyond these limits.
-        counts = detector if detector.bins is None else detector.bins.size()
-        flat_counts = counts.flatten(to='dummy')
-        with_counts = flat_counts[flat_counts.data > sc.scalar(0.0, unit=counts.unit)]
-        det_coord = with_counts.coords[dim]
-
-    lo = det_coord.nanmin()
-    hi = det_coord.nanmax()
+        # Mask zero count bins, which are an artifact from the rectangular 2-D binning.
+        # The wavelength of those bins must be excluded when determining the integration
+        # range.
+        counts = (
+            detector.copy(deep=False) if detector.bins is None else detector.bins.size()
+        )
+        counts.masks[uuid4().hex] = counts.data == sc.scalar(0.0, unit=counts.unit)
+        det_coord = detector.coords[dim]
+        edge_dims = {
+            dim: size == det_coord.sizes[dim] + 1 for dim, size in counts.sizes.items()
+        }
+        if len(edge_dims) != 1:
+            raise sc.CoordError(
+                f"Cannot determine edges of coordinate '{dim}' in detector data."
+            )
+        edge_dim = next(iter(edge_dims))
+        lo = counts.assign(det_coord[edge_dim, :-1]).nanmin().data
+        hi = counts.assign(det_coord[edge_dim, 1:]).nanmax().data
+
     if monitor.coords[dim].min() > lo or monitor.coords[dim].max() < hi:
         raise ValueError(
             "Cannot normalize by monitor: The wavelength range of the monitor is "

src/ess/powder/grouping.py

@@ -21,7 +21,7 @@
 def _reconstruct_wavelength(
     dspacing_bins: DspacingBins, two_theta_bins: TwoThetaBins
 ) -> sc.Variable:
-    dspacing = sc.midpoints(dspacing_bins)
+    dspacing = dspacing_bins
     two_theta = sc.midpoints(two_theta_bins)
     return (2 * dspacing * sc.sin(two_theta / 2)).to(unit='angstrom')
 

tests/powder/correction_test.py

@@ -376,8 +376,6 @@ def test_normalize_by_monitor_integrated_expected_results():
         sc.arange('wavelength', 1, 4, unit='counts'),
         coords={'wavelength': sc.arange('wavelength', 3.0, unit='Å')},
     ).bin(wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3], unit='Å'))
-    # "midpoints" at bounds to ensure we include that range.
-    detector.coords['wavelength'] = detector.coords['wavelength'][::2]
     monitor = sc.DataArray(
         sc.array(dims=['wavelength'], values=[4.0, 5.0, 6.0], unit='counts'),
         coords={
@@ -391,19 +389,32 @@ def test_normalize_by_monitor_integrated_expected_results():
         monitor=WavelengthMonitor[SampleRun, CaveMonitor](monitor),
         uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail,
     )
+    # Last event is at 2, so the monitor bin with value 6.0 is not used.
     expected = ScaledCountsDspacing[SampleRun](
-        detector / sc.scalar(4 * 0.5 + 5 * 1.5 + 6 * 1, unit='counts * Å')
+        detector / sc.scalar(4 * 0.5 + 5 * 1.5, unit='counts * Å')
     )
     sc.testing.assert_identical(normalized, expected)
 
 
-def test_normalize_by_monitor_integrated_ignores_monitor_values_out_of_range():
+@pytest.mark.parametrize('event_coord', [True, False])
+def test_normalize_by_monitor_integrated_ignores_monitor_values_out_of_range(
+    event_coord: bool,
+):
     detector = sc.DataArray(
-        sc.arange('wavelength', 3, unit='counts'),
-        coords={'wavelength': sc.arange('wavelength', 3.0, unit='Å')},
-    ).bin(wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3], unit='Å'))
-    # "midpoints" at bounds to ensure we include that range.
-    detector.coords['wavelength'] = detector.coords['wavelength'][::2]
+        sc.arange('wavelength', 4, unit='counts'),
+        coords={'wavelength': sc.arange('wavelength', 4.0, unit='Å')},
+    )
+    if event_coord:
+        # Make sure event at 3 is included
+        detector = detector.bin(
+            wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3.1], unit='Å')
+        )
+        del detector.coords['wavelength']
+    else:
+        detector = detector.bin(
+            wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3], unit='Å')
+        )
+        del detector.bins.coords['wavelength']
     monitor = sc.DataArray(
         sc.array(dims=['wavelength'], values=[4.0, 10.0], unit='counts'),
         coords={
@@ -421,13 +432,25 @@ def test_normalize_by_monitor_integrated_ignores_monitor_values_out_of_range():
     sc.testing.assert_identical(normalized, expected)
 
 
-def test_normalize_by_monitor_integrated_uses_monitor_values_at_boundary():
+@pytest.mark.parametrize('event_coord', [True, False])
+def test_normalize_by_monitor_integrated_uses_monitor_values_at_boundary(
+    event_coord: bool,
+):
     detector = sc.DataArray(
-        sc.arange('wavelength', 3, unit='counts'),
-        coords={'wavelength': sc.arange('wavelength', 3.0, unit='Å')},
-    ).bin(wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3], unit='Å'))
-    # "midpoints" at bounds to ensure we include that range.
-    detector.coords['wavelength'] = detector.coords['wavelength'][::2]
+        sc.arange('wavelength', 4, unit='counts'),
+        coords={'wavelength': sc.arange('wavelength', 4.0, unit='Å')},
+    )
+    if event_coord:
+        # Make sure event at 3 is included
+        detector = detector.bin(
+            wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3.1], unit='Å')
+        )
+        del detector.coords['wavelength']
+    else:
+        detector = detector.bin(
+            wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3], unit='Å')
+        )
+        del detector.bins.coords['wavelength']
     monitor = sc.DataArray(
         sc.array(dims=['wavelength'], values=[4.0, 10.0], unit='counts'),
         coords={
@@ -450,8 +473,6 @@ def test_normalize_by_monitor_integrated_raises_if_monitor_range_too_narrow():
         sc.arange('wavelength', 3, unit='counts'),
         coords={'wavelength': sc.arange('wavelength', 3.0, unit='Å')},
     ).bin(wavelength=sc.array(dims=['wavelength'], values=[0.0, 2, 3], unit='Å'))
-    # "midpoints" at bounds to ensure we include that range.
-    detector.coords['wavelength'] = detector.coords['wavelength'][::2]
     monitor = sc.DataArray(
         sc.array(dims=['wavelength'], values=[4.0, 10.0], unit='counts'),
         coords={