Merge pull request #67 from scipp/masking-with-calibration

Compute d-spacing from calibration

docs/user-guide/dream/dream-data-reduction.ipynb

@@ -45,6 +45,7 @@
     "    Filename[SampleRun]: dream.data.simulated_diamond_sample(),\n",
     "    Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(),\n",
     "    Filename[EmptyCanRun]: dream.data.simulated_empty_can(),\n",
+    "    CalibrationFilename: None,\n",
     "    NeXusDetectorName: \"mantle\",\n",
     "    # The upper bounds mode is not yet implemented.\n",
     "    UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,\n",
@@ -286,7 +287,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,

docs/user-guide/sns-instruments/POWGEN_data_reduction.ipynb

@@ -386,7 +386,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,

src/ess/dream/io/geant4.py

@@ -7,6 +7,8 @@
 import sciline
 import scipp as sc
 from ess.powder.types import (
+    CalibrationData,
+    CalibrationFilename,
     Filename,
     NeXusDetectorDimensions,
     NeXusDetectorName,
@@ -196,6 +198,17 @@ def geant4_detector_dimensions(
     return NeXusDetectorDimensions[NeXusDetectorName](data.sizes)
 
 
+def geant4_load_calibration(
+    filename: CalibrationFilename,
+) -> CalibrationData:
+    if filename is not None:
+        # Needed to build a complete pipeline.
+        raise NotImplementedError(
+            "Calibration data loading is not implemented for DREAM GEANT4 data."
+        )
+    return CalibrationFilename(None)
+
+
 providers = (
     extract_geant4_detector,
     extract_geant4_detector_data,
@@ -204,5 +217,6 @@ def geant4_detector_dimensions(
     get_source_position,
     patch_detector_data,
     geant4_detector_dimensions,
+    geant4_load_calibration,
 )
 """Geant4-providers for Sciline pipelines."""

src/ess/powder/__init__.py

@@ -25,7 +25,7 @@
 del importlib
 
 providers = (
-    *conversion.providers_with_positions,
+    *conversion.providers,
     *correction.providers,
     *filtering.providers,
     *grouping.providers,

src/ess/powder/conversion.py

@@ -14,6 +14,7 @@
     DataWithScatteringCoordinates,
     DspacingData,
     ElasticCoordTransformGraph,
+    MaskedData,
     NormalizedByProtonCharge,
     RunType,
 )
@@ -83,9 +84,9 @@ def _consume_positions(position, sample_position, source_position):
 
 
 def to_dspacing_with_calibration(
-    data: NormalizedByProtonCharge[RunType],
-    calibration: CalibrationData,
-) -> DspacingData[RunType]:
+    data: sc.DataArray,
+    calibration: sc.Dataset,
+) -> sc.DataArray:
     """
     Transform coordinates to d-spacing from calibration parameters.
 
@@ -177,7 +178,7 @@ def add_scattering_coordinates_from_positions(
     data: NormalizedByProtonCharge[RunType], graph: ElasticCoordTransformGraph
 ) -> DataWithScatteringCoordinates[RunType]:
     """
-    Add ``wavelength``, ``two_theta``, and ``dspacing`` coordinates to the data.
+    Add ``wavelength`` and ``two_theta`` 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.
@@ -190,15 +191,28 @@ def add_scattering_coordinates_from_positions(
         Coordinate transformation graph.
     """
     out = data.transform_coords(
-        ["two_theta", "wavelength", "dspacing"], graph=graph, keep_intermediate=False
+        ["two_theta", "wavelength"], graph=graph, keep_intermediate=False
     )
     return DataWithScatteringCoordinates[RunType](out)
 
 
-providers_with_calibration = (to_dspacing_with_calibration,)
-"""Sciline providers for coordinate transformations."""
+def convert_to_dspacing(
+    data: MaskedData[RunType],
+    graph: ElasticCoordTransformGraph,
+    calibration: CalibrationData,
+) -> DspacingData[RunType]:
+    if calibration is None:
+        out = data.transform_coords(["dspacing"], graph=graph, keep_intermediate=False)
+        return DspacingData[RunType](out)
+    out = to_dspacing_with_calibration(data, calibration=calibration)
+    for key in ('wavelength', 'two_theta'):
+        if key in out.coords.keys():
+            out.coords.set_aligned(key, False)
+    return DspacingData[RunType](out)
+
 
-providers_with_positions = (
+providers = (
     powder_coordinate_transformation_graph,
     add_scattering_coordinates_from_positions,
+    convert_to_dspacing,
 )

src/ess/powder/correction.py

@@ -179,11 +179,11 @@ def merge_calibration(*, into: sc.DataArray, calibration: sc.Dataset) -> sc.Data
     --------
     ess.powder.load_calibration
     """
-    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."
-        )
+    for name, coord in calibration.coords.items():
+        if not sc.identical(into.coords[name], coord):
+            raise ValueError(
+                f"Coordinate {name} of calibration and target dataset do not agree."
+            )
     out = into.copy(deep=False)
     for name in ("difa", "difc", "tzero"):
         if name in out.coords:

src/ess/powder/external/powgen/beamline.py

@@ -7,12 +7,9 @@
 import scipp as sc
 
 from ...types import (
-    CalibrationData,
     NeXusDetectorDimensions,
     NeXusDetectorName,
-    RawCalibrationData,
 )
-from .types import DetectorInfo
 
 DETECTOR_BANK_SIZES = {"powgen_detector": {"bank": 23, "column": 154, "row": 7}}
 
@@ -56,18 +53,6 @@ def map_detector_to_spectrum(
     return out.rename_dims({"detector": "spectrum"})
 
 
-def preprocess_calibration_data(
-    data: RawCalibrationData, detector_info: DetectorInfo
-) -> CalibrationData:
-    """Convert calibration data to a format that can be used by Scipp.
-
-    The raw calibration data is encoded in terms of a `'detector'` coordinate.
-    This needs to be converted to a `'spectrum'` coordinate to align
-    if with sample data.
-    """
-    return CalibrationData(map_detector_to_spectrum(data, detector_info=detector_info))
-
-
 def powgen_detector_dimensions(
     detector_name: NeXusDetectorName,
 ) -> NeXusDetectorDimensions[NeXusDetectorName]:
@@ -77,5 +62,5 @@ def powgen_detector_dimensions(
     )
 
 
-providers = (preprocess_calibration_data, powgen_detector_dimensions)
+providers = (powgen_detector_dimensions,)
 """Sciline providers for POWGEN beamline processing."""

src/ess/powder/external/powgen/data.py

@@ -7,18 +7,16 @@
 
 from ...types import (
     AccumulatedProtonCharge,
+    CalibrationData,
     CalibrationFilename,
     Filename,
     NeXusDetectorDimensions,
     NeXusDetectorName,
     ProtonCharge,
-    RawCalibrationData,
     RawDataAndMetadata,
     ReducibleDetectorData,
     RunType,
-    SampleRun,
 )
-from .types import DetectorInfo
 
 _version = "1"
 
@@ -41,6 +39,7 @@ def _make_pooch():
             "PG3_4844_event.zip": "md5:a644c74f5e740385469b67431b690a3e",
             "PG3_4866_event.zip": "md5:5bc49def987f0faeb212a406b92b548e",
             "PG3_FERNS_d4832_2011_08_24.zip": "md5:0fef4ed5f61465eaaa3f87a18f5bb80d",
+            "PG3_FERNS_d4832_2011_08_24_spectrum.h5": "md5:7aee0b40deee22d57e21558baa7a6a1a",  # noqa: E501
         },
     )
 
@@ -89,7 +88,7 @@ def powgen_tutorial_vanadium_file() -> str:
 
 
 def powgen_tutorial_calibration_file() -> str:
-    return _get_path("PG3_FERNS_d4832_2011_08_24.zip")
+    return _get_path("PG3_FERNS_d4832_2011_08_24_spectrum.h5")
 
 
 def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]:
@@ -103,9 +102,21 @@ def pooch_load(filename: Filename[RunType]) -> RawDataAndMetadata[RunType]:
     return RawDataAndMetadata[RunType](sc.io.load_hdf5(filename))
 
 
-def pooch_load_calibration(filename: CalibrationFilename) -> RawCalibrationData:
+def pooch_load_calibration(
+    filename: CalibrationFilename,
+    detector_dimensions: NeXusDetectorDimensions[NeXusDetectorName],
+) -> CalibrationData:
     """Load the calibration data for the POWGEN test data."""
-    return RawCalibrationData(sc.io.load_hdf5(filename))
+    if filename is None:
+        return CalibrationFilename(None)
+    ds = sc.io.load_hdf5(filename)
+    ds = sc.Dataset(
+        {
+            key: da.fold(dim='spectrum', sizes=detector_dimensions)
+            for key, da in ds.items()
+        }
+    )
+    return CalibrationData(ds)
 
 
 def extract_raw_data(
@@ -120,11 +131,6 @@ def extract_raw_data(
     return ReducibleDetectorData[RunType](out)
 
 
-def extract_detector_info(dg: RawDataAndMetadata[SampleRun]) -> DetectorInfo:
-    """Return the detector info from a loaded data group."""
-    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"])
@@ -141,7 +147,6 @@ def extract_accumulated_proton_charge(
     pooch_load,
     pooch_load_calibration,
     extract_accumulated_proton_charge,
-    extract_detector_info,
     extract_proton_charge,
     extract_raw_data,
 )

src/ess/powder/grouping.py

@@ -4,24 +4,24 @@
 
 from .types import (
     DspacingBins,
+    DspacingData,
     FocussedDataDspacing,
     FocussedDataDspacingTwoTheta,
-    MaskedData,
     RunType,
     TwoThetaBins,
 )
 
 
 def focus_data_dspacing(
-    data: MaskedData[RunType],
+    data: DspacingData[RunType],
     dspacing_bins: DspacingBins,
 ) -> FocussedDataDspacing[RunType]:
     out = data.bins.concat().bin({dspacing_bins.dim: dspacing_bins})
     return FocussedDataDspacing[RunType](out)
 
 
 def focus_data_dspacing_and_two_theta(
-    data: MaskedData[RunType],
+    data: DspacingData[RunType],
     dspacing_bins: DspacingBins,
     twotheta_bins: TwoThetaBins,
 ) -> FocussedDataDspacingTwoTheta[RunType]:

src/ess/powder/types.py

@@ -30,7 +30,7 @@
 
 # 2 Workflow parameters
 
-CalibrationFilename = NewType("CalibrationFilename", str)
+CalibrationFilename = NewType("CalibrationFilename", str | None)
 """Filename of the instrument calibration file."""
 
 
@@ -78,7 +78,7 @@ def __init__(self, *args: Any, **kwargs: Any) -> None:
         super().__init__(*args, **kwargs)
 
 
-CalibrationData = NewType("CalibrationData", sc.Dataset)
+CalibrationData = NewType("CalibrationData", sc.Dataset | None)
 """Detector calibration data."""
 
 DataFolder = NewType("DataFolder", str)
@@ -155,10 +155,6 @@ class ProtonCharge(sciline.Scope[RunType, sc.DataArray], sc.DataArray):
     """Time-dependent proton charge."""
 
 
-RawCalibrationData = NewType("RawCalibrationData", sc.Dataset)
-"""Calibration data as loaded from file, needs preprocessing before using."""
-
-
 class RawDataAndMetadata(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup):
     """Raw data and associated metadata."""
 

tests/dream/geant4_reduction_test.py

@@ -7,6 +7,7 @@
 from ess import powder
 from ess.powder.types import (
     AccumulatedProtonCharge,
+    CalibrationFilename,
     DspacingBins,
     EmptyCanRun,
     Filename,
@@ -49,6 +50,7 @@ def params(request):
         Filename[SampleRun]: dream.data.simulated_diamond_sample(),
         Filename[VanadiumRun]: dream.data.simulated_vanadium_sample(),
         Filename[EmptyCanRun]: dream.data.simulated_empty_can(),
+        CalibrationFilename: None,
         UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop,
         DspacingBins: sc.linspace('dspacing', 0.0, 2.3434, 201, unit='angstrom'),
         TofMask: lambda x: (x < sc.scalar(0.0, unit='ns'))

tests/powder/conversion_test.py

@@ -162,4 +162,3 @@ def test_add_scattering_coordinates_from_positions():
 
     assert 'wavelength' in result.coords
     assert 'two_theta' in result.coords
-    assert 'dspacing' in result.coords