From ab04d1d826fe6c08a06f6dc2057ae17214b41a0d Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Mon, 6 May 2024 16:35:30 +0200 Subject: [PATCH 01/26] add types for masking, fix types in pooch paths and add initial providers for masking in pixels, wavelength and theta --- src/ess/powder/__init__.py | 14 +++- src/ess/powder/conversion.py | 96 ++++++++++++++++++++------ src/ess/powder/external/powgen/data.py | 86 +++++++++++++---------- src/ess/powder/filtering.py | 28 ++++---- src/ess/powder/masking.py | 88 +++++++++++++++++++++++ src/ess/powder/types.py | 93 +++++++++++++++++++------ 6 files changed, 310 insertions(+), 95 deletions(-) create mode 100644 src/ess/powder/masking.py diff --git a/src/ess/powder/__init__.py b/src/ess/powder/__init__.py index e2c69c75..bdf9af1c 100644 --- a/src/ess/powder/__init__.py +++ b/src/ess/powder/__init__.py @@ -1,13 +1,20 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -# flake8: noqa: F401 """ Components for powder diffraction experiments. """ import importlib.metadata -from . import conversion, correction, filtering, grouping, smoothing, uncertainty +from . import ( + conversion, + correction, + filtering, + grouping, + masking, + smoothing, + uncertainty, +) try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -17,9 +24,10 @@ del importlib providers = ( - *conversion.providers, + *conversion.providers_with_calibration, *correction.providers, *filtering.providers, *grouping.providers, + *masking.providers, ) """Sciline providers for powder diffraction.""" diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index 4ab63ee3..e0e222a4 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -14,10 +14,16 @@ from .types import ( CalibrationData, DspacingData, + ElasticCoordTransformGraph, NormalizedByProtonCharge, + PixelMaskedData, RawSample, RawSource, RunType, + TwoThetaData, + TwoThetaMaskedData, + WavelengthData, + WavelengthMaskedData, ) @@ -130,26 +136,67 @@ def to_dspacing_with_calibration( out = _restore_tof_if_in_wavelength(out) graph = { - 'dspacing': _dspacing_from_diff_calibration, + "dspacing": _dspacing_from_diff_calibration, } # `_dspacing_from_diff_calibration` does not need positions but conceptually, # the conversion maps from positions to d-spacing. # The mechanism with `_tag_positions_consumed` is meant to ensure that, # if positions are present, they are consumed (mad unaligned or dropped) # by the coordinate transform similarly to `to_dspacing_with_positions`. - if 'position' in out.coords or ( - out.bins is not None and 'position' in out.bins.coords + if "position" in out.coords or ( + out.bins is not None and "position" in out.bins.coords ): - graph['_tag_positions_consumed'] = _consume_positions + graph["_tag_positions_consumed"] = _consume_positions else: - graph['_tag_positions_consumed'] = lambda: sc.scalar(0) - out = out.transform_coords('dspacing', graph=graph, keep_intermediate=False) - out.coords.pop('_tag_positions_consumed', None) + graph["_tag_positions_consumed"] = lambda: sc.scalar(0) + out = out.transform_coords("dspacing", graph=graph, keep_intermediate=False) + out.coords.pop("_tag_positions_consumed", None) return DspacingData[RunType](out) +def powder_coordinate_transformation_graph() -> ElasticCoordTransformGraph: + """ + Generate a coordinate transformation graph for powder diffraction. + + Returns + ------- + : + A dictionary with the graph for the transformation. + """ + return ElasticCoordTransformGraph( + { + **scn.conversion.graph.beamline.beamline(scatter=True), + **scn.conversion.graph.tof.elastic("tof"), + } + ) + + +def to_wavelength_with_positions( + data: PixelMaskedData[RunType], + graph: ElasticCoordTransformGraph, +) -> WavelengthData[RunType]: + """ + Transform coordinates to wavelength using detector positions. + """ + return WavelengthData[RunType]( + data.transform_coords("wavelength", graph=graph, keep_intermediate=False) + ) + + +def to_twotheta_with_positions( + data: WavelengthMaskedData[RunType], + graph: ElasticCoordTransformGraph, +) -> TwoThetaData[RunType]: + """ + Transform coordinates to two-theta using detector positions. + """ + return TwoThetaData[RunType]( + data.transform_coords("two_theta", graph=graph, keep_intermediate=False) + ) + + def to_dspacing_with_positions( - data: NormalizedByProtonCharge[RunType], + data: TwoThetaMaskedData[RunType], *, sample: Optional[RawSample[RunType]] = None, source: Optional[RawSource] = None, @@ -190,42 +237,49 @@ def to_dspacing_with_positions( """ graph = { **scn.conversion.graph.beamline.beamline(scatter=True), - **scn.conversion.graph.tof.elastic_dspacing('tof'), + **scn.conversion.graph.tof.elastic_dspacing("tof"), } if sample is not None: - graph['sample_position'] = lambda: sample['position'] + graph["sample_position"] = lambda: sample["position"] if source is not None: - graph['source_position'] = lambda: source['position'] + graph["source_position"] = lambda: source["position"] out = _restore_tof_if_in_wavelength(data) - out = out.transform_coords('dspacing', graph=graph, keep_intermediate=False) + out = out.transform_coords("dspacing", graph=graph, keep_intermediate=False) # Add coords to ensure the result is the same whether sample or source are # coords in the input or separate function arguments. if sample is not None: - out.coords['sample_position'] = sample['position'] - out.coords.set_aligned('sample_position', False) + out.coords["sample_position"] = sample["position"] + out.coords.set_aligned("sample_position", False) if source is not None: - out.coords['source_position'] = source['position'] - out.coords.set_aligned('source_position', False) + out.coords["source_position"] = source["position"] + out.coords.set_aligned("source_position", False) return DspacingData[RunType](out) def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: out = data.copy(deep=False) - outer = out.coords.pop('wavelength', None) + outer = out.coords.pop("wavelength", None) if out.bins is not None: - binned = out.bins.coords.pop('wavelength', None) + binned = out.bins.coords.pop("wavelength", None) else: binned = None if outer is not None or binned is not None: get_logger().info("Discarded coordinate 'wavelength' in favor of 'tof'.") - if 'wavelength' in out.dims: - out = out.rename_dims(wavelength='tof') + if "wavelength" in out.dims: + out = out.rename_dims(wavelength="tof") return out -providers = (to_dspacing_with_calibration,) +providers_with_calibration = (to_dspacing_with_calibration,) """Sciline providers for coordinate transformations.""" + +providers_with_positions = ( + powder_coordinate_transformation_graph, + to_wavelength_with_positions, + to_twotheta_with_positions, + to_dspacing_with_positions, +) diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 13197bc8..7a298c09 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -9,6 +9,8 @@ AccumulatedProtonCharge, CalibrationFilename, Filename, + FilenameType, + FilePath, ProtonCharge, RawCalibrationData, RawDataAndMetadata, @@ -18,29 +20,29 @@ ) from .types import DetectorInfo -_version = '1' +_version = "1" -__all__ = ['get_path'] +__all__ = ["_get_path"] def _make_pooch(): import pooch return pooch.create( - path=pooch.os_cache('ess/powgen'), - env='ESS_DATA_DIR', - base_url='https://public.esss.dk/groups/scipp/ess/powgen/{version}/', + path=pooch.os_cache("ess/powgen"), + env="ESS_DATA_DIR", + base_url="https://public.esss.dk/groups/scipp/ess/powgen/{version}/", version=_version, registry={ # Files loadable with Mantid - 'PG3_4844_event.nxs': 'md5:d5ae38871d0a09a28ae01f85d969de1e', - 'PG3_4866_event.nxs': 'md5:3d543bc6a646e622b3f4542bc3435e7e', - 'PG3_5226_event.nxs': 'md5:58b386ebdfeb728d34fd3ba00a2d4f1e', - 'PG3_FERNS_d4832_2011_08_24.cal': 'md5:c181221ebef9fcf30114954268c7a6b6', + "PG3_4844_event.nxs": "md5:d5ae38871d0a09a28ae01f85d969de1e", + "PG3_4866_event.nxs": "md5:3d543bc6a646e622b3f4542bc3435e7e", + "PG3_5226_event.nxs": "md5:58b386ebdfeb728d34fd3ba00a2d4f1e", + "PG3_FERNS_d4832_2011_08_24.cal": "md5:c181221ebef9fcf30114954268c7a6b6", # Zipped Scipp HDF5 files - 'PG3_4844_event.zip': 'md5:a644c74f5e740385469b67431b690a3e', - 'PG3_4866_event.zip': 'md5:5bc49def987f0faeb212a406b92b548e', - 'PG3_FERNS_d4832_2011_08_24.zip': 'md5:0fef4ed5f61465eaaa3f87a18f5bb80d', + "PG3_4844_event.zip": "md5:a644c74f5e740385469b67431b690a3e", + "PG3_4866_event.zip": "md5:5bc49def987f0faeb212a406b92b548e", + "PG3_FERNS_d4832_2011_08_24.zip": "md5:0fef4ed5f61465eaaa3f87a18f5bb80d", }, ) @@ -48,7 +50,7 @@ def _make_pooch(): _pooch = _make_pooch() -def get_path(name: str, unzip: bool = False) -> str: +def _get_path(name: str, unzip: bool = False) -> str: """ Return the path to a data file bundled with scippneutron. @@ -61,37 +63,50 @@ def get_path(name: str, unzip: bool = False) -> str: def mantid_sample_file() -> str: - return get_path('PG3_4844_event.nxs') + return _get_path("PG3_4844_event.nxs") def mantid_vanadium_file() -> str: - return get_path('PG3_4866_event.nxs') + return _get_path("PG3_4866_event.nxs") def mantid_empty_instrument_file() -> str: - return get_path('PG3_5226_event.nxs') + return _get_path("PG3_5226_event.nxs") def mantid_calibration_file() -> str: - return get_path('PG3_FERNS_d4832_2011_08_24.cal') + return _get_path("PG3_FERNS_d4832_2011_08_24.cal") def sample_file() -> str: - (path,) = get_path('PG3_4844_event.zip', unzip=True) + (path,) = _get_path("PG3_4844_event.zip", unzip=True) return path def vanadium_file() -> str: - (path,) = get_path('PG3_4866_event.zip', unzip=True) + (path,) = _get_path("PG3_4866_event.zip", unzip=True) return path def calibration_file() -> str: - (path,) = get_path('PG3_FERNS_d4832_2011_08_24.zip', unzip=True) + (path,) = _get_path("PG3_FERNS_d4832_2011_08_24.zip", unzip=True) return path -def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]: +def get_path(filename: FilenameType) -> FilePath[FilenameType]: + """Translate any filename to a path to the file obtained from pooch registry.""" + if filename.endswith(".zip"): + (path,) = _get_path(filename, unzip=True) + else: + path = _get_path(filename) + return FilePath[FilenameType](path) + + +def _load_hdf5(filename: str) -> sc.DataArray: + return sc.io.load_hdf5(filename) + + +def pooch_load(filename: FilePath[Filename[RunType]]) -> RawDataAndMetadata[RunType]: """Load a file with pooch. If the file is a zip archive, it is extracted and a path to the contained @@ -99,45 +114,44 @@ def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]: The loaded data holds both the events and any metadata from the file. """ - if filename.endswith('.zip'): - (path,) = get_path(filename, unzip=True) - else: - path = get_path(filename) - return RawDataAndMetadata[RunType](sc.io.load_hdf5(path)) + return RawDataAndMetadata[RunType](_load_hdf5(filename)) -def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData: +def pooch_load_calibration( + filename: FilePath[CalibrationFilename], +) -> RawCalibrationData: """Load the calibration data for the POWGEN test data.""" - if filename.endswith('.zip'): - (path,) = get_path(filename, unzip=True) - else: - path = get_path(filename) - return RawCalibrationData(sc.io.load_hdf5(path)) + # if filename.endswith(".zip"): + # (path,) = _get_path(filename, unzip=True) + # else: + # path = _get_path(filename) + return RawCalibrationData(_load_hdf5(filename)) def extract_raw_data(dg: RawDataAndMetadata[RunType]) -> RawDetectorData[RunType]: """Return the events from a loaded data group.""" - return RawDetectorData[RunType](dg['data']) + return RawDetectorData[RunType](dg["data"]) def extract_detector_info(dg: RawDataAndMetadata[SampleRun]) -> DetectorInfo: """Return the detector info from a loaded data group.""" - return DetectorInfo(dg['detector_info']) + return DetectorInfo(dg["detector_info"]) def extract_proton_charge(dg: RawDataAndMetadata[RunType]) -> ProtonCharge[RunType]: """Return the proton charge from a loaded data group.""" - return ProtonCharge[RunType](dg['proton_charge']) + return ProtonCharge[RunType](dg["proton_charge"]) def extract_accumulated_proton_charge( data: RawDetectorData[RunType], ) -> AccumulatedProtonCharge[RunType]: """Return the stored accumulated proton charge from a loaded data group.""" - return AccumulatedProtonCharge[RunType](data.coords['gd_prtn_chrg']) + return AccumulatedProtonCharge[RunType](data.coords["gd_prtn_chrg"]) providers = ( + get_path, pooch_load, pooch_load_calibration, extract_accumulated_proton_charge, diff --git a/src/ess/powder/filtering.py b/src/ess/powder/filtering.py index 45f9e247..3f8ba022 100644 --- a/src/ess/powder/filtering.py +++ b/src/ess/powder/filtering.py @@ -17,10 +17,10 @@ def _equivalent_bin_indices(a, b) -> bool: - a_begin = a.bins.constituents['begin'].flatten(to='') - a_end = a.bins.constituents['end'].flatten(to='') - b_begin = b.bins.constituents['begin'].flatten(to='') - b_end = b.bins.constituents['end'].flatten(to='') + a_begin = a.bins.constituents["begin"].flatten(to="") + a_end = a.bins.constituents["end"].flatten(to="") + b_begin = b.bins.constituents["begin"].flatten(to="") + b_end = b.bins.constituents["end"].flatten(to="") non_empty = a_begin != a_end return ( sc.all((a_begin == b_begin)[non_empty]).value @@ -33,10 +33,10 @@ def _temporary_bin_coord(data: sc.DataArray, name: str, coord: sc.Variable) -> N if not _equivalent_bin_indices(data, coord): raise ValueError("data and coord do not have equivalent bin indices") coord = sc.bins( - data=coord.bins.constituents['data'], - begin=data.bins.coords['pulse_time'].bins.constituents['begin'], - end=data.bins.coords['pulse_time'].bins.constituents['end'], - dim=coord.bins.constituents['dim'], + data=coord.bins.constituents["data"], + begin=data.bins.coords["pulse_time"].bins.constituents["begin"], + end=data.bins.coords["pulse_time"].bins.constituents["end"], + dim=coord.bins.constituents["dim"], ) data.bins.coords[name] = coord yield @@ -46,7 +46,7 @@ def _temporary_bin_coord(data: sc.DataArray, name: str, coord: sc.Variable) -> N # TODO non-monotonic proton charge -> raise? def _with_pulse_time_edges(da: sc.DataArray, dim: str) -> sc.DataArray: pulse_time = da.coords[dim] - one = sc.scalar(1, dtype='int64', unit=pulse_time.unit) + one = sc.scalar(1, dtype="int64", unit=pulse_time.unit) lo = pulse_time[0] - one hi = pulse_time[-1] + one mid = sc.midpoints(pulse_time) @@ -64,12 +64,12 @@ def remove_bad_pulses( good_pulse = _with_pulse_time_edges(proton_charge >= min_charge, proton_charge.dim) with _temporary_bin_coord( data, - 'good_pulse', + "good_pulse", sc.lookup(good_pulse, good_pulse.dim)[data.bins.coords[good_pulse.dim]], ): - filtered = data.group(sc.array(dims=['good_pulse'], values=[True])) - filtered = filtered.squeeze('good_pulse').copy(deep=False) - del filtered.coords['good_pulse'] + filtered = data.group(sc.array(dims=["good_pulse"], values=[True])) + filtered = filtered.squeeze("good_pulse").copy(deep=False) + del filtered.coords["good_pulse"] return filtered @@ -92,7 +92,7 @@ def crop_tof( : Cropped data. """ - tof = event_or_outer_coord(data, 'tof') + tof = event_or_outer_coord(data, "tof") tof_unit = elem_unit(tof) tof_dtype = elem_dtype(tof) return TofCroppedData[RunType]( diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py new file mode 100644 index 00000000..22211e1a --- /dev/null +++ b/src/ess/powder/masking.py @@ -0,0 +1,88 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) +""" +Masking functions for the powder workflow. +""" +import numpy as np +import sciline +import scipp as sc + +from ess.reduce.masking import mask_range + +from .types import ( + FilePath, + MaskedDetectorIDs, + NormalizedByProtonCharge, + PixelMaskedData, + PixelMaskFilename, + RunType, + TwoThetaData, + TwoThetaMask, + TwoThetaMaskedData, + WavelengthData, + WavelengthMask, + WavelengthMaskedData, +) + + +def read_pixel_masks( + filename: FilePath[PixelMaskFilename], +) -> MaskedDetectorIDs: + """Read a pixel mask from a Scipp hdf5 file. + + Parameters + ---------- + filename: + Path to the hdf5 file. + """ + return MaskedDetectorIDs(sc.io.load_hdf5(filename)) + + +def apply_pixel_masks( + data: NormalizedByProtonCharge[RunType], + masked_ids: sciline.Series[PixelMaskFilename, MaskedDetectorIDs], +) -> PixelMaskedData[RunType]: + """Apply pixel-specific masks to raw data. + The masks are based on detector IDs stored in XML files. + + Parameters + ---------- + data: + Raw data with configured component positions. + masks: + A series of masks. + """ + masked_ids = {"pix_mask": sc.arange("spectrum", 1, 101)} + if len(masked_ids) > 0: + data = data.copy(deep=False) + key = ( + set(data.coords.keys()) & {"detector_number", "detector_id", "spectrum"} + ).pop() + ids = data.coords[key] + for name, masked in masked_ids.items(): + mask = sc.zeros(sizes=ids.sizes, dtype="bool") + mask.values[np.isin(ids.values, masked.values)] = True + data.masks[name] = mask + return PixelMaskedData[RunType](data) + + +def apply_wavelength_masks( + da: WavelengthData[RunType], mask: WavelengthMask +) -> WavelengthMaskedData[RunType]: + if "wavelength" in da.coords and da.coords["wavelength"].ndim > 1: + da = da.bin(wavelength=1) + return WavelengthMaskedData[RunType](mask_range(da, mask=mask)) + + +def apply_twotheta_masks( + da: TwoThetaData[RunType], mask: TwoThetaMask +) -> TwoThetaMaskedData[RunType]: + return TwoThetaMaskedData[RunType](mask_range(da, mask=mask)) + + +providers = ( + read_pixel_masks, + apply_pixel_masks, + apply_wavelength_masks, + apply_twotheta_masks, +) diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 5a4310c5..20582817 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -17,21 +17,23 @@ # 1 TypeVars used to parametrize the generic parts of the workflow # 1.1 Run types -EmptyCanRun = NewType('EmptyCanRun', int) +EmptyCanRun = NewType("EmptyCanRun", int) """Empty sample can run.""" -EmptyInstrumentRun = NewType('EmptyInstrumentRun', int) +EmptyInstrumentRun = NewType("EmptyInstrumentRun", int) """Empty instrument run.""" -SampleRun = NewType('SampleRun', int) +SampleRun = NewType("SampleRun", int) """Sample run.""" -VanadiumRun = NewType('VanadiumRun', int) +VanadiumRun = NewType("VanadiumRun", int) """Vanadium run.""" -RunType = TypeVar('RunType', EmptyInstrumentRun, SampleRun, VanadiumRun) +RunType = TypeVar("RunType", EmptyInstrumentRun, SampleRun, VanadiumRun) """TypeVar used for specifying the run.""" +FilenameType = TypeVar("FilenameType", bound=str) + # 2 Workflow parameters -CalibrationFilename = NewType('CalibrationFilename', str) +CalibrationFilename = NewType("CalibrationFilename", str) """Filename of the instrument calibration file.""" @@ -39,13 +41,13 @@ class DetectorName(str, Enum): """Name of a detector.""" - mantle = 'mantle' - high_resolution = 'high_resolution' - endcap_backward = 'endcap_backward' - endcap_forward = 'endcap_forward' + mantle = "mantle" + high_resolution = "high_resolution" + endcap_backward = "endcap_backward" + endcap_forward = "endcap_forward" -DspacingBins = NewType('DSpacingBins', sc.Variable) +DspacingBins = NewType("DSpacingBins", sc.Variable) """Bin edges for d-spacing.""" @@ -53,14 +55,14 @@ class Filename(sciline.Scope[RunType, str], str): """Name of an input file.""" -class FilePath(sciline.Scope[RunType, Path], Path): +class FilePath(sciline.Scope[FilenameType, Path], Path): """Path to an input file on disk.""" -OutFilename = NewType('OutFilename', str) +OutFilename = NewType("OutFilename", str) """Filename of the output.""" -TwoThetaBins = NewType('TwoThetaBins', sc.Variable) +TwoThetaBins = NewType("TwoThetaBins", sc.Variable) """Bin edges for grouping in 2theta. This is used by an alternative focussing step that groups detector @@ -68,14 +70,14 @@ class FilePath(sciline.Scope[RunType, Path], Path): """ UncertaintyBroadcastMode = Enum( - 'UncertaintyBroadcastMode', ['drop', 'upper_bound', 'fail'] + "UncertaintyBroadcastMode", ["drop", "upper_bound", "fail"] ) """Mode for broadcasting uncertainties. See https://doi.org/10.3233/JNR-220049 for context. """ -ValidTofRange = NewType('ValidTofRange', sc.Variable) +ValidTofRange = NewType("ValidTofRange", sc.Variable) """Min and max tof value of the instrument.""" # 3 Workflow (intermediate) results @@ -90,9 +92,11 @@ def __init__(self, *args: Any, **kwargs: Any) -> None: super().__init__(*args, **kwargs) -CalibrationData = NewType('CalibrationData', sc.Dataset) +CalibrationData = NewType("CalibrationData", sc.Dataset) """Detector calibration data.""" +DataFolder = NewType("DataFolder", str) + class DetectorDimensions(sciline.Scope[DetectorName, tuple[str, ...]], tuple[str, ...]): """Logical detector dimensions.""" @@ -106,9 +110,12 @@ class DspacingDataWithoutVariances(sciline.Scope[RunType, sc.DataArray], sc.Data """Data converted to d-spacing where variances where removed.""" -DspacingHistogram = NewType('DspacingHistogram', sc.DataArray) +DspacingHistogram = NewType("DspacingHistogram", sc.DataArray) """Histogrammed intensity vs d-spacing.""" +ElasticCoordTransformGraph = NewType("ElasticCoordTransformGraph", dict) +"""Graph for transforming coordinates in elastic scattering.""" + class FilteredData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Raw data without invalid events.""" @@ -118,19 +125,31 @@ class FocussedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Intensity vs d-spacing after focussing pixels.""" +MaskedDetectorIDs = NewType("MaskedDetectorIDs", sc.Variable) +"""1-D variable listing all masked detector IDs.""" + + class NormalizedByProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data that has been normalized by proton charge.""" -NormalizedByVanadium = NewType('NormalizedByVanadium', sc.DataArray) +NormalizedByVanadium = NewType("NormalizedByVanadium", sc.DataArray) """Data that has been normalized by a vanadium run.""" +class PixelMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with masked pixels.""" + + +PixelMaskFilename = NewType("PixelMaskFilename", str) +"""Filename of a pixel mask.""" + + class ProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Time-dependent proton charge.""" -RawCalibrationData = NewType('RawCalibrationData', sc.Dataset) +RawCalibrationData = NewType("RawCalibrationData", sc.Dataset) """Calibration data as loaded from file, needs preprocessing before using.""" @@ -150,7 +169,7 @@ class RawSample(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Raw data from a loaded sample.""" -RawSource = NewType('RawSource', sc.DataGroup) +RawSource = NewType("RawSource", sc.DataGroup) """Raw data from a loaded neutron source.""" @@ -158,4 +177,36 @@ class TofCroppedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Raw data cropped to the valid TOF range.""" +class TwoThetaData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data converted to 2theta.""" + + +TwoThetaMask = NewType("TwoThetaMask", sc.DataArray) +"""A data array defining a mask to be applied in 2theta. +Only one-dimensional masks are supported. +The data array should contain a bin-edge coordinate which represents +the edges of the ranges to be masked. The values of the data array represent the +mask values (``True`` or ``False``) inside each range defined by the coordinate.""" + + +class TwoThetaMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with masks in 2theta.""" + + +class WavelengthData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data converted to wavelength.""" + + +WavelengthMask = NewType("WavelengthMask", sc.DataArray) +"""A data array defining a mask to be applied in wavelength. +Only one-dimensional masks are supported. +The data array should contain a bin-edge coordinate which represents +the edges of the ranges to be masked. The values of the data array represent the +mask values (``True`` or ``False``) inside each range defined by the coordinate.""" + + +class WavelengthMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with masks in wavelength.""" + + del sc, sciline, NewType, TypeVar From 1b9e13a58a1d529b9ce38ae42473b2f4dd44024d Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 15 May 2024 16:28:38 +0200 Subject: [PATCH 02/26] work on the raw events --- src/ess/powder/conversion.py | 17 ++++-- src/ess/powder/external/powgen/data.py | 72 ++++++++++---------------- src/ess/powder/filtering.py | 3 +- src/ess/powder/masking.py | 53 ++++++++++++++----- src/ess/powder/types.py | 36 ++++++------- 5 files changed, 102 insertions(+), 79 deletions(-) diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index e0e222a4..a37b3574 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -13,6 +13,7 @@ from .logging import get_logger from .types import ( CalibrationData, + DataWithScatteringCoordinates, DspacingData, ElasticCoordTransformGraph, NormalizedByProtonCharge, @@ -274,12 +275,22 @@ def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: return out +def add_scattering_coordinates( + data: PixelMaskedData[RunType], graph: ElasticCoordTransformGraph +) -> DataWithScatteringCoordinates[RunType]: + out = data.transform_coords( + ["two_theta", "wavelength", "dspacing"], graph=graph, keep_intermediate=False + ) + return DataWithScatteringCoordinates[RunType](out) + + providers_with_calibration = (to_dspacing_with_calibration,) """Sciline providers for coordinate transformations.""" providers_with_positions = ( powder_coordinate_transformation_graph, - to_wavelength_with_positions, - to_twotheta_with_positions, - to_dspacing_with_positions, + add_scattering_coordinates, + # to_wavelength_with_positions, + # to_twotheta_with_positions, + # to_dspacing_with_positions, ) diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 7a298c09..3d26e1a7 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -9,8 +9,6 @@ AccumulatedProtonCharge, CalibrationFilename, Filename, - FilenameType, - FilePath, ProtonCharge, RawCalibrationData, RawDataAndMetadata, @@ -22,8 +20,6 @@ _version = "1" -__all__ = ["_get_path"] - def _make_pooch(): import pooch @@ -50,7 +46,7 @@ def _make_pooch(): _pooch = _make_pooch() -def _get_path(name: str, unzip: bool = False) -> str: +def _get_path(name: str) -> str: """ Return the path to a data file bundled with scippneutron. @@ -59,54 +55,42 @@ def _get_path(name: str, unzip: bool = False) -> str: """ import pooch - return _pooch.fetch(name, processor=pooch.Unzip() if unzip else None) + if name.endswith(".zip"): + (path,) = _pooch.fetch(name, processor=pooch.Unzip()) + else: + path = _pooch.fetch(name) + return path -def mantid_sample_file() -> str: +def powgen_tutorial_mantid_sample_file() -> str: return _get_path("PG3_4844_event.nxs") -def mantid_vanadium_file() -> str: +def powgen_tutorial_mantid_vanadium_file() -> str: return _get_path("PG3_4866_event.nxs") -def mantid_empty_instrument_file() -> str: +def powgen_tutorial_mantid_empty_instrument_file() -> str: return _get_path("PG3_5226_event.nxs") -def mantid_calibration_file() -> str: +def powgen_tutorial_mantid_calibration_file() -> str: return _get_path("PG3_FERNS_d4832_2011_08_24.cal") -def sample_file() -> str: - (path,) = _get_path("PG3_4844_event.zip", unzip=True) - return path - +def powgen_tutorial_sample_file() -> str: + return _get_path("PG3_4844_event.zip") -def vanadium_file() -> str: - (path,) = _get_path("PG3_4866_event.zip", unzip=True) - return path - -def calibration_file() -> str: - (path,) = _get_path("PG3_FERNS_d4832_2011_08_24.zip", unzip=True) - return path - - -def get_path(filename: FilenameType) -> FilePath[FilenameType]: - """Translate any filename to a path to the file obtained from pooch registry.""" - if filename.endswith(".zip"): - (path,) = _get_path(filename, unzip=True) - else: - path = _get_path(filename) - return FilePath[FilenameType](path) +def powgen_tutorial_vanadium_file() -> str: + return _get_path("PG3_4866_event.zip") -def _load_hdf5(filename: str) -> sc.DataArray: - return sc.io.load_hdf5(filename) +def powgen_tutorial_calibration_file() -> str: + return _get_path("PG3_FERNS_d4832_2011_08_24.zip") -def pooch_load(filename: FilePath[Filename[RunType]]) -> RawDataAndMetadata[RunType]: +def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]: """Load a file with pooch. If the file is a zip archive, it is extracted and a path to the contained @@ -114,23 +98,24 @@ def pooch_load(filename: FilePath[Filename[RunType]]) -> RawDataAndMetadata[RunT The loaded data holds both the events and any metadata from the file. """ - return RawDataAndMetadata[RunType](_load_hdf5(filename)) + return RawDataAndMetadata[RunType](sc.io.load_hdf5(filename)) -def pooch_load_calibration( - filename: FilePath[CalibrationFilename], -) -> RawCalibrationData: +def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData: """Load the calibration data for the POWGEN test data.""" - # if filename.endswith(".zip"): - # (path,) = _get_path(filename, unzip=True) - # else: - # path = _get_path(filename) - return RawCalibrationData(_load_hdf5(filename)) + return RawCalibrationData(sc.io.load_hdf5(filename)) def extract_raw_data(dg: RawDataAndMetadata[RunType]) -> RawDetectorData[RunType]: """Return the events from a loaded data group.""" - return RawDetectorData[RunType](dg["data"]) + out = dg["data"].squeeze() + del out.coords["tof"] + out.bins.coords["position"] = sc.bins_like(out, out.coords["position"]) + out.bins.coords["spectrum"] = sc.bins_like(out, out.coords["spectrum"]) + raw_events = out.bins.constituents["data"].copy() + for c in ("gd_prtn_chrg", "sample_position", "source_position"): + raw_events.coords[c] = out.coords[c] + return RawDetectorData[RunType](raw_events) def extract_detector_info(dg: RawDataAndMetadata[SampleRun]) -> DetectorInfo: @@ -151,7 +136,6 @@ def extract_accumulated_proton_charge( providers = ( - get_path, pooch_load, pooch_load_calibration, extract_accumulated_proton_charge, diff --git a/src/ess/powder/filtering.py b/src/ess/powder/filtering.py index 3f8ba022..f7112326 100644 --- a/src/ess/powder/filtering.py +++ b/src/ess/powder/filtering.py @@ -100,7 +100,8 @@ def crop_tof( ) -def filter_events(data: TofCroppedData[RunType]) -> FilteredData[RunType]: +# def filter_events(data: TofCroppedData[RunType]) -> FilteredData[RunType]: +def filter_events(data: RawDetectorData[RunType]) -> FilteredData[RunType]: """Remove bad events. Attention diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py index 22211e1a..8c164694 100644 --- a/src/ess/powder/masking.py +++ b/src/ess/powder/masking.py @@ -3,6 +3,8 @@ """ Masking functions for the powder workflow. """ +from typing import Union, Optional + import numpy as np import sciline import scipp as sc @@ -10,12 +12,14 @@ from ess.reduce.masking import mask_range from .types import ( - FilePath, + DataWithScatteringCoordinates, MaskedDetectorIDs, NormalizedByProtonCharge, PixelMaskedData, PixelMaskFilename, RunType, + TofMask, + TofMaskedData, TwoThetaData, TwoThetaMask, TwoThetaMaskedData, @@ -26,7 +30,7 @@ def read_pixel_masks( - filename: FilePath[PixelMaskFilename], + filename: PixelMaskFilename, ) -> MaskedDetectorIDs: """Read a pixel mask from a Scipp hdf5 file. @@ -66,23 +70,48 @@ def apply_pixel_masks( return PixelMaskedData[RunType](data) -def apply_wavelength_masks( - da: WavelengthData[RunType], mask: WavelengthMask +# def apply_wavelength_masks( +# da: WavelengthData[RunType], mask: WavelengthMask +# ) -> WavelengthMaskedData[RunType]: +# if "wavelength" in da.coords and da.coords["wavelength"].ndim > 1: +# da = da.bin(wavelength=1) +# return WavelengthMaskedData[RunType](mask_range(da, mask=mask)) + + +# def apply_twotheta_masks( +# da: TwoThetaData[RunType], mask: TwoThetaMask +# ) -> TwoThetaMaskedData[RunType]: +# return TwoThetaMaskedData[RunType](mask_range(da, mask=mask)) + + +def apply_tof_masking( + da: DataWithScatteringCoordinates[RunType], mask_func: TofMask +) -> TofMaskedData[RunType]: + out = da.copy(deep=False) + out.masks["tof"] = mask_func(da.coords["tof"]) + return TofMaskedData[RunType](out) + + +def apply_wavelength_masking( + da: TofMaskedData[RunType], mask_func: WavelengthMask ) -> WavelengthMaskedData[RunType]: - if "wavelength" in da.coords and da.coords["wavelength"].ndim > 1: - da = da.bin(wavelength=1) - return WavelengthMaskedData[RunType](mask_range(da, mask=mask)) + out = da.copy(deep=False) + out.masks["wavelength"] = mask_func(da.coords["wavelength"]) + return WavelengthMaskedData[RunType](out) -def apply_twotheta_masks( - da: TwoThetaData[RunType], mask: TwoThetaMask +def apply_twotheta_masking( + da: WavelengthMaskedData[RunType], mask_func: TwoThetaMask ) -> TwoThetaMaskedData[RunType]: - return TwoThetaMaskedData[RunType](mask_range(da, mask=mask)) + out = da.copy(deep=False) + out.masks["two_theta"] = mask_func(da.coords["two_theta"]) + return TwoThetaMaskedData[RunType](out) providers = ( read_pixel_masks, apply_pixel_masks, - apply_wavelength_masks, - apply_twotheta_masks, + apply_tof_masking, + apply_twotheta_masking, + apply_wavelength_masking, ) diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 20582817..8b5cca92 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -9,7 +9,7 @@ from enum import Enum from pathlib import Path -from typing import Any, NewType, TypeVar +from typing import Any, Callable, NewType, TypeVar import sciline import scipp as sc @@ -28,8 +28,6 @@ RunType = TypeVar("RunType", EmptyInstrumentRun, SampleRun, VanadiumRun) """TypeVar used for specifying the run.""" -FilenameType = TypeVar("FilenameType", bound=str) - # 2 Workflow parameters @@ -55,10 +53,6 @@ class Filename(sciline.Scope[RunType, str], str): """Name of an input file.""" -class FilePath(sciline.Scope[FilenameType, Path], Path): - """Path to an input file on disk.""" - - OutFilename = NewType("OutFilename", str) """Filename of the output.""" @@ -98,6 +92,11 @@ def __init__(self, *args: Any, **kwargs: Any) -> None: DataFolder = NewType("DataFolder", str) +class DataWithScatteringCoordinates(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with scattering coordinates computed for all events: wavelength, 2theta, + d-spacing.""" + + class DetectorDimensions(sciline.Scope[DetectorName, tuple[str, ...]], tuple[str, ...]): """Logical detector dimensions.""" @@ -172,6 +171,13 @@ class RawSample(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): RawSource = NewType("RawSource", sc.DataGroup) """Raw data from a loaded neutron source.""" +TofMask = NewType("TofMask", Callable) +"""""" + + +class TofMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with masked TOF values.""" + class TofCroppedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Raw data cropped to the valid TOF range.""" @@ -181,12 +187,8 @@ class TwoThetaData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data converted to 2theta.""" -TwoThetaMask = NewType("TwoThetaMask", sc.DataArray) -"""A data array defining a mask to be applied in 2theta. -Only one-dimensional masks are supported. -The data array should contain a bin-edge coordinate which represents -the edges of the ranges to be masked. The values of the data array represent the -mask values (``True`` or ``False``) inside each range defined by the coordinate.""" +TwoThetaMask = NewType("TwoThetaMask", Callable) +"""""" class TwoThetaMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): @@ -197,12 +199,8 @@ class WavelengthData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data converted to wavelength.""" -WavelengthMask = NewType("WavelengthMask", sc.DataArray) -"""A data array defining a mask to be applied in wavelength. -Only one-dimensional masks are supported. -The data array should contain a bin-edge coordinate which represents -the edges of the ranges to be masked. The values of the data array represent the -mask values (``True`` or ``False``) inside each range defined by the coordinate.""" +WavelengthMask = NewType("WavelengthMask", Callable) +"""""" class WavelengthMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): From f8ab74878da0cbc31d8e2da69ca58cc829a5f5fb Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Thu, 16 May 2024 15:40:14 +0200 Subject: [PATCH 03/26] now use data binned in logical dimensions --- src/ess/powder/conversion.py | 2 +- src/ess/powder/external/powgen/data.py | 14 +-- src/ess/powder/masking.py | 115 +++++++++++++------------ src/ess/powder/types.py | 5 ++ 4 files changed, 75 insertions(+), 61 deletions(-) diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index a37b3574..f9d48c36 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -276,7 +276,7 @@ def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: def add_scattering_coordinates( - data: PixelMaskedData[RunType], graph: ElasticCoordTransformGraph + data: NormalizedByProtonCharge[RunType], graph: ElasticCoordTransformGraph ) -> DataWithScatteringCoordinates[RunType]: out = data.transform_coords( ["two_theta", "wavelength", "dspacing"], graph=graph, keep_intermediate=False diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 3d26e1a7..7d2e8ead 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -110,12 +110,14 @@ def extract_raw_data(dg: RawDataAndMetadata[RunType]) -> RawDetectorData[RunType """Return the events from a loaded data group.""" out = dg["data"].squeeze() del out.coords["tof"] - out.bins.coords["position"] = sc.bins_like(out, out.coords["position"]) - out.bins.coords["spectrum"] = sc.bins_like(out, out.coords["spectrum"]) - raw_events = out.bins.constituents["data"].copy() - for c in ("gd_prtn_chrg", "sample_position", "source_position"): - raw_events.coords[c] = out.coords[c] - return RawDetectorData[RunType](raw_events) + # out = out.fold(dim="spectrum", sizes={"column": 154, "row": 7, "bank": 23}) + out = out.fold(dim="spectrum", sizes={"bank": 23, "column": 154, "row": 7}) + # out.bins.coords["position"] = sc.bins_like(out, out.coords["position"]) + # out.bins.coords["spectrum"] = sc.bins_like(out, out.coords["spectrum"]) + # raw_events = out.bins.constituents["data"].copy() + # for c in ("gd_prtn_chrg", "sample_position", "source_position"): + # raw_events.coords[c] = out.coords[c] + return RawDetectorData[RunType](out) def extract_detector_info(dg: RawDataAndMetadata[SampleRun]) -> DetectorInfo: diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py index 8c164694..e592fff7 100644 --- a/src/ess/powder/masking.py +++ b/src/ess/powder/masking.py @@ -13,6 +13,7 @@ from .types import ( DataWithScatteringCoordinates, + MaskedData, MaskedDetectorIDs, NormalizedByProtonCharge, PixelMaskedData, @@ -42,76 +43,82 @@ def read_pixel_masks( return MaskedDetectorIDs(sc.io.load_hdf5(filename)) -def apply_pixel_masks( - data: NormalizedByProtonCharge[RunType], - masked_ids: sciline.Series[PixelMaskFilename, MaskedDetectorIDs], -) -> PixelMaskedData[RunType]: - """Apply pixel-specific masks to raw data. - The masks are based on detector IDs stored in XML files. - - Parameters - ---------- - data: - Raw data with configured component positions. - masks: - A series of masks. - """ - masked_ids = {"pix_mask": sc.arange("spectrum", 1, 101)} - if len(masked_ids) > 0: - data = data.copy(deep=False) +def apply_masks( + data: DataWithScatteringCoordinates[RunType], + masked_pixel_ids: sciline.Series[PixelMaskFilename, MaskedDetectorIDs], + tof_mask_func: Optional[TofMask] = None, + wavelength_mask_func: Optional[WavelengthMask] = None, + two_theta_mask_func: Optional[TwoThetaMask] = None, +) -> MaskedData[RunType]: + """ """ + out = data.copy(deep=False) + masked_pixel_ids = {"pix_mask": sc.arange("spectrum", 1, 101)} + if len(masked_pixel_ids) > 0: key = ( - set(data.coords.keys()) & {"detector_number", "detector_id", "spectrum"} + set(out.coords.keys()) & {"detector_number", "detector_id", "spectrum"} ).pop() - ids = data.coords[key] - for name, masked in masked_ids.items(): + ids = out.coords[key] + for name, masked in masked_pixel_ids.items(): mask = sc.zeros(sizes=ids.sizes, dtype="bool") mask.values[np.isin(ids.values, masked.values)] = True - data.masks[name] = mask - return PixelMaskedData[RunType](data) + out.masks[name] = mask + for dim, mask in { + "tof": tof_mask_func, + "wavelength": wavelength_mask_func, + "two_theta": two_theta_mask_func, + }.items(): + if mask is not None: + if dim in out.bins.coords: + out.bins.masks[dim] = mask(out.bins.coords[dim]) + else: + out.masks[dim] = mask(out.coords[dim]) -# def apply_wavelength_masks( -# da: WavelengthData[RunType], mask: WavelengthMask -# ) -> WavelengthMaskedData[RunType]: -# if "wavelength" in da.coords and da.coords["wavelength"].ndim > 1: -# da = da.bin(wavelength=1) -# return WavelengthMaskedData[RunType](mask_range(da, mask=mask)) + return MaskedData[RunType](out) -# def apply_twotheta_masks( -# da: TwoThetaData[RunType], mask: TwoThetaMask -# ) -> TwoThetaMaskedData[RunType]: -# return TwoThetaMaskedData[RunType](mask_range(da, mask=mask)) +# # def apply_wavelength_masks( +# # da: WavelengthData[RunType], mask: WavelengthMask +# # ) -> WavelengthMaskedData[RunType]: +# # if "wavelength" in da.coords and da.coords["wavelength"].ndim > 1: +# # da = da.bin(wavelength=1) +# # return WavelengthMaskedData[RunType](mask_range(da, mask=mask)) + +# # def apply_twotheta_masks( +# # da: TwoThetaData[RunType], mask: TwoThetaMask +# # ) -> TwoThetaMaskedData[RunType]: +# # return TwoThetaMaskedData[RunType](mask_range(da, mask=mask)) -def apply_tof_masking( - da: DataWithScatteringCoordinates[RunType], mask_func: TofMask -) -> TofMaskedData[RunType]: - out = da.copy(deep=False) - out.masks["tof"] = mask_func(da.coords["tof"]) - return TofMaskedData[RunType](out) +# def apply_tof_masking( +# da: DataWithScatteringCoordinates[RunType], mask_func: TofMask +# ) -> TofMaskedData[RunType]: +# out = da.copy(deep=False) +# out.masks["tof"] = mask_func(da.coords["tof"]) +# return TofMaskedData[RunType](out) -def apply_wavelength_masking( - da: TofMaskedData[RunType], mask_func: WavelengthMask -) -> WavelengthMaskedData[RunType]: - out = da.copy(deep=False) - out.masks["wavelength"] = mask_func(da.coords["wavelength"]) - return WavelengthMaskedData[RunType](out) +# def apply_wavelength_masking( +# da: TofMaskedData[RunType], mask_func: WavelengthMask +# ) -> WavelengthMaskedData[RunType]: +# out = da.copy(deep=False) +# out.masks["wavelength"] = mask_func(da.coords["wavelength"]) +# return WavelengthMaskedData[RunType](out) -def apply_twotheta_masking( - da: WavelengthMaskedData[RunType], mask_func: TwoThetaMask -) -> TwoThetaMaskedData[RunType]: - out = da.copy(deep=False) - out.masks["two_theta"] = mask_func(da.coords["two_theta"]) - return TwoThetaMaskedData[RunType](out) + +# def apply_twotheta_masking( +# da: WavelengthMaskedData[RunType], mask_func: TwoThetaMask +# ) -> TwoThetaMaskedData[RunType]: +# out = da.copy(deep=False) +# out.masks["two_theta"] = mask_func(da.coords["two_theta"]) +# return TwoThetaMaskedData[RunType](out) providers = ( read_pixel_masks, - apply_pixel_masks, - apply_tof_masking, - apply_twotheta_masking, - apply_wavelength_masking, + apply_masks, + # apply_tof_masking, + # apply_twotheta_masking, + # apply_wavelength_masking, ) diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 8b5cca92..9cab812d 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -124,6 +124,11 @@ class FocussedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Intensity vs d-spacing after focussing pixels.""" +class MaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data with masked pixels, tof regions, wavelength regions, 2theta regions, or + dspacing regions.""" + + MaskedDetectorIDs = NewType("MaskedDetectorIDs", sc.Variable) """1-D variable listing all masked detector IDs.""" From 6500f9986d9e76479670c6acb480c927eb71707e Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Fri, 17 May 2024 17:02:51 +0200 Subject: [PATCH 04/26] update Powgen notebook --- .../POWGEN_data_reduction.ipynb | 133 ++++++---------- src/ess/dream/data.py | 50 +++--- src/ess/dream/io/geant4.py | 50 +++--- src/ess/powder/__init__.py | 2 +- src/ess/powder/conversion.py | 149 +++++++----------- src/ess/powder/correction.py | 32 ++-- src/ess/powder/external/powgen/beamline.py | 16 +- src/ess/powder/external/powgen/data.py | 7 +- src/ess/powder/grouping.py | 64 +++----- src/ess/powder/types.py | 9 +- 10 files changed, 208 insertions(+), 304 deletions(-) diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index f7853190..4b904c93 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -23,6 +23,7 @@ "outputs": [], "source": [ "import scipp as sc\n", + "import plopp as pp\n", "import scippneutron as scn\n", "import sciline\n", "\n", @@ -51,9 +52,9 @@ "source": [ "params = {\n", " # Input data\n", - " Filename[SampleRun]: 'PG3_4844_event.zip',\n", - " Filename[VanadiumRun]: 'PG3_4866_event.zip',\n", - " CalibrationFilename: 'PG3_FERNS_d4832_2011_08_24.zip',\n", + " Filename[SampleRun]: powgen.data.powgen_tutorial_sample_file(),\n", + " Filename[VanadiumRun]: powgen.data.powgen_tutorial_vanadium_file(),\n", + " CalibrationFilename: powgen.data.powgen_tutorial_calibration_file(),\n", "\n", " # The upper bounds mode is not yet implemented.\n", " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", @@ -62,8 +63,15 @@ " ValidTofRange: sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us'),\n", "\n", " # Edges for binning in d-spacing\n", - " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 200, unit='angstrom'),\n", - "}" + " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'),\n", + "}\n", + "\n", + "# Mask in time-of-flight\n", + "params[TofMask] = lambda x: (x < sc.scalar(0., unit='us')) | (x > sc.scalar(16666.67, unit='us'))\n", + "# Mask in wavelength\n", + "params[WavelengthMask] = lambda x: (x > sc.scalar(0.18, unit='angstrom')) & (x < sc.scalar(0.21, unit='angstrom'))\n", + "# Pixel masks\n", + "pixel_masks = []" ] }, { @@ -83,14 +91,13 @@ "metadata": {}, "outputs": [], "source": [ - "providers = [\n", - " *powder.providers,\n", - " *powgen.providers,\n", - "]\n", + "providers = [*powder.providers, *powgen.providers]\n", "pipeline = sciline.Pipeline(\n", " providers,\n", - " params=params,\n", - ")" + " params=params\n", + ")\n", + "\n", + "pipeline.set_param_series(PixelMaskFilename, pixel_masks)" ] }, { @@ -108,21 +115,11 @@ { "cell_type": "code", "execution_count": null, - "id": "7d72a6ad-b43d-4e05-a70f-7f334ca06b7a", + "id": "d3069da8-de3c-4cc3-bd14-ff8fcbf58e91", "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(DspacingHistogram, graph_attr={'rankdir': 'LR'})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d09bc8d5-4251-483f-bb45-aa8cf67f69b2", - "metadata": {}, - "outputs": [], - "source": [ - "dspacing_histogram = pipeline.get(DspacingHistogram)" + "pipeline.visualize(NormalizedByVanadium, graph_attr={'rankdir': 'LR'})" ] }, { @@ -140,7 +137,7 @@ "metadata": {}, "outputs": [], "source": [ - "result = dspacing_histogram.compute()\n", + "result = pipeline.compute(NormalizedByVanadium)\n", "result" ] }, @@ -151,7 +148,8 @@ "metadata": {}, "outputs": [], "source": [ - "result.plot()" + "dspacing_histogram = result.hist()\n", + "dspacing_histogram.plot()" ] }, { @@ -176,10 +174,10 @@ "metadata": {}, "outputs": [], "source": [ - "def save_xye(reduced_data: DspacingHistogram,\n", + "def save_xye(reduced_data: NormalizedByVanadium,\n", " out_filename: OutFilename,\n", ") -> None:\n", - " data = reduced_data.copy(deep=False)\n", + " data = reduced_data.hist()\n", " data.coords['dspacing'] = sc.midpoints(data.coords['dspacing'])\n", " scn.io.save_xye(out_filename, data, coord='dspacing')" ] @@ -243,6 +241,7 @@ "source": [ "results = pipeline.compute((\n", " RawDetectorData[SampleRun],\n", + " MaskedData[SampleRun],\n", " FilteredData[SampleRun],\n", " FilteredData[VanadiumRun],\n", "))" @@ -265,7 +264,7 @@ "metadata": {}, "outputs": [], "source": [ - "scn.instrument_view(results[RawDetectorData[SampleRun]].hist())" + "scn.instrument_view(results[MaskedData[SampleRun]].hist(), pixel_size=0.05)" ] }, { @@ -276,8 +275,8 @@ "outputs": [], "source": [ "tof_data = sc.DataGroup(\n", - " sample=results[FilteredData[SampleRun]].bins.concat('spectrum'),\n", - " vanadium=results[FilteredData[VanadiumRun]].bins.concat('spectrum'),\n", + " sample=results[FilteredData[SampleRun]].bins.concat(),\n", + " vanadium=results[FilteredData[VanadiumRun]].bins.concat(),\n", ")\n", "tof_data.hist(tof=100).plot()" ] @@ -290,104 +289,68 @@ "## Group by scattering angle\n", "\n", "The above pipeline focuses the data by merging all instrument pixels to produce a 1d d-spacing curve.\n", - "If instead we want to group into $2\\theta$ bins, we can alter the pipeline by replacing the focussing step:" + "If instead we want to group into $2\\theta$ bins, we can alter the pipeline parameters by adding some binning in $2\\theta$:" ] }, { "cell_type": "code", "execution_count": null, - "id": "d24c8bf8-ad59-4211-ae6a-3ee29b0556a3", + "id": "a4b68853-a70b-42d6-a8cc-58c77e83eaec", "metadata": {}, "outputs": [], "source": [ - "from ess.powder.grouping import group_by_two_theta, merge_all_pixels\n", - "\n", - "grouping_providers = list(providers)\n", - "grouping_providers.remove(merge_all_pixels)\n", - "grouping_providers = (*grouping_providers, group_by_two_theta)" + "pipeline[TwoThetaBins] = sc.linspace(dim='two_theta', unit='deg', start=25.0, stop=90.0, num=17).to(unit='rad')" ] }, { "cell_type": "markdown", - "id": "3e83ed6c-46fc-4396-8331-3654993def94", - "metadata": {}, - "source": [ - "We also need to specify the grouping with a new parameter:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a4b68853-a70b-42d6-a8cc-58c77e83eaec", + "id": "2a8a0853-2829-49ab-b343-31deca3de31c", "metadata": {}, - "outputs": [], "source": [ - "params[TwoThetaBins] = sc.linspace(dim='two_theta', unit='deg', start=25.0, stop=90.0, num=16)" + "Compute and plot the result:" ] }, { "cell_type": "code", "execution_count": null, - "id": "ea65b6bd-b9fa-496d-bfa2-f459e33f75fc", + "id": "a63e037d-3cfc-4ac9-963c-1e2d0881d482", "metadata": {}, "outputs": [], "source": [ - "grouping_pipeline = sciline.Pipeline(\n", - " grouping_providers,\n", - " params=params\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "59ae5b45-8e11-4de1-bef5-c4b58df75804", - "metadata": {}, - "source": [ - "Inspect the graph to check that the new provider has been inserted:" + "grouped_dspacing = pipeline.compute(NormalizedByVanadium)\n", + "grouped_dspacing" ] }, { "cell_type": "code", "execution_count": null, - "id": "a63e037d-3cfc-4ac9-963c-1e2d0881d482", + "id": "7ca04947-b019-4814-9467-4d0519d8d384", "metadata": {}, "outputs": [], "source": [ - "grouped_dspacing = grouping_pipeline.get(DspacingHistogram)\n", - "grouped_dspacing.visualize(graph_attr={'rankdir': 'LR'})" + "angle = sc.midpoints(grouped_dspacing.coords['two_theta'])\n", + "sc.plot({\n", + " f'{angle[group].value:.3f} {angle[group].unit}': grouped_dspacing['two_theta', group].hist()\n", + " for group in range(2, 6)\n", + "})" ] }, { "cell_type": "markdown", - "id": "2a8a0853-2829-49ab-b343-31deca3de31c", + "id": "3de8061a-7c28-48e2-b569-a08fa09f1295", "metadata": {}, "source": [ - "Compute and plot the result:" + "Or we can view it as a 2D plot, which should display powder peaks as vertical bright lines:" ] }, { "cell_type": "code", "execution_count": null, - "id": "bc73276a-dadb-4c21-9954-e69bebc2ff90", + "id": "4c9ee765-0699-4ca2-9bd8-f8dea27a261c", "metadata": {}, "outputs": [], "source": [ - "grouped_result = grouped_dspacing.compute()\n", - "grouped_result" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7ca04947-b019-4814-9467-4d0519d8d384", - "metadata": {}, - "outputs": [], - "source": [ - "angle = sc.midpoints(grouped_result.coords['two_theta'])\n", - "sc.plot({\n", - " f'{angle[group].value:.3f} {angle[group].unit}': grouped_result['two_theta', group]\n", - " for group in range(2, 6)\n", - "})" + "grouped_dspacing.hist().plot()" ] } ], @@ -407,7 +370,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/src/ess/dream/data.py b/src/ess/dream/data.py index fe1b694d..e4cc8f3e 100644 --- a/src/ess/dream/data.py +++ b/src/ess/dream/data.py @@ -2,30 +2,30 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Data for tests and documentation with DREAM.""" -from pathlib import Path +# from pathlib import Path -_version = '1' +_version = "1" -__all__ = ['get_path'] +__all__ = ["get_path"] def _make_pooch(): import pooch return pooch.create( - path=pooch.os_cache('ess/dream'), - env='ESS_DATA_DIR', - base_url='https://public.esss.dk/groups/scipp/ess/dream/{version}/', + path=pooch.os_cache("ess/dream"), + env="ESS_DATA_DIR", + base_url="https://public.esss.dk/groups/scipp/ess/dream/{version}/", version=_version, registry={ - 'data_dream_with_sectors.csv.zip': 'md5:52ae6eb3705e5e54306a001bc0ae85d8', - 'data_dream0_new_hkl_Si_pwd.csv.zip': 'md5:d0ae518dc1b943bb817b3d18c354ed01', # noqa: E501 - 'DREAM_nexus_sorted-2023-12-07.nxs': 'md5:22824e14f6eb950d24a720b2a0e2cb66', - 'DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union.csv.zip': 'md5:33302d0506b36aab74003b8aed4664cc', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union_run2.csv.zip': 'md5:c7758682f978d162dcb91e47c79abb83', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip': 'md5:1e22917b2bb68b5cacfb506b72700a4d', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip': 'md5:e5addfc06768140c76533946433fa2ec', # noqa: E501 - 'DREAM_simple_pwd_workflow/data_dream_vanadium_inc_coh.csv.zip': 'md5:39d1a44e248b12966b26f7c2f6c602a2', # noqa: E501 + "data_dream_with_sectors.csv.zip": "md5:52ae6eb3705e5e54306a001bc0ae85d8", + "data_dream0_new_hkl_Si_pwd.csv.zip": "md5:d0ae518dc1b943bb817b3d18c354ed01", # noqa: E501 + "DREAM_nexus_sorted-2023-12-07.nxs": "md5:22824e14f6eb950d24a720b2a0e2cb66", + "DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union.csv.zip": "md5:33302d0506b36aab74003b8aed4664cc", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_diamond_vana_container_sample_union_run2.csv.zip": "md5:c7758682f978d162dcb91e47c79abb83", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip": "md5:1e22917b2bb68b5cacfb506b72700a4d", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip": "md5:e5addfc06768140c76533946433fa2ec", # noqa: E501 + "DREAM_simple_pwd_workflow/data_dream_vanadium_inc_coh.csv.zip": "md5:39d1a44e248b12966b26f7c2f6c602a2", # noqa: E501 }, ) @@ -33,7 +33,7 @@ def _make_pooch(): _pooch = _make_pooch() -def get_path(name: str, unzip: bool = False) -> Path: +def get_path(name: str, unzip: bool = False) -> str: """ Return the path to a data file bundled with ess.dream. @@ -42,10 +42,10 @@ def get_path(name: str, unzip: bool = False) -> Path: """ import pooch - return Path(_pooch.fetch(name, processor=pooch.Unzip() if unzip else None)) + return _pooch.fetch(name, processor=pooch.Unzip() if unzip else None) -def simulated_diamond_sample() -> Path: +def simulated_diamond_sample() -> str: """Path to a GEANT4 CSV file for a diamond sample. **Sample**: @@ -70,12 +70,12 @@ def simulated_diamond_sample() -> Path: absorption of 36.73 1/m """ return get_path( - 'DREAM_simple_pwd_workflow/' - 'data_dream_diamond_vana_container_sample_union.csv.zip' + "DREAM_simple_pwd_workflow/" + "data_dream_diamond_vana_container_sample_union.csv.zip" ) -def simulated_vanadium_sample() -> Path: +def simulated_vanadium_sample() -> str: """Path to a GEANT4 CSV file for a vanadium sample. Contains both coherent and incoherent scattering. @@ -90,10 +90,10 @@ def simulated_vanadium_sample() -> Path: unit cell volume=27.66 angstrom^3 absorption of 36.73 1/m """ - return get_path('DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip') + return get_path("DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip") -def simulated_vanadium_sample_incoherent() -> Path: +def simulated_vanadium_sample_incoherent() -> str: """Path to a GEANT4 CSV file for a vanadium sample with only incoherent scattering. **Sample**: @@ -102,10 +102,10 @@ def simulated_vanadium_sample_incoherent() -> Path: - vertical dimension of sample (along y)=0.01 m - packing factor=1 """ - return get_path('DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip') + return get_path("DREAM_simple_pwd_workflow/data_dream_vanadium.csv.zip") -def simulated_empty_can() -> Path: +def simulated_empty_can() -> str: """Path to a GEANT4 CSV file for an empty can measurement. **Container**: @@ -119,5 +119,5 @@ def simulated_empty_can() -> Path: absorption of 36.73 1/m """ return get_path( - 'DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip' + "DREAM_simple_pwd_workflow/data_dream_vana_container_sample_union.csv.zip" ) diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index eb602a2c..860393fe 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -10,7 +10,7 @@ import scipp as sc from ess.powder.types import ( DetectorName, - FilePath, + Filename, RawDetector, RawDetectorData, RunType, @@ -27,9 +27,7 @@ class AllRawDetectors(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Raw data for all detectors.""" -def load_geant4_csv( - file_path: Union[FilePath[RunType], str, StringIO, BytesIO], -) -> AllRawDetectors[RunType]: +def load_geant4_csv(file_path: Filename[RunType]) -> AllRawDetectors[RunType]: """Load a GEANT4 CSV file for DREAM. Parameters @@ -54,7 +52,7 @@ def load_geant4_csv( detectors = _group(detectors) return AllRawDetectors[RunType]( - sc.DataGroup({'instrument': sc.DataGroup(detectors)}) + sc.DataGroup({"instrument": sc.DataGroup(detectors)}) ) @@ -62,7 +60,7 @@ def extract_geant4_detector( detectors: AllRawDetectors[RunType], detector_name: DetectorName ) -> RawDetector[RunType]: """Extract a single detector from a loaded GEANT4 simulation.""" - return RawDetector[RunType](detectors['instrument'][detector_name.name]) + return RawDetector[RunType](detectors["instrument"][detector_name.name]) def extract_geant4_detector_data( @@ -72,42 +70,40 @@ def extract_geant4_detector_data( return RawDetectorData[RunType](extract_detector_data(detector)) -def _load_raw_events( - file_path: Union[str, os.PathLike, StringIO, BytesIO], -) -> sc.DataArray: +def _load_raw_events(file_path: str) -> sc.DataArray: table = sc.io.load_csv( - file_path, sep='\t', header_parser='bracket', data_columns=[] + file_path, sep="\t", header_parser="bracket", data_columns=[] ) - table = table.rename_dims(row='event') + table = table.rename_dims(row="event") return sc.DataArray( - sc.ones(sizes=table.sizes, with_variances=True, unit='counts'), + sc.ones(sizes=table.sizes, with_variances=True, unit="counts"), coords=table.coords, ) def _adjust_coords(da: sc.DataArray) -> None: - da.coords['wavelength'] = da.coords.pop('lambda') - da.coords['position'] = sc.spatial.as_vectors( - da.coords.pop('x_pos'), da.coords.pop('y_pos'), da.coords.pop('z_pos') + da.coords["wavelength"] = da.coords.pop("lambda") + da.coords["position"] = sc.spatial.as_vectors( + da.coords.pop("x_pos"), da.coords.pop("y_pos"), da.coords.pop("z_pos") ) def _group(detectors: Dict[str, sc.DataArray]) -> Dict[str, sc.DataGroup]: - elements = ('module', 'segment', 'counter', 'wire', 'strip') + elements = ("module", "segment", "counter", "wire", "strip") def group(key: str, da: sc.DataArray) -> sc.DataArray: - if key in ['high_resolution', 'sans']: + if key in ["high_resolution", "sans"]: # Only the HR and SANS detectors have sectors. - return da.group('sector', *elements) + return da.group("sector", *elements) res = da.group(*elements) - res.bins.coords.pop('sector', None) + res.bins.coords.pop("sector", None) return res return {key: sc.DataGroup(events=group(key, da)) for key, da in detectors.items()} def _split_detectors( - data: sc.DataArray, detector_id_name: str = 'det ID' + data: sc.DataArray, detector_id_name: str = "det ID" ) -> Dict[str, sc.DataArray]: groups = data.group( sc.concat( @@ -124,15 +120,15 @@ def _split_detectors( if ( mantle := _extract_detector(groups, detector_id_name, MANTLE_DETECTOR_ID) ) is not None: - detectors['mantle'] = mantle.copy() + detectors["mantle"] = mantle.copy() if ( high_res := _extract_detector(groups, detector_id_name, HIGH_RES_DETECTOR_ID) ) is not None: - detectors['high_resolution'] = high_res.copy() + detectors["high_resolution"] = high_res.copy() if ( sans := _extract_detector(groups, detector_id_name, SANS_DETECTOR_ID) ) is not None: - detectors['sans'] = sans.copy() + detectors["sans"] = sans.copy() endcaps_list = [ det @@ -143,13 +139,13 @@ def _split_detectors( endcaps = sc.concat(endcaps_list, data.dim) endcaps = endcaps.bin( z_pos=sc.array( - dims=['z_pos'], + dims=["z_pos"], values=[-np.inf, 0.0, np.inf], - unit=endcaps.coords['z_pos'].unit, + unit=endcaps.coords["z_pos"].unit, ) ) - detectors['endcap_backward'] = endcaps[0].bins.concat().value.copy() - detectors['endcap_forward'] = endcaps[1].bins.concat().value.copy() + detectors["endcap_backward"] = endcaps[0].bins.concat().value.copy() + detectors["endcap_forward"] = endcaps[1].bins.concat().value.copy() return detectors diff --git a/src/ess/powder/__init__.py b/src/ess/powder/__init__.py index bdf9af1c..37a34832 100644 --- a/src/ess/powder/__init__.py +++ b/src/ess/powder/__init__.py @@ -24,7 +24,7 @@ del importlib providers = ( - *conversion.providers_with_calibration, + *conversion.providers_with_positions, *correction.providers, *filtering.providers, *grouping.providers, diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index f9d48c36..a9bed1d2 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -172,91 +172,67 @@ def powder_coordinate_transformation_graph() -> ElasticCoordTransformGraph: ) -def to_wavelength_with_positions( - data: PixelMaskedData[RunType], - graph: ElasticCoordTransformGraph, -) -> WavelengthData[RunType]: - """ - Transform coordinates to wavelength using detector positions. - """ - return WavelengthData[RunType]( - data.transform_coords("wavelength", graph=graph, keep_intermediate=False) - ) - - -def to_twotheta_with_positions( - data: WavelengthMaskedData[RunType], - graph: ElasticCoordTransformGraph, -) -> TwoThetaData[RunType]: - """ - Transform coordinates to two-theta using detector positions. - """ - return TwoThetaData[RunType]( - data.transform_coords("two_theta", graph=graph, keep_intermediate=False) - ) - - -def to_dspacing_with_positions( - data: TwoThetaMaskedData[RunType], - *, - sample: Optional[RawSample[RunType]] = None, - source: Optional[RawSource] = None, -) -> DspacingData[RunType]: - """ - Transform coordinates to d-spacing using detector positions. - - Computes d-spacing from time-of-flight stored in `data`. - - Attention - --------- - `data` may have a wavelength coordinate and dimension, - but those are discarded. - Only the stored time-of-flight is used, that is, any modifications to - the wavelength coordinate after it was computed from time-of-flight are lost. - - Raises - ------ - KeyError - If `data` does not contain a 'tof' coordinate. - - Parameters - ---------- - data: - Input data in tof or wavelength dimension. - Must have a tof coordinate. - sample: - Sample data with a position. - If not given, ``data`` must contain a 'sample_position' coordinate. - source: - Source data with a position. - If not given, ``data`` must contain a 'source_position' coordinate. - - Returns - ------- - : - A DataArray with the same data as the input and a 'dspacing' coordinate. - """ - graph = { - **scn.conversion.graph.beamline.beamline(scatter=True), - **scn.conversion.graph.tof.elastic_dspacing("tof"), - } - if sample is not None: - graph["sample_position"] = lambda: sample["position"] - if source is not None: - graph["source_position"] = lambda: source["position"] - - out = _restore_tof_if_in_wavelength(data) - out = out.transform_coords("dspacing", graph=graph, keep_intermediate=False) - # Add coords to ensure the result is the same whether sample or source are - # coords in the input or separate function arguments. - if sample is not None: - out.coords["sample_position"] = sample["position"] - out.coords.set_aligned("sample_position", False) - if source is not None: - out.coords["source_position"] = source["position"] - out.coords.set_aligned("source_position", False) - - return DspacingData[RunType](out) +# def to_dspacing_with_positions( +# data: TwoThetaMaskedData[RunType], +# *, +# sample: Optional[RawSample[RunType]] = None, +# source: Optional[RawSource] = None, +# ) -> DspacingData[RunType]: +# """ +# Transform coordinates to d-spacing using detector positions. + +# Computes d-spacing from time-of-flight stored in `data`. + +# Attention +# --------- +# `data` may have a wavelength coordinate and dimension, +# but those are discarded. +# Only the stored time-of-flight is used, that is, any modifications to +# the wavelength coordinate after it was computed from time-of-flight are lost. + +# Raises +# ------ +# KeyError +# If `data` does not contain a 'tof' coordinate. + +# Parameters +# ---------- +# data: +# Input data in tof or wavelength dimension. +# Must have a tof coordinate. +# sample: +# Sample data with a position. +# If not given, ``data`` must contain a 'sample_position' coordinate. +# source: +# Source data with a position. +# If not given, ``data`` must contain a 'source_position' coordinate. + +# Returns +# ------- +# : +# A DataArray with the same data as the input and a 'dspacing' coordinate. +# """ +# graph = { +# **scn.conversion.graph.beamline.beamline(scatter=True), +# **scn.conversion.graph.tof.elastic_dspacing("tof"), +# } +# if sample is not None: +# graph["sample_position"] = lambda: sample["position"] +# if source is not None: +# graph["source_position"] = lambda: source["position"] + +# out = _restore_tof_if_in_wavelength(data) +# out = out.transform_coords("dspacing", graph=graph, keep_intermediate=False) +# # Add coords to ensure the result is the same whether sample or source are +# # coords in the input or separate function arguments. +# if sample is not None: +# out.coords["sample_position"] = sample["position"] +# out.coords.set_aligned("sample_position", False) +# if source is not None: +# out.coords["source_position"] = source["position"] +# out.coords.set_aligned("source_position", False) + +# return DspacingData[RunType](out) def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: @@ -290,7 +266,4 @@ def add_scattering_coordinates( providers_with_positions = ( powder_coordinate_transformation_graph, add_scattering_coordinates, - # to_wavelength_with_positions, - # to_twotheta_with_positions, - # to_dspacing_with_positions, ) diff --git a/src/ess/powder/correction.py b/src/ess/powder/correction.py index d950f993..6ccbbe89 100644 --- a/src/ess/powder/correction.py +++ b/src/ess/powder/correction.py @@ -56,9 +56,9 @@ def normalize_by_monitor( : `data` normalized by a monitor. """ - if 'wavelength' not in monitor.coords: + if "wavelength" not in monitor.coords: monitor = monitor.transform_coords( - 'wavelength', + "wavelength", graph={**beamline.beamline(scatter=False), **tof.elastic("tof")}, keep_inputs=False, keep_intermediate=False, @@ -73,14 +73,13 @@ def normalize_by_monitor( "ess.powder.smoothing.lowpass with %s.", smooth_args, ) - monitor = lowpass(monitor, dim='wavelength', **smooth_args) - return data.bins / sc.lookup(func=monitor, dim='wavelength') + monitor = lowpass(monitor, dim="wavelength", **smooth_args) + return data.bins / sc.lookup(func=monitor, dim="wavelength") def normalize_by_vanadium( data: FocussedData[SampleRun], vanadium: FocussedData[VanadiumRun], - edges: DspacingBins, uncertainty_broadcast_mode: Optional[UncertaintyBroadcastMode] = None, ) -> NormalizedByVanadium: """ @@ -92,24 +91,17 @@ def normalize_by_vanadium( Sample data. vanadium: Vanadium data. - edges: - `vanadium` is histogrammed into these bins before dividing the data by it. uncertainty_broadcast_mode: Choose how uncertainties of vanadium are broadcast to the sample data. Defaults to ``UncertaintyBroadcastMode.fail``. - - Returns - ------- - : - `data` normalized by `vanadium`. """ vanadium = broadcast_uncertainties(vanadium, uncertainty_broadcast_mode) - norm = sc.lookup(vanadium.hist({edges.dim: edges}), dim=edges.dim) + norm = vanadium.hist() # Converting to unit 'one' because the division might produce a unit # with a large scale if the proton charges in data and vanadium were # measured with different units. - return (data.bins / norm).to(unit='one', copy=False) + return NormalizedByVanadium((data / norm).to(unit="one", copy=False)) def normalize_by_proton_charge( @@ -156,22 +148,22 @@ def merge_calibration(*, into: sc.DataArray, calibration: sc.Dataset) -> sc.Data dim = calibration.dim if not sc.identical(into.coords[dim], calibration.coords[dim]): raise ValueError( - f'Coordinate {dim} of calibration and target dataset do not agree.' + f"Coordinate {dim} of calibration and target dataset do not agree." ) out = into.copy(deep=False) - for name in ('difa', 'difc', 'tzero'): + for name in ("difa", "difc", "tzero"): if name in out.coords: raise ValueError( f"Cannot add calibration parameter '{name}' to data, " "there already is metadata with the same name." ) out.coords[name] = calibration[name].data - if 'calibration' in out.masks: + if "calibration" in out.masks: raise ValueError( "Cannot add calibration mask 'calibration' tp data, " "there already is a mask with the same name." ) - out.masks['calibration'] = calibration['mask'].data + out.masks["calibration"] = calibration["mask"].data return out @@ -207,8 +199,8 @@ def apply_lorentz_correction(da: sc.DataArray) -> sc.DataArray: # The implementation is optimized under the assumption that two_theta # is small and dspacing and the data are large. out = _shallow_copy(da) - dspacing = event_or_outer_coord(da, 'dspacing') - two_theta = event_or_outer_coord(da, 'two_theta') + dspacing = event_or_outer_coord(da, "dspacing") + two_theta = event_or_outer_coord(da, "two_theta") theta = 0.5 * two_theta d4 = dspacing.broadcast(sizes=out.sizes) ** 4 diff --git a/src/ess/powder/external/powgen/beamline.py b/src/ess/powder/external/powgen/beamline.py index b17c1359..b8dac938 100644 --- a/src/ess/powder/external/powgen/beamline.py +++ b/src/ess/powder/external/powgen/beamline.py @@ -29,24 +29,24 @@ def map_detector_to_spectrum( `data` with 'detector' coord and dim replaced by 'spectrum'. """ if not sc.identical( - data.coords['detector'].to( - dtype=detector_info.coords['detector'].dtype, copy=False + data.coords["detector"].to( + dtype=detector_info.coords["detector"].dtype, copy=False ), - detector_info.coords['detector'], + detector_info.coords["detector"], ): raise sc.CoordError( "The 'detector' coords of `data` and `detector_info` do not match." ) out = data.copy(deep=False) - del out.coords['detector'] + del out.coords["detector"] # Add 1 because spectrum numbers in the data start at 1 but # detector_info contains spectrum indices which start at 0. - out.coords['spectrum'] = detector_info.coords['spectrum'] + sc.index( - 1, dtype=detector_info.coords['spectrum'].dtype + out.coords["spectrum"] = detector_info.coords["spectrum"] + sc.index( + 1, dtype=detector_info.coords["spectrum"].dtype ) - return out.rename_dims({'detector': 'spectrum'}) + return out.rename_dims({"detector": "spectrum"}) def preprocess_calibration_data( @@ -63,7 +63,7 @@ def preprocess_calibration_data( def powgen_detector_dimensions() -> DetectorDimensions: """Dimensions used by POWGEN detectors.""" - return DetectorDimensions(('spectrum',)) + return DetectorDimensions({"bank": 23, "column": 154, "row": 7}) providers = (preprocess_calibration_data, powgen_detector_dimensions) diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 7d2e8ead..7407aed6 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -8,6 +8,7 @@ from ...types import ( AccumulatedProtonCharge, CalibrationFilename, + DetectorDimensions, Filename, ProtonCharge, RawCalibrationData, @@ -106,12 +107,14 @@ def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData: return RawCalibrationData(sc.io.load_hdf5(filename)) -def extract_raw_data(dg: RawDataAndMetadata[RunType]) -> RawDetectorData[RunType]: +def extract_raw_data( + dg: RawDataAndMetadata[RunType], sizes: DetectorDimensions +) -> RawDetectorData[RunType]: """Return the events from a loaded data group.""" out = dg["data"].squeeze() del out.coords["tof"] # out = out.fold(dim="spectrum", sizes={"column": 154, "row": 7, "bank": 23}) - out = out.fold(dim="spectrum", sizes={"bank": 23, "column": 154, "row": 7}) + out = out.fold(dim="spectrum", sizes=sizes) # out.bins.coords["position"] = sc.bins_like(out, out.coords["position"]) # out.bins.coords["spectrum"] = sc.bins_like(out, out.coords["spectrum"]) # raw_events = out.bins.constituents["data"].copy() diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index 94d459a4..e71a1cf4 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -1,6 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Grouping and merging of pixels / voxels.""" +from typing import Optional from scippneutron.conversion.graph import beamline @@ -10,59 +11,32 @@ DspacingData, DspacingHistogram, FocussedData, + MaskedData, NormalizedByVanadium, RunType, TwoThetaBins, ) -def group_by_two_theta( - data: DspacingData[RunType], edges: TwoThetaBins +def focus_data( + data: MaskedData[RunType], + detector_dims: DetectorDimensions, + dspacing_bins: DspacingBins, + twotheta_bins: Optional[TwoThetaBins] = None, ) -> FocussedData[RunType]: - """ - Group data into two_theta bins. - - Parameters - ---------- - data: - Input data array with events. Must contain a coord called 'two_theta' - or coords that can be used to compute it. - edges: - Bin edges in two_theta. `data` is grouped into those bins. - - Returns - ------- - : - `data` grouped into two_theta bins. - """ - out = data.transform_coords('two_theta', graph=beamline.beamline(scatter=True)) - return FocussedData[RunType]( - out.bin(two_theta=edges.to(unit=out.coords['two_theta'].unit, copy=False)) - ) + bins = {} + if twotheta_bins is not None: + bins["two_theta"] = twotheta_bins + bins[dspacing_bins.dim] = dspacing_bins + if twotheta_bins is None: + # In this case merge data from all pixels + # Put the dims into the same order as in the data. + # See https://github.com/scipp/scipp/issues/3408 + to_concat = tuple(dim for dim in data.dims if dim in detector_dims) + data = data.bins.concat(to_concat) -def merge_all_pixels( - data: DspacingData[RunType], dims: DetectorDimensions -) -> FocussedData[RunType]: - """Combine all pixels (spectra) of the detector. - - Parameters - ---------- - data: - Input data with pixel dimensions. - dims: - The dimensions to reduce over. - Corresponds to the pixel dimensions of the detector. - - Returns - ------- - : - The input without pixel dimensions. - """ - # Put the dims into the same order as in the data. - # See https://github.com/scipp/scipp/issues/3408 - to_concat = tuple(dim for dim in data.dims if dim in dims) - return FocussedData(data.bins.concat(to_concat)) + return FocussedData[RunType](data.bin(**bins)) def finalize_histogram( @@ -87,5 +61,5 @@ def finalize_histogram( return DspacingHistogram(data.hist(dspacing=edges)) -providers = (merge_all_pixels, finalize_histogram) +providers = (finalize_histogram, focus_data) """Sciline providers for grouping pixels.""" diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 9cab812d..e755e13b 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -9,7 +9,7 @@ from enum import Enum from pathlib import Path -from typing import Any, Callable, NewType, TypeVar +from typing import Any, Callable, Dict, NewType, TypeVar import sciline import scipp as sc @@ -97,8 +97,11 @@ class DataWithScatteringCoordinates(sciline.Scope[RunType, sc.DataArray], sc.Dat d-spacing.""" -class DetectorDimensions(sciline.Scope[DetectorName, tuple[str, ...]], tuple[str, ...]): - """Logical detector dimensions.""" +# class DetectorDimensions(sciline.Scope[DetectorName, tuple[str, ...]], tuple[str, ...]): +# """Logical detector dimensions.""" + +DetectorDimensions = NewType("DetectorDimensions", Dict[str, int]) +"""Logical detector dimensions.""" class DspacingData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): From f3924a13fa54f18d92012e31867bb900d995127c Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Fri, 17 May 2024 19:51:29 +0200 Subject: [PATCH 05/26] also make dream workflow work --- .../dream/dream-data-reduction.ipynb | 98 +++++++++++++------ src/ess/dream/instrument_view.py | 32 +++--- src/ess/dream/io/geant4.py | 38 ++++++- src/ess/powder/external/powgen/data.py | 13 +-- src/ess/powder/filtering.py | 11 ++- src/ess/powder/grouping.py | 2 +- src/ess/powder/masking.py | 1 - src/ess/powder/types.py | 17 +++- 8 files changed, 151 insertions(+), 61 deletions(-) diff --git a/docs/user-guide/dream/dream-data-reduction.ipynb b/docs/user-guide/dream/dream-data-reduction.ipynb index 06cd8586..5281dd41 100644 --- a/docs/user-guide/dream/dream-data-reduction.ipynb +++ b/docs/user-guide/dream/dream-data-reduction.ipynb @@ -36,27 +36,37 @@ " *dream.providers,\n", ")\n", "params = {\n", - " FilePath[SampleRun]: dream.data.simulated_diamond_sample(),\n", - " FilePath[VanadiumRun]: dream.data.simulated_vanadium_sample(),\n", - " FilePath[EmptyCanRun]: dream.data.simulated_empty_can(),\n", + " Filename[SampleRun]: dream.data.simulated_diamond_sample(),\n", + " Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(),\n", + " Filename[EmptyCanRun]: dream.data.simulated_empty_can(),\n", " DetectorName: DetectorName('mantle'),\n", "\n", " # The upper bounds mode is not yet implemented.\n", " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", "\n", - " # Crop data to this range in time-of-flight\n", - " ValidTofRange: sc.array(dims=['tof'], values=[0.0, 86_000.0], unit='us'),\n", - "\n", " # Edges for binning in d-spacing\n", - " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 200, unit='angstrom'),\n", + " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'),\n", + "\n", + " # Mask in time-of-flight\n", + " TofMask: lambda x: (x < sc.scalar(0., unit='ns')) | (x > sc.scalar(86e6, unit='ns')),\n", + " # Mask in wavelength\n", + " WavelengthMask: lambda x: (x > sc.scalar(2.0, unit='angstrom')) & (x < sc.scalar(2.7, unit='angstrom'))\n", + "}\n", + "\n", + "# Not available in simulated data\n", + "sample = sc.DataGroup(position=sc.vector([0., 0., 0.], unit='mm'))\n", + "params[RawSample[SampleRun]] = sample\n", + "params[RawSample[VanadiumRun]] = sample\n", + "\n", + "source = sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit='mm'))\n", + "params[RawSource[SampleRun]] = source\n", + "params[RawSource[VanadiumRun]] = source\n", "\n", - " # Not available in simulated data\n", - " RawSample[SampleRun]: sc.DataGroup(position=sc.vector([0., 0., 0.], unit='mm')),\n", - " RawSample[VanadiumRun]: sc.DataGroup(position=sc.vector([0., 0., 0.], unit='mm')),\n", - " RawSource: sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit='mm')),\n", - " AccumulatedProtonCharge[SampleRun]: sc.scalar(1.0, unit='µAh'),\n", - " AccumulatedProtonCharge[VanadiumRun]: sc.scalar(1.0, unit='µAh'),\n", - "}" + "charge = sc.scalar(1.0, unit='µAh')\n", + "params[AccumulatedProtonCharge[SampleRun]] = charge\n", + "params[AccumulatedProtonCharge[VanadiumRun]] = charge\n", + "\n", + "pixel_masks = []" ] }, { @@ -66,7 +76,9 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline = sciline.Pipeline(providers, params=params)" + "pipeline = sciline.Pipeline(providers, params=params)\n", + "\n", + "pipeline.set_param_series(PixelMaskFilename, pixel_masks)" ] }, { @@ -81,63 +93,89 @@ { "cell_type": "code", "execution_count": null, - "id": "3dd37911-7173-413f-a30b-362da8e7f326", + "id": "93f2ba0d-f256-4b64-b8fa-42cd1081cc41", "metadata": {}, "outputs": [], "source": [ - "from ess.powder.conversion import to_dspacing_with_positions\n", - "pipeline.insert(to_dspacing_with_positions)" + "pipeline.visualize(NormalizedByVanadium, graph_attr={'rankdir': 'LR'})" ] }, { "cell_type": "code", "execution_count": null, - "id": "93f2ba0d-f256-4b64-b8fa-42cd1081cc41", + "id": "6a17cf1d-0407-41dd-8a84-0975371ac70a", "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(DspacingHistogram, graph_attr={'rankdir': 'LR'})" + "result = pipeline.compute(NormalizedByVanadium)\n", + "result" ] }, { "cell_type": "code", "execution_count": null, - "id": "6a17cf1d-0407-41dd-8a84-0975371ac70a", + "id": "527623d2-0eee-456d-a4dd-764df8618ee4", "metadata": {}, "outputs": [], "source": [ - "result = pipeline.compute(DspacingHistogram)" + "dspacing_histogram = result.hist()\n", + "dspacing_histogram.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "39dc0ea3-b30d-4a88-878c-86cf0e8e76e9", + "metadata": {}, + "source": [ + "We can now save the result to disk:" ] }, { "cell_type": "code", "execution_count": null, - "id": "c05859ef-753b-4ea6-8d90-a17aa4a58c8f", + "id": "f37eba41-6178-4bcb-bf3c-fa177e6112b0", "metadata": {}, "outputs": [], "source": [ - "result" + "dspacing_histogram.coords['dspacing'] = sc.midpoints(dspacing_histogram.coords['dspacing'])\n", + "scn.io.save_xye('dspacing.xye', dspacing_histogram)" + ] + }, + { + "cell_type": "markdown", + "id": "0898265e-7edd-4093-9dfb-6e6b2f05c766", + "metadata": {}, + "source": [ + "### Compute intermediate results\n", + "\n", + "For inspection and debugging purposes, we can also compute intermediate results.\n", + "To avoid repeated computation (including costly loading of files), we can request multiple results at once, including the final result, if desired.\n", + "For example:" ] }, { "cell_type": "code", "execution_count": null, - "id": "4b2c884c-776e-4bb0-b517-10f5ddbf8aba", + "id": "4bd5a7d5-581f-4d04-bf11-955e55646199", "metadata": {}, "outputs": [], "source": [ - "result.plot()" + "intermediates = pipeline.compute((\n", + " DataWithScatteringCoordinates[SampleRun],\n", + " MaskedData[SampleRun],\n", + "))\n", + "\n", + "intermediates[DataWithScatteringCoordinates[SampleRun]].bins.coords.keys()" ] }, { "cell_type": "code", "execution_count": null, - "id": "27793efc-ecf2-4e1a-a73f-940ab652007c", + "id": "185cb6ee-8e55-43d0-9505-ebe2761785ed", "metadata": {}, "outputs": [], "source": [ - "result.coords['dspacing'] = sc.midpoints(result.coords['dspacing'])\n", - "scn.io.save_xye('dspacing.xye', result)" + "intermediates[MaskedData[SampleRun]].bins.concat().hist(two_theta=300, wavelength=300).plot(norm='log')" ] } ], @@ -157,7 +195,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.13" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/src/ess/dream/instrument_view.py b/src/ess/dream/instrument_view.py index fedb9bde..c48ffddd 100644 --- a/src/ess/dream/instrument_view.py +++ b/src/ess/dream/instrument_view.py @@ -50,7 +50,7 @@ def instrument_view( def _to_data_group(data: Union[sc.DataArray, sc.DataGroup, dict]) -> sc.DataGroup: if isinstance(data, sc.DataArray): - data = sc.DataGroup({data.name or 'data': data}) + data = sc.DataGroup({data.name or "data": data}) elif isinstance(data, dict): data = sc.DataGroup(data) return data @@ -61,8 +61,8 @@ def _pre_process(da: sc.DataArray, dim: str) -> sc.DataArray: dims = list(da.dims) if dim is not None: dims.remove(dim) - out = da.flatten(dims=dims, to='pixel') - sel = sc.isfinite(out.coords['position']) + out = da.flatten(dims=dims, to="pixel") + sel = sc.isfinite(out.coords["position"]) return out[sel] @@ -87,7 +87,7 @@ def __init__( if dim is not None: self.slider = SliceWidget(next(iter(self.data.values())), dims=[dim]) - self.slider.controls[dim]['slider'].layout = {'width': '600px'} + self.slider.controls[dim]["slider"].layout = {"width": "600px"} self.slider_node = pp.widget_node(self.slider) self.slice_nodes = { key: slice_dims(n, self.slider_node) @@ -101,8 +101,8 @@ def __init__( self.scatter = pp.scatter3d( to_scatter, - pos='position', - pixel_size=1.0 * sc.Unit('cm') if pixel_size is None else pixel_size, + pos="position", + pixel_size=1.0 * sc.Unit("cm") if pixel_size is None else pixel_size, **kwargs, ) @@ -114,10 +114,10 @@ def __init__( def _add_module_control(self): import ipywidgets as ipw - self.fig = self.scatter[0] - self.cutting_tool = self.scatter[1] + # self.fig = self.scatter[0] + self.cutting_tool = self.scatter.bottom_bar[0] self.artist_mapping = dict( - zip(self.data.keys(), self.fig.artists.keys(), strict=True) + zip(self.data.keys(), self.scatter.artists.keys(), strict=True) ) self.checkboxes = { key: ipw.Checkbox( @@ -137,10 +137,10 @@ def _add_module_control(self): ) for key, ch in self.checkboxes.items(): ch.key = key - ch.observe(self._check_visibility, names='value') - self.cutting_tool.cut_x.button.observe(self._check_visibility, names="value") - self.cutting_tool.cut_y.button.observe(self._check_visibility, names="value") - self.cutting_tool.cut_z.button.observe(self._check_visibility, names="value") + ch.observe(self._check_visibility, names="value") + # self.cutting_tool.cut_x.button.observe(self._check_visibility, names="value") + # self.cutting_tool.cut_y.button.observe(self._check_visibility, names="value") + # self.cutting_tool.cut_z.button.observe(self._check_visibility, names="value") self.children.insert(0, self.modules_widget) def _check_visibility(self, _): @@ -151,8 +151,8 @@ def _check_visibility(self, _): for name, ch in self.checkboxes.items(): key = self.artist_mapping[name] val = ch.value - self.fig.artists[key].points.visible = val + self.scatter.artists[key].points.visible = val for c in "xyz": - cut_nodes = getattr(self.cutting_tool, f'cut_{c}').select_nodes + cut_nodes = getattr(self.cutting_tool, f"cut_{c}").select_nodes if key in cut_nodes: - self.fig.artists[cut_nodes[key].id].points.visible = val + self.scatter.artists[cut_nodes[key].id].points.visible = val diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index 860393fe..de764313 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -13,7 +13,12 @@ Filename, RawDetector, RawDetectorData, + RawSample, + RawSource, + ReducibleDetectorData, RunType, + SamplePosition, + SourcePosition, ) from ess.reduce.nexus import extract_detector_data @@ -83,6 +88,7 @@ def _load_raw_events(file_path: str) -> sc.DataArray: def _adjust_coords(da: sc.DataArray) -> None: da.coords["wavelength"] = da.coords.pop("lambda") + da.coords["wavelength"].unit = "angstrom" da.coords["position"] = sc.spatial.as_vectors( da.coords.pop("x_pos"), da.coords.pop("y_pos"), da.coords.pop("z_pos") ) @@ -159,5 +165,35 @@ def _extract_detector( return events -providers = (extract_geant4_detector, extract_geant4_detector_data, load_geant4_csv) +def get_source_position( + raw_source: RawSource[RunType], +) -> SourcePosition[RunType]: + return SourcePosition[RunType](raw_source["position"]) + + +def get_sample_position( + raw_sample: RawSample[RunType], +) -> SamplePosition[RunType]: + return SamplePosition[RunType](raw_sample["position"]) + + +def patch_detector_data( + detector_data: RawDetectorData[RunType], + source_position: SourcePosition[RunType], + sample_position: SamplePosition[RunType], +) -> ReducibleDetectorData[RunType]: + out = detector_data.copy(deep=False) + out.coords["source_position"] = source_position + out.coords["sample_position"] = sample_position + return ReducibleDetectorData[RunType](out) + + +providers = ( + extract_geant4_detector, + extract_geant4_detector_data, + load_geant4_csv, + get_sample_position, + get_source_position, + patch_detector_data, +) """Geant4-providers for Sciline pipelines.""" diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 7407aed6..78caf38f 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -14,6 +14,7 @@ RawCalibrationData, RawDataAndMetadata, RawDetectorData, + ReducibleDetectorData, RunType, SampleRun, ) @@ -109,18 +110,14 @@ def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData: def extract_raw_data( dg: RawDataAndMetadata[RunType], sizes: DetectorDimensions -) -> RawDetectorData[RunType]: +) -> ReducibleDetectorData[RunType]: """Return the events from a loaded data group.""" + # Remove the tof binning and dimension, as it is not needed and it gets in the way + # of masking. out = dg["data"].squeeze() del out.coords["tof"] - # out = out.fold(dim="spectrum", sizes={"column": 154, "row": 7, "bank": 23}) out = out.fold(dim="spectrum", sizes=sizes) - # out.bins.coords["position"] = sc.bins_like(out, out.coords["position"]) - # out.bins.coords["spectrum"] = sc.bins_like(out, out.coords["spectrum"]) - # raw_events = out.bins.constituents["data"].copy() - # for c in ("gd_prtn_chrg", "sample_position", "source_position"): - # raw_events.coords[c] = out.coords[c] - return RawDetectorData[RunType](out) + return ReducibleDetectorData[RunType](out) def extract_detector_info(dg: RawDataAndMetadata[SampleRun]) -> DetectorInfo: diff --git a/src/ess/powder/filtering.py b/src/ess/powder/filtering.py index f7112326..eacac467 100644 --- a/src/ess/powder/filtering.py +++ b/src/ess/powder/filtering.py @@ -13,7 +13,14 @@ import scipp as sc from ._util import elem_dtype, elem_unit, event_or_outer_coord -from .types import FilteredData, RawDetectorData, RunType, TofCroppedData, ValidTofRange +from .types import ( + FilteredData, + RawDetectorData, + ReducibleDetectorData, + RunType, + TofCroppedData, + ValidTofRange, +) def _equivalent_bin_indices(a, b) -> bool: @@ -101,7 +108,7 @@ def crop_tof( # def filter_events(data: TofCroppedData[RunType]) -> FilteredData[RunType]: -def filter_events(data: RawDetectorData[RunType]) -> FilteredData[RunType]: +def filter_events(data: ReducibleDetectorData[RunType]) -> FilteredData[RunType]: """Remove bad events. Attention diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index e71a1cf4..9ee9492b 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -29,7 +29,7 @@ def focus_data( bins["two_theta"] = twotheta_bins bins[dspacing_bins.dim] = dspacing_bins - if twotheta_bins is None: + if (twotheta_bins is None) or ("two_theta" in data.bins.coords): # In this case merge data from all pixels # Put the dims into the same order as in the data. # See https://github.com/scipp/scipp/issues/3408 diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py index e592fff7..8e6f39b3 100644 --- a/src/ess/powder/masking.py +++ b/src/ess/powder/masking.py @@ -52,7 +52,6 @@ def apply_masks( ) -> MaskedData[RunType]: """ """ out = data.copy(deep=False) - masked_pixel_ids = {"pix_mask": sc.arange("spectrum", 1, 101)} if len(masked_pixel_ids) > 0: key = ( set(out.coords.keys()) & {"detector_number", "detector_id", "spectrum"} diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index e755e13b..d9083199 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -176,8 +176,21 @@ class RawSample(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Raw data from a loaded sample.""" -RawSource = NewType("RawSource", sc.DataGroup) -"""Raw data from a loaded neutron source.""" +class RawSource(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): + """Raw data from a loaded neutron source.""" + + +class ReducibleDetectorData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Data that is in a state ready for reduction.""" + + +class SamplePosition(sciline.Scope[RunType, sc.Variable], sc.Variable): + """Sample position""" + + +class SourcePosition(sciline.Scope[RunType, sc.Variable], sc.Variable): + """Source position""" + TofMask = NewType("TofMask", Callable) """""" From b54d91f9c8ef57bfc7fcb607d8427413742193f7 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Tue, 21 May 2024 14:35:40 +0200 Subject: [PATCH 06/26] remove param tables --- docs/conf.py | 92 +++++------ .../dream/dream-data-reduction.ipynb | 148 ++++++++++++++---- .../dream/dream-instrument-view.ipynb | 26 ++- .../POWGEN_data_reduction.ipynb | 85 +++++----- .../sns-instruments/preprocess_files.ipynb | 54 ++++--- src/ess/dream/beamline.py | 74 ++++----- src/ess/dream/io/geant4.py | 15 +- src/ess/powder/conversion.py | 20 ++- src/ess/powder/external/powgen/beamline.py | 18 ++- src/ess/powder/external/powgen/data.py | 7 +- src/ess/powder/filtering.py | 43 +---- src/ess/powder/grouping.py | 8 +- src/ess/powder/masking.py | 67 +------- src/ess/powder/types.py | 63 +++----- 14 files changed, 359 insertions(+), 361 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 17b45932..37ffb8fa 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -5,35 +5,35 @@ import sys from importlib.metadata import version as get_version -sys.path.insert(0, os.path.abspath('.')) +sys.path.insert(0, os.path.abspath(".")) # General information about the project. -project = 'ESSdiffraction' -copyright = '2024 Scipp contributors' -author = 'Scipp contributors' +project = "ESSdiffraction" +copyright = "2024 Scipp contributors" +author = "Scipp contributors" html_show_sourcelink = True extensions = [ - 'sphinx.ext.autodoc', - 'sphinx.ext.autosummary', - 'sphinx.ext.doctest', - 'sphinx.ext.githubpages', - 'sphinx.ext.intersphinx', - 'sphinx.ext.mathjax', - 'sphinx.ext.napoleon', - 'sphinx.ext.viewcode', - 'sphinx_autodoc_typehints', - 'sphinx_copybutton', - 'sphinx_design', - 'nbsphinx', - 'myst_parser', + "sphinx.ext.autodoc", + "sphinx.ext.autosummary", + "sphinx.ext.doctest", + "sphinx.ext.githubpages", + "sphinx.ext.intersphinx", + "sphinx.ext.mathjax", + "sphinx.ext.napoleon", + "sphinx.ext.viewcode", + "sphinx_autodoc_typehints", + "sphinx_copybutton", + "sphinx_design", + "nbsphinx", + "myst_parser", ] try: import sciline.sphinxext.domain_types # noqa: F401 - extensions.append('sciline.sphinxext.domain_types') + extensions.append("sciline.sphinxext.domain_types") except ModuleNotFoundError: pass @@ -56,13 +56,13 @@ myst_heading_anchors = 3 autodoc_type_aliases = { - 'array_like': 'array_like', + "array_like": "array_like", } intersphinx_mapping = { - 'python': ('https://docs.python.org/3', None), - 'numpy': ('https://numpy.org/doc/stable/', None), - 'scipp': ('https://scipp.github.io/', None), + "python": ("https://docs.python.org/3", None), + "numpy": ("https://numpy.org/doc/stable/", None), + "scipp": ("https://scipp.github.io/", None), } # autodocs includes everything, even irrelevant API internals. autosummary @@ -79,32 +79,32 @@ # objects without namespace: numpy "ndarray": "~numpy.ndarray", } -typehints_defaults = 'comma' +typehints_defaults = "comma" typehints_use_rtype = False -sciline_domain_types_prefix = 'ess.diffraction' +sciline_domain_types_prefix = "ess.diffraction" sciline_domain_types_aliases = { - 'scipp._scipp.core.DataArray': 'scipp.DataArray', - 'scipp._scipp.core.Dataset': 'scipp.Dataset', - 'scipp._scipp.core.DType': 'scipp.DType', - 'scipp._scipp.core.Unit': 'scipp.Unit', - 'scipp._scipp.core.Variable': 'scipp.Variable', - 'scipp.core.data_group.DataGroup': 'scipp.DataGroup', + "scipp._scipp.core.DataArray": "scipp.DataArray", + "scipp._scipp.core.Dataset": "scipp.Dataset", + "scipp._scipp.core.DType": "scipp.DType", + "scipp._scipp.core.Unit": "scipp.Unit", + "scipp._scipp.core.Variable": "scipp.Variable", + "scipp.core.data_group.DataGroup": "scipp.DataGroup", } # Add any paths that contain templates here, relative to this directory. -templates_path = ['_templates'] +templates_path = ["_templates"] # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: # -source_suffix = ['.rst', '.md'] -html_sourcelink_suffix = '' # Avoid .ipynb.txt extensions in sources +source_suffix = [".rst", ".md"] +html_sourcelink_suffix = "" # Avoid .ipynb.txt extensions in sources # The master toctree document. -master_doc = 'index' +master_doc = "index" # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -112,7 +112,7 @@ # release = get_version("essdiffraction") -version = ".".join(release.split('.')[:3]) # CalVer +version = ".".join(release.split(".")[:3]) # CalVer warning_is_error = True @@ -127,15 +127,15 @@ # directories to ignore when looking for source files. # This patterns also effect to html_static_path and html_extra_path exclude_patterns = [ - '_build', - 'Thumbs.db', - '.DS_Store', - '**.ipynb_checkpoints', - 'user-guide/sns-instruments/preprocess_files.ipynb', + "_build", + "Thumbs.db", + ".DS_Store", + "**.ipynb_checkpoints", + "user-guide/sns-instruments/preprocess_files.ipynb", ] # The name of the Pygments (syntax highlighting) style to use. -pygments_style = 'sphinx' +pygments_style = "sphinx" # If true, `todo` and `todoList` produce output, else they produce nothing. todo_include_todos = False @@ -200,14 +200,14 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['_static'] +html_static_path = ["_static"] html_css_files = [] html_js_files = ["anaconda-icon.js"] # -- Options for HTMLHelp output ------------------------------------------ # Output file base name for HTML help builder. -htmlhelp_basename = 'essdiffractiondoc' +htmlhelp_basename = "essdiffractiondoc" # -- Options for Matplotlib in notebooks ---------------------------------- @@ -223,7 +223,7 @@ # In addition, there is no need to make plots in doctest as the documentation # build already tests if those plots can be made. # So we simply disable plots in doctests. -doctest_global_setup = ''' +doctest_global_setup = """ import numpy as np try: @@ -240,7 +240,7 @@ def do_not_plot(*args, **kwargs): except ImportError: # Scipp is not needed by docs if it is not installed. pass -''' +""" # Using normalize whitespace because many __str__ functions in scipp produce # extraneous empty lines and it would look strange to include them in the docs. @@ -255,5 +255,5 @@ def do_not_plot(*args, **kwargs): linkcheck_ignore = [ # Specific lines in Github blobs cannot be found by linkcheck. - r'https?://github\.com/.*?/blob/[a-f0-9]+/.+?#', + r"https?://github\.com/.*?/blob/[a-f0-9]+/.+?#", ] diff --git a/docs/user-guide/dream/dream-data-reduction.ipynb b/docs/user-guide/dream/dream-data-reduction.ipynb index 5281dd41..2455e7f0 100644 --- a/docs/user-guide/dream/dream-data-reduction.ipynb +++ b/docs/user-guide/dream/dream-data-reduction.ipynb @@ -23,6 +23,17 @@ "from ess.dream.io.geant4 import providers as geant4_providers" ] }, + { + "cell_type": "markdown", + "id": "1252feab-12d2-46ac-bf74-70b32344473d", + "metadata": {}, + "source": [ + "## Define reduction parameters\n", + "\n", + "We define a dictionary containing the reduction parameters.\n", + "The keys are types defined in [essdiffraction.types](../generated/modules/ess.diffraction.types.rst)." + ] + }, { "cell_type": "code", "execution_count": null, @@ -30,43 +41,45 @@ "metadata": {}, "outputs": [], "source": [ - "providers = (\n", - " *geant4_providers,\n", - " *powder.providers,\n", - " *dream.providers,\n", - ")\n", "params = {\n", " Filename[SampleRun]: dream.data.simulated_diamond_sample(),\n", " Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(),\n", " Filename[EmptyCanRun]: dream.data.simulated_empty_can(),\n", - " DetectorName: DetectorName('mantle'),\n", - "\n", + " NeXusDetectorName: \"mantle\",\n", " # The upper bounds mode is not yet implemented.\n", " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", - "\n", " # Edges for binning in d-spacing\n", - " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'),\n", - "\n", + " DspacingBins: sc.linspace(\"dspacing\", 0.0, 2.3434, 201, unit=\"angstrom\"),\n", " # Mask in time-of-flight\n", - " TofMask: lambda x: (x < sc.scalar(0., unit='ns')) | (x > sc.scalar(86e6, unit='ns')),\n", + " TofMask: lambda x: (x < sc.scalar(0.0, unit=\"ns\"))\n", + " | (x > sc.scalar(86e6, unit=\"ns\")),\n", " # Mask in wavelength\n", - " WavelengthMask: lambda x: (x > sc.scalar(2.0, unit='angstrom')) & (x < sc.scalar(2.7, unit='angstrom'))\n", + " WavelengthMask: lambda x: (x > sc.scalar(2.0, unit=\"angstrom\"))\n", + " & (x < sc.scalar(2.7, unit=\"angstrom\")),\n", "}\n", "\n", "# Not available in simulated data\n", - "sample = sc.DataGroup(position=sc.vector([0., 0., 0.], unit='mm'))\n", + "sample = sc.DataGroup(position=sc.vector([0.0, 0.0, 0.0], unit=\"mm\"))\n", "params[RawSample[SampleRun]] = sample\n", "params[RawSample[VanadiumRun]] = sample\n", "\n", - "source = sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit='mm'))\n", + "source = sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit=\"mm\"))\n", "params[RawSource[SampleRun]] = source\n", "params[RawSource[VanadiumRun]] = source\n", "\n", - "charge = sc.scalar(1.0, unit='µAh')\n", + "charge = sc.scalar(1.0, unit=\"µAh\")\n", "params[AccumulatedProtonCharge[SampleRun]] = charge\n", - "params[AccumulatedProtonCharge[VanadiumRun]] = charge\n", + "params[AccumulatedProtonCharge[VanadiumRun]] = charge" + ] + }, + { + "cell_type": "markdown", + "id": "21cb87f2-4ff7-436e-b603-cc8f60c73e7a", + "metadata": {}, + "source": [ + "## Create pipeline using Sciline\n", "\n", - "pixel_masks = []" + "We use the basic providers available in `essdiffraction` as well as the specialised `powder` and `geant4` providers." ] }, { @@ -76,9 +89,12 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline = sciline.Pipeline(providers, params=params)\n", + "providers = (\n", + " *geant4_providers,\n", + " *powder.providers,\n", + ")\n", "\n", - "pipeline.set_param_series(PixelMaskFilename, pixel_masks)" + "pipeline = sciline.Pipeline(providers, params=params)" ] }, { @@ -86,8 +102,7 @@ "id": "21fb4492-e836-41d3-a2d4-9678df43b9f9", "metadata": {}, "source": [ - "We don't have calibration data yet.\n", - "So convert to d-spacing without calibration." + "We can visualize the graph for computing the final normalized result for intensity as a function of d-spacing:" ] }, { @@ -97,7 +112,15 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(NormalizedByVanadium, graph_attr={'rankdir': 'LR'})" + "pipeline.visualize(NormalizedByVanadium, graph_attr={\"rankdir\": \"LR\"})" + ] + }, + { + "cell_type": "markdown", + "id": "478f580e-273c-4a46-8dfc-9a2ed31b0cd3", + "metadata": {}, + "source": [ + "We then call `compute()` to compute the result:" ] }, { @@ -137,8 +160,10 @@ "metadata": {}, "outputs": [], "source": [ - "dspacing_histogram.coords['dspacing'] = sc.midpoints(dspacing_histogram.coords['dspacing'])\n", - "scn.io.save_xye('dspacing.xye', dspacing_histogram)" + "dspacing_histogram.coords[\"dspacing\"] = sc.midpoints(\n", + " dspacing_histogram.coords[\"dspacing\"]\n", + ")\n", + "scn.io.save_xye(\"dspacing.xye\", dspacing_histogram)" ] }, { @@ -146,7 +171,7 @@ "id": "0898265e-7edd-4093-9dfb-6e6b2f05c766", "metadata": {}, "source": [ - "### Compute intermediate results\n", + "## Compute intermediate results\n", "\n", "For inspection and debugging purposes, we can also compute intermediate results.\n", "To avoid repeated computation (including costly loading of files), we can request multiple results at once, including the final result, if desired.\n", @@ -160,10 +185,12 @@ "metadata": {}, "outputs": [], "source": [ - "intermediates = pipeline.compute((\n", - " DataWithScatteringCoordinates[SampleRun],\n", - " MaskedData[SampleRun],\n", - "))\n", + "intermediates = pipeline.compute(\n", + " (\n", + " DataWithScatteringCoordinates[SampleRun],\n", + " MaskedData[SampleRun],\n", + " )\n", + ")\n", "\n", "intermediates[DataWithScatteringCoordinates[SampleRun]].bins.coords.keys()" ] @@ -175,7 +202,68 @@ "metadata": {}, "outputs": [], "source": [ - "intermediates[MaskedData[SampleRun]].bins.concat().hist(two_theta=300, wavelength=300).plot(norm='log')" + "intermediates[MaskedData[SampleRun]].bins.concat().hist(\n", + " two_theta=300, wavelength=300\n", + ").plot(norm=\"log\")" + ] + }, + { + "cell_type": "markdown", + "id": "7f866ad4-5a0d-4c98-b5ab-a2436f97074d", + "metadata": {}, + "source": [ + "## Grouping by scattering angle" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc516311-1ab4-4eb7-9a0b-d83ff61e036b", + "metadata": {}, + "outputs": [], + "source": [ + "pipeline[TwoThetaBins] = sc.linspace(\n", + " dim=\"two_theta\", unit=\"rad\", start=0.8, stop=2.4, num=17\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b19ab368-e0c4-4a03-ac20-61a442e8826e", + "metadata": {}, + "outputs": [], + "source": [ + "grouped_dspacing = pipeline.compute(NormalizedByVanadium)\n", + "grouped_dspacing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc25ee9c-9395-481f-bd31-3ad139859686", + "metadata": {}, + "outputs": [], + "source": [ + "angle = sc.midpoints(grouped_dspacing.coords[\"two_theta\"])\n", + "sc.plot(\n", + " {\n", + " f\"{angle[group].value:.3f} {angle[group].unit}\": grouped_dspacing[\n", + " \"two_theta\", group\n", + " ].hist()\n", + " for group in range(2, 6)\n", + " }\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e414419-ccf4-420b-899d-a1fa4a6ce4ca", + "metadata": {}, + "outputs": [], + "source": [ + "grouped_dspacing.hist().plot(norm=\"log\")" ] } ], diff --git a/docs/user-guide/dream/dream-instrument-view.ipynb b/docs/user-guide/dream/dream-instrument-view.ipynb index 3a3765a2..e1bc9cff 100644 --- a/docs/user-guide/dream/dream-instrument-view.ipynb +++ b/docs/user-guide/dream/dream-instrument-view.ipynb @@ -43,19 +43,15 @@ "metadata": {}, "outputs": [], "source": [ - "dg = dream.io.load_geant4_csv(\n", - " dream.data.get_path('data_dream0_new_hkl_Si_pwd.csv.zip')\n", - ")\n", - "dg = dg['instrument'] # Extract the instrument data\n", + "dg = dream.io.load_geant4_csv(dream.data.get_path(\"data_dream0_new_hkl_Si_pwd.csv.zip\"))\n", + "dg = dg[\"instrument\"] # Extract the instrument data\n", "\n", "# Extract the events from nested data groups\n", - "dg = sc.DataGroup({\n", - " key: detector['events'] for key, detector in dg.items()\n", - "})\n", + "dg = sc.DataGroup({key: detector[\"events\"] for key, detector in dg.items()})\n", "\n", "# Construct the pixel positions from event positions\n", "for da in dg.values():\n", - " da.coords['position'] = da.bins.coords['position'].bins.mean()\n", + " da.coords[\"position\"] = da.bins.coords[\"position\"].bins.mean()\n", "\n", "dg" ] @@ -75,7 +71,7 @@ "outputs": [], "source": [ "# Only plot half of the pixels to reduce html docs size\n", - "dg = dg['counter', 0]" + "dg = dg[\"counter\", 0]" ] }, { @@ -108,7 +104,7 @@ }, "outputs": [], "source": [ - "tof_edges = sc.linspace('tof', 1.0e7, 1.0e8, 51, unit='ns', dtype=int)\n", + "tof_edges = sc.linspace(\"tof\", 1.0e7, 1.0e8, 51, unit=\"ns\", dtype=int)\n", "data = dg.hist(tof=tof_edges)" ] }, @@ -140,7 +136,7 @@ }, "outputs": [], "source": [ - "full_view = dream.instrument_view(data, dim='tof')\n", + "full_view = dream.instrument_view(data, dim=\"tof\")\n", "full_view" ] }, @@ -158,7 +154,7 @@ }, "outputs": [], "source": [ - "full_view[2].controls['tof']['slider'].value = 35" + "full_view[2].controls[\"tof\"][\"slider\"].value = 35" ] }, { @@ -192,7 +188,7 @@ }, "outputs": [], "source": [ - "mantle_view = dream.instrument_view(dg['mantle'].hist(wavelength=50), dim='wavelength')\n", + "mantle_view = dream.instrument_view(dg[\"mantle\"].hist(wavelength=50), dim=\"wavelength\")\n", "mantle_view" ] }, @@ -210,7 +206,7 @@ }, "outputs": [], "source": [ - "mantle_view[1].controls['wavelength']['slider'].value = 43" + "mantle_view[1].controls[\"wavelength\"][\"slider\"].value = 43" ] }, { @@ -231,7 +227,7 @@ "metadata": {}, "outputs": [], "source": [ - "dream.instrument_view(dg['endcap_backward']['module', 0].hist(tof=1))" + "dream.instrument_view(dg[\"endcap_backward\"][\"module\", 0].hist(tof=1))" ] } ], diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index 4b904c93..0742c9de 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -51,27 +51,22 @@ "outputs": [], "source": [ "params = {\n", + " NeXusDetectorName: \"powgen_detector\",\n", " # Input data\n", " Filename[SampleRun]: powgen.data.powgen_tutorial_sample_file(),\n", " Filename[VanadiumRun]: powgen.data.powgen_tutorial_vanadium_file(),\n", " CalibrationFilename: powgen.data.powgen_tutorial_calibration_file(),\n", - "\n", " # The upper bounds mode is not yet implemented.\n", " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", - "\n", - " # Crop data to this range in time-of-flight\n", - " ValidTofRange: sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us'),\n", - "\n", " # Edges for binning in d-spacing\n", - " DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'),\n", - "}\n", - "\n", - "# Mask in time-of-flight\n", - "params[TofMask] = lambda x: (x < sc.scalar(0., unit='us')) | (x > sc.scalar(16666.67, unit='us'))\n", - "# Mask in wavelength\n", - "params[WavelengthMask] = lambda x: (x > sc.scalar(0.18, unit='angstrom')) & (x < sc.scalar(0.21, unit='angstrom'))\n", - "# Pixel masks\n", - "pixel_masks = []" + " DspacingBins: sc.linspace(\"dspacing\", 0.0, 2.3434, 201, unit=\"angstrom\"),\n", + " # Mask in time-of-flight\n", + " TofMask: lambda x: (x < sc.scalar(0.0, unit=\"us\"))\n", + " | (x > sc.scalar(16666.67, unit=\"us\")),\n", + " # Mask in wavelength\n", + " WavelengthMask: lambda x: (x > sc.scalar(0.18, unit=\"angstrom\"))\n", + " & (x < sc.scalar(0.21, unit=\"angstrom\")),\n", + "}" ] }, { @@ -92,12 +87,7 @@ "outputs": [], "source": [ "providers = [*powder.providers, *powgen.providers]\n", - "pipeline = sciline.Pipeline(\n", - " providers,\n", - " params=params\n", - ")\n", - "\n", - "pipeline.set_param_series(PixelMaskFilename, pixel_masks)" + "pipeline = sciline.Pipeline(providers, params=params)" ] }, { @@ -119,7 +109,7 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(NormalizedByVanadium, graph_attr={'rankdir': 'LR'})" + "pipeline.visualize(NormalizedByVanadium, graph_attr={\"rankdir\": \"LR\"})" ] }, { @@ -174,12 +164,13 @@ "metadata": {}, "outputs": [], "source": [ - "def save_xye(reduced_data: NormalizedByVanadium,\n", - " out_filename: OutFilename,\n", + "def save_xye(\n", + " reduced_data: NormalizedByVanadium,\n", + " out_filename: OutFilename,\n", ") -> None:\n", " data = reduced_data.hist()\n", - " data.coords['dspacing'] = sc.midpoints(data.coords['dspacing'])\n", - " scn.io.save_xye(out_filename, data, coord='dspacing')" + " data.coords[\"dspacing\"] = sc.midpoints(data.coords[\"dspacing\"])\n", + " scn.io.save_xye(out_filename, data, coord=\"dspacing\")" ] }, { @@ -198,7 +189,7 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline[OutFilename] = 'reduced.xye'" + "pipeline[OutFilename] = \"reduced.xye\"" ] }, { @@ -239,12 +230,14 @@ "metadata": {}, "outputs": [], "source": [ - "results = pipeline.compute((\n", - " RawDetectorData[SampleRun],\n", - " MaskedData[SampleRun],\n", - " FilteredData[SampleRun],\n", - " FilteredData[VanadiumRun],\n", - "))" + "results = pipeline.compute(\n", + " (\n", + " RawDetectorData[SampleRun],\n", + " MaskedData[SampleRun],\n", + " FilteredData[SampleRun],\n", + " FilteredData[VanadiumRun],\n", + " )\n", + ")" ] }, { @@ -299,7 +292,19 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline[TwoThetaBins] = sc.linspace(dim='two_theta', unit='deg', start=25.0, stop=90.0, num=17).to(unit='rad')" + "pipeline[TwoThetaBins] = sc.linspace(\n", + " dim=\"two_theta\", unit=\"deg\", start=25.0, stop=90.0, num=17\n", + ").to(unit=\"rad\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e88eb36e-9377-4415-8123-a8916fd14793", + "metadata": {}, + "outputs": [], + "source": [ + "pipeline.visualize(NormalizedByVanadium, graph_attr={\"rankdir\": \"LR\"})" ] }, { @@ -328,11 +333,15 @@ "metadata": {}, "outputs": [], "source": [ - "angle = sc.midpoints(grouped_dspacing.coords['two_theta'])\n", - "sc.plot({\n", - " f'{angle[group].value:.3f} {angle[group].unit}': grouped_dspacing['two_theta', group].hist()\n", - " for group in range(2, 6)\n", - "})" + "angle = sc.midpoints(grouped_dspacing.coords[\"two_theta\"])\n", + "sc.plot(\n", + " {\n", + " f\"{angle[group].value:.3f} {angle[group].unit}\": grouped_dspacing[\n", + " \"two_theta\", group\n", + " ].hist()\n", + " for group in range(2, 6)\n", + " }\n", + ")" ] }, { diff --git a/docs/user-guide/sns-instruments/preprocess_files.ipynb b/docs/user-guide/sns-instruments/preprocess_files.ipynb index 4cd44b42..e1a8971c 100644 --- a/docs/user-guide/sns-instruments/preprocess_files.ipynb +++ b/docs/user-guide/sns-instruments/preprocess_files.ipynb @@ -44,7 +44,8 @@ " powgen.data.mantid_sample_file(),\n", " advanced_geometry=True,\n", " load_pulse_times=False,\n", - " mantid_args={\"LoadMonitors\": False})" + " mantid_args={\"LoadMonitors\": False},\n", + ")" ] }, { @@ -64,13 +65,15 @@ "metadata": {}, "outputs": [], "source": [ - "sample_data = sample['data']\n", - "sample_data.coords['gd_prtn_chrg'] = sample['gd_prtn_chrg']\n", - "sample_data.coords.set_aligned('gd_prtn_chrg', False)\n", - "sample_dg = sc.DataGroup({\n", - " 'data': sample_data,\n", - " 'detector_info': sample_data.coords.pop('detector_info').value\n", - "})" + "sample_data = sample[\"data\"]\n", + "sample_data.coords[\"gd_prtn_chrg\"] = sample[\"gd_prtn_chrg\"]\n", + "sample_data.coords.set_aligned(\"gd_prtn_chrg\", False)\n", + "sample_dg = sc.DataGroup(\n", + " {\n", + " \"data\": sample_data,\n", + " \"detector_info\": sample_data.coords.pop(\"detector_info\").value,\n", + " }\n", + ")" ] }, { @@ -90,7 +93,7 @@ "metadata": {}, "outputs": [], "source": [ - "sample_dg.save_hdf5('PG3_4844_event.h5')" + "sample_dg.save_hdf5(\"PG3_4844_event.h5\")" ] }, { @@ -100,7 +103,7 @@ "metadata": {}, "outputs": [], "source": [ - "sc.io.load_hdf5('PG3_4844_event.h5')" + "sc.io.load_hdf5(\"PG3_4844_event.h5\")" ] }, { @@ -122,7 +125,8 @@ " powgen.data.mantid_vanadium_file(),\n", " advanced_geometry=False,\n", " load_pulse_times=True,\n", - " mantid_args={\"LoadMonitors\": False})" + " mantid_args={\"LoadMonitors\": False},\n", + ")" ] }, { @@ -132,13 +136,15 @@ "metadata": {}, "outputs": [], "source": [ - "vana_data = vana['data']\n", - "vana_data.coords['gd_prtn_chrg'] = vana['gd_prtn_chrg']\n", - "vana_data.coords.set_aligned('gd_prtn_chrg', False)\n", - "vana_dg = sc.DataGroup({\n", - " 'data': vana_data,\n", - " 'proton_charge': vana['proton_charge'].rename(time='pulse_time')\n", - "})" + "vana_data = vana[\"data\"]\n", + "vana_data.coords[\"gd_prtn_chrg\"] = vana[\"gd_prtn_chrg\"]\n", + "vana_data.coords.set_aligned(\"gd_prtn_chrg\", False)\n", + "vana_dg = sc.DataGroup(\n", + " {\n", + " \"data\": vana_data,\n", + " \"proton_charge\": vana[\"proton_charge\"].rename(time=\"pulse_time\"),\n", + " }\n", + ")" ] }, { @@ -158,7 +164,7 @@ "metadata": {}, "outputs": [], "source": [ - "vana_dg['data'].bins.constituents['data']" + "vana_dg[\"data\"].bins.constituents[\"data\"]" ] }, { @@ -168,7 +174,7 @@ "metadata": {}, "outputs": [], "source": [ - "vana_dg.save_hdf5('PG3_4866_event.h5')" + "vana_dg.save_hdf5(\"PG3_4866_event.h5\")" ] }, { @@ -178,7 +184,7 @@ "metadata": {}, "outputs": [], "source": [ - "sc.io.load_hdf5('PG3_4866_event.h5')" + "sc.io.load_hdf5(\"PG3_4866_event.h5\")" ] }, { @@ -198,7 +204,7 @@ "source": [ "cal = load_calibration(\n", " powgen.data.mantid_calibration_file(),\n", - " instrument_filename='POWGEN_Definition_2011-02-25.xml',\n", + " instrument_filename=\"POWGEN_Definition_2011-02-25.xml\",\n", ")" ] }, @@ -219,7 +225,7 @@ "metadata": {}, "outputs": [], "source": [ - "cal.save_hdf5('PG3_FERNS_d4832_2011_08_24.h5')" + "cal.save_hdf5(\"PG3_FERNS_d4832_2011_08_24.h5\")" ] }, { @@ -229,7 +235,7 @@ "metadata": {}, "outputs": [], "source": [ - "sc.io.load_hdf5('PG3_FERNS_d4832_2011_08_24.h5')" + "sc.io.load_hdf5(\"PG3_FERNS_d4832_2011_08_24.h5\")" ] } ], diff --git a/src/ess/dream/beamline.py b/src/ess/dream/beamline.py index 2b0a5bb1..44f0424d 100644 --- a/src/ess/dream/beamline.py +++ b/src/ess/dream/beamline.py @@ -3,48 +3,34 @@ """Beamline tools for DREAM.""" -from ess.powder.types import DetectorDimensions, RawDetectorData, SampleRun - - -def dream_detector_dimensions(data: RawDetectorData[SampleRun]) -> DetectorDimensions: - """Logical dimensions used by a DREAM detector. - - The dimensions differ between simulated data loaded from GEANT4 CSV files - and measured data loaded from NeXus files. - The dimensions returned by this function match the dimensions found - in the ``data`` argument. - - Parameters - ---------- - data: - Dimensions are deduced based on this data. - - Returns - ------- - : - Logical dimensions used by the given DREAM detector. - """ - geant4_dims = { - 'module', - 'segment', - 'counter', - 'wire', - 'strip', - 'sector', - } - nexus_dims = { - 'wire', - 'mounting_unit', - 'cassette', - 'counter', - 'strip', - 'sector', - 'sumo_cass_ctr', - 'other', - } - dims = (geant4_dims | nexus_dims) & set(data.dims) - return DetectorDimensions(tuple(dims)) - - -providers = (dream_detector_dimensions,) +from ess.powder.types import NeXusDetectorDimensions, NeXusDetectorName + +DETECTOR_BANK_SHAPES_DAY1 = { + "endcap_backward": { + "strip": 16, + "wire": 16, + "module": 11, + "segment": 28, + "counter": 2, + }, + "endcap_forward": { + "strip": 16, + "wire": 16, + "module": 5, + "segment": 28, + "counter": 2, + }, + "mantle": {"wire": 32, "module": 5, "segment": 6, "strip": 256, "counter": 2}, + # TODO: missing "high_resolution" and "sans" detectors +} + + +def dream_detector_dimensions_day1( + detector_name: NeXusDetectorName, +) -> NeXusDetectorDimensions[NeXusDetectorName]: + """Logical dimensions of a NeXus DREAM detector for the day 1 configuration.""" + return NeXusDetectorDimensions(DETECTOR_BANK_SHAPES_DAY1[detector_name]) + + +providers = (dream_detector_dimensions_day1,) """Sciline providers for DREAM detector handling.""" diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index de764313..0801af4a 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -9,8 +9,9 @@ import sciline import scipp as sc from ess.powder.types import ( - DetectorName, Filename, + NeXusDetectorDimensions, + NeXusDetectorName, RawDetector, RawDetectorData, RawSample, @@ -18,6 +19,7 @@ ReducibleDetectorData, RunType, SamplePosition, + SampleRun, SourcePosition, ) from ess.reduce.nexus import extract_detector_data @@ -62,10 +64,10 @@ def load_geant4_csv(file_path: Filename[RunType]) -> AllRawDetectors[RunType]: def extract_geant4_detector( - detectors: AllRawDetectors[RunType], detector_name: DetectorName + detectors: AllRawDetectors[RunType], detector_name: NeXusDetectorName ) -> RawDetector[RunType]: """Extract a single detector from a loaded GEANT4 simulation.""" - return RawDetector[RunType](detectors["instrument"][detector_name.name]) + return RawDetector[RunType](detectors["instrument"][detector_name]) def extract_geant4_detector_data( @@ -188,6 +190,12 @@ def patch_detector_data( return ReducibleDetectorData[RunType](out) +def geant4_detector_dimensions( + data: RawDetectorData[SampleRun], +) -> NeXusDetectorDimensions[NeXusDetectorName]: + return NeXusDetectorDimensions[NeXusDetectorName](data.sizes) + + providers = ( extract_geant4_detector, extract_geant4_detector_data, @@ -195,5 +203,6 @@ def patch_detector_data( get_sample_position, get_source_position, patch_detector_data, + geant4_detector_dimensions, ) """Geant4-providers for Sciline pipelines.""" diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index a9bed1d2..0d130830 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -17,14 +17,7 @@ DspacingData, ElasticCoordTransformGraph, NormalizedByProtonCharge, - PixelMaskedData, - RawSample, - RawSource, RunType, - TwoThetaData, - TwoThetaMaskedData, - WavelengthData, - WavelengthMaskedData, ) @@ -254,6 +247,19 @@ def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: def add_scattering_coordinates( data: NormalizedByProtonCharge[RunType], graph: ElasticCoordTransformGraph ) -> DataWithScatteringCoordinates[RunType]: + """ + Add ``wavelength``, ``two_theta``, and ``dspacing`` coordinates to the data. + The input ``data`` must have a ``tof`` coordinate, as well as the necessary + positions of the beamline components (source, sample, detectors) to compute + the scattering coordinates. + + Parameters + ---------- + data: + Input data with a ``tof`` coordinate. + graph: + Coordinate transformation graph. + """ out = data.transform_coords( ["two_theta", "wavelength", "dspacing"], graph=graph, keep_intermediate=False ) diff --git a/src/ess/powder/external/powgen/beamline.py b/src/ess/powder/external/powgen/beamline.py index b8dac938..56426c5b 100644 --- a/src/ess/powder/external/powgen/beamline.py +++ b/src/ess/powder/external/powgen/beamline.py @@ -6,10 +6,18 @@ import scipp as sc -from ...types import CalibrationData, DetectorDimensions, RawCalibrationData +from ...types import ( + CalibrationData, + NeXusDetectorDimensions, + NeXusDetectorName, + RawCalibrationData, +) from .types import DetectorInfo +DETECTOR_BANK_SHAPES = {"powgen_detector": {"bank": 23, "column": 154, "row": 7}} + + def map_detector_to_spectrum( data: sc.Dataset, *, detector_info: sc.Dataset ) -> sc.Dataset: @@ -61,9 +69,13 @@ def preprocess_calibration_data( return CalibrationData(map_detector_to_spectrum(data, detector_info=detector_info)) -def powgen_detector_dimensions() -> DetectorDimensions: +def powgen_detector_dimensions( + detector_name: NeXusDetectorName, +) -> NeXusDetectorDimensions[NeXusDetectorName]: """Dimensions used by POWGEN detectors.""" - return DetectorDimensions({"bank": 23, "column": 154, "row": 7}) + return NeXusDetectorDimensions[NeXusDetectorName]( + DETECTOR_BANK_SHAPES[detector_name] + ) providers = (preprocess_calibration_data, powgen_detector_dimensions) diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 78caf38f..70ac38ce 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -8,8 +8,9 @@ from ...types import ( AccumulatedProtonCharge, CalibrationFilename, - DetectorDimensions, Filename, + NeXusDetectorDimensions, + NeXusDetectorName, ProtonCharge, RawCalibrationData, RawDataAndMetadata, @@ -109,7 +110,7 @@ def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData: def extract_raw_data( - dg: RawDataAndMetadata[RunType], sizes: DetectorDimensions + dg: RawDataAndMetadata[RunType], sizes: NeXusDetectorDimensions[NeXusDetectorName] ) -> ReducibleDetectorData[RunType]: """Return the events from a loaded data group.""" # Remove the tof binning and dimension, as it is not needed and it gets in the way @@ -131,7 +132,7 @@ def extract_proton_charge(dg: RawDataAndMetadata[RunType]) -> ProtonCharge[RunTy def extract_accumulated_proton_charge( - data: RawDetectorData[RunType], + data: ReducibleDetectorData[RunType], ) -> AccumulatedProtonCharge[RunType]: """Return the stored accumulated proton charge from a loaded data group.""" return AccumulatedProtonCharge[RunType](data.coords["gd_prtn_chrg"]) diff --git a/src/ess/powder/filtering.py b/src/ess/powder/filtering.py index eacac467..f49ab736 100644 --- a/src/ess/powder/filtering.py +++ b/src/ess/powder/filtering.py @@ -12,15 +12,7 @@ import scipp as sc -from ._util import elem_dtype, elem_unit, event_or_outer_coord -from .types import ( - FilteredData, - RawDetectorData, - ReducibleDetectorData, - RunType, - TofCroppedData, - ValidTofRange, -) +from .types import FilteredData, ReducibleDetectorData, RunType def _equivalent_bin_indices(a, b) -> bool: @@ -80,34 +72,6 @@ def remove_bad_pulses( return filtered -def crop_tof( - data: RawDetectorData[RunType], tof_range: ValidTofRange -) -> TofCroppedData[RunType]: - """Remove events outside the specified TOF range. - - Parameters - ---------- - data: - Data to be cropped. - Expected to have a coordinate called `'tof'`. - tof_range: - 1d, len-2 variable containing the lower and upper bounds for - time-of-flight. - - Returns - ------- - : - Cropped data. - """ - tof = event_or_outer_coord(data, "tof") - tof_unit = elem_unit(tof) - tof_dtype = elem_dtype(tof) - return TofCroppedData[RunType]( - data.bin(tof=tof_range.to(unit=tof_unit, dtype=tof_dtype)) - ) - - -# def filter_events(data: TofCroppedData[RunType]) -> FilteredData[RunType]: def filter_events(data: ReducibleDetectorData[RunType]) -> FilteredData[RunType]: """Remove bad events. @@ -132,8 +96,5 @@ def filter_events(data: ReducibleDetectorData[RunType]) -> FilteredData[RunType] return FilteredData[RunType](data) -providers = ( - crop_tof, - filter_events, -) +providers = (filter_events,) """Sciline providers for event filtering.""" diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index 9ee9492b..b4b452cd 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -3,15 +3,13 @@ """Grouping and merging of pixels / voxels.""" from typing import Optional -from scippneutron.conversion.graph import beamline - from .types import ( - DetectorDimensions, DspacingBins, - DspacingData, DspacingHistogram, FocussedData, MaskedData, + NeXusDetectorDimensions, + NeXusDetectorName, NormalizedByVanadium, RunType, TwoThetaBins, @@ -20,7 +18,7 @@ def focus_data( data: MaskedData[RunType], - detector_dims: DetectorDimensions, + detector_dims: NeXusDetectorDimensions[NeXusDetectorName], dspacing_bins: DspacingBins, twotheta_bins: Optional[TwoThetaBins] = None, ) -> FocussedData[RunType]: diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py index 8e6f39b3..0b485119 100644 --- a/src/ess/powder/masking.py +++ b/src/ess/powder/masking.py @@ -3,35 +3,25 @@ """ Masking functions for the powder workflow. """ -from typing import Union, Optional +from typing import Optional import numpy as np -import sciline import scipp as sc -from ess.reduce.masking import mask_range - from .types import ( DataWithScatteringCoordinates, MaskedData, MaskedDetectorIDs, - NormalizedByProtonCharge, - PixelMaskedData, PixelMaskFilename, RunType, TofMask, - TofMaskedData, - TwoThetaData, TwoThetaMask, - TwoThetaMaskedData, - WavelengthData, WavelengthMask, - WavelengthMaskedData, ) def read_pixel_masks( - filename: PixelMaskFilename, + filename: Optional[PixelMaskFilename] = None, ) -> MaskedDetectorIDs: """Read a pixel mask from a Scipp hdf5 file. @@ -40,12 +30,15 @@ def read_pixel_masks( filename: Path to the hdf5 file. """ - return MaskedDetectorIDs(sc.io.load_hdf5(filename)) + masked_ids = {} + if filename is not None: + masked_ids = {filename: sc.io.load_hdf5(filename)} + return MaskedDetectorIDs(masked_ids) def apply_masks( data: DataWithScatteringCoordinates[RunType], - masked_pixel_ids: sciline.Series[PixelMaskFilename, MaskedDetectorIDs], + masked_pixel_ids: MaskedDetectorIDs, tof_mask_func: Optional[TofMask] = None, wavelength_mask_func: Optional[WavelengthMask] = None, two_theta_mask_func: Optional[TwoThetaMask] = None, @@ -76,48 +69,4 @@ def apply_masks( return MaskedData[RunType](out) -# # def apply_wavelength_masks( -# # da: WavelengthData[RunType], mask: WavelengthMask -# # ) -> WavelengthMaskedData[RunType]: -# # if "wavelength" in da.coords and da.coords["wavelength"].ndim > 1: -# # da = da.bin(wavelength=1) -# # return WavelengthMaskedData[RunType](mask_range(da, mask=mask)) - - -# # def apply_twotheta_masks( -# # da: TwoThetaData[RunType], mask: TwoThetaMask -# # ) -> TwoThetaMaskedData[RunType]: -# # return TwoThetaMaskedData[RunType](mask_range(da, mask=mask)) - - -# def apply_tof_masking( -# da: DataWithScatteringCoordinates[RunType], mask_func: TofMask -# ) -> TofMaskedData[RunType]: -# out = da.copy(deep=False) -# out.masks["tof"] = mask_func(da.coords["tof"]) -# return TofMaskedData[RunType](out) - - -# def apply_wavelength_masking( -# da: TofMaskedData[RunType], mask_func: WavelengthMask -# ) -> WavelengthMaskedData[RunType]: -# out = da.copy(deep=False) -# out.masks["wavelength"] = mask_func(da.coords["wavelength"]) -# return WavelengthMaskedData[RunType](out) - - -# def apply_twotheta_masking( -# da: WavelengthMaskedData[RunType], mask_func: TwoThetaMask -# ) -> TwoThetaMaskedData[RunType]: -# out = da.copy(deep=False) -# out.masks["two_theta"] = mask_func(da.coords["two_theta"]) -# return TwoThetaMaskedData[RunType](out) - - -providers = ( - read_pixel_masks, - apply_masks, - # apply_tof_masking, - # apply_twotheta_masking, - # apply_wavelength_masking, -) +providers = (read_pixel_masks, apply_masks) diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index d9083199..4b9fd967 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -35,15 +35,17 @@ """Filename of the instrument calibration file.""" -# In Python 3.11, this can be replaced with a StrEnum -class DetectorName(str, Enum): - """Name of a detector.""" +# # In Python 3.11, this can be replaced with a StrEnum +# class DetectorName(str, Enum): +# """Name of a detector.""" - mantle = "mantle" - high_resolution = "high_resolution" - endcap_backward = "endcap_backward" - endcap_forward = "endcap_forward" +# mantle = "mantle" +# high_resolution = "high_resolution" +# endcap_backward = "endcap_backward" +# endcap_forward = "endcap_forward" +NeXusDetectorName = NewType("NeXusDetectorName", str) +"""Name of detector entry in NeXus file""" DspacingBins = NewType("DSpacingBins", sc.Variable) """Bin edges for d-spacing.""" @@ -97,11 +99,14 @@ class DataWithScatteringCoordinates(sciline.Scope[RunType, sc.DataArray], sc.Dat d-spacing.""" -# class DetectorDimensions(sciline.Scope[DetectorName, tuple[str, ...]], tuple[str, ...]): -# """Logical detector dimensions.""" +class NeXusDetectorDimensions( + sciline.Scope[NeXusDetectorName, Dict[str, int]], Dict[str, int] +): + """Logical detector dimensions.""" -DetectorDimensions = NewType("DetectorDimensions", Dict[str, int]) -"""Logical detector dimensions.""" + +# DetectorDimensions = NewType("DetectorDimensions", Dict[str, int]) +# """Logical detector dimensions.""" class DspacingData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): @@ -132,7 +137,7 @@ class MaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): dspacing regions.""" -MaskedDetectorIDs = NewType("MaskedDetectorIDs", sc.Variable) +MaskedDetectorIDs = NewType("MaskedDetectorIDs", Dict[str, sc.Variable]) """1-D variable listing all masked detector IDs.""" @@ -144,10 +149,6 @@ class NormalizedByProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArra """Data that has been normalized by a vanadium run.""" -class PixelMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Data with masked pixels.""" - - PixelMaskFilename = NewType("PixelMaskFilename", str) """Filename of a pixel mask.""" @@ -193,39 +194,15 @@ class SourcePosition(sciline.Scope[RunType, sc.Variable], sc.Variable): TofMask = NewType("TofMask", Callable) -"""""" - - -class TofMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Data with masked TOF values.""" - - -class TofCroppedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Raw data cropped to the valid TOF range.""" - - -class TwoThetaData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Data converted to 2theta.""" +"""TofMask is a callable that returns a mask for a given TofData.""" TwoThetaMask = NewType("TwoThetaMask", Callable) -"""""" - - -class TwoThetaMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Data with masks in 2theta.""" - - -class WavelengthData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Data converted to wavelength.""" +"""TwoThetaMask is a callable that returns a mask for a given TwoThetaData.""" WavelengthMask = NewType("WavelengthMask", Callable) -"""""" - - -class WavelengthMaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Data with masks in wavelength.""" +"""WavelengthMask is a callable that returns a mask for a given WavelengthData.""" del sc, sciline, NewType, TypeVar From f2163ffaef015c6d6167607aa002dc9b546e8c06 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Tue, 21 May 2024 14:49:54 +0200 Subject: [PATCH 07/26] formatting --- .../POWGEN_data_reduction.ipynb | 1 - src/ess/dream/instrument_view.py | 14 ++-- src/ess/dream/io/geant4.py | 6 +- src/ess/powder/__init__.py | 11 +++ src/ess/powder/conversion.py | 67 +------------------ src/ess/powder/correction.py | 1 - src/ess/powder/external/powgen/beamline.py | 1 - src/ess/powder/external/powgen/data.py | 1 - src/ess/powder/grouping.py | 1 + src/ess/powder/masking.py | 1 + src/ess/powder/types.py | 14 ---- 11 files changed, 25 insertions(+), 93 deletions(-) diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index 0742c9de..762d34b5 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -23,7 +23,6 @@ "outputs": [], "source": [ "import scipp as sc\n", - "import plopp as pp\n", "import scippneutron as scn\n", "import sciline\n", "\n", diff --git a/src/ess/dream/instrument_view.py b/src/ess/dream/instrument_view.py index c48ffddd..7a53a4d6 100644 --- a/src/ess/dream/instrument_view.py +++ b/src/ess/dream/instrument_view.py @@ -114,10 +114,10 @@ def __init__( def _add_module_control(self): import ipywidgets as ipw - # self.fig = self.scatter[0] + self.fig = self.scatter[0] self.cutting_tool = self.scatter.bottom_bar[0] self.artist_mapping = dict( - zip(self.data.keys(), self.scatter.artists.keys(), strict=True) + zip(self.data.keys(), self.fig.artists.keys(), strict=True) ) self.checkboxes = { key: ipw.Checkbox( @@ -138,9 +138,9 @@ def _add_module_control(self): for key, ch in self.checkboxes.items(): ch.key = key ch.observe(self._check_visibility, names="value") - # self.cutting_tool.cut_x.button.observe(self._check_visibility, names="value") - # self.cutting_tool.cut_y.button.observe(self._check_visibility, names="value") - # self.cutting_tool.cut_z.button.observe(self._check_visibility, names="value") + self.cutting_tool.cut_x.button.observe(self._check_visibility, names="value") + self.cutting_tool.cut_y.button.observe(self._check_visibility, names="value") + self.cutting_tool.cut_z.button.observe(self._check_visibility, names="value") self.children.insert(0, self.modules_widget) def _check_visibility(self, _): @@ -151,8 +151,8 @@ def _check_visibility(self, _): for name, ch in self.checkboxes.items(): key = self.artist_mapping[name] val = ch.value - self.scatter.artists[key].points.visible = val + self.fig.artists[key].points.visible = val for c in "xyz": cut_nodes = getattr(self.cutting_tool, f"cut_{c}").select_nodes if key in cut_nodes: - self.scatter.artists[cut_nodes[key].id].points.visible = val + self.fig.artists[cut_nodes[key].id].points.visible = val diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index 0801af4a..e9ff9206 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -1,9 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -import os -from io import BytesIO, StringIO -from typing import Dict, Optional, Union +from typing import Dict, Optional import numpy as np import sciline @@ -193,6 +191,8 @@ def patch_detector_data( def geant4_detector_dimensions( data: RawDetectorData[SampleRun], ) -> NeXusDetectorDimensions[NeXusDetectorName]: + # Fir geant4 data, we group by detector identifier, so the data already has + # logical dimensions, so we simply return the dimensions of the detector. return NeXusDetectorDimensions[NeXusDetectorName](data.sizes) diff --git a/src/ess/powder/__init__.py b/src/ess/powder/__init__.py index 37a34832..240592fc 100644 --- a/src/ess/powder/__init__.py +++ b/src/ess/powder/__init__.py @@ -31,3 +31,14 @@ *masking.providers, ) """Sciline providers for powder diffraction.""" + +__all__ = [ + "conversion", + "correction", + "filtering", + "grouping", + "masking", + "providers", + "smoothing", + "uncertainty", +] diff --git a/src/ess/powder/conversion.py b/src/ess/powder/conversion.py index 0d130830..982ecf59 100644 --- a/src/ess/powder/conversion.py +++ b/src/ess/powder/conversion.py @@ -165,69 +165,6 @@ def powder_coordinate_transformation_graph() -> ElasticCoordTransformGraph: ) -# def to_dspacing_with_positions( -# data: TwoThetaMaskedData[RunType], -# *, -# sample: Optional[RawSample[RunType]] = None, -# source: Optional[RawSource] = None, -# ) -> DspacingData[RunType]: -# """ -# Transform coordinates to d-spacing using detector positions. - -# Computes d-spacing from time-of-flight stored in `data`. - -# Attention -# --------- -# `data` may have a wavelength coordinate and dimension, -# but those are discarded. -# Only the stored time-of-flight is used, that is, any modifications to -# the wavelength coordinate after it was computed from time-of-flight are lost. - -# Raises -# ------ -# KeyError -# If `data` does not contain a 'tof' coordinate. - -# Parameters -# ---------- -# data: -# Input data in tof or wavelength dimension. -# Must have a tof coordinate. -# sample: -# Sample data with a position. -# If not given, ``data`` must contain a 'sample_position' coordinate. -# source: -# Source data with a position. -# If not given, ``data`` must contain a 'source_position' coordinate. - -# Returns -# ------- -# : -# A DataArray with the same data as the input and a 'dspacing' coordinate. -# """ -# graph = { -# **scn.conversion.graph.beamline.beamline(scatter=True), -# **scn.conversion.graph.tof.elastic_dspacing("tof"), -# } -# if sample is not None: -# graph["sample_position"] = lambda: sample["position"] -# if source is not None: -# graph["source_position"] = lambda: source["position"] - -# out = _restore_tof_if_in_wavelength(data) -# out = out.transform_coords("dspacing", graph=graph, keep_intermediate=False) -# # Add coords to ensure the result is the same whether sample or source are -# # coords in the input or separate function arguments. -# if sample is not None: -# out.coords["sample_position"] = sample["position"] -# out.coords.set_aligned("sample_position", False) -# if source is not None: -# out.coords["source_position"] = source["position"] -# out.coords.set_aligned("source_position", False) - -# return DspacingData[RunType](out) - - def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: out = data.copy(deep=False) outer = out.coords.pop("wavelength", None) @@ -244,7 +181,7 @@ def _restore_tof_if_in_wavelength(data: sc.DataArray) -> sc.DataArray: return out -def add_scattering_coordinates( +def add_scattering_coordinates_from_positions( data: NormalizedByProtonCharge[RunType], graph: ElasticCoordTransformGraph ) -> DataWithScatteringCoordinates[RunType]: """ @@ -271,5 +208,5 @@ def add_scattering_coordinates( providers_with_positions = ( powder_coordinate_transformation_graph, - add_scattering_coordinates, + add_scattering_coordinates_from_positions, ) diff --git a/src/ess/powder/correction.py b/src/ess/powder/correction.py index 6ccbbe89..687dbbb6 100644 --- a/src/ess/powder/correction.py +++ b/src/ess/powder/correction.py @@ -12,7 +12,6 @@ from .smoothing import lowpass from .types import ( AccumulatedProtonCharge, - DspacingBins, FilteredData, FocussedData, NormalizedByProtonCharge, diff --git a/src/ess/powder/external/powgen/beamline.py b/src/ess/powder/external/powgen/beamline.py index 56426c5b..c20ae358 100644 --- a/src/ess/powder/external/powgen/beamline.py +++ b/src/ess/powder/external/powgen/beamline.py @@ -14,7 +14,6 @@ ) from .types import DetectorInfo - DETECTOR_BANK_SHAPES = {"powgen_detector": {"bank": 23, "column": 154, "row": 7}} diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index 70ac38ce..b4f7ae46 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -14,7 +14,6 @@ ProtonCharge, RawCalibrationData, RawDataAndMetadata, - RawDetectorData, ReducibleDetectorData, RunType, SampleRun, diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index b4b452cd..40c0425a 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -1,6 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Grouping and merging of pixels / voxels.""" + from typing import Optional from .types import ( diff --git a/src/ess/powder/masking.py b/src/ess/powder/masking.py index 0b485119..cd51ded8 100644 --- a/src/ess/powder/masking.py +++ b/src/ess/powder/masking.py @@ -3,6 +3,7 @@ """ Masking functions for the powder workflow. """ + from typing import Optional import numpy as np diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 4b9fd967..45a4e920 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -8,7 +8,6 @@ """ from enum import Enum -from pathlib import Path from typing import Any, Callable, Dict, NewType, TypeVar import sciline @@ -35,15 +34,6 @@ """Filename of the instrument calibration file.""" -# # In Python 3.11, this can be replaced with a StrEnum -# class DetectorName(str, Enum): -# """Name of a detector.""" - -# mantle = "mantle" -# high_resolution = "high_resolution" -# endcap_backward = "endcap_backward" -# endcap_forward = "endcap_forward" - NeXusDetectorName = NewType("NeXusDetectorName", str) """Name of detector entry in NeXus file""" @@ -105,10 +95,6 @@ class NeXusDetectorDimensions( """Logical detector dimensions.""" -# DetectorDimensions = NewType("DetectorDimensions", Dict[str, int]) -# """Logical detector dimensions.""" - - class DspacingData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data converted to d-spacing.""" From edbfa0918c7fb6e29f26fc54f454dc4cbf2f63eb Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Tue, 21 May 2024 22:54:42 +0200 Subject: [PATCH 08/26] fix tests --- tests/dream/io/geant4_test.py | 169 +++++++++--------- tests/powder/conversion_test.py | 25 ++- .../external/powgen/powgen_reduction_test.py | 25 +-- 3 files changed, 110 insertions(+), 109 deletions(-) diff --git a/tests/dream/io/geant4_test.py b/tests/dream/io/geant4_test.py index 5c082f80..15178da6 100644 --- a/tests/dream/io/geant4_test.py +++ b/tests/dream/io/geant4_test.py @@ -11,30 +11,30 @@ import scipp as sc import scipp.testing from ess.dream import data, load_geant4_csv -from ess.powder.types import DetectorName, FilePath, RawDetectorData, SampleRun +from ess.powder.types import NeXusDetectorName, Filename, RawDetectorData, SampleRun -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def file_path(): - return data.get_path('data_dream0_new_hkl_Si_pwd.csv.zip') + return data.get_path("data_dream0_new_hkl_Si_pwd.csv.zip") -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def file_path_without_sans(): - return data.get_path('data_dream_with_sectors.csv.zip') + return data.get_path("data_dream_with_sectors.csv.zip") # Load file into memory only once -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def load_file(file_path): - with zipfile.ZipFile(file_path, 'r') as archive: + with zipfile.ZipFile(file_path, "r") as archive: return archive.read(archive.namelist()[0]) # Load file into memory only once -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def load_file_without_sans(file_path_without_sans): - with zipfile.ZipFile(file_path_without_sans, 'r') as archive: + with zipfile.ZipFile(file_path_without_sans, "r") as archive: return archive.read(archive.namelist()[0]) @@ -53,7 +53,7 @@ def assert_index_coord( ) -> None: assert coord.ndim == 1 assert coord.unit is None - assert coord.dtype == 'int64' + assert coord.dtype == "int64" if values is not None: assert set(np.unique(coord.values)) == values @@ -61,120 +61,120 @@ def assert_index_coord( def test_load_geant4_csv_loads_expected_structure(file): loaded = load_geant4_csv(file) assert isinstance(loaded, sc.DataGroup) - assert loaded.keys() == {'instrument'} + assert loaded.keys() == {"instrument"} - instrument = loaded['instrument'] + instrument = loaded["instrument"] assert isinstance(instrument, sc.DataGroup) assert instrument.keys() == { - 'mantle', - 'high_resolution', - 'sans', - 'endcap_forward', - 'endcap_backward', + "mantle", + "high_resolution", + "sans", + "endcap_forward", + "endcap_backward", } def test_load_geant4_csv_loads_expected_structure_without_sans(file_without_sans): loaded = load_geant4_csv(file_without_sans) assert isinstance(loaded, sc.DataGroup) - assert loaded.keys() == {'instrument'} + assert loaded.keys() == {"instrument"} - instrument = loaded['instrument'] + instrument = loaded["instrument"] assert isinstance(instrument, sc.DataGroup) assert instrument.keys() == { - 'mantle', - 'high_resolution', - 'endcap_forward', - 'endcap_backward', + "mantle", + "high_resolution", + "endcap_forward", + "endcap_backward", } @pytest.mark.parametrize( - 'key', ['mantle', 'high_resolution', 'endcap_forward', 'endcap_backward'] + "key", ["mantle", "high_resolution", "endcap_forward", "endcap_backward"] ) def test_load_gean4_csv_set_weights_to_one(file, key): - detector = load_geant4_csv(file)['instrument'][key]['events'] - events = detector.bins.constituents['data'].data + detector = load_geant4_csv(file)["instrument"][key]["events"] + events = detector.bins.constituents["data"].data sc.testing.assert_identical( - events, sc.ones(sizes=events.sizes, with_variances=True, unit='counts') + events, sc.ones(sizes=events.sizes, with_variances=True, unit="counts") ) def test_load_geant4_csv_mantle_has_expected_coords(file): # Only testing ranges that will not change in the future - mantle = load_geant4_csv(file)['instrument']['mantle']['events'] - assert_index_coord(mantle.coords['module']) - assert_index_coord(mantle.coords['segment']) - assert_index_coord(mantle.coords['counter']) - assert_index_coord(mantle.coords['wire'], values=set(range(1, 33))) - assert_index_coord(mantle.coords['strip'], values=set(range(1, 257))) - assert 'sector' not in mantle.coords + mantle = load_geant4_csv(file)["instrument"]["mantle"]["events"] + assert_index_coord(mantle.coords["module"]) + assert_index_coord(mantle.coords["segment"]) + assert_index_coord(mantle.coords["counter"]) + assert_index_coord(mantle.coords["wire"], values=set(range(1, 33))) + assert_index_coord(mantle.coords["strip"], values=set(range(1, 257))) + assert "sector" not in mantle.coords - assert 'sector' not in mantle.bins.coords - assert 'tof' in mantle.bins.coords - assert 'wavelength' in mantle.bins.coords - assert 'position' in mantle.bins.coords + assert "sector" not in mantle.bins.coords + assert "tof" in mantle.bins.coords + assert "wavelength" in mantle.bins.coords + assert "position" in mantle.bins.coords def test_load_geant4_csv_endcap_backward_has_expected_coords(file): - endcap = load_geant4_csv(file)['instrument']['endcap_backward']['events'] - assert_index_coord(endcap.coords['module']) - assert_index_coord(endcap.coords['segment']) - assert_index_coord(endcap.coords['counter']) - assert_index_coord(endcap.coords['wire'], values=set(range(1, 17))) - assert_index_coord(endcap.coords['strip'], values=set(range(1, 17))) - assert 'sector' not in endcap.coords + endcap = load_geant4_csv(file)["instrument"]["endcap_backward"]["events"] + assert_index_coord(endcap.coords["module"]) + assert_index_coord(endcap.coords["segment"]) + assert_index_coord(endcap.coords["counter"]) + assert_index_coord(endcap.coords["wire"], values=set(range(1, 17))) + assert_index_coord(endcap.coords["strip"], values=set(range(1, 17))) + assert "sector" not in endcap.coords - assert 'sector' not in endcap.bins.coords - assert 'tof' in endcap.bins.coords - assert 'wavelength' in endcap.bins.coords - assert 'position' in endcap.bins.coords + assert "sector" not in endcap.bins.coords + assert "tof" in endcap.bins.coords + assert "wavelength" in endcap.bins.coords + assert "position" in endcap.bins.coords def test_load_geant4_csv_endcap_forward_has_expected_coords(file): - endcap = load_geant4_csv(file)['instrument']['endcap_forward']['events'] - assert_index_coord(endcap.coords['module']) - assert_index_coord(endcap.coords['segment']) - assert_index_coord(endcap.coords['counter']) - assert_index_coord(endcap.coords['wire'], values=set(range(1, 17))) - assert_index_coord(endcap.coords['strip'], values=set(range(1, 17))) - assert 'sector' not in endcap.coords + endcap = load_geant4_csv(file)["instrument"]["endcap_forward"]["events"] + assert_index_coord(endcap.coords["module"]) + assert_index_coord(endcap.coords["segment"]) + assert_index_coord(endcap.coords["counter"]) + assert_index_coord(endcap.coords["wire"], values=set(range(1, 17))) + assert_index_coord(endcap.coords["strip"], values=set(range(1, 17))) + assert "sector" not in endcap.coords - assert 'sector' not in endcap.bins.coords - assert 'tof' in endcap.bins.coords - assert 'wavelength' in endcap.bins.coords - assert 'position' in endcap.bins.coords + assert "sector" not in endcap.bins.coords + assert "tof" in endcap.bins.coords + assert "wavelength" in endcap.bins.coords + assert "position" in endcap.bins.coords def test_load_geant4_csv_high_resolution_has_expected_coords(file): - hr = load_geant4_csv(file)['instrument']['high_resolution']['events'] - assert_index_coord(hr.coords['module']) - assert_index_coord(hr.coords['segment']) - assert_index_coord(hr.coords['counter']) - assert_index_coord(hr.coords['wire'], values=set(range(1, 17))) - assert_index_coord(hr.coords['strip'], values=set(range(1, 33))) - assert_index_coord(hr.coords['sector'], values=set(range(1, 5))) + hr = load_geant4_csv(file)["instrument"]["high_resolution"]["events"] + assert_index_coord(hr.coords["module"]) + assert_index_coord(hr.coords["segment"]) + assert_index_coord(hr.coords["counter"]) + assert_index_coord(hr.coords["wire"], values=set(range(1, 17))) + assert_index_coord(hr.coords["strip"], values=set(range(1, 33))) + assert_index_coord(hr.coords["sector"], values=set(range(1, 5))) - assert 'tof' in hr.bins.coords - assert 'wavelength' in hr.bins.coords - assert 'position' in hr.bins.coords + assert "tof" in hr.bins.coords + assert "wavelength" in hr.bins.coords + assert "position" in hr.bins.coords def test_load_geant4_csv_sans_has_expected_coords(file): - sans = load_geant4_csv(file)['instrument']['sans']['events'] - assert_index_coord(sans.coords['module']) - assert_index_coord(sans.coords['segment']) - assert_index_coord(sans.coords['counter']) + sans = load_geant4_csv(file)["instrument"]["sans"]["events"] + assert_index_coord(sans.coords["module"]) + assert_index_coord(sans.coords["segment"]) + assert_index_coord(sans.coords["counter"]) # check ranges only if csv file contains events from SANS detectors - if len(sans.coords['module'].values) > 0: - assert_index_coord(sans.coords['wire'], values=set(range(1, 17))) - assert_index_coord(sans.coords['strip'], values=set(range(1, 33))) - assert_index_coord(sans.coords['sector'], values=set(range(1, 5))) + if len(sans.coords["module"].values) > 0: + assert_index_coord(sans.coords["wire"], values=set(range(1, 17))) + assert_index_coord(sans.coords["strip"], values=set(range(1, 33))) + assert_index_coord(sans.coords["sector"], values=set(range(1, 5))) - assert 'tof' in sans.bins.coords - assert 'wavelength' in sans.bins.coords - assert 'position' in sans.bins.coords + assert "tof" in sans.bins.coords + assert "wavelength" in sans.bins.coords + assert "position" in sans.bins.coords def test_geant4_in_pipeline(file_path, file): @@ -182,8 +182,11 @@ def test_geant4_in_pipeline(file_path, file): pipeline = sciline.Pipeline( providers, - params={FilePath[SampleRun]: file_path, DetectorName: DetectorName('mantle')}, + params={ + Filename[SampleRun]: file_path, + NeXusDetectorName: NeXusDetectorName("mantle"), + }, ) detector = pipeline.compute(RawDetectorData[SampleRun]) - expected = load_geant4_csv(file)['instrument']['mantle']['events'] + expected = load_geant4_csv(file)["instrument"]["mantle"]["events"] sc.testing.assert_identical(detector, expected) diff --git a/tests/powder/conversion_test.py b/tests/powder/conversion_test.py index 08af8a77..281f37f3 100644 --- a/tests/powder/conversion_test.py +++ b/tests/powder/conversion_test.py @@ -6,8 +6,8 @@ import scipp.testing import scippneutron as scn from ess.powder.conversion import ( + add_scattering_coordinates_from_positions, to_dspacing_with_calibration, - to_dspacing_with_positions, ) @@ -138,7 +138,7 @@ def test_dspacing_with_calibration_does_not_use_positions(calibration): ) -def test_dspacing_with_positions(): +def test_add_scattering_coordinates_from_positions(): position = sc.vectors( dims=['spectrum'], values=np.arange(14 * 3).reshape((14, 3)), unit='m' ) @@ -149,20 +149,17 @@ def test_dspacing_with_positions(): coords={ 'position': position, 'tof': sc.linspace('tof', 1.0, 1000.0, 27, unit='us'), + 'sample_position': sample_position, + 'source_position': source_position, }, ) - dspacing = to_dspacing_with_positions( - tof, - sample=sc.DataGroup(position=sample_position), - source=sc.DataGroup(position=source_position), - ) - graph = { **scn.conversion.graph.beamline.beamline(scatter=True), - **scn.conversion.graph.tof.elastic_dspacing('tof'), + **scn.conversion.graph.tof.elastic('tof'), } - tof = tof.assign_coords( - {'sample_position': sample_position, 'source_position': source_position} - ) - expected = tof.transform_coords('dspacing', graph=graph, keep_intermediate=False) - sc.testing.assert_identical(dspacing, expected) + + result = add_scattering_coordinates_from_positions(tof, graph) + + assert 'wavelength' in result.coords + assert 'two_theta' in result.coords + assert 'dspacing' in result.coords diff --git a/tests/powder/external/powgen/powgen_reduction_test.py b/tests/powder/external/powgen/powgen_reduction_test.py index 64bba88f..7885f030 100644 --- a/tests/powder/external/powgen/powgen_reduction_test.py +++ b/tests/powder/external/powgen/powgen_reduction_test.py @@ -9,11 +9,12 @@ DspacingBins, DspacingHistogram, Filename, + NeXusDetectorName, NormalizedByProtonCharge, SampleRun, + TofMask, TwoThetaBins, UncertaintyBroadcastMode, - ValidTofRange, VanadiumRun, ) @@ -28,13 +29,17 @@ def providers(): @pytest.fixture() def params(): + from ess.powder.external import powgen + return { - Filename[SampleRun]: 'PG3_4844_event.zip', - Filename[VanadiumRun]: 'PG3_4866_event.zip', - CalibrationFilename: 'PG3_FERNS_d4832_2011_08_24.zip', + NeXusDetectorName: "powgen_detector", + Filename[SampleRun]: powgen.data.powgen_tutorial_sample_file(), + Filename[VanadiumRun]: powgen.data.powgen_tutorial_vanadium_file(), + CalibrationFilename: powgen.data.powgen_tutorial_calibration_file(), UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, - ValidTofRange: sc.array(dims=['tof'], values=[0.0, 16666.67], unit='us'), DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 200, unit='angstrom'), + TofMask: lambda x: (x < sc.scalar(0.0, unit="us")) + | (x > sc.scalar(16666.67, unit="us")), } @@ -64,17 +69,13 @@ def test_workflow_is_deterministic(providers, params): def test_pipeline_can_compute_intermediate_results(providers, params): pipeline = sciline.Pipeline(providers, params=params) result = pipeline.compute(NormalizedByProtonCharge[SampleRun]) - assert set(result.dims) == {'spectrum', 'tof'} + assert set(result.dims) == {'bank', 'column', 'row'} def test_pipeline_group_by_two_theta(providers, params): - from ess.powder.grouping import group_by_two_theta, merge_all_pixels - - providers.remove(merge_all_pixels) - providers.append(group_by_two_theta) params[TwoThetaBins] = sc.linspace( dim='two_theta', unit='deg', start=25.0, stop=90.0, num=16 - ) + ).to(unit='rad') pipeline = sciline.Pipeline(providers, params=params) result = pipeline.compute(DspacingHistogram) assert result.sizes == { @@ -82,4 +83,4 @@ def test_pipeline_group_by_two_theta(providers, params): 'dspacing': len(params[DspacingBins]) - 1, } assert sc.identical(result.coords['dspacing'], params[DspacingBins]) - assert sc.allclose(result.coords['two_theta'].to(unit='deg'), params[TwoThetaBins]) + assert sc.allclose(result.coords['two_theta'], params[TwoThetaBins]) From 6b6ac2633974f5b26e0fb41b4af379e0893276d7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Tue, 21 May 2024 20:55:18 +0000 Subject: [PATCH 09/26] Apply automatic formatting --- tests/dream/io/geant4_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/dream/io/geant4_test.py b/tests/dream/io/geant4_test.py index 15178da6..5b78cd25 100644 --- a/tests/dream/io/geant4_test.py +++ b/tests/dream/io/geant4_test.py @@ -11,7 +11,7 @@ import scipp as sc import scipp.testing from ess.dream import data, load_geant4_csv -from ess.powder.types import NeXusDetectorName, Filename, RawDetectorData, SampleRun +from ess.powder.types import Filename, NeXusDetectorName, RawDetectorData, SampleRun @pytest.fixture(scope="module") From ad54240252f1dc7779dc00b6ea1e8fa4840d2d7e Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 22 May 2024 15:48:59 +0200 Subject: [PATCH 10/26] add some tests --- src/ess/dream/io/nexus.py | 339 ++++++++++++------ src/ess/powder/grouping.py | 9 +- src/ess/powder/types.py | 5 + tests/dream/geant4_reduction_test.py | 137 +++++++ .../external/powgen/powgen_reduction_test.py | 40 +++ 5 files changed, 410 insertions(+), 120 deletions(-) create mode 100644 tests/dream/geant4_reduction_test.py diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index dd485e61..d1e38ad0 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -7,129 +7,244 @@ from typing import Any, Dict, Union import scipp as sc -import scippnexus as snx - - -def load_nexus( - path: Union[str, os.PathLike], - *, - load_pixel_shape: bool = False, - entry: str = 'entry', - fold_detectors: bool = True, -) -> sc.DataGroup: - """ - Load an unprocessed DREAM NeXus file. - - See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 - and the ICD DREAM interface specification for details. - - Notes (2023-12-06): - - - Mounting-unit, cassette, and counter roughly correspond to the azimuthal angle - in the mantle detector. However, counter is reversed in the current files. This - may also be the case in the other detectors. - - The endcap detectors have a irregular structure that cannot be fully folded. - This is not a problem but note again that the counter may be reversed. It is - currently not clear if this is a bug in the file. - - The high-resolution detector has a very odd numbering scheme. The SANS detector - is using the same, but is not populated in the current files. The scheme - attempts to follows some sort of physical ordering in space (x,y,z), but it - looks partially messed up. - - Parameters - ---------- - path: - Path to the NeXus file. - load_pixel_shape: - If True, load the pixel shape from the file's NXoff_geometry group. This is - often unused by would slow down file loading. Default is False. - entry: - Name of the entry to load. Default is "entry". - fold_detectors: - If True, fold the detector data to (partially) mimic the logical detector - structure. Default is True. - - Returns - ------- - : - A data group with the loaded file contents. - """ - definitions = snx.base_definitions() - if not load_pixel_shape: - definitions["NXdetector"] = FilteredDetector - - with snx.File(path, definitions=definitions) as f: - dg = f[entry][()] - dg = snx.compute_positions(dg) - return fold_nexus_detectors(dg) if fold_detectors else dg - - -def fold_nexus_detectors(dg: sc.DataGroup) -> sc.DataGroup: - """ - Fold the detector data in a DREAM NeXus file. - - The detector banks in the returned data group will have a multi-dimensional shape, - following the logical structure as far as possible. Note that the full structure - cannot be folded, as some dimensions are irregular. - """ - dg = dg.copy() - dg['instrument'] = dg['instrument'].copy() - instrument = dg['instrument'] - mantle = instrument['mantle_detector'] - mantle['mantle_event_data'] = mantle['mantle_event_data'].fold( +from ess.powder.types import ( + Filename, + LoadedNeXusDetector, + NeXusDetectorName, + RawSample, + RawSource, + RunType, + SamplePosition, + SourcePosition, +) +from ess.reduce import nexus + +DETECTOR_BANK_RESHAPING = { + 'endcap_backward_detector': lambda x: x.fold( dim='detector_number', sizes={ - 'wire': 32, - 'mounting_unit': 5, - 'cassette': 6, - 'counter': 2, - 'strip': 256, + "strip": 16, + "wire": 16, + "module": 11, + "segment": 28, + "counter": 2, }, - ) - for direction in ('backward', 'forward'): - endcap = instrument[f'endcap_{direction}_detector'] - endcap[f'endcap_{direction}_event_data'] = endcap[ - f'endcap_{direction}_event_data' - ].fold( - dim='detector_number', - sizes={ - 'strip': 16, - 'wire': 16, - 'sector': 5 if direction == 'forward' else 11, - 'sumo_cass_ctr': -1, # sumo*cassette*counter, irregular, cannot fold - }, - ) - high_resolution = instrument['high_resolution_detector'] - high_resolution['high_resolution_event_data'] = high_resolution[ - 'high_resolution_event_data' - ].fold( + ), + 'endcap_forward_detector': lambda x: x.fold( dim='detector_number', sizes={ - 'strip': 32, - 'other': -1, # some random order that is hard to follow + "strip": 16, + "wire": 16, + "module": 5, + "segment": 28, + "counter": 2, }, - ) - sans = instrument['sans_detector'] - sans['sans_event_data'] = sans['sans_event_data'].fold( + ), + 'mantle_detector': lambda x: x.fold( dim='detector_number', sizes={ - 'strip': 32, - 'other': -1, # some random order that is hard to follow + "wire": 32, + "module": 5, + "segment": 6, + "strip": 256, + "counter": 2, }, + ), + 'high_resolution_detector': lambda x: x.fold( + dim='detector_number', + sizes={ + "strip": 32, + "other": -1, + }, + ), + 'sans_detector': lambda x: x.fold( + dim='detector_number', + sizes={ + "strip": 32, + "other": -1, + }, + ), +} + + +def load_nexus_sample(file_path: Filename[RunType]) -> RawSample[RunType]: + return RawSample[RunType](nexus.load_sample(file_path)) + + +def dummy_load_sample(file_path: Filename[RunType]) -> RawSample[RunType]: + return RawSample[RunType]( + sc.DataGroup({'position': sc.vector(value=[0, 0, 0], unit='m')}) ) - return dg -def _skip(name: str, obj: Union[snx.Field, snx.Group]) -> bool: - skip_classes = (snx.NXoff_geometry,) - return isinstance(obj, snx.Group) and (obj.nx_class in skip_classes) +def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: + return RawSource[RunType](nexus.load_source(file_path)) + + +def load_nexus_detector( + file_path: Filename[RunType], detector_name: NeXusDetectorName +) -> LoadedNeXusDetector[RunType]: + return LoadedNeXusDetector[RunType]( + nexus.load_detector(file_path=file_path, detector_name=detector_name) + ) + + +def get_source_position( + raw_source: RawSource[RunType], +) -> SourcePosition[RunType]: + return SourcePosition[RunType](raw_source['position']) + + +def get_sample_position( + raw_sample: RawSample[RunType], +) -> SamplePosition[RunType]: + return SamplePosition[RunType](raw_sample['position']) + + +def get_detector_data( + detector: LoadedNeXusDetector[RunType], + detector_name: NeXusDetectorName, +) -> RawData[RunType]: + da = nexus.extract_detector_data(detector) + if detector_name in DETECTOR_BANK_RESHAPING: + da = DETECTOR_BANK_RESHAPING[detector_name](da) + return RawData[RunType](da) + + +def patch_detector_data( + detector_data: RawData[ScatteringRunType], + source_position: SourcePosition[ScatteringRunType], + sample_position: SamplePosition[ScatteringRunType], +) -> ConfiguredReducibleDataData[ScatteringRunType]: + return ConfiguredReducibleDataData[ScatteringRunType]( + _add_variances_and_coordinates( + da=detector_data, + source_position=source_position, + sample_position=sample_position, + ) + ) -class FilteredDetector(snx.NXdetector): - def __init__( - self, attrs: Dict[str, Any], children: Dict[str, Union[snx.Field, snx.Group]] - ): - children = { - name: child for name, child in children.items() if not _skip(name, child) - } - super().__init__(attrs=attrs, children=children) +# def load_nexus( +# path: Union[str, os.PathLike], +# *, +# load_pixel_shape: bool = False, +# entry: str = 'entry', +# fold_detectors: bool = True, +# ) -> sc.DataGroup: +# """ +# Load an unprocessed DREAM NeXus file. + +# See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 +# and the ICD DREAM interface specification for details. + +# Notes (2023-12-06): + +# - Mounting-unit, cassette, and counter roughly correspond to the azimuthal angle +# in the mantle detector. However, counter is reversed in the current files. This +# may also be the case in the other detectors. +# - The endcap detectors have a irregular structure that cannot be fully folded. +# This is not a problem but note again that the counter may be reversed. It is +# currently not clear if this is a bug in the file. +# - The high-resolution detector has a very odd numbering scheme. The SANS detector +# is using the same, but is not populated in the current files. The scheme +# attempts to follows some sort of physical ordering in space (x,y,z), but it +# looks partially messed up. + +# Parameters +# ---------- +# path: +# Path to the NeXus file. +# load_pixel_shape: +# If True, load the pixel shape from the file's NXoff_geometry group. This is +# often unused by would slow down file loading. Default is False. +# entry: +# Name of the entry to load. Default is "entry". +# fold_detectors: +# If True, fold the detector data to (partially) mimic the logical detector +# structure. Default is True. + +# Returns +# ------- +# : +# A data group with the loaded file contents. +# """ +# definitions = snx.base_definitions() +# if not load_pixel_shape: +# definitions["NXdetector"] = FilteredDetector + +# with snx.File(path, definitions=definitions) as f: +# dg = f[entry][()] +# dg = snx.compute_positions(dg) +# return fold_nexus_detectors(dg) if fold_detectors else dg + + +# def fold_nexus_detectors(dg: sc.DataGroup) -> sc.DataGroup: +# """ +# Fold the detector data in a DREAM NeXus file. + +# The detector banks in the returned data group will have a multi-dimensional shape, +# following the logical structure as far as possible. Note that the full structure +# cannot be folded, as some dimensions are irregular. +# """ +# dg = dg.copy() +# dg['instrument'] = dg['instrument'].copy() +# instrument = dg['instrument'] +# mantle = instrument['mantle_detector'] +# mantle['mantle_event_data'] = mantle['mantle_event_data'].fold( +# dim='detector_number', +# sizes={ +# 'wire': 32, +# 'mounting_unit': 5, +# 'cassette': 6, +# 'counter': 2, +# 'strip': 256, +# }, +# ) +# for direction in ('backward', 'forward'): +# endcap = instrument[f'endcap_{direction}_detector'] +# endcap[f'endcap_{direction}_event_data'] = endcap[ +# f'endcap_{direction}_event_data' +# ].fold( +# dim='detector_number', +# sizes={ +# 'strip': 16, +# 'wire': 16, +# 'sector': 5 if direction == 'forward' else 11, +# 'sumo_cass_ctr': -1, # sumo*cassette*counter, irregular, cannot fold +# }, +# ) +# high_resolution = instrument['high_resolution_detector'] +# high_resolution['high_resolution_event_data'] = high_resolution[ +# 'high_resolution_event_data' +# ].fold( +# dim='detector_number', +# sizes={ +# 'strip': 32, +# 'other': -1, # some random order that is hard to follow +# }, +# ) +# sans = instrument['sans_detector'] +# sans['sans_event_data'] = sans['sans_event_data'].fold( +# dim='detector_number', +# sizes={ +# 'strip': 32, +# 'other': -1, # some random order that is hard to follow +# }, +# ) +# return dg + + +# def _skip(name: str, obj: Union[snx.Field, snx.Group]) -> bool: +# skip_classes = (snx.NXoff_geometry,) +# return isinstance(obj, snx.Group) and (obj.nx_class in skip_classes) + + +# class FilteredDetector(snx.NXdetector): +# def __init__( +# self, attrs: Dict[str, Any], children: Dict[str, Union[snx.Field, snx.Group]] +# ): +# children = { +# name: child for name, child in children.items() if not _skip(name, child) +# } +# super().__init__(attrs=attrs, children=children) diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index 40c0425a..974e14e8 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -9,8 +9,6 @@ DspacingHistogram, FocussedData, MaskedData, - NeXusDetectorDimensions, - NeXusDetectorName, NormalizedByVanadium, RunType, TwoThetaBins, @@ -19,7 +17,6 @@ def focus_data( data: MaskedData[RunType], - detector_dims: NeXusDetectorDimensions[NeXusDetectorName], dspacing_bins: DspacingBins, twotheta_bins: Optional[TwoThetaBins] = None, ) -> FocussedData[RunType]: @@ -29,11 +26,7 @@ def focus_data( bins[dspacing_bins.dim] = dspacing_bins if (twotheta_bins is None) or ("two_theta" in data.bins.coords): - # In this case merge data from all pixels - # Put the dims into the same order as in the data. - # See https://github.com/scipp/scipp/issues/3408 - to_concat = tuple(dim for dim in data.dims if dim in detector_dims) - data = data.bins.concat(to_concat) + data = data.bins.concat() return FocussedData[RunType](data.bin(**bins)) diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 45a4e920..89c2226c 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -118,6 +118,11 @@ class FocussedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Intensity vs d-spacing after focussing pixels.""" +class LoadedNeXusDetector(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): + """Detector data, loaded from a NeXus file, containing not only neutron events + but also pixel shape information, transformations, ...""" + + class MaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data with masked pixels, tof regions, wavelength regions, 2theta regions, or dspacing regions.""" diff --git a/tests/dream/geant4_reduction_test.py b/tests/dream/geant4_reduction_test.py new file mode 100644 index 00000000..57c48e13 --- /dev/null +++ b/tests/dream/geant4_reduction_test.py @@ -0,0 +1,137 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +import pytest +import sciline +import scipp as sc +from ess.powder.types import ( + AccumulatedProtonCharge, + DspacingBins, + DspacingHistogram, + EmptyCanRun, + Filename, + MaskedData, + NeXusDetectorName, + NormalizedByProtonCharge, + RawSample, + RawSource, + SampleRun, + TofMask, + TwoThetaBins, + TwoThetaMask, + UncertaintyBroadcastMode, + VanadiumRun, + WavelengthMask, +) + + +@pytest.fixture() +def providers(): + from ess import powder + from ess.dream.io.geant4 import providers as geant4_providers + + return [*powder.providers, *geant4_providers] + + +@pytest.fixture(params=["mantle", "endcap_backward", "endcap_forward"]) +def params(request): + from ess import dream + + # Not available in simulated data + sample = sc.DataGroup(position=sc.vector([0.0, 0.0, 0.0], unit='mm')) + source = sc.DataGroup(position=sc.vector([-3.478, 0.0, -76550], unit='mm')) + charge = sc.scalar(1.0, unit='µAh') + + return { + NeXusDetectorName: request.param, + Filename[SampleRun]: dream.data.simulated_diamond_sample(), + Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(), + Filename[EmptyCanRun]: dream.data.simulated_empty_can(), + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, + DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'), + TofMask: lambda x: (x < sc.scalar(0.0, unit='ns')) + | (x > sc.scalar(86e6, unit='ns')), + RawSample[SampleRun]: sample, + RawSample[VanadiumRun]: sample, + RawSource[SampleRun]: source, + RawSource[VanadiumRun]: source, + AccumulatedProtonCharge[SampleRun]: charge, + AccumulatedProtonCharge[VanadiumRun]: charge, + } + + +def test_can_create_pipeline(providers, params): + sciline.Pipeline(providers, params=params) + + +def test_pipeline_can_compute_dspacing_histogram(providers, params): + pipeline = sciline.Pipeline(providers, params=params) + result = pipeline.compute(DspacingHistogram) + assert result.sizes == { + 'dspacing': len(params[DspacingBins]) - 1, + } + assert sc.identical(result.coords['dspacing'], params[DspacingBins]) + + +def test_workflow_is_deterministic(providers, params): + pipeline = sciline.Pipeline(providers, params=params) + # This is Sciline's default scheduler, but we want to be explicit here + scheduler = sciline.scheduler.DaskScheduler() + graph = pipeline.get(DspacingHistogram, scheduler=scheduler) + reference = graph.compute().data + result = graph.compute().data + assert sc.identical(sc.values(result), sc.values(reference)) + + +def test_pipeline_can_compute_intermediate_results(providers, params): + pipeline = sciline.Pipeline(providers, params=params) + result = pipeline.compute(NormalizedByProtonCharge[SampleRun]) + assert set(result.dims) == {'segment', 'wire', 'counter', 'strip', 'module'} + + +def test_pipeline_group_by_two_theta(providers, params): + params[TwoThetaBins] = sc.linspace( + dim='two_theta', unit='rad', start=0.8, stop=2.4, num=17 + ) + pipeline = sciline.Pipeline(providers, params=params) + result = pipeline.compute(DspacingHistogram) + assert result.sizes == { + 'two_theta': 16, + 'dspacing': len(params[DspacingBins]) - 1, + } + assert sc.identical(result.coords['dspacing'], params[DspacingBins]) + assert sc.allclose(result.coords['two_theta'], params[TwoThetaBins]) + + +def test_pipeline_wavelength_masking(providers, params): + wmin = sc.scalar(0.18, unit="angstrom") + wmax = sc.scalar(0.21, unit="angstrom") + params[WavelengthMask] = lambda x: (x > wmin) & (x < wmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'wavelength' in masked_sample.bins.masks + sum_in_masked_region = ( + masked_sample.bin(wavelength=sc.concat([wmin, wmax], dim='wavelength')) + .sum() + .data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) + + +def test_pipeline_two_theta_masking(providers, params): + tmin = sc.scalar(1.0, unit="rad") + tmax = sc.scalar(1.2, unit="rad") + params[TwoThetaMask] = lambda x: (x > tmin) & (x < tmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'two_theta' in masked_sample.bins.masks + sum_in_masked_region = ( + masked_sample.bin(two_theta=sc.concat([tmin, tmax], dim='two_theta')).sum().data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) diff --git a/tests/powder/external/powgen/powgen_reduction_test.py b/tests/powder/external/powgen/powgen_reduction_test.py index 7885f030..e9564010 100644 --- a/tests/powder/external/powgen/powgen_reduction_test.py +++ b/tests/powder/external/powgen/powgen_reduction_test.py @@ -9,13 +9,16 @@ DspacingBins, DspacingHistogram, Filename, + MaskedData, NeXusDetectorName, NormalizedByProtonCharge, SampleRun, TofMask, TwoThetaBins, + TwoThetaMask, UncertaintyBroadcastMode, VanadiumRun, + WavelengthMask, ) @@ -84,3 +87,40 @@ def test_pipeline_group_by_two_theta(providers, params): } assert sc.identical(result.coords['dspacing'], params[DspacingBins]) assert sc.allclose(result.coords['two_theta'], params[TwoThetaBins]) + + +def test_pipeline_wavelength_masking(providers, params): + wmin = sc.scalar(0.18, unit="angstrom") + wmax = sc.scalar(0.21, unit="angstrom") + params[WavelengthMask] = lambda x: (x > wmin) & (x < wmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'wavelength' in masked_sample.bins.masks + sum_in_masked_region = ( + masked_sample.bin(wavelength=sc.concat([wmin, wmax], dim='wavelength')) + .sum() + .data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) + + +def test_pipeline_two_theta_masking(providers, params): + tmin = sc.scalar(0.8, unit="rad") + tmax = sc.scalar(1.0, unit="rad") + params[TwoThetaMask] = lambda x: (x > tmin) & (x < tmax) + pipeline = sciline.Pipeline(providers, params=params) + masked_sample = pipeline.compute(MaskedData[SampleRun]) + assert 'two_theta' in masked_sample.masks + sum_in_masked_region = ( + masked_sample.flatten(to='pixel') + .hist(two_theta=sc.concat([tmin, tmax], dim='two_theta')) + .sum() + .data + ) + assert sc.allclose( + sum_in_masked_region, + sc.scalar(0.0, unit=sum_in_masked_region.unit), + ) From 6134e3e19aea952c7e1aa63067b3e4ba47342573 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 22 May 2024 16:47:42 +0200 Subject: [PATCH 11/26] copy style of nexus utilities from esssans --- src/ess/dream/beamline.py | 36 ------------------ src/ess/dream/io/nexus.py | 80 +++++++++++++++++++++------------------ 2 files changed, 44 insertions(+), 72 deletions(-) delete mode 100644 src/ess/dream/beamline.py diff --git a/src/ess/dream/beamline.py b/src/ess/dream/beamline.py deleted file mode 100644 index 44f0424d..00000000 --- a/src/ess/dream/beamline.py +++ /dev/null @@ -1,36 +0,0 @@ -# SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) - -"""Beamline tools for DREAM.""" - -from ess.powder.types import NeXusDetectorDimensions, NeXusDetectorName - -DETECTOR_BANK_SHAPES_DAY1 = { - "endcap_backward": { - "strip": 16, - "wire": 16, - "module": 11, - "segment": 28, - "counter": 2, - }, - "endcap_forward": { - "strip": 16, - "wire": 16, - "module": 5, - "segment": 28, - "counter": 2, - }, - "mantle": {"wire": 32, "module": 5, "segment": 6, "strip": 256, "counter": 2}, - # TODO: missing "high_resolution" and "sans" detectors -} - - -def dream_detector_dimensions_day1( - detector_name: NeXusDetectorName, -) -> NeXusDetectorDimensions[NeXusDetectorName]: - """Logical dimensions of a NeXus DREAM detector for the day 1 configuration.""" - return NeXusDetectorDimensions(DETECTOR_BANK_SHAPES_DAY1[detector_name]) - - -providers = (dream_detector_dimensions_day1,) -"""Sciline providers for DREAM detector handling.""" diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index d1e38ad0..8b31e499 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -1,18 +1,27 @@ # SPDX-License-Identifier: BSD-3-Clause -# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) -"""NeXus input/output for DREAM.""" +"""NeXus input/output for DREAM. + +Notes on the detector dimensions (2024-05-22): + +See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 +and the ICD DREAM interface specification for details. + + - The high-resolution and SANS detectors have a very odd numbering scheme. + The scheme attempts to follows some sort of physical ordering in space (x,y,z), + but it is not possible to reshape the data into all the logical dimensions. +""" -import os -from typing import Any, Dict, Union -import scipp as sc from ess.powder.types import ( Filename, LoadedNeXusDetector, NeXusDetectorName, + RawDetectorData, RawSample, RawSource, + ReducibleDetectorData, RunType, SamplePosition, SourcePosition, @@ -20,8 +29,8 @@ from ess.reduce import nexus DETECTOR_BANK_RESHAPING = { - 'endcap_backward_detector': lambda x: x.fold( - dim='detector_number', + "endcap_backward_detector": lambda x: x.fold( + dim="detector_number", sizes={ "strip": 16, "wire": 16, @@ -30,8 +39,8 @@ "counter": 2, }, ), - 'endcap_forward_detector': lambda x: x.fold( - dim='detector_number', + "endcap_forward_detector": lambda x: x.fold( + dim="detector_number", sizes={ "strip": 16, "wire": 16, @@ -40,8 +49,8 @@ "counter": 2, }, ), - 'mantle_detector': lambda x: x.fold( - dim='detector_number', + "mantle_detector": lambda x: x.fold( + dim="detector_number", sizes={ "wire": 32, "module": 5, @@ -50,15 +59,15 @@ "counter": 2, }, ), - 'high_resolution_detector': lambda x: x.fold( - dim='detector_number', + "high_resolution_detector": lambda x: x.fold( + dim="detector_number", sizes={ "strip": 32, "other": -1, }, ), - 'sans_detector': lambda x: x.fold( - dim='detector_number', + "sans_detector": lambda x: x.fold( + dim="detector_number", sizes={ "strip": 32, "other": -1, @@ -71,12 +80,6 @@ def load_nexus_sample(file_path: Filename[RunType]) -> RawSample[RunType]: return RawSample[RunType](nexus.load_sample(file_path)) -def dummy_load_sample(file_path: Filename[RunType]) -> RawSample[RunType]: - return RawSample[RunType]( - sc.DataGroup({'position': sc.vector(value=[0, 0, 0], unit='m')}) - ) - - def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: return RawSource[RunType](nexus.load_source(file_path)) @@ -92,37 +95,42 @@ def load_nexus_detector( def get_source_position( raw_source: RawSource[RunType], ) -> SourcePosition[RunType]: - return SourcePosition[RunType](raw_source['position']) + return SourcePosition[RunType](raw_source["position"]) def get_sample_position( raw_sample: RawSample[RunType], ) -> SamplePosition[RunType]: - return SamplePosition[RunType](raw_sample['position']) + return SamplePosition[RunType](raw_sample["position"]) def get_detector_data( detector: LoadedNeXusDetector[RunType], detector_name: NeXusDetectorName, -) -> RawData[RunType]: +) -> RawDetectorData[RunType]: da = nexus.extract_detector_data(detector) if detector_name in DETECTOR_BANK_RESHAPING: da = DETECTOR_BANK_RESHAPING[detector_name](da) - return RawData[RunType](da) + return RawDetectorData[RunType](da) def patch_detector_data( - detector_data: RawData[ScatteringRunType], - source_position: SourcePosition[ScatteringRunType], - sample_position: SamplePosition[ScatteringRunType], -) -> ConfiguredReducibleDataData[ScatteringRunType]: - return ConfiguredReducibleDataData[ScatteringRunType]( - _add_variances_and_coordinates( - da=detector_data, - source_position=source_position, - sample_position=sample_position, - ) - ) + detector_data: RawDetectorData[RunType], + source_position: SourcePosition[RunType], + sample_position: SamplePosition[RunType], +) -> ReducibleDetectorData[RunType]: + """ + Patch a detector data object with source and sample positions. + Also adds variances to the event data if they are missing. + """ + out = detector_data.copy(deep=False) + if out.bins is not None: + content = out.bins.constituents["data"] + if content.variances is None: + content.variances = content.values + out.coords["sample_position"] = sample_position + out.coords["source_position"] = source_position + return ReducibleDetectorData[RunType](out) # def load_nexus( From 90ae8bbb5a7096dfbae456a53be546ae659af7fb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Wed, 22 May 2024 14:48:21 +0000 Subject: [PATCH 12/26] Apply automatic formatting --- src/ess/dream/io/nexus.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index 8b31e499..46b2d6c7 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -13,7 +13,6 @@ but it is not possible to reshape the data into all the logical dimensions. """ - from ess.powder.types import ( Filename, LoadedNeXusDetector, From 5e5f5ecd0cea907774145b52b7aa513dacc56aca Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 22 May 2024 16:51:36 +0200 Subject: [PATCH 13/26] fix notebook --- .../sns-instruments/POWGEN_data_reduction.ipynb | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index 762d34b5..0a2ba67b 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -231,7 +231,7 @@ "source": [ "results = pipeline.compute(\n", " (\n", - " RawDetectorData[SampleRun],\n", + " ReducibleDetectorData[SampleRun],\n", " MaskedData[SampleRun],\n", " FilteredData[SampleRun],\n", " FilteredData[VanadiumRun],\n", @@ -246,7 +246,7 @@ "metadata": {}, "outputs": [], "source": [ - "results[RawDetectorData[SampleRun]]" + "results[ReducibleDetectorData[SampleRun]]" ] }, { @@ -256,7 +256,7 @@ "metadata": {}, "outputs": [], "source": [ - "scn.instrument_view(results[MaskedData[SampleRun]].hist(), pixel_size=0.05)" + "results[MaskedData[SampleRun]].bins.concat().hist(wavelength=300).plot()" ] }, { @@ -377,8 +377,7 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.12" + "pygments_lexer": "ipython3" } }, "nbformat": 4, From af1a8977fc3e332023b9a0166e0b0eaf631f4d4d Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 22 May 2024 17:03:34 +0200 Subject: [PATCH 14/26] fix imports for nexus and add nexus providers --- src/ess/dream/__init__.py | 9 +- src/ess/dream/io/__init__.py | 4 +- src/ess/dream/io/nexus.py | 136 +++---------------------- src/ess/powder/external/powgen/data.py | 2 +- 4 files changed, 19 insertions(+), 132 deletions(-) diff --git a/src/ess/dream/__init__.py b/src/ess/dream/__init__.py index bb50dc01..020bf577 100644 --- a/src/ess/dream/__init__.py +++ b/src/ess/dream/__init__.py @@ -7,9 +7,9 @@ import importlib.metadata -from . import beamline, data +from . import data from .instrument_view import instrument_view -from .io import fold_nexus_detectors, load_geant4_csv, load_nexus +from .io import load_geant4_csv, nexus try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -18,13 +18,12 @@ del importlib -providers = (*beamline.providers,) +providers = (*nexus.providers,) __all__ = [ 'data', 'beamline', - 'fold_nexus_detectors', 'instrument_view', 'load_geant4_csv', - 'load_nexus', + 'nexus', ] diff --git a/src/ess/dream/io/__init__.py b/src/ess/dream/io/__init__.py index 605738ed..2224a078 100644 --- a/src/ess/dream/io/__init__.py +++ b/src/ess/dream/io/__init__.py @@ -3,7 +3,7 @@ """Input/output for DREAM.""" +from . import nexus from .geant4 import load_geant4_csv -from .nexus import fold_nexus_detectors, load_nexus -__all__ = ["fold_nexus_detectors", "load_geant4_csv", "load_nexus"] +__all__ = ["nexus", "load_geant4_csv"] diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index 46b2d6c7..7a16ac33 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -131,127 +131,15 @@ def patch_detector_data( out.coords["source_position"] = source_position return ReducibleDetectorData[RunType](out) - -# def load_nexus( -# path: Union[str, os.PathLike], -# *, -# load_pixel_shape: bool = False, -# entry: str = 'entry', -# fold_detectors: bool = True, -# ) -> sc.DataGroup: -# """ -# Load an unprocessed DREAM NeXus file. - -# See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 -# and the ICD DREAM interface specification for details. - -# Notes (2023-12-06): - -# - Mounting-unit, cassette, and counter roughly correspond to the azimuthal angle -# in the mantle detector. However, counter is reversed in the current files. This -# may also be the case in the other detectors. -# - The endcap detectors have a irregular structure that cannot be fully folded. -# This is not a problem but note again that the counter may be reversed. It is -# currently not clear if this is a bug in the file. -# - The high-resolution detector has a very odd numbering scheme. The SANS detector -# is using the same, but is not populated in the current files. The scheme -# attempts to follows some sort of physical ordering in space (x,y,z), but it -# looks partially messed up. - -# Parameters -# ---------- -# path: -# Path to the NeXus file. -# load_pixel_shape: -# If True, load the pixel shape from the file's NXoff_geometry group. This is -# often unused by would slow down file loading. Default is False. -# entry: -# Name of the entry to load. Default is "entry". -# fold_detectors: -# If True, fold the detector data to (partially) mimic the logical detector -# structure. Default is True. - -# Returns -# ------- -# : -# A data group with the loaded file contents. -# """ -# definitions = snx.base_definitions() -# if not load_pixel_shape: -# definitions["NXdetector"] = FilteredDetector - -# with snx.File(path, definitions=definitions) as f: -# dg = f[entry][()] -# dg = snx.compute_positions(dg) -# return fold_nexus_detectors(dg) if fold_detectors else dg - - -# def fold_nexus_detectors(dg: sc.DataGroup) -> sc.DataGroup: -# """ -# Fold the detector data in a DREAM NeXus file. - -# The detector banks in the returned data group will have a multi-dimensional shape, -# following the logical structure as far as possible. Note that the full structure -# cannot be folded, as some dimensions are irregular. -# """ -# dg = dg.copy() -# dg['instrument'] = dg['instrument'].copy() -# instrument = dg['instrument'] -# mantle = instrument['mantle_detector'] -# mantle['mantle_event_data'] = mantle['mantle_event_data'].fold( -# dim='detector_number', -# sizes={ -# 'wire': 32, -# 'mounting_unit': 5, -# 'cassette': 6, -# 'counter': 2, -# 'strip': 256, -# }, -# ) -# for direction in ('backward', 'forward'): -# endcap = instrument[f'endcap_{direction}_detector'] -# endcap[f'endcap_{direction}_event_data'] = endcap[ -# f'endcap_{direction}_event_data' -# ].fold( -# dim='detector_number', -# sizes={ -# 'strip': 16, -# 'wire': 16, -# 'sector': 5 if direction == 'forward' else 11, -# 'sumo_cass_ctr': -1, # sumo*cassette*counter, irregular, cannot fold -# }, -# ) -# high_resolution = instrument['high_resolution_detector'] -# high_resolution['high_resolution_event_data'] = high_resolution[ -# 'high_resolution_event_data' -# ].fold( -# dim='detector_number', -# sizes={ -# 'strip': 32, -# 'other': -1, # some random order that is hard to follow -# }, -# ) -# sans = instrument['sans_detector'] -# sans['sans_event_data'] = sans['sans_event_data'].fold( -# dim='detector_number', -# sizes={ -# 'strip': 32, -# 'other': -1, # some random order that is hard to follow -# }, -# ) -# return dg - - -# def _skip(name: str, obj: Union[snx.Field, snx.Group]) -> bool: -# skip_classes = (snx.NXoff_geometry,) -# return isinstance(obj, snx.Group) and (obj.nx_class in skip_classes) - - -# class FilteredDetector(snx.NXdetector): -# def __init__( -# self, attrs: Dict[str, Any], children: Dict[str, Union[snx.Field, snx.Group]] -# ): -# children = { -# name: child for name, child in children.items() if not _skip(name, child) -# } -# super().__init__(attrs=attrs, children=children) +providers = ( + load_nexus_sample, + load_nexus_source, + load_nexus_detector, + get_source_position, + get_sample_position, + get_detector_data, + patch_detector_data, +) +""" +Providers for loading and processing DREAM NeXus data. +""" diff --git a/src/ess/powder/external/powgen/data.py b/src/ess/powder/external/powgen/data.py index b4f7ae46..8be5b5c0 100644 --- a/src/ess/powder/external/powgen/data.py +++ b/src/ess/powder/external/powgen/data.py @@ -115,7 +115,7 @@ def extract_raw_data( # Remove the tof binning and dimension, as it is not needed and it gets in the way # of masking. out = dg["data"].squeeze() - del out.coords["tof"] + out.coords.pop("tof", None) out = out.fold(dim="spectrum", sizes=sizes) return ReducibleDetectorData[RunType](out) From 60dbca6ed479ba6ca94a8f046197f8f2d2933665 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci-lite[bot]" <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com> Date: Wed, 22 May 2024 15:04:15 +0000 Subject: [PATCH 15/26] Apply automatic formatting --- src/ess/dream/io/nexus.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index 7a16ac33..07d95206 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -131,6 +131,7 @@ def patch_detector_data( out.coords["source_position"] = source_position return ReducibleDetectorData[RunType](out) + providers = ( load_nexus_sample, load_nexus_source, From 80745c28aba7abb54a763002b55fb25ba55d480e Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 22 May 2024 17:55:25 +0200 Subject: [PATCH 16/26] update nexus tests --- src/ess/dream/io/nexus.py | 17 ++++- tests/dream/io/nexus_test.py | 132 +++++++++++++++++++---------------- 2 files changed, 85 insertions(+), 64 deletions(-) diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index 7a16ac33..beb9096c 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -13,6 +13,7 @@ but it is not possible to reshape the data into all the logical dimensions. """ +import scipp as sc from ess.powder.types import ( Filename, LoadedNeXusDetector, @@ -79,6 +80,15 @@ def load_nexus_sample(file_path: Filename[RunType]) -> RawSample[RunType]: return RawSample[RunType](nexus.load_sample(file_path)) +def dummy_load_sample(file_path: Filename[RunType]) -> RawSample[RunType]: + """ + In test files there is not always a sample, so we need a dummy. + """ + return RawSample[RunType]( + sc.DataGroup({'position': sc.vector(value=[0, 0, 0], unit='m')}) + ) + + def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: return RawSource[RunType](nexus.load_source(file_path)) @@ -86,9 +96,9 @@ def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: def load_nexus_detector( file_path: Filename[RunType], detector_name: NeXusDetectorName ) -> LoadedNeXusDetector[RunType]: - return LoadedNeXusDetector[RunType]( - nexus.load_detector(file_path=file_path, detector_name=detector_name) - ) + out = nexus.load_detector(file_path=file_path, detector_name=detector_name) + out.pop("pixel_shape", None) + return LoadedNeXusDetector[RunType](out) def get_source_position( @@ -131,6 +141,7 @@ def patch_detector_data( out.coords["source_position"] = source_position return ReducibleDetectorData[RunType](out) + providers = ( load_nexus_sample, load_nexus_source, diff --git a/tests/dream/io/nexus_test.py b/tests/dream/io/nexus_test.py index 6f7fe5f9..126197ac 100644 --- a/tests/dream/io/nexus_test.py +++ b/tests/dream/io/nexus_test.py @@ -1,76 +1,86 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) import pytest +import sciline from ess import dream +from ess.dream import nexus +from ess.powder.types import ( + Filename, + NeXusDetectorName, + RawDetectorData, + ReducibleDetectorData, + SampleRun, +) - -@pytest.fixture() -def filename(): - return dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs') +bank_dims = {'wire', 'module', 'segment', 'strip', 'counter'} +hr_sans_dims = {'strip', 'other'} -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_loads_file(filename): - dg = dream.load_nexus(filename) - assert 'instrument' in dg - instr = dg['instrument'] - for name in ( - 'mantle', - 'endcap_backward', - 'endcap_forward', - 'high_resolution', - 'sans', - ): - assert f'{name}_detector' in instr - det = instr[f'{name}_detector'] - assert 'pixel_shape' not in det +@pytest.fixture() +def providers(): + return ( + nexus.dummy_load_sample, + nexus.load_nexus_source, + nexus.load_nexus_detector, + nexus.get_sample_position, + nexus.get_source_position, + nexus.get_detector_data, + nexus.patch_detector_data, + ) -def test_load_nexus_fails_if_entry_not_found(filename): - with pytest.raises(KeyError): - dream.load_nexus(filename, entry='foo') +@pytest.fixture( + params=[ + 'mantle_detector', + 'endcap_backward_detector', + 'endcap_forward_detector', + 'high_resolution_detector', + # TODO: the 'sans_detector' is strange in the current files + ] +) +def params(request): + params = { + Filename[SampleRun]: dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs'), + NeXusDetectorName: request.param, + } + return params -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_folds_detectors_by_default(filename): - dg = dream.load_nexus(filename) - instr = dg['instrument'] - # sans_detector is not populated in the current files - for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): - det = instr[f'{name}_detector'] - # There may be other dims, but some are irregular and this may be subject to - # change - assert 'strip' in det.dims +def test_can_load_nexus_detector_data(providers, params): + pipeline = sciline.Pipeline(params=params, providers=providers) + result = pipeline.compute(RawDetectorData[SampleRun]) + assert ( + set(result.dims) == hr_sans_dims + if params[NeXusDetectorName] + in ( + 'high_resolution_detector', + 'sans_detector', + ) + else bank_dims + ) -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_with_disabled_fold(filename): - dg = dream.load_nexus(filename, fold_detectors=False) - instr = dg['instrument'] - for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): - det = instr[f'{name}_detector'] - assert det.dims == ('detector_number',) +def test_load_fails_with_bad_detector_name(providers): + params = { + Filename[SampleRun]: dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs'), + NeXusDetectorName: 'bad_detector', + } + pipeline = sciline.Pipeline(params=params, providers=providers) + with pytest.raises(KeyError, match='bad_detector'): + pipeline.compute(RawDetectorData[SampleRun]) -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") -@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") -def test_load_nexus_with_pixel_shape(filename): - dg = dream.load_nexus(filename, load_pixel_shape=True) - assert 'instrument' in dg - instr = dg['instrument'] - # sans_detector is not populated in the current files - for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): - assert f'{name}_detector' in instr - det = instr[f'{name}_detector'] - assert 'pixel_shape' in det +def test_patch_nexus_detector_data(providers, params): + pipeline = sciline.Pipeline(params=params, providers=providers) + result = pipeline.compute(ReducibleDetectorData[SampleRun]) + assert ( + set(result.dims) == hr_sans_dims + if params[NeXusDetectorName] + in ( + 'high_resolution_detector', + 'sans_detector', + ) + else bank_dims + ) + assert "source_position" in result.coords + assert "sample_position" in result.coords From bfa37113128b6bc5f6382651c256bea3aab89639 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Thu, 23 May 2024 11:09:33 +0200 Subject: [PATCH 17/26] fix instrument_view --- src/ess/dream/instrument_view.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ess/dream/instrument_view.py b/src/ess/dream/instrument_view.py index 7a53a4d6..6b371444 100644 --- a/src/ess/dream/instrument_view.py +++ b/src/ess/dream/instrument_view.py @@ -115,7 +115,7 @@ def _add_module_control(self): import ipywidgets as ipw self.fig = self.scatter[0] - self.cutting_tool = self.scatter.bottom_bar[0] + self.cutting_tool = self.scatter[1] self.artist_mapping = dict( zip(self.data.keys(), self.fig.artists.keys(), strict=True) ) From 7a941cc9dd0860e5d33657a8b2c6b94e3f1ad194 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Thu, 23 May 2024 11:21:22 +0200 Subject: [PATCH 18/26] fix dream reduction notebook --- docs/user-guide/dream/dream-data-reduction.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user-guide/dream/dream-data-reduction.ipynb b/docs/user-guide/dream/dream-data-reduction.ipynb index 2455e7f0..f063f687 100644 --- a/docs/user-guide/dream/dream-data-reduction.ipynb +++ b/docs/user-guide/dream/dream-data-reduction.ipynb @@ -192,7 +192,7 @@ " )\n", ")\n", "\n", - "intermediates[DataWithScatteringCoordinates[SampleRun]].bins.coords.keys()" + "intermediates[DataWithScatteringCoordinates[SampleRun]]" ] }, { From ecca9ceca1ccffd15016e157f08ef8272df62a09 Mon Sep 17 00:00:00 2001 From: Neil Vaytet <39047984+nvaytet@users.noreply.github.com> Date: Wed, 29 May 2024 10:48:57 +0200 Subject: [PATCH 19/26] Update src/ess/dream/io/geant4.py Co-authored-by: Jan-Lukas Wynen --- src/ess/dream/io/geant4.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index e9ff9206..09bb5a4b 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -191,7 +191,7 @@ def patch_detector_data( def geant4_detector_dimensions( data: RawDetectorData[SampleRun], ) -> NeXusDetectorDimensions[NeXusDetectorName]: - # Fir geant4 data, we group by detector identifier, so the data already has + # For geant4 data, we group by detector identifier, so the data already has # logical dimensions, so we simply return the dimensions of the detector. return NeXusDetectorDimensions[NeXusDetectorName](data.sizes) From da001b804fea9166d15a8fad837a1941566bf9da Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 15:01:58 +0200 Subject: [PATCH 20/26] remove wavelength mask in dream notebook and also remove mention of basic providers --- docs/user-guide/dream/dream-data-reduction.ipynb | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/docs/user-guide/dream/dream-data-reduction.ipynb b/docs/user-guide/dream/dream-data-reduction.ipynb index f063f687..7aa35c2d 100644 --- a/docs/user-guide/dream/dream-data-reduction.ipynb +++ b/docs/user-guide/dream/dream-data-reduction.ipynb @@ -50,12 +50,9 @@ " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", " # Edges for binning in d-spacing\n", " DspacingBins: sc.linspace(\"dspacing\", 0.0, 2.3434, 201, unit=\"angstrom\"),\n", - " # Mask in time-of-flight\n", + " # Mask in time-of-flight to crop to valid range\n", " TofMask: lambda x: (x < sc.scalar(0.0, unit=\"ns\"))\n", " | (x > sc.scalar(86e6, unit=\"ns\")),\n", - " # Mask in wavelength\n", - " WavelengthMask: lambda x: (x > sc.scalar(2.0, unit=\"angstrom\"))\n", - " & (x < sc.scalar(2.7, unit=\"angstrom\")),\n", "}\n", "\n", "# Not available in simulated data\n", @@ -79,7 +76,7 @@ "source": [ "## Create pipeline using Sciline\n", "\n", - "We use the basic providers available in `essdiffraction` as well as the specialised `powder` and `geant4` providers." + "We use the `powder` and `geant4` providers to build our pipeline." ] }, { @@ -282,8 +279,7 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.12" + "pygments_lexer": "ipython3" } }, "nbformat": 4, From 5febee809cd7f47899017eb8bde5b31610ab33b5 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 15:04:59 +0200 Subject: [PATCH 21/26] remove wavelength mask from powgen notebook --- docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index 0a2ba67b..6d31fce2 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -59,12 +59,9 @@ " UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n", " # Edges for binning in d-spacing\n", " DspacingBins: sc.linspace(\"dspacing\", 0.0, 2.3434, 201, unit=\"angstrom\"),\n", - " # Mask in time-of-flight\n", + " # Mask in time-of-flight to crop to valid range\n", " TofMask: lambda x: (x < sc.scalar(0.0, unit=\"us\"))\n", " | (x > sc.scalar(16666.67, unit=\"us\")),\n", - " # Mask in wavelength\n", - " WavelengthMask: lambda x: (x > sc.scalar(0.18, unit=\"angstrom\"))\n", - " & (x < sc.scalar(0.21, unit=\"angstrom\")),\n", "}" ] }, From 0ab2596a2fbfed25332017661991545e047873c6 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 15:06:11 +0200 Subject: [PATCH 22/26] fix indentation --- src/ess/dream/io/nexus.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index beb9096c..d380d7d5 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -8,9 +8,9 @@ See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 and the ICD DREAM interface specification for details. - - The high-resolution and SANS detectors have a very odd numbering scheme. - The scheme attempts to follows some sort of physical ordering in space (x,y,z), - but it is not possible to reshape the data into all the logical dimensions. +- The high-resolution and SANS detectors have a very odd numbering scheme. + The scheme attempts to follows some sort of physical ordering in space (x,y,z), + but it is not possible to reshape the data into all the logical dimensions. """ import scipp as sc From b83ebb8fc15b892fa2f54f73eb77aca1fc19bb27 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 15:11:59 +0200 Subject: [PATCH 23/26] DETECTOR_BANK_* -> DETECTOR_BANK_SIZES --- src/ess/dream/io/nexus.py | 68 +++++++++------------- src/ess/powder/external/powgen/beamline.py | 4 +- 2 files changed, 30 insertions(+), 42 deletions(-) diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index d380d7d5..3de86b88 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -28,44 +28,32 @@ ) from ess.reduce import nexus -DETECTOR_BANK_RESHAPING = { - "endcap_backward_detector": lambda x: x.fold( - dim="detector_number", - sizes={ - "strip": 16, - "wire": 16, - "module": 11, - "segment": 28, - "counter": 2, - }, - ), - "endcap_forward_detector": lambda x: x.fold( - dim="detector_number", - sizes={ - "strip": 16, - "wire": 16, - "module": 5, - "segment": 28, - "counter": 2, - }, - ), - "mantle_detector": lambda x: x.fold( - dim="detector_number", - sizes={ - "wire": 32, - "module": 5, - "segment": 6, - "strip": 256, - "counter": 2, - }, - ), - "high_resolution_detector": lambda x: x.fold( - dim="detector_number", - sizes={ - "strip": 32, - "other": -1, - }, - ), +DETECTOR_BANK_SIZES = { + "endcap_backward_detector": { + "strip": 16, + "wire": 16, + "module": 11, + "segment": 28, + "counter": 2, + }, + "endcap_forward_detector": { + "strip": 16, + "wire": 16, + "module": 5, + "segment": 28, + "counter": 2, + }, + "mantle_detector": { + "wire": 32, + "module": 5, + "segment": 6, + "strip": 256, + "counter": 2, + }, + "high_resolution_detector": { + "strip": 32, + "other": -1, + }, "sans_detector": lambda x: x.fold( dim="detector_number", sizes={ @@ -118,8 +106,8 @@ def get_detector_data( detector_name: NeXusDetectorName, ) -> RawDetectorData[RunType]: da = nexus.extract_detector_data(detector) - if detector_name in DETECTOR_BANK_RESHAPING: - da = DETECTOR_BANK_RESHAPING[detector_name](da) + if detector_name in DETECTOR_BANK_SIZES: + da = da.fold(dim="detector_number", sizes=DETECTOR_BANK_SIZES[detector_name]) return RawDetectorData[RunType](da) diff --git a/src/ess/powder/external/powgen/beamline.py b/src/ess/powder/external/powgen/beamline.py index c20ae358..3c4e109f 100644 --- a/src/ess/powder/external/powgen/beamline.py +++ b/src/ess/powder/external/powgen/beamline.py @@ -14,7 +14,7 @@ ) from .types import DetectorInfo -DETECTOR_BANK_SHAPES = {"powgen_detector": {"bank": 23, "column": 154, "row": 7}} +DETECTOR_BANK_SIZES = {"powgen_detector": {"bank": 23, "column": 154, "row": 7}} def map_detector_to_spectrum( @@ -73,7 +73,7 @@ def powgen_detector_dimensions( ) -> NeXusDetectorDimensions[NeXusDetectorName]: """Dimensions used by POWGEN detectors.""" return NeXusDetectorDimensions[NeXusDetectorName]( - DETECTOR_BANK_SHAPES[detector_name] + DETECTOR_BANK_SIZES[detector_name] ) From 77e738e9b7092240ab6354d421515995afcbfda7 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 15:12:39 +0200 Subject: [PATCH 24/26] remove commented code --- src/ess/dream/data.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/ess/dream/data.py b/src/ess/dream/data.py index e4cc8f3e..00ae7e88 100644 --- a/src/ess/dream/data.py +++ b/src/ess/dream/data.py @@ -2,8 +2,6 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Data for tests and documentation with DREAM.""" -# from pathlib import Path - _version = "1" __all__ = ["get_path"] From 3d977a992bea074cbb4b160b5108c87b76b7591d Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 15:44:01 +0200 Subject: [PATCH 25/26] use different type to request data binned in dspacing and 2theta instead of using Optional[TwoThetaBins] --- .../dream/dream-data-reduction.ipynb | 11 +-- .../POWGEN_data_reduction.ipynb | 18 +++-- src/ess/powder/correction.py | 67 +++++++++++++++---- src/ess/powder/grouping.py | 54 +++++---------- src/ess/powder/types.py | 17 +++-- tests/dream/geant4_reduction_test.py | 11 +-- .../external/powgen/powgen_reduction_test.py | 15 +++-- 7 files changed, 115 insertions(+), 78 deletions(-) diff --git a/docs/user-guide/dream/dream-data-reduction.ipynb b/docs/user-guide/dream/dream-data-reduction.ipynb index 7aa35c2d..92bad4de 100644 --- a/docs/user-guide/dream/dream-data-reduction.ipynb +++ b/docs/user-guide/dream/dream-data-reduction.ipynb @@ -109,7 +109,7 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(NormalizedByVanadium, graph_attr={\"rankdir\": \"LR\"})" + "pipeline.visualize(IofDspacing, graph_attr={\"rankdir\": \"LR\"})" ] }, { @@ -127,7 +127,7 @@ "metadata": {}, "outputs": [], "source": [ - "result = pipeline.compute(NormalizedByVanadium)\n", + "result = pipeline.compute(IofDspacing)\n", "result" ] }, @@ -209,7 +209,10 @@ "id": "7f866ad4-5a0d-4c98-b5ab-a2436f97074d", "metadata": {}, "source": [ - "## Grouping by scattering angle" + "## Grouping by scattering angle\n", + "\n", + "The above pipeline focuses the data by merging all instrument pixels to produce a 1d d-spacing curve.\n", + "If instead we want to group into $2\\theta$ bins, we can alter the pipeline parameters by adding some binning in $2\\theta$:" ] }, { @@ -231,7 +234,7 @@ "metadata": {}, "outputs": [], "source": [ - "grouped_dspacing = pipeline.compute(NormalizedByVanadium)\n", + "grouped_dspacing = pipeline.compute(IofDspacingTwoTheta)\n", "grouped_dspacing" ] }, diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index 6d31fce2..477b6c7b 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -95,7 +95,7 @@ "\n", "### Compute final result\n", "\n", - "We can get the graph for computing the final d-spacing histogram:" + "We can get the graph for computing the final intensity as a function of d-spacing:" ] }, { @@ -105,7 +105,7 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(NormalizedByVanadium, graph_attr={\"rankdir\": \"LR\"})" + "pipeline.visualize(IofDspacing, graph_attr={\"rankdir\": \"LR\"})" ] }, { @@ -123,7 +123,7 @@ "metadata": {}, "outputs": [], "source": [ - "result = pipeline.compute(NormalizedByVanadium)\n", + "result = pipeline.compute(IofDspacing)\n", "result" ] }, @@ -293,6 +293,14 @@ ").to(unit=\"rad\")" ] }, + { + "cell_type": "markdown", + "id": "c8377091-43b5-467a-8f86-b58a96a3dc89", + "metadata": {}, + "source": [ + "We then have to request a final result that depends on both d-spacing and $2\\theta$:" + ] + }, { "cell_type": "code", "execution_count": null, @@ -300,7 +308,7 @@ "metadata": {}, "outputs": [], "source": [ - "pipeline.visualize(NormalizedByVanadium, graph_attr={\"rankdir\": \"LR\"})" + "pipeline.visualize(IofDspacingTwoTheta, graph_attr={\"rankdir\": \"LR\"})" ] }, { @@ -318,7 +326,7 @@ "metadata": {}, "outputs": [], "source": [ - "grouped_dspacing = pipeline.compute(NormalizedByVanadium)\n", + "grouped_dspacing = pipeline.compute(IofDspacingTwoTheta)\n", "grouped_dspacing" ] }, diff --git a/src/ess/powder/correction.py b/src/ess/powder/correction.py index 687dbbb6..46363e7b 100644 --- a/src/ess/powder/correction.py +++ b/src/ess/powder/correction.py @@ -13,9 +13,11 @@ from .types import ( AccumulatedProtonCharge, FilteredData, - FocussedData, + FocussedDataDspacing, + FocussedDataDspacingTwoTheta, + IofDspacing, + IofDspacingTwoTheta, NormalizedByProtonCharge, - NormalizedByVanadium, RunType, SampleRun, UncertaintyBroadcastMode, @@ -76,13 +78,26 @@ def normalize_by_monitor( return data.bins / sc.lookup(func=monitor, dim="wavelength") -def normalize_by_vanadium( - data: FocussedData[SampleRun], - vanadium: FocussedData[VanadiumRun], +def _normalize_by_vanadium( + data: sc.DataArray, + vanadium: sc.DataArray, + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> sc.DataArray: + vanadium = broadcast_uncertainties(vanadium, uncertainty_broadcast_mode) + norm = vanadium.hist() + # Converting to unit 'one' because the division might produce a unit + # with a large scale if the proton charges in data and vanadium were + # measured with different units. + return (data / norm).to(unit="one", copy=False) + + +def normalize_by_vanadium_dspacing( + data: FocussedDataDspacing[SampleRun], + vanadium: FocussedDataDspacing[VanadiumRun], uncertainty_broadcast_mode: Optional[UncertaintyBroadcastMode] = None, -) -> NormalizedByVanadium: +) -> IofDspacing: """ - Normalize sample data by a vanadium measurement. + Normalize sample data by a vanadium measurement and return intensity vs d-spacing. Parameters ---------- @@ -94,13 +109,33 @@ def normalize_by_vanadium( Choose how uncertainties of vanadium are broadcast to the sample data. Defaults to ``UncertaintyBroadcastMode.fail``. """ - vanadium = broadcast_uncertainties(vanadium, uncertainty_broadcast_mode) + return IofDspacing( + _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode) + ) - norm = vanadium.hist() - # Converting to unit 'one' because the division might produce a unit - # with a large scale if the proton charges in data and vanadium were - # measured with different units. - return NormalizedByVanadium((data / norm).to(unit="one", copy=False)) + +def normalize_by_vanadium_dspacing_and_two_theta( + data: FocussedDataDspacingTwoTheta[SampleRun], + vanadium: FocussedDataDspacingTwoTheta[VanadiumRun], + uncertainty_broadcast_mode: Optional[UncertaintyBroadcastMode] = None, +) -> IofDspacingTwoTheta: + """ + Normalize sample data by a vanadium measurement and return intensity vs + (d-spacing, 2theta). + + Parameters + ---------- + data: + Sample data. + vanadium: + Vanadium data. + uncertainty_broadcast_mode: + Choose how uncertainties of vanadium are broadcast to the sample data. + Defaults to ``UncertaintyBroadcastMode.fail``. + """ + return IofDspacingTwoTheta( + _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode) + ) def normalize_by_proton_charge( @@ -222,5 +257,9 @@ def _shallow_copy(da: sc.DataArray) -> sc.DataArray: return out -providers = (normalize_by_proton_charge, normalize_by_vanadium) +providers = ( + normalize_by_proton_charge, + normalize_by_vanadium_dspacing, + normalize_by_vanadium_dspacing_and_two_theta, +) """Sciline providers for powder diffraction corrections.""" diff --git a/src/ess/powder/grouping.py b/src/ess/powder/grouping.py index 974e14e8..ecd74092 100644 --- a/src/ess/powder/grouping.py +++ b/src/ess/powder/grouping.py @@ -2,56 +2,34 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) """Grouping and merging of pixels / voxels.""" -from typing import Optional - from .types import ( DspacingBins, - DspacingHistogram, - FocussedData, + FocussedDataDspacing, + FocussedDataDspacingTwoTheta, MaskedData, - NormalizedByVanadium, RunType, TwoThetaBins, ) -def focus_data( +def focus_data_dspacing( data: MaskedData[RunType], dspacing_bins: DspacingBins, - twotheta_bins: Optional[TwoThetaBins] = None, -) -> FocussedData[RunType]: - bins = {} - if twotheta_bins is not None: - bins["two_theta"] = twotheta_bins - bins[dspacing_bins.dim] = dspacing_bins - - if (twotheta_bins is None) or ("two_theta" in data.bins.coords): - data = data.bins.concat() - - return FocussedData[RunType](data.bin(**bins)) - +) -> FocussedDataDspacing[RunType]: + out = data.bins.concat().bin({dspacing_bins.dim: dspacing_bins}) + return FocussedDataDspacing[RunType](out) -def finalize_histogram( - data: NormalizedByVanadium, edges: DspacingBins -) -> DspacingHistogram: - """Finalize the d-spacing histogram. - Histograms the input data into the given d-spacing bins. - - Parameters - ---------- - data: - Data to be histogrammed. - edges: - Bin edges in d-spacing. - - Returns - ------- - : - Histogrammed data. - """ - return DspacingHistogram(data.hist(dspacing=edges)) +def focus_data_dspacing_and_two_theta( + data: MaskedData[RunType], + dspacing_bins: DspacingBins, + twotheta_bins: TwoThetaBins, +) -> FocussedDataDspacingTwoTheta[RunType]: + bins = {twotheta_bins.dim: twotheta_bins, dspacing_bins.dim: dspacing_bins} + if "two_theta" in data.bins.coords: + data = data.bins.concat() + return FocussedDataDspacingTwoTheta[RunType](data.bin(**bins)) -providers = (finalize_histogram, focus_data) +providers = (focus_data_dspacing, focus_data_dspacing_and_two_theta) """Sciline providers for grouping pixels.""" diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 89c2226c..8887af7e 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -114,10 +114,21 @@ class FilteredData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Raw data without invalid events.""" -class FocussedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): +class FocussedDataDspacing(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Intensity vs d-spacing after focussing pixels.""" +class FocussedDataDspacingTwoTheta(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Intensity vs (d-spacing, 2theta) after focussing pixels.""" + + +IofDspacing = NewType("IofDspacing", sc.DataArray) +"""Data that has been normalized by a vanadium run.""" + +IofDspacingTwoTheta = NewType("IofDspacingTwoTheta", sc.DataArray) +"""Data that has been normalized by a vanadium run, and grouped into 2theta bins.""" + + class LoadedNeXusDetector(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Detector data, loaded from a NeXus file, containing not only neutron events but also pixel shape information, transformations, ...""" @@ -136,10 +147,6 @@ class NormalizedByProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArra """Data that has been normalized by proton charge.""" -NormalizedByVanadium = NewType("NormalizedByVanadium", sc.DataArray) -"""Data that has been normalized by a vanadium run.""" - - PixelMaskFilename = NewType("PixelMaskFilename", str) """Filename of a pixel mask.""" diff --git a/tests/dream/geant4_reduction_test.py b/tests/dream/geant4_reduction_test.py index 57c48e13..085fa5a0 100644 --- a/tests/dream/geant4_reduction_test.py +++ b/tests/dream/geant4_reduction_test.py @@ -7,9 +7,10 @@ from ess.powder.types import ( AccumulatedProtonCharge, DspacingBins, - DspacingHistogram, EmptyCanRun, Filename, + IofDspacing, + IofDspacingTwoTheta, MaskedData, NeXusDetectorName, NormalizedByProtonCharge, @@ -64,9 +65,9 @@ def test_can_create_pipeline(providers, params): sciline.Pipeline(providers, params=params) -def test_pipeline_can_compute_dspacing_histogram(providers, params): +def test_pipeline_can_compute_dspacing_result(providers, params): pipeline = sciline.Pipeline(providers, params=params) - result = pipeline.compute(DspacingHistogram) + result = pipeline.compute(IofDspacing) assert result.sizes == { 'dspacing': len(params[DspacingBins]) - 1, } @@ -77,7 +78,7 @@ def test_workflow_is_deterministic(providers, params): pipeline = sciline.Pipeline(providers, params=params) # This is Sciline's default scheduler, but we want to be explicit here scheduler = sciline.scheduler.DaskScheduler() - graph = pipeline.get(DspacingHistogram, scheduler=scheduler) + graph = pipeline.get(IofDspacing, scheduler=scheduler) reference = graph.compute().data result = graph.compute().data assert sc.identical(sc.values(result), sc.values(reference)) @@ -94,7 +95,7 @@ def test_pipeline_group_by_two_theta(providers, params): dim='two_theta', unit='rad', start=0.8, stop=2.4, num=17 ) pipeline = sciline.Pipeline(providers, params=params) - result = pipeline.compute(DspacingHistogram) + result = pipeline.compute(IofDspacingTwoTheta) assert result.sizes == { 'two_theta': 16, 'dspacing': len(params[DspacingBins]) - 1, diff --git a/tests/powder/external/powgen/powgen_reduction_test.py b/tests/powder/external/powgen/powgen_reduction_test.py index e9564010..a723a0e4 100644 --- a/tests/powder/external/powgen/powgen_reduction_test.py +++ b/tests/powder/external/powgen/powgen_reduction_test.py @@ -7,8 +7,9 @@ from ess.powder.types import ( CalibrationFilename, DspacingBins, - DspacingHistogram, Filename, + IofDspacing, + IofDspacingTwoTheta, MaskedData, NeXusDetectorName, NormalizedByProtonCharge, @@ -50,9 +51,9 @@ def test_can_create_pipeline(providers, params): sciline.Pipeline(providers, params=params) -def test_pipeline_can_compute_dspacing_histogram(providers, params): +def test_pipeline_can_compute_dspacing_result(providers, params): pipeline = sciline.Pipeline(providers, params=params) - result = pipeline.compute(DspacingHistogram) + result = pipeline.compute(IofDspacing) assert result.sizes == { 'dspacing': len(params[DspacingBins]) - 1, } @@ -63,9 +64,9 @@ def test_workflow_is_deterministic(providers, params): pipeline = sciline.Pipeline(providers, params=params) # This is Sciline's default scheduler, but we want to be explicit here scheduler = sciline.scheduler.DaskScheduler() - graph = pipeline.get(DspacingHistogram, scheduler=scheduler) - reference = graph.compute().data - result = graph.compute().data + graph = pipeline.get(IofDspacing, scheduler=scheduler) + reference = graph.compute().hist().data + result = graph.compute().hist().data assert sc.identical(sc.values(result), sc.values(reference)) @@ -80,7 +81,7 @@ def test_pipeline_group_by_two_theta(providers, params): dim='two_theta', unit='deg', start=25.0, stop=90.0, num=16 ).to(unit='rad') pipeline = sciline.Pipeline(providers, params=params) - result = pipeline.compute(DspacingHistogram) + result = pipeline.compute(IofDspacingTwoTheta) assert result.sizes == { 'two_theta': 15, 'dspacing': len(params[DspacingBins]) - 1, From ea5f681e38fe4784b2e644e88b52096cbb2c9dd3 Mon Sep 17 00:00:00 2001 From: Neil Vaytet Date: Wed, 29 May 2024 16:01:54 +0200 Subject: [PATCH 26/26] fix docs --- docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb index 477b6c7b..c6c55efb 100644 --- a/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb +++ b/docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb @@ -161,7 +161,7 @@ "outputs": [], "source": [ "def save_xye(\n", - " reduced_data: NormalizedByVanadium,\n", + " reduced_data: IofDspacing,\n", " out_filename: OutFilename,\n", ") -> None:\n", " data = reduced_data.hist()\n",