Merge pull request #144 from scipp/subtract-empty-can

Subtract sample run by empty can

Author: SimonHeybrock

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

@@ -56,7 +56,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "workflow = dream.DreamGeant4Workflow(run_norm=powder.RunNormalization.monitor_histogram)"
+    "workflow = dream.DreamGeant4Workflow(\n",
+    "    run_norm=powder.RunNormalization.monitor_histogram,\n",
+    "    subtract_empty_can=True,\n",
+    ")"
    ]
   },
   {

src/ess/dream/data.py

@@ -29,9 +29,10 @@ def _make_pooch():
             "DREAM_simple_pwd_workflow/Cave_TOF_Monitor_vana_inc_coh.dat": "md5:701d66792f20eb283a4ce76bae0c8f8f",  # noqa: E501
             "DREAM-high-flux-tof-lookup-table.h5": "md5:1b95a359fa7b0d8b4277806ece9bf279",  # noqa: E501
             # Smaller files for unit tests
-            "DREAM_simple_pwd_workflow/TEST_data_dream_diamond_vana_container_sample_union.csv.zip": "md5:018a87e0934c1dd0f07a708e9d497891",  # noqa: E501
-            "DREAM_simple_pwd_workflow/TEST_data_dream_vana_container_sample_union.csv.zip": "md5:d244126cd8012f9ed186f4e08a19f88d",  # noqa: E501
-            "DREAM_simple_pwd_workflow/TEST_data_dream_vanadium.csv.zip": "md5:178f9bea9f35dbdef693e38ff893c258",  # noqa: E501
+            "DREAM_simple_pwd_workflow/TEST_data_dream_diamond_vana_container_sample_union.csv.zip": "md5:405df9b5ade9d61ab71fe8d8c19bb51b",  # noqa: E501
+            "DREAM_simple_pwd_workflow/TEST_data_dream_vana_container_sample_union.csv.zip": "md5:20186119d1debfb0c2352f9db384cd0a",  # noqa: E501
+            "DREAM_simple_pwd_workflow/TEST_data_dream_vanadium.csv.zip": "md5:e5624a973f30dfe440642319c5bf0da6",  # noqa: E501
+            "DREAM_simple_pwd_workflow/TEST_data_dream_vanadium_inc_coh.csv.zip": "md5:a2e7756b264f65ab5ed7ff7fb5eca280",  # noqa: E501
             "TEST_data_dream0_new_hkl_Si_pwd.csv.zip": "md5:df6c41f4b7b21e129915808f625828f6",  # noqa: E501
             "TEST_data_dream_with_sectors.csv.zip": "md5:2a6b5e40e6b67f6c71b25373bf4b11a1",  # noqa: E501
             # The TEST_DREAM_nexus_sorted-2023-12-07.nxs file was created using the
@@ -147,7 +148,9 @@ def simulated_vanadium_sample(*, small: bool = False) -> str:
         ```
     """
     prefix = "TEST_" if small else ""
-    return get_path(f"DREAM_simple_pwd_workflow/{prefix}data_dream_vanadium.csv.zip")
+    return get_path(
+        f"DREAM_simple_pwd_workflow/{prefix}data_dream_vanadium_inc_coh.csv.zip"
+    )
 
 
 def simulated_vanadium_sample_incoherent(*, small: bool = False) -> str:

src/ess/dream/io/geant4.py

@@ -8,6 +8,7 @@
 from scippneutron.metadata import ESS_SOURCE
 
 from ess.powder.types import (
+    BackgroundRun,
     Beamline,
     CalibratedBeamline,
     CalibratedDetector,
@@ -320,7 +321,7 @@ def LoadGeant4Workflow() -> sciline.Pipeline:
     Workflow for loading NeXus data.
     """
     wf = GenericTofWorkflow(
-        run_types=[SampleRun, VanadiumRun], monitor_types=[CaveMonitor]
+        run_types=[SampleRun, VanadiumRun, BackgroundRun], monitor_types=[CaveMonitor]
     )
     wf.insert(extract_geant4_detector)
     wf.insert(load_geant4_csv)

src/ess/dream/workflow.py

@@ -15,10 +15,12 @@
 from ess.powder.conversion import convert_monitor_to_wavelength
 from ess.powder.correction import (
     RunNormalization,
+    add_empty_can_subtraction,
     insert_run_normalization,
 )
 from ess.powder.types import (
     AccumulatedProtonCharge,
+    BackgroundRun,
     CaveMonitorPosition,  # Should this be a DREAM-only parameter?
     MonitorType,
     PixelMaskFilename,
@@ -75,10 +77,13 @@ def default_parameters() -> dict:
     return {
         Position[snx.NXsample, SampleRun]: sample_position,
         Position[snx.NXsample, VanadiumRun]: sample_position,
+        Position[snx.NXsample, BackgroundRun]: sample_position,
         Position[snx.NXsource, SampleRun]: source_position,
         Position[snx.NXsource, VanadiumRun]: source_position,
+        Position[snx.NXsource, BackgroundRun]: source_position,
         AccumulatedProtonCharge[SampleRun]: charge,
         AccumulatedProtonCharge[VanadiumRun]: charge,
+        AccumulatedProtonCharge[BackgroundRun]: charge,
         TofMask: None,
         WavelengthMask: None,
         TwoThetaMask: None,
@@ -110,16 +115,34 @@ def convert_dream_monitor_to_wavelength(
     return convert_monitor_to_wavelength(monitor)
 
 
-def DreamGeant4Workflow(*, run_norm: RunNormalization) -> sciline.Pipeline:
+def DreamGeant4Workflow(
+    *, run_norm: RunNormalization, subtract_empty_can: bool = False
+) -> sciline.Pipeline:
     """
     Workflow with default parameters for the Dream Geant4 simulation.
+
+    Parameters
+    ----------
+    run_norm:
+        Select how to normalize each run (sample, vanadium, etc.).
+    subtract_empty_can:
+        If ``True``, subtract the same data by an empty can / empty instrument
+        measurement before normalizing by vanadium.
+        This requires specifying a filename parameter for the empty can run.
+
+    Returns
+    -------
+    :
+        A workflow object for DREAM.
     """
     wf = LoadGeant4Workflow()
     for provider in itertools.chain(powder_providers, _dream_providers):
         wf.insert(provider)
     wf.insert(convert_dream_monitor_to_wavelength)
     wf.insert(resample_monitor_time_of_flight_data)
     insert_run_normalization(wf, run_norm)
+    if subtract_empty_can:
+        add_empty_can_subtraction(wf)
     for key, value in itertools.chain(
         default_parameters().items(), time_of_flight.default_parameters().items()
     ):

src/ess/powder/correction.py

@@ -12,6 +12,9 @@
 from ._util import event_or_outer_coord
 from .types import (
     AccumulatedProtonCharge,
+    BackgroundRun,
+    BackgroundSubtractedData,
+    BackgroundSubtractedDataTwoTheta,
     CaveMonitor,
     DataWithScatteringCoordinates,
     FocussedDataDspacing,
@@ -154,13 +157,49 @@ def _normalize_by_vanadium(
     return normed
 
 
+def normalize_by_vanadium_dspacing_removed_empty_can(
+    data: BackgroundSubtractedData[SampleRun],
+    vanadium: FocussedDataDspacing[VanadiumRun],
+    uncertainty_broadcast_mode: UncertaintyBroadcastMode,
+) -> IofDspacing:
+    """
+    Normalize sample data by a vanadium measurement and return intensity vs. d-spacing.
+
+    Parameters
+    ----------
+    data:
+        Sample data where events from an empty can / empty instrument measurement
+        have been subtracted.
+    vanadium:
+        Vanadium data.
+    uncertainty_broadcast_mode:
+        Choose how uncertainties of vanadium are broadcast to the sample data.
+        Defaults to ``UncertaintyBroadcastMode.fail``.
+
+    Returns
+    -------
+    :
+        ``data / vanadium``.
+        May contain a mask "zero_vanadium" which is ``True``
+        for bins where vanadium is zero.
+
+    See also
+    --------
+    normalize_by_vanadium_dspacing:
+        The same function but using data where no empty can has been subtracted.
+    """
+    return IofDspacing(
+        _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode)
+    )
+
+
 def normalize_by_vanadium_dspacing(
     data: FocussedDataDspacing[SampleRun],
     vanadium: FocussedDataDspacing[VanadiumRun],
     uncertainty_broadcast_mode: UncertaintyBroadcastMode,
 ) -> IofDspacing:
     """
-    Normalize sample data by a vanadium measurement and return intensity vs d-spacing.
+    Normalize sample data by a vanadium measurement and return intensity vs. d-spacing.
 
     Parameters
     ----------
@@ -178,19 +217,62 @@ def normalize_by_vanadium_dspacing(
         ``data / vanadium``.
         May contain a mask "zero_vanadium" which is ``True``
         for bins where vanadium is zero.
+
+    See also
+    --------
+    normalize_by_vanadium_dspacing_removed_empty_can:
+        The same function but using data where events from an empty can /
+        empty instrument measurement have been subtracted.
     """
     return IofDspacing(
         _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode)
     )
 
 
+def normalize_by_vanadium_dspacing_and_two_theta_removed_empty_can(
+    data: BackgroundSubtractedDataTwoTheta[SampleRun],
+    vanadium: FocussedDataDspacingTwoTheta[VanadiumRun],
+    uncertainty_broadcast_mode: UncertaintyBroadcastMode,
+) -> IofDspacingTwoTheta:
+    """
+    Normalize sample data by a vanadium measurement and return intensity vs.
+    (d-spacing, 2theta).
+
+    Parameters
+    ----------
+    data:
+        Sample data where events from an empty can / empty instrument measurement
+        have been subtracted.
+    vanadium:
+        Vanadium data.
+    uncertainty_broadcast_mode:
+        Choose how uncertainties of vanadium are broadcast to the sample data.
+        Defaults to ``UncertaintyBroadcastMode.fail``.
+
+    Returns
+    -------
+    :
+        ``data / vanadium``.
+        May contain a mask "zero_vanadium" which is ``True``
+        for bins where vanadium is zero.
+
+    See also
+    --------
+    normalize_by_vanadium_dspacing_and_two_theta:
+        The same function but using data where no empty can has been subtracted.
+    """
+    return IofDspacingTwoTheta(
+        _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode)
+    )
+
+
 def normalize_by_vanadium_dspacing_and_two_theta(
     data: FocussedDataDspacingTwoTheta[SampleRun],
     vanadium: FocussedDataDspacingTwoTheta[VanadiumRun],
     uncertainty_broadcast_mode: UncertaintyBroadcastMode,
 ) -> IofDspacingTwoTheta:
     """
-    Normalize sample data by a vanadium measurement and return intensity vs
+    Normalize sample data by a vanadium measurement and return intensity vs.
     (d-spacing, 2theta).
 
     Parameters
@@ -209,6 +291,12 @@ def normalize_by_vanadium_dspacing_and_two_theta(
         ``data / vanadium``.
         May contain a mask "zero_vanadium" which is ``True``
         for bins where vanadium is zero.
+
+    See also
+    --------
+    normalize_by_vanadium_dspacing_and_two_theta_removed_empty_can:
+        The same function but using data where events from an empty can /
+        empty instrument measurement have been subtracted.
     """
     return IofDspacingTwoTheta(
         _normalize_by_vanadium(data, vanadium, uncertainty_broadcast_mode)
@@ -335,6 +423,22 @@ def _shallow_copy(da: sc.DataArray) -> sc.DataArray:
     return out
 
 
+def subtract_background(
+    data: FocussedDataDspacing[SampleRun],
+    background: FocussedDataDspacing[BackgroundRun],
+) -> BackgroundSubtractedData[SampleRun]:
+    return BackgroundSubtractedData[SampleRun](data.bins.concatenate(-background))
+
+
+def subtract_background_two_theta(
+    data: FocussedDataDspacingTwoTheta[SampleRun],
+    background: FocussedDataDspacingTwoTheta[BackgroundRun],
+) -> BackgroundSubtractedDataTwoTheta[SampleRun]:
+    return BackgroundSubtractedDataTwoTheta[SampleRun](
+        data.bins.concatenate(-background)
+    )
+
+
 class RunNormalization(enum.Enum):
     """Type of normalization applied to each run."""
 
@@ -356,7 +460,15 @@ def insert_run_normalization(
             workflow.insert(normalize_by_proton_charge)
 
 
+def add_empty_can_subtraction(workflow: sciline.Pipeline) -> None:
+    """Insert providers to subtract empty can events from sample data."""
+    workflow.insert(normalize_by_vanadium_dspacing_removed_empty_can)
+    workflow.insert(normalize_by_vanadium_dspacing_and_two_theta_removed_empty_can)
+
+
 providers = (
+    subtract_background,
+    subtract_background_two_theta,
     normalize_by_proton_charge,
     normalize_by_vanadium_dspacing,
     normalize_by_vanadium_dspacing_and_two_theta,

src/ess/powder/types.py

@@ -58,7 +58,7 @@
 TimeOfFlightLookupTableFilename = tof_t.TimeOfFlightLookupTableFilename
 SimulationResults = tof_t.SimulationResults
 
-RunType = TypeVar("RunType", SampleRun, VanadiumRun)
+RunType = TypeVar("RunType", SampleRun, VanadiumRun, BackgroundRun)
 MonitorType = TypeVar("MonitorType", CaveMonitor, BunkerMonitor)
 
 
@@ -175,6 +175,16 @@ class NormalizedRunData(sciline.Scope[RunType, sc.DataArray], sc.DataArray):
     """Data that has been normalized by proton charge."""
 
 
+class BackgroundSubtractedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray):
+    """Data where background has been subtracted."""
+
+
+class BackgroundSubtractedDataTwoTheta(
+    sciline.Scope[RunType, sc.DataArray], sc.DataArray
+):
+    """Data with 2theta bins where background has been subtracted."""
+
+
 PixelMaskFilename = NewType("PixelMaskFilename", str)
 """Filename of a pixel mask."""
 

tests/dream/geant4_reduction_test.py

@@ -101,10 +101,16 @@ def params_for_det(request):
 
 
 @pytest.fixture
-def workflow(params_for_det):
+def workflow_no_empy_can(params_for_det):
     return make_workflow(params_for_det, run_norm=powder.RunNormalization.proton_charge)
 
 
+@pytest.fixture
+def workflow(workflow_no_empy_can):
+    powder.correction.add_empty_can_subtraction(workflow_no_empy_can)
+    return workflow_no_empy_can
+
+
 def make_workflow(params_for_det, *, run_norm):
     wf = dream.DreamGeant4Workflow(run_norm=run_norm)
     for key, value in params_for_det.items():
@@ -119,6 +125,15 @@ def test_pipeline_can_compute_dspacing_result(workflow):
     assert sc.identical(result.coords['dspacing'], params[DspacingBins])
 
 
+def test_pipeline_can_compute_dspacing_result_without_empty_can(workflow_no_empy_can):
+    workflow_no_empy_can[Filename[BackgroundRun]] = None
+    workflow_no_empy_can[MonitorFilename[BackgroundRun]] = None
+    workflow_no_empy_can = powder.with_pixel_mask_filenames(workflow_no_empy_can, [])
+    result = workflow_no_empy_can.compute(IofDspacing)
+    assert result.sizes == {'dspacing': len(params[DspacingBins]) - 1}
+    assert sc.identical(result.coords['dspacing'], params[DspacingBins])
+
+
 def test_pipeline_can_compute_dspacing_result_using_lookup_table_filename(workflow):
     workflow = powder.with_pixel_mask_filenames(workflow, [])
     workflow[TimeOfFlightLookupTableFilename] = dream.data.tof_lookup_table_high_flux()