Bin by pixel

src/ess/beer/clustering.py

@@ -12,20 +12,25 @@ def cluster_events_by_streak(da: RawDetector[RunType]) -> StreakClusteredData[Ru
     da = da.copy(deep=False)
 
     t = time_of_arrival(
-        da.coords['event_time_offset'],
-        da.coords['tc'].to(unit=da.coords['event_time_offset'].unit),
+        da.bins.coords['event_time_offset'],
+        da.coords['tc'].to(unit=da.bins.coords['event_time_offset'].unit),
     )
     approximate_t0 = t0_estimate(
         da.coords['wavelength_estimate'], da.coords['L0'], da.coords['Ltotal']
     ).to(unit=t.unit)
 
-    da.coords['coarse_d'] = dspacing_from_tof(
+    da.bins.coords['coarse_d'] = dspacing_from_tof(
         tof=t - approximate_t0,
         Ltotal=da.coords['Ltotal'],
         two_theta=da.coords['two_theta'],
     ).to(unit='angstrom')
 
-    h = da.hist(coarse_d=1000)
+    # We need to keep these coordinates after binning,
+    # adding them to the binned data coords achieves this.
+    for coord in ('two_theta', 'Ltotal'):
+        da.bins.coords[coord] = sc.bins_like(da, da.coords[coord])
+
+    h = da.bins.concat().hist(coarse_d=1000)
     i_peaks, _ = find_peaks(
         h.data.values, height=medfilt(h.values, kernel_size=99), distance=3
     )
@@ -52,8 +57,10 @@ def cluster_events_by_streak(da: RawDetector[RunType]) -> StreakClusteredData[Ru
         )
     ]
     has_peak = peaks.bin(coarse_d=filtered_valleys).bins.size().data
-    b = da.bin(coarse_d=filtered_valleys).assign_masks(
-        no_peak=has_peak != sc.scalar(1, unit=None)
+    b = (
+        da.bins.concat()
+        .bin(coarse_d=filtered_valleys)
+        .assign_masks(no_peak=has_peak != sc.scalar(1, unit=None))
     )
     b = b.drop_coords(('coarse_d',))
     b = b.bins.drop_coords(('coarse_d',))

src/ess/beer/conversions.py

@@ -92,10 +92,10 @@ def _compute_d(
     sinth = sc.sin(theta)
     t = time_of_arrival
 
-    d = sc.empty(dims=sinth.dims, shape=sinth.shape, unit=dhkl_list[0].unit)
-    d[:] = sc.scalar(float('nan'), unit=dhkl_list[0].unit)
-    dtfound = sc.empty(dims=sinth.dims, shape=sinth.shape, dtype='float64', unit=t.unit)
-    dtfound[:] = sc.scalar(float('nan'), unit=t.unit)
+    d = sc.full_like(
+        time_of_arrival, value=float('nan'), unit=dhkl_list[0].unit, dtype='float64'
+    )
+    dtfound = sc.full_like(time_of_arrival, value=float('nan'), dtype='float64')
 
     const = (2 * sinth * L0 / (scipp.constants.h / scipp.constants.m_n)).to(
         unit=f'{time_of_arrival.unit}/angstrom'
@@ -145,9 +145,8 @@ def _tof_from_dhkl(
     c = (-2 * 1.0 / (scipp.constants.h / scipp.constants.m_n)).to(
         unit=f'{time_of_arrival.unit}/m/angstrom'
     )
-    out = sc.sin(theta)
-    out *= c
-    out *= coarse_dhkl
+    out = c * coarse_dhkl
+    out *= sc.sin(theta)
     out *= Ltotal
     out += time_of_arrival
     out -= time0

src/ess/beer/io.py

@@ -183,28 +183,38 @@ def _load_beer_mcstas(f, bank=1):
         list(map(float, da.coords.pop('position').value.split(' '))), unit='m'
     )
 
+    da.coords.pop('n')
     da.coords['x'].unit = 'm'
     da.coords['y'].unit = 'm'
     da.coords['t'].unit = 's'
 
-    z = sc.norm(da.coords['detector_position'] - da.coords['sample_position'])
+    da = da.bin(
+        y=sc.linspace('y', -0.5, 0.5, 500, unit='m'),
+        x=sc.linspace('x', -0.5, 0.5, 500, unit='m'),
+    )
+    da.coords['position'] = (
+        da.coords['detector_position']
+        + sc.spatial.as_vectors(
+            sc.scalar(0.0, unit='m'),
+            sc.midpoints(da.coords['y']),
+            sc.midpoints(-da.coords['x'] if bank == 1 else da.coords['x']),
+        )
+    ).transpose(da.dims)
+
     L1 = sc.norm(da.coords['sample_position'] - da.coords['chopper_position'])
-    L2 = sc.sqrt(da.coords['x'] ** 2 + da.coords['y'] ** 2 + z**2)
+    L2 = sc.norm(da.coords['position'] - da.coords['sample_position'])
 
     # Source is assumed to be at the origin
     da.coords['L0'] = L1 + L2 + sc.norm(da.coords['chopper_position'])
     da.coords['Ltotal'] = L1 + L2
     da.coords['two_theta'] = sc.acos(
-        (-da.coords['x'] if bank == 1 else da.coords['x']) / L2
+        (da.coords['position'] - da.coords['sample_position']).fields.z / L2
     )
 
-    # Save some space
-    da.coords.pop('x')
-    da.coords.pop('y')
-    da.coords.pop('n')
-
-    t = da.coords.pop('t')
-    da.coords['event_time_offset'] = t % sc.scalar(1 / 14, unit='s').to(unit=t.unit)
+    t = da.bins.coords['t']
+    da.bins.coords['event_time_offset'] = t % sc.scalar(1 / 14, unit='s').to(
+        unit=t.unit
+    )
     da.coords["tc"] = (
         sc.constants.m_n
         / sc.constants.h

tests/beer/mcstas_reduction_test.py

@@ -22,7 +22,7 @@ def test_can_reduce_using_known_peaks_workflow():
     assert 'bank1' in da
     assert 'bank2' in da
     da = da['bank1']
-    assert 'tof' in da.coords
+    assert 'tof' in da.bins.coords
     # assert dataarray has all coords required to compute dspacing
     da = da.transform_coords(
         ('dspacing',),
@@ -64,7 +64,7 @@ def test_pulse_shaping_workflow():
     assert 'bank1' in da
     assert 'bank2' in da
     da = da['bank1']
-    assert 'tof' in da.coords
+    assert 'tof' in da.bins.coords
     # assert dataarray has all coords required to compute dspacing
     da = da.transform_coords(
         ('dspacing',),