fix: load chopper settings in all modes - truncate time to limits of event_time_offset

Author: jokasimr

src/ess/beer/__init__.py

@@ -9,6 +9,7 @@
 
 from .io import load_beer_mcstas
 from .workflow import (
+    BeerMcStasWorkflowPulseShaping,
     BeerModMcStasWorkflow,
     BeerModMcStasWorkflowKnownPeaks,
     default_parameters,
@@ -22,6 +23,7 @@
 del importlib
 
 __all__ = [
+    'BeerMcStasWorkflowPulseShaping',
     'BeerModMcStasWorkflow',
     'BeerModMcStasWorkflowKnownPeaks',
     '__version__',

src/ess/beer/clustering.py

@@ -2,6 +2,7 @@
 from scippneutron.conversion.tof import dspacing_from_tof
 from scipy.signal import find_peaks, medfilt
 
+from .conversions import _t0_estimate, _time_of_arrival
 from .types import DetectorData, RunType, StreakClusteredData
 
 
@@ -10,12 +11,17 @@ def cluster_events_by_streak(da: DetectorData[RunType]) -> StreakClusteredData[R
         return sc.DataGroup({k: cluster_events_by_streak(v) for k, v in da.items()})
     da = da.copy(deep=False)
 
-    # TODO: approximate t0 should depend on chopper information
-    approximate_t0 = sc.scalar(6e-3, unit='s')
+    t = _time_of_arrival(
+        da.coords['event_time_offset'],
+        da.coords['tc'].to(unit=da.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(
-        tof=da.coords['event_time_offset'] - approximate_t0,
-        Ltotal=da.coords['L0'],
+        tof=t - approximate_t0,
+        Ltotal=da.coords['Ltotal'],
         two_theta=da.coords['two_theta'],
     ).to(unit='angstrom')
 

src/ess/beer/conversions.py

@@ -38,7 +38,10 @@ def compute_tof_in_each_cluster(
 
     max_distance_from_streak_line = mod_period / 3
     sin_theta_L = sc.sin(da.bins.coords['two_theta'] / 2) * da.bins.coords['Ltotal']
-    t = da.bins.coords['event_time_offset']
+    t = _time_of_arrival(
+        da.bins.coords['event_time_offset'],
+        da.coords['tc'].to(unit=da.bins.coords['event_time_offset'].unit),
+    )
     for _ in range(15):
         s, t0 = _linear_regression_by_bin(sin_theta_L, t, da.data)
 
@@ -77,7 +80,7 @@ def _linear_regression_by_bin(
 
 
 def _compute_d(
-    event_time_offset: sc.Variable,
+    time_of_arrival: sc.Variable,
     theta: sc.Variable,
     dhkl_list: sc.Variable,
     pulse_length: sc.Variable,
@@ -87,15 +90,15 @@ def _compute_d(
     given a list of known peaks."""
     # Source: https://www2.mcstas.org/download/components/3.4/contrib/NPI_tof_dhkl_detector.comp
     sinth = sc.sin(theta)
-    t = event_time_offset
+    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)
 
     const = (2 * sinth * L0 / (scipp.constants.h / scipp.constants.m_n)).to(
-        unit=f'{event_time_offset.unit}/angstrom'
+        unit=f'{time_of_arrival.unit}/angstrom'
     )
 
     for dhkl in dhkl_list:
@@ -112,8 +115,18 @@ def _compute_d(
     return d
 
 
-def _tof_from_dhkl(
+def _time_of_arrival(
     event_time_offset: sc.Variable,
+    tc: sc.Variable,
+):
+    _eto = event_time_offset
+    T = sc.scalar(1 / 14, unit='s').to(unit=_eto.unit)
+    tc = tc.to(unit=_eto.unit)
+    return sc.where(_eto >= tc % T, _eto, _eto + T)
+
+
+def _tof_from_dhkl(
+    time_of_arrival: sc.Variable,
     theta: sc.Variable,
     coarse_dhkl: sc.Variable,
     Ltotal: sc.Variable,
@@ -123,25 +136,25 @@ def _tof_from_dhkl(
     '''Computes tof for BEER given the dhkl peak that the event belongs to'''
     # Source: https://www2.mcstas.org/download/components/3.4/contrib/NPI_tof_dhkl_detector.comp
     # tref = 2 * d_hkl * sin(theta) / hm * Ltotal
-    # tc = event_time_zero - time0 - tref
+    # tc = time_of_arrival - time0 - tref
     # dt = floor(tc / mod_period + 0.5) * mod_period + time0
-    # tof = event_time_offset - dt
+    # tof = time_of_arrival - dt
     c = (-2 * 1.0 / (scipp.constants.h / scipp.constants.m_n)).to(
-        unit=f'{event_time_offset.unit}/m/angstrom'
+        unit=f'{time_of_arrival.unit}/m/angstrom'
     )
     out = sc.sin(theta)
     out *= c
     out *= coarse_dhkl
     out *= Ltotal
-    out += event_time_offset
+    out += time_of_arrival
     out -= time0
     out /= mod_period
     out += 0.5
     sc.floor(out, out=out)
     out *= mod_period
     out += time0
     out *= -1
-    out += event_time_offset
+    out += time_of_arrival
     return out
 
 
@@ -152,10 +165,10 @@ def tof_from_known_dhkl_graph(
     dhkl_list: DHKLList,
 ) -> TofCoordTransformGraph:
     def _compute_coarse_dspacing(
-        event_time_offset,
+        time_of_arrival: sc.Variable,
         theta: sc.Variable,
         pulse_length: sc.Variable,
-        L0,
+        L0: sc.Variable,
     ):
         '''To capture dhkl_list, otherwise it causes an error when
         ``.transform_coords`` is called unless it is called with
@@ -164,7 +177,7 @@ def _compute_coarse_dspacing(
         with dimensions not present on the data.
         '''
         return _compute_d(
-            event_time_offset=event_time_offset,
+            time_of_arrival=time_of_arrival,
             theta=theta,
             pulse_length=pulse_length,
             L0=L0,
@@ -176,19 +189,52 @@ def _compute_coarse_dspacing(
         'mod_period': lambda: mod_period,
         'time0': lambda: time0,
         'tof': _tof_from_dhkl,
+        'time_of_arrival': _time_of_arrival,
         'coarse_dhkl': _compute_coarse_dspacing,
         'theta': lambda two_theta: two_theta / 2,
     }
 
 
-def compute_tof_from_known_peaks(
+def _t0_estimate(
+    wavelength_estimate: sc.Variable,
+    L0: sc.Variable,
+    Ltotal: sc.Variable,
+) -> sc.Variable:
+    return (
+        sc.constants.m_n
+        / sc.constants.h
+        * wavelength_estimate
+        * (L0 - Ltotal).to(unit=wavelength_estimate.unit)
+    ).to(unit='s')
+
+
+def _tof_from_t0_estimate(
+    time_of_arrival: sc.Variable,
+    t0_estimate: sc.Variable,
+) -> sc.Variable:
+    return time_of_arrival - t0_estimate
+
+
+def tof_from_t0_estimate() -> TofCoordTransformGraph:
+    return {
+        't0_estimate': _t0_estimate,
+        'tof': _tof_from_t0_estimate,
+        'time_of_arrival': _time_of_arrival,
+    }
+
+
+def compute_tof(
     da: DetectorData[RunType], graph: TofCoordTransformGraph
 ) -> DetectorTofData[RunType]:
     return da.transform_coords(('tof',), graph=graph)
 
 
 convert_from_known_peaks_providers = (
     tof_from_known_dhkl_graph,
-    compute_tof_from_known_peaks,
+    compute_tof,
+)
+convert_pulse_shaping = (
+    tof_from_t0_estimate,
+    compute_tof,
 )
 providers = (compute_tof_in_each_cluster,)

src/ess/beer/io.py

@@ -48,22 +48,32 @@ def _load_beer_mcstas(f, bank=1):
     positions = {
         name: f'/entry1/instrument/components/{key}/Position'
         for key in f['/entry1/instrument/components']
-        for name in ['sampleMantid', 'PSC1', 'PSC2', 'PSC3', 'MCA', 'MCC']
+        for name in ['sampleMantid', 'PSC1', 'PSC2', 'PSC3', 'MCA', 'MCB', 'MCC']
         if name in key
     }
-    data, events, params, sample_pos, psc1_pos, psc2_pos, psc3_pos, mca_pos, mcc_pos = (
-        _load_h5(
-            f,
-            f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t',
-            f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t/events',
-            'NXentry/simulation/Param',
-            positions['sampleMantid'],
-            positions['PSC1'],
-            positions['PSC2'],
-            positions['PSC3'],
-            positions['MCA'],
-            positions['MCC'],
-        )
+    (
+        data,
+        events,
+        params,
+        sample_pos,
+        psc1_pos,
+        psc2_pos,
+        psc3_pos,
+        mca_pos,
+        mcb_pos,
+        mcc_pos,
+    ) = _load_h5(
+        f,
+        f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t',
+        f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t/events',
+        'NXentry/simulation/Param',
+        positions['sampleMantid'],
+        positions['PSC1'],
+        positions['PSC2'],
+        positions['PSC3'],
+        positions['MCA'],
+        positions['MCB'],
+        positions['MCC'],
     )
 
     events = events[()]
@@ -86,7 +96,10 @@ def _load_beer_mcstas(f, bank=1):
         if k in ('mode', 'sample_filename', 'lambda'):
             da.coords[k] = sc.scalar(v)
 
-    if da.coords['lambda'].value == '0':
+    if 'lambda' in da.coords:
+        da.coords['wavelength_estimate'] = da.coords.pop('lambda')
+
+    if da.coords['wavelength_estimate'].value == '0':
         if da.coords['mode'].value in [
             '0',
             '3',
@@ -100,28 +113,28 @@ def _load_beer_mcstas(f, bank=1):
             '15',
             '16',
         ]:
-            da.coords['lambda'] = sc.scalar(2.1, unit='angstrom')
+            da.coords['wavelength_estimate'] = sc.scalar(2.1, unit='angstrom')
         elif da.coords['mode'].value in ['1', '2']:
-            da.coords['lambda'] = sc.scalar(3.1, unit='angstrom')
+            da.coords['wavelength_estimate'] = sc.scalar(3.1, unit='angstrom')
         elif da.coords['mode'].value == '11':
-            da.coords['lambda'] = sc.scalar(3.0, unit='angstrom')
+            da.coords['wavelength_estimate'] = sc.scalar(3.0, unit='angstrom')
         elif da.coords['mode'].value == '12':
-            da.coords['lambda'] = sc.scalar(3.5, unit='angstrom')
+            da.coords['wavelength_estimate'] = sc.scalar(3.5, unit='angstrom')
         elif da.coords['mode'].value == '13':
-            da.coords['lambda'] = sc.scalar(6.0, unit='angstrom')
+            da.coords['wavelength_estimate'] = sc.scalar(6.0, unit='angstrom')
         elif da.coords['mode'].value == '14':
-            da.coords['lambda'] = sc.scalar(4.0, unit='angstrom')
+            da.coords['wavelength_estimate'] = sc.scalar(4.0, unit='angstrom')
     else:
-        da.coords['lambda'] = sc.scalar(
-            float(da.coords['lambda'].value), unit='angstrom'
+        da.coords['wavelength_estimate'] = sc.scalar(
+            float(da.coords['wavelength_estimate'].value), unit='angstrom'
         )
 
     da.coords['sample_position'] = sc.vector(sample_pos[:], unit='m')
     da.coords['detector_position'] = sc.vector(
         list(map(float, da.coords.pop('position').value.split(' '))), unit='m'
     )
 
-    if da.coords['mode'].value in ['0', '1', '2', '11', '16']:
+    if da.coords['mode'].value in ['0', '1', '2', '11']:
         da.coords['chopper_position'] = sc.vector([0.0, 0.0, 0.0], unit='m')
     elif da.coords['mode'].value in ['3', '4', '12', '13', '15']:
         da.coords['chopper_position'] = sc.vector(
@@ -132,9 +145,13 @@ def _load_beer_mcstas(f, bank=1):
             0.5 * (psc1_pos[:] + psc2_pos[:]), unit='m'
         )
     elif da.coords['mode'].value in ['7', '8', '9', '10']:
-        da.coords['chopper_position'] = sc.vector(0.5 * mca_pos[:], unit='m')
+        da.coords['chopper_position'] = sc.vector(mca_pos[:], unit='m')
     elif da.coords['mode'].value == '14':
         da.coords['chopper_position'] = sc.vector(mcc_pos[:], unit='m')
+    elif da.coords['mode'].value == '16':
+        da.coords['chopper_position'] = sc.vector(
+            0.5 * (mca_pos[:] + mcb_pos[:]), unit='m'
+        )
     else:
         raise ValueError(f'Unkonwn chopper mode {da.coords["mode"].value}.')
 
@@ -143,7 +160,6 @@ def _load_beer_mcstas(f, bank=1):
     da.coords['t'].unit = 's'
 
     z = sc.norm(da.coords['detector_position'] - da.coords['sample_position'])
-    # z = da.coords['position'].fields.z - da.coords['sample_position'].fields.z
     L1 = sc.norm(da.coords['sample_position'] - da.coords['chopper_position'])
     L2 = sc.sqrt(da.coords['x'] ** 2 + da.coords['y'] ** 2 + z**2)
 
@@ -159,17 +175,15 @@ def _load_beer_mcstas(f, bank=1):
     da.coords.pop('y')
     da.coords.pop('n')
 
-    # expression of temporal offset delta_t checked for pulse shaping mode: 4,5,6
-    delta_t = (
+    t = da.coords.pop('t')
+    da.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
-        * da.coords['lambda']
-        * sc.norm(da.coords['chopper_position']).to(unit='angstrom')
-    ).to(unit='s')
+        * da.coords['wavelength_estimate']
+        * da.coords['L0'].min().to(unit='angstrom')
+    ).to(unit='s') - sc.scalar(1 / 14, unit='s') / 2
 
-    t = da.coords.pop('t')
-    da.coords['event_time_offset'] = t % sc.scalar(1 / 14, unit=t.unit)
-    da.coords["approximate_tof"] = da.coords['event_time_offset'] - delta_t
     return da
 
 
@@ -218,10 +232,8 @@ def mcstas_chopper_delay_from_mode(
     da: DetectorData[SampleRun],
 ) -> WavelengthDefinitionChopperDelay:
     mode = next(iter(d.coords['mode'] for d in da.values())).value
-    if mode in ('7', '8'):
-        return sc.scalar(0.00245635, unit='s')
-    if mode in ('9', '10'):
-        return sc.scalar(0.0033730158730158727, unit='s')
+    if mode in ('7', '8', '9', '10'):
+        return sc.scalar(0.0024730158730158727, unit='s')
     if mode == '16':
         return sc.scalar(0.000876984126984127, unit='s')
     raise ValueError(f'Mode {mode} is not known.')

src/ess/beer/workflow.py

@@ -4,7 +4,7 @@
 import scipp as sc
 
 from .clustering import providers as clustering_providers
-from .conversions import convert_from_known_peaks_providers
+from .conversions import convert_from_known_peaks_providers, convert_pulse_shaping
 from .conversions import providers as conversion_providers
 from .io import mcstas_providers
 from .types import (
@@ -41,3 +41,12 @@ def BeerModMcStasWorkflowKnownPeaks():
         params=default_parameters,
         constraints={RunType: (SampleRun,)},
     )
+
+
+def BeerMcStasWorkflowPulseShaping():
+    '''Workflow to process BEER (pulse shaping modes) McStas files'''
+    return sl.Pipeline(
+        (*mcstas_providers, *convert_pulse_shaping),
+        params=default_parameters,
+        constraints={RunType: (SampleRun,)},
+    )