diff --git a/src/ess/beer/__init__.py b/src/ess/beer/__init__.py index fee0b817..829dab28 100644 --- a/src/ess/beer/__init__.py +++ b/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__', diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index 59f6aeae..40c848a4 100644 --- a/src/ess/beer/clustering.py +++ b/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 RawDetector, RunType, StreakClusteredData @@ -10,12 +11,17 @@ def cluster_events_by_streak(da: RawDetector[RunType]) -> StreakClusteredData[Ru 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') diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index be08c8b6..0475b101 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -18,7 +18,7 @@ def compute_tof_in_each_cluster( da: StreakClusteredData[RunType], mod_period: ModulationPeriod, ) -> TofDetector[RunType]: - '''Fits a line through each cluster, the intercept of the line is t0. + """Fits a line through each cluster, the intercept of the line is t0. The line is fitted using linear regression with an outlier removal procedure. The algorithm is: @@ -30,7 +30,7 @@ def compute_tof_in_each_cluster( of the points in the cluster, and probably should belong to another cluster or are part of the background. 3. Go back to 1) and iterate until convergence. A few iterations should be enough. - ''' + """ if isinstance(da, sc.DataGroup): return sc.DataGroup( {k: compute_tof_in_each_cluster(v, mod_period) for k, v in da.items()} @@ -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) @@ -57,10 +60,10 @@ def compute_tof_in_each_cluster( def _linear_regression_by_bin( x: sc.Variable, y: sc.Variable, w: sc.Variable ) -> tuple[sc.Variable, sc.Variable]: - '''Performs a weighted linear regression of the points + """Performs a weighted linear regression of the points in the binned variables ``x`` and ``y`` weighted by ``w``. Returns ``b1`` and ``b0`` such that ``y = b1 * x + b0``. - ''' + """ w = sc.values(w) tot_w = w.bins.sum() @@ -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,7 +90,7 @@ 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) @@ -95,7 +98,7 @@ def _compute_d( 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,28 +115,41 @@ def _compute_d( return d -def _tof_from_dhkl( +def time_of_arrival( event_time_offset: sc.Variable, + tc: sc.Variable, +): + """Does frame unwrapping for pulse shaping chopper modes. + + Events before the "cutoff time" `tc` are assumed to come from the previous pulse.""" + _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, _eto, _eto + T) + + +def _tof_from_dhkl( + time_of_arrival: sc.Variable, theta: sc.Variable, coarse_dhkl: sc.Variable, Ltotal: sc.Variable, mod_period: sc.Variable, time0: sc.Variable, ) -> sc.Variable: - '''Computes tof for BEER given the dhkl peak that the event belongs to''' + """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 @@ -141,7 +157,7 @@ def _tof_from_dhkl( out *= mod_period out += time0 out *= -1 - out += event_time_offset + out += time_of_arrival return out @@ -151,20 +167,23 @@ def tof_from_known_dhkl_graph( time0: WavelengthDefinitionChopperDelay, dhkl_list: DHKLList, ) -> TofCoordTransformGraph: + """Graph computing ``tof`` in modulation chopper modes using + list of peak positions.""" + 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 + """To capture dhkl_list, otherwise it causes an error when ``.transform_coords`` is called unless it is called with ``keep_indermediates=False``. The error happens because data arrays do not allow coordinates 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 +195,56 @@ 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: + """Estimates the time-at-chopper by assuming the wavelength.""" + return ( + sc.constants.m_n + / sc.constants.h + * wavelength_estimate + * (L0 - Ltotal).to(unit=wavelength_estimate.unit) + ).to(unit='s') + + +def _tof_from_t0( + time_of_arrival: sc.Variable, + t0: sc.Variable, +) -> sc.Variable: + """Computes time-of-flight by subtracting a start time.""" + return time_of_arrival - t0 + + +def tof_from_t0_estimate_graph() -> TofCoordTransformGraph: + """Graph for computing ``tof`` in pulse shaping chopper modes.""" + return { + 't0': t0_estimate, + 'tof': _tof_from_t0, + 'time_of_arrival': time_of_arrival, + } + + +def compute_tof( da: RawDetector[RunType], graph: TofCoordTransformGraph ) -> TofDetector[RunType]: + """Uses the transformation graph to compute ``tof``.""" 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_graph, + compute_tof, ) providers = (compute_tof_in_each_cluster,) diff --git a/src/ess/beer/data.py b/src/ess/beer/data.py index 70bd4391..e9e7e15e 100644 --- a/src/ess/beer/data.py +++ b/src/ess/beer/data.py @@ -28,6 +28,15 @@ "silicon-dhkl.tab": "md5:59ee9ed57a7c039ce416c8df886da9cc", "duplex-dhkl.tab": "md5:b4c6c2fcd66466ad291f306b2d6b346e", "dhkl_quartz_nc.tab": "md5:40887d736e3acf859e44488bfd9a9213", + # Simulations from new model with corrected(?) L0. + # For correct reduction you need to use + # beer.io.mcstas_chopper_delay_from_mode_new_simulations + # to obtain the correct WavelengthDefinitionChopperDelay for these files. + "silicon-mode10-new-model.h5": "md5:98500830f27700fc719634e1acd49944", + "silicon-mode16-new-model.h5": "md5:393f9287e7d3f97ceedbe64343918413", + "silicon-mode7-new-model.h5": "md5:d2070d3132722bb551d99b243c62752f", + "silicon-mode8-new-model.h5": "md5:d6dfdf7e87eccedf4f83c67ec552ca22", + "silicon-mode9-new-model.h5": "md5:694a17fb616b7f1c20e94d9da113d201", }, ) @@ -54,6 +63,13 @@ def mcstas_silicon_medium_resolution() -> Path: return _registry.get_path('silicon-mode09.h5') +def mcstas_silicon_new_model(mode: int) -> Path: + """ + Simulated intensity from duplex sample with ``mode`` chopper configuration. + """ + return _registry.get_path(f'silicon-mode{mode}-new-model.h5') + + def duplex_peaks() -> Path: return _registry.get_path('duplex-dhkl.tab') diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index e800e659..96a9bed0 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -4,6 +4,7 @@ import h5py import scipp as sc +import scipp.constants from .types import ( Filename, @@ -43,21 +44,95 @@ def _unique_child_group_h5( return out -def _load_beer_mcstas(f, bank=1): - for key in f['/entry1/instrument/components']: - if 'sampleMantid' in key: - sample_position_path = f'/entry1/instrument/components/{key}/Position' - break +def _wavelength_estimate_from_mode(mode, value_in_file): + if value_in_file != '0': + return sc.scalar(float(value_in_file), unit='angstrom') + if mode in [ + '0', + '3', + '4', + '5', + '6', + '7', + '8', + '9', + '10', + '15', + '16', + ]: + return sc.scalar(2.1, unit='angstrom') + elif mode in ['1', '2']: + return sc.scalar(3.1, unit='angstrom') + elif mode == '11': + return sc.scalar(3.0, unit='angstrom') + elif mode == '12': + return sc.scalar(3.5, unit='angstrom') + elif mode == '13': + return sc.scalar(6.0, unit='angstrom') + elif mode == '14': + return sc.scalar(4.0, unit='angstrom') + else: + raise ValueError(f'Unkonwn chopper mode {mode}.') + + +def _chopper_position_from_mode( + mode, + *, + psc1_pos, + psc2_pos, + psc3_pos, + mca_pos, + mcb_pos, + mcc_pos, +): + if mode in ['0', '1', '2', '11']: + return sc.vector([0.0, 0.0, 0.0], unit='m') + elif mode in ['3', '4', '12', '13', '15']: + return sc.vector(0.5 * (psc1_pos + psc3_pos), unit='m') + elif mode in ['5', '6']: + return sc.vector(0.5 * (psc1_pos + psc2_pos), unit='m') + elif mode in ['7', '8', '9', '10']: + return sc.vector(mca_pos, unit='m') + elif mode == '14': + return sc.vector(mcc_pos, unit='m') + elif mode == '16': + return sc.vector(0.5 * (mca_pos + mcb_pos), unit='m') else: - raise ValueError('Sample position entry not found in file.') - data, events, params, sample_pos, chopper_pos = _load_h5( + raise ValueError(f'Unkonwn chopper mode {mode}.') + + +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', 'MCB', 'MCC'] + if name in key + } + ( + 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', - sample_position_path, - '/entry1/instrument/components/0017_cMCA/Position', + positions['sampleMantid'], + positions['PSC1'], + positions['PSC2'], + positions['PSC3'], + positions['MCA'], + positions['MCB'], + positions['MCC'], ) + events = events[()] da = sc.DataArray( sc.array(dims=['events'], values=events[:, 0], variances=events[:, 0] ** 2), @@ -75,14 +150,29 @@ def _load_beer_mcstas(f, bank=1): v = v[0] if isinstance(v, bytes): v = v.decode() - if k in ('mode', 'sample_filename'): + if k in ('mode', 'sample_filename', 'lambda'): da.coords[k] = sc.scalar(v) + da.coords['wavelength_estimate'] = _wavelength_estimate_from_mode( + da.coords['mode'].value, + da.coords.pop('lambda').value if 'lambda' in da.coords else '0', + ) + + da.coords['chopper_position'] = _chopper_position_from_mode( + da.coords['mode'].value, + psc1_pos=psc1_pos[:], + psc2_pos=psc2_pos[:], + psc3_pos=psc3_pos[:], + mca_pos=mca_pos[:], + mcb_pos=mcb_pos[:], + mcc_pos=mcc_pos[:], + ) + 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' ) - da.coords['chopper_position'] = sc.vector(chopper_pos[:], unit='m') + da.coords['x'].unit = 'm' da.coords['y'].unit = 'm' da.coords['t'].unit = 's' @@ -90,6 +180,7 @@ def _load_beer_mcstas(f, bank=1): z = sc.norm(da.coords['detector_position'] - da.coords['sample_position']) L1 = sc.norm(da.coords['sample_position'] - da.coords['chopper_position']) L2 = sc.sqrt(da.coords['x'] ** 2 + da.coords['y'] ** 2 + z**2) + # Source is assumed to be at the origin da.coords['L0'] = L1 + L2 + sc.norm(da.coords['chopper_position']) da.coords['Ltotal'] = L1 + L2 @@ -102,7 +193,15 @@ def _load_beer_mcstas(f, bank=1): da.coords.pop('y') da.coords.pop('n') - da.coords['event_time_offset'] = da.coords.pop('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['wavelength_estimate'] + * da.coords['L0'].min().to(unit='angstrom') + ).to(unit='s') - sc.scalar(1 / 14, unit='s') / 2 + return da @@ -150,16 +249,37 @@ def load_beer_mcstas_provider( def mcstas_chopper_delay_from_mode( da: RawDetector[SampleRun], ) -> WavelengthDefinitionChopperDelay: + '''These settings are good for the set of McStas runs that we + use in the docs currently. + Eventually we will want to determine this from the chopper information + in the files, but that information is not in the simulation output.''' 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.') +def mcstas_chopper_delay_from_mode_new_simulations( + da: RawDetector[SampleRun], +) -> WavelengthDefinitionChopperDelay: + '''Celine has a new simulation with some changes to the chopper placement(?). + For those simulations we need to adapt the chopper delay values.''' + mode = next(iter(d.coords['mode'] for d in da.values())).value + if mode == '7': + return sc.scalar(0.001370158730158727, unit='s') + if mode == '8': + return sc.scalar(0.001370158730158727, unit='s') + if mode == '9': + return sc.scalar(0.0022630158730158727, unit='s') + if mode == '10': + return sc.scalar(0.0022630158730158727, unit='s') + if mode == '16': + return sc.scalar(0.000476984126984127, unit='s') + raise ValueError(f'Mode {mode} is not known.') + + def mcstas_modulation_period_from_mode(da: RawDetector[SampleRun]) -> ModulationPeriod: mode = next(iter(d.coords['mode'] for d in da.values())).value if mode in ('7', '8'): diff --git a/src/ess/beer/types.py b/src/ess/beer/types.py index dd90f586..d0caafb1 100644 --- a/src/ess/beer/types.py +++ b/src/ess/beer/types.py @@ -31,17 +31,17 @@ class StreakClusteredData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): TofCoordTransformGraph = NewType("TofCoordTransformGraph", dict) PulseLength = NewType('PulseLength', sc.Variable) -'''Length of the neutron source pulse in time.''' +"""Length of the neutron source pulse in time.""" ModulationPeriod = NewType('ModulationPeriod', sc.Variable) -'''The effective period of the modulating chopper: +"""The effective period of the modulating chopper: ``1 / (K * F)`` where ``K`` is the number of chopper openings and -``F`` is the chopper frequency.''' +``F`` is the chopper frequency.""" WavelengthDefinitionChopperDelay = NewType( 'WavelengthDefinitionChopperDelay', sc.Variable ) -'''Wavelength definition chopper time delay relative to source pulse.''' +"""Wavelength definition chopper time delay relative to source pulse.""" DHKLList = NewType('DHKLList', sc.Variable) -'''List of peak position estimates.''' +"""List of peak position estimates.""" diff --git a/src/ess/beer/workflow.py b/src/ess/beer/workflow.py index f47fdcb0..6545267b 100644 --- a/src/ess/beer/workflow.py +++ b/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 ( @@ -24,8 +24,8 @@ def BeerModMcStasWorkflow(): - '''Workflow to process BEER (modulation regime) McStas files without a list - of estimated peak positions.''' + """Workflow to process BEER (modulation regime) McStas files without a list + of estimated peak positions.""" return sl.Pipeline( (*mcstas_providers, *clustering_providers, *conversion_providers), params=default_parameters, @@ -34,10 +34,19 @@ def BeerModMcStasWorkflow(): def BeerModMcStasWorkflowKnownPeaks(): - '''Workflow to process BEER (modulation regime) McStas files using a list - of estimated peak positions.''' + """Workflow to process BEER (modulation regime) McStas files using a list + of estimated peak positions.""" return sl.Pipeline( (*mcstas_providers, *convert_from_known_peaks_providers), 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,)}, + )