Merge pull request #212 from scipp/fix-run-mapping

Fix run mapping functionality

Author: SimonHeybrock

src/ess/isissans/__init__.py

@@ -4,7 +4,13 @@
 import importlib.metadata
 
 from . import general, io, sans2d, zoom
-from .general import DetectorBankOffset, MonitorOffset, SampleOffset, default_parameters
+from .general import (
+    DetectorBankOffset,
+    MonitorOffset,
+    MonitorSpectrumNumber,
+    SampleOffset,
+    default_parameters,
+)
 from .io import CalibrationFilename
 from .visualization import plot_flat_detector_xy
 
@@ -21,6 +27,7 @@
     'CalibrationFilename',
     'DetectorBankOffset',
     'MonitorOffset',
+    'MonitorSpectrumNumber',
     'SampleOffset',
     'default_parameters',
     'io',

src/ess/isissans/general.py

@@ -4,7 +4,8 @@
 Providers for the ISIS instruments.
 """
 
-from typing import NewType
+from dataclasses import dataclass
+from typing import Generic, NewType
 
 import sciline
 import scipp as sc
@@ -45,6 +46,17 @@ class MonitorOffset(sciline.Scope[MonitorType, sc.Variable], sc.Variable):
     """
 
 
+@dataclass
+class MonitorSpectrumNumber(Generic[MonitorType]):
+    """
+    Spectrum number for monitor data.
+
+    This is used to identify the monitor in ISIS histogram data files.
+    """
+
+    value: int
+
+
 DetectorBankOffset = NewType('DetectorBankOffset', sc.Variable)
 SampleOffset = NewType('SampleOffset', sc.Variable)
 
@@ -61,6 +73,10 @@ def default_parameters() -> dict:
         SampleOffset: SampleOffset(sc.vector([0, 0, 0], unit='m')),
         NonBackgroundWavelengthRange: None,
         Period: None,
+        # Not used for event files, setting default so the workflow works. If histogram
+        # data is used, the user should set this to the correct value.
+        MonitorSpectrumNumber[Incident]: MonitorSpectrumNumber[Incident](-1),
+        MonitorSpectrumNumber[Transmission]: MonitorSpectrumNumber[Transmission](-1),
     }
 
 
@@ -132,12 +148,43 @@ def get_calibrated_isis_detector(
 
 
 def get_monitor_data(
-    dg: LoadedFileContents[RunType], nexus_name: NeXusMonitorName[MonitorType]
+    dg: LoadedFileContents[RunType],
+    nexus_name: NeXusMonitorName[MonitorType],
+    spectrum_number: MonitorSpectrumNumber[MonitorType],
 ) -> NeXusComponent[MonitorType, RunType]:
+    """
+    Extract monitor data that was loaded together with detector data.
+
+    If the raw file is histogram data, Mantid stores this as a Workspace2D, where some
+    or all spectra correspond to monitors.
+
+    Parameters
+    ----------
+    dg:
+        Data loaded with Mantid and converted to Scipp.
+    nexus_name:
+        Name of the monitor in the NeXus file, e.g. 'incident_monitor' or
+        'transmission_monitor'. Used when raw data is an EventWorkspace, where
+        monitors are stored in a group with this name.
+    spectrum_number:
+        Spectrum number of the monitor in the NeXus file, e.g., 3 for incident monitor
+        or 4 for transmission monitor. This is used when the raw data is a
+        Workspace2D, where the monitor data is stored in a spectrum with this number.
+
+    Returns
+    -------
+    :
+        Monitor data extracted from the loaded file.
+    """
     # The generic NeXus workflow will try to extract 'data' from this, which is exactly
     # what we also have in the Mantid data. We use the generic workflow since it also
     # applies offsets, etc.
-    monitor = dg['monitors'][nexus_name]['data']
+    if 'monitors' in dg:
+        # From EventWorkspace
+        monitor = dg['monitors'][nexus_name]['data']
+    else:
+        # From Workspace2D
+        monitor = sc.values(dg["data"]["spectrum", sc.index(spectrum_number.value)])
     return NeXusComponent[MonitorType, RunType](
         sc.DataGroup(data=monitor, position=monitor.coords['position'])
     )

src/ess/isissans/mantidio.py

@@ -4,6 +4,7 @@
 File loading functions for ISIS data using Mantid.
 """
 
+import threading
 from typing import NewType, NoReturn
 
 import sciline
@@ -99,12 +100,28 @@ def from_data_workspace(
 
 
 def load_run(filename: Filename[RunType], period: Period) -> DataWorkspace[RunType]:
-    loaded = _mantid_simpleapi.Load(
-        Filename=str(filename), LoadMonitors=True, StoreInADS=False
-    )
+    # Loading many small files with Mantid is, for some reason, very slow when using
+    # the default number of threads in the Dask threaded scheduler (1 thread worked
+    # best, 2 is a bit slower but still fast). We can either limit that thread count,
+    # or add a lock here, which is more specific.
+    with load_run.lock:
+        try:
+            loaded = _mantid_simpleapi.Load(
+                Filename=str(filename), LoadMonitors=True, StoreInADS=False
+            )
+        except TypeError:
+            # Not loaded using LoadEventNexus, so LoadMonitor option is not available.
+            loaded = _mantid_simpleapi.Load(Filename=str(filename), StoreInADS=False)
     if isinstance(loaded, _mantid_api.Workspace):
         # A single workspace
         data_ws = loaded
+        if isinstance(data_ws, _mantid_api.WorkspaceGroup):
+            if period is None:
+                raise ValueError(
+                    f'Needs {Period} to be set to know what '
+                    'section of the event data to load'
+                )
+            data_ws = data_ws.getItem(period)
     else:
         # Separate data and monitor workspaces
         data_ws = loaded.OutputWorkspace
@@ -121,6 +138,8 @@ def load_run(filename: Filename[RunType], period: Period) -> DataWorkspace[RunTy
     return DataWorkspace[RunType](data_ws)
 
 
+load_run.lock = threading.Lock()
+
 providers = (
     from_data_workspace,
     load_calibration,

src/ess/sans/workflow.py

@@ -12,7 +12,6 @@
 from . import common, conversions, i_of_q, masking, normalization
 from .types import (
     BackgroundRun,
-    CleanSummedQ,
     CorrectForGravity,
     Denominator,
     DetectorBankSizes,
@@ -30,6 +29,7 @@
     TransmissionRun,
     WavelengthBands,
     WavelengthMask,
+    WavelengthScaledQ,
 )
 
 
@@ -99,8 +99,8 @@ def _set_runs(
     pipeline = pipeline.copy()
     runs = pd.DataFrame({Filename[key]: runs}).rename_axis(axis_name)
     for part in (Numerator, Denominator):
-        pipeline[CleanSummedQ[key, part]] = (
-            pipeline[CleanSummedQ[key, part]]
+        pipeline[WavelengthScaledQ[key, part]] = (
+            pipeline[WavelengthScaledQ[key, part]]
             .map(runs)
             .reduce(index=axis_name, func=merge_contributions)
         )

tests/loki/iofq_test.py

@@ -16,7 +16,6 @@
     BackgroundSubtractedIofQ,
     BackgroundSubtractedIofQxy,
     BeamCenter,
-    CleanSummedQ,
     CleanWavelength,
     CorrectForGravity,
     Denominator,
@@ -33,6 +32,7 @@
     UncertaintyBroadcastMode,
     WavelengthBands,
     WavelengthBins,
+    WavelengthScaledQ,
 )
 
 sys.path.insert(0, str(Path(__file__).resolve().parent))
@@ -252,17 +252,21 @@ def test_pipeline_IofQ_merging_events_yields_consistent_results():
     assert all(sc.variances(iofq1.data) > sc.variances(iofq3.data))
     assert sc.allclose(
         sc.values(
-            pipeline_single.compute(CleanSummedQ[SampleRun, Numerator]).hist().data
+            pipeline_single.compute(WavelengthScaledQ[SampleRun, Numerator]).hist().data
         )
         * N,
         sc.values(
-            pipeline_triple.compute(CleanSummedQ[SampleRun, Numerator]).hist().data
+            pipeline_triple.compute(WavelengthScaledQ[SampleRun, Numerator]).hist().data
         ),
     )
     assert sc.allclose(
-        sc.values(pipeline_single.compute(CleanSummedQ[SampleRun, Denominator]).data)
+        sc.values(
+            pipeline_single.compute(WavelengthScaledQ[SampleRun, Denominator]).data
+        )
         * N,
-        sc.values(pipeline_triple.compute(CleanSummedQ[SampleRun, Denominator]).data),
+        sc.values(
+            pipeline_triple.compute(WavelengthScaledQ[SampleRun, Denominator]).data
+        ),
     )