[essimaging] Orca normalization by proton charge per frame

packages/essimaging/src/ess/imaging/normalization.py

@@ -10,16 +10,31 @@
     BackgroundSubtractedDetector,
     DarkBackgroundRun,
     FluxNormalizedDetector,
+    MeanDarkFrame,
     NormalizedImage,
     OpenBeamRun,
     SampleRun,
     UncertaintyBroadcastMode,
 )
 
 
+def average_dark_frames(
+    dark_frames: FluxNormalizedDetector[DarkBackgroundRun],
+) -> MeanDarkFrame:
+    """
+    Average the dark frames (background runs) over time to obtain a mean dark frame.
+
+    Parameters
+    ----------
+    dark_frames:
+        Dark frames to average.
+    """
+    return MeanDarkFrame(dark_frames.mean('time'))
+
+
 def subtract_background_sample(
     data: FluxNormalizedDetector[SampleRun],
-    background: FluxNormalizedDetector[DarkBackgroundRun],
+    background: MeanDarkFrame,
     uncertainties: UncertaintyBroadcastMode,
 ) -> BackgroundSubtractedDetector[SampleRun]:
     """
@@ -42,7 +57,8 @@ def subtract_background_sample(
 
 def subtract_background_openbeam(
     data: FluxNormalizedDetector[OpenBeamRun],
-    background: FluxNormalizedDetector[DarkBackgroundRun],
+    background: MeanDarkFrame,
+    uncertainties: UncertaintyBroadcastMode,
 ) -> BackgroundSubtractedDetector[OpenBeamRun]:
     """
     Subtract background from proton-charge normalized open beam data.
@@ -53,8 +69,13 @@ def subtract_background_openbeam(
         Open beam data to process (normalized to proton charge).
     background:
         Background (dark frame) data to subtract (normalized to proton charge).
+    uncertainties:
+        Mode to use when broadcasting uncertainties from background to data (data
+        typically has multiple frames, while background usually has only one frame).
     """
-    return BackgroundSubtractedDetector[OpenBeamRun](data - background)
+    return BackgroundSubtractedDetector[OpenBeamRun](
+        data - broadcast_uncertainties(background, prototype=data, mode=uncertainties)
+    )
 
 
 def sample_over_openbeam(
@@ -78,11 +99,14 @@ def sample_over_openbeam(
     """
     return NormalizedImage(
         sample
-        / broadcast_uncertainties(open_beam, prototype=sample, mode=uncertainties)
+        / broadcast_uncertainties(
+            open_beam.mean('time'), prototype=sample, mode=uncertainties
+        )
     )
 
 
 providers = (
+    average_dark_frames,
     subtract_background_sample,
     subtract_background_openbeam,
     sample_over_openbeam,

packages/essimaging/src/ess/imaging/types.py

@@ -87,6 +87,10 @@ class BackgroundSubtractedDetector(sciline.Scope[RunType, sc.DataArray], sc.Data
     """Detector counts with dark background subtracted."""
 
 
+MeanDarkFrame = NewType("MeanDarkFrame", sc.DataArray)
+"""Mean dark frame: average of background (dark) frames."""
+
+
 NormalizedImage = NewType("NormalizedImage", sc.DataArray)
 """Final image: background-subtracted sample run divided by background-subtracted open
 beam run."""

packages/essimaging/src/ess/tbl/data.py

@@ -8,16 +8,16 @@
 
 _registry = make_registry(
     'ess/tbl',
-    version="2",
+    version="3",
     files={
         "tbl_sample_data_2025-03.hdf": "md5:12db6bc06721278b3abe47992eac3e77",
         "TBL-wavelength-lookup-table-no-choppers-5m-35m.h5": "md5:e28793b7e1c12986ee63a1df68723268",  # noqa: E501
         'tbl-orca-focussing.hdf.zip': Entry(
             alg='md5', chk='f365acd9ea45dd205c0b9398d163cfa4', unzip=True
         ),
-        "ymir_lego_dark_run.hdf": "md5:c0ed089dd7663986042e29fb47514130",
-        "ymir_lego_openbeam_run.hdf": "md5:00375becd54d2ed3be096413dc30883c",
-        "ymir_lego_sample_run.hdf": "md5:ae56a335cf3d4e87ef090ec4e51da69c",
+        "ymir_lego_dark_run.hdf": "md5:2d9ae7d6f1d0502c17b0030bcb3dad1d",
+        "ymir_lego_openbeam_run.hdf": "md5:f55e6064a2602de404c535f72beb25a3",
+        "ymir_lego_sample_run.hdf": "md5:d3eb65d9b5f8b90777059ec2638297ab",
     },
 )
 

packages/essimaging/src/ess/tbl/orca.py

@@ -7,13 +7,8 @@
 import sciline as sl
 import scipp as sc
 
-from ess.reduce.nexus import GenericNeXusWorkflow, load_from_path
-from ess.reduce.nexus.types import (
-    NeXusDetectorName,
-    NeXusFileSpec,
-    NeXusLocationSpec,
-    NeXusName,
-)
+from ess.reduce.nexus import GenericNeXusWorkflow
+from ess.reduce.nexus.types import NeXusDetectorName, NeXusName
 
 from .. import imaging
 from ..imaging.types import (
@@ -28,110 +23,39 @@
 )
 
 
-def load_exposure_time(
-    file: NeXusFileSpec[RunType], path: NeXusName[ExposureTime[RunType]]
-) -> ExposureTime[RunType]:
-    # Note that putting '/value' at the end of the 'path' in the default_parameters
-    # yields different results as it can return a Variable instead of a DataArray,
-    # depending on the contents of the NeXus file.
-    return ExposureTime[RunType](
-        load_from_path(NeXusLocationSpec(filename=file.value, component_name=path))[
-            "value"
-        ].squeeze()
-    )
-
-
-def _compute_proton_charge_per_exposure(
-    data: sc.DataArray, proton_charge: sc.DataArray, exposure_time: sc.DataArray
-) -> sc.DataArray:
-    # A note on timings:
-    # We want to sum the proton charge inside each time bin (defined by the duration of
-    # each frame). However, the time dimension of the data recorded at the detector is
-    # not the same time as the proton charge (which is when the protons hit the
-    # target). We need to shift the proton charge time to account for the time it takes
-    # for neutrons to travel from the target to the detector. Does this mean we cannot
-    # do the normalization without computing time of flight?
-
-    t = data.coords['time']
-    exp = exposure_time.data.to(unit=t.unit)
-    # The following assumes that the different between successive frames (time stamps)
-    # is larger than the exposure time. We need to check that this is indeed the case.
-    if (t[1:] - t[:-1]).min() < exp:
-        raise ValueError(
-            "normalize_by_proton_charge_orca: Exposure time is larger than the "
-            "smallest time between successive frames."
-        )
-
-    bins = sc.sort(sc.concat([t, t + exp], dim=t.dim), t.dim)
-    # We select every second bin, as the odd bins lie between the end of the exposure
-    # of one frame and the start of the next frame.
-    return proton_charge.hist(time=bins).data[::2]
-
-
 def normalize_by_proton_charge_orca(
-    data: CorrectedDetector[RunType],
-    proton_charge: ProtonCharge[RunType],
-    exposure_time: ExposureTime[RunType],
+    data: CorrectedDetector[RunType], proton_charge: ProtonCharge[RunType]
 ) -> FluxNormalizedDetector[RunType]:
     """
-    Normalize detector data by the proton charge (dark and open beam runs).
-    We find the time stamps for the data, which mark the start of an exposure.
-    We sum the proton charge within the intervals timestamp to timestamp +
-    exposure_time.
-    We then sum the data along the time dimension, and divide by the sum of the proton
-    charge accumulated during each exposure.
+    Normalize detector data by the proton charge.
+    We find the time stamps for the data, which mark the start of an exposure,
+    find the corresponding proton charge for each time stamp, and divide the data by
+    the proton charge.
 
     Parameters
     ----------
     data:
         Corrected detector data to be normalized.
     proton_charge:
         Proton charge data for normalization.
-    exposure_time:
-        Exposure time for each image in the data.
     """
 
-    charge_per_frame = _compute_proton_charge_per_exposure(
-        data, proton_charge, exposure_time
-    )
-
-    return FluxNormalizedDetector[RunType](
-        data.sum('time') / charge_per_frame.sum('time')
-    )
-
+    # A note on timings:
+    # We want to sum the proton charge inside each time bin (defined by the duration of
+    # each frame). However, the time dimension of the data recorded at the detector is
+    # not the same time as the proton charge (which is when the protons hit the
+    # target). We need to shift the proton charge time to account for the time it takes
+    # for neutrons to travel from the target to the detector. Does this mean we cannot
+    # do the normalization without computing time of flight?
 
-def normalize_by_proton_charge_orca_sample(
-    data: CorrectedDetector[SampleRun],
-    proton_charge: ProtonCharge[SampleRun],
-    exposure_time: ExposureTime[SampleRun],
-) -> FluxNormalizedDetector[SampleRun]:
-    """
-    Normalize sample run detector data by the proton charge.
-    The handling of the SampleRun is different from the other runs:
-    The sample may have a time dimension representing multiple frames.
-    In this case, we need to sum the proton charge for each frame separately.
+    proton_charge_lookup = sc.lookup(proton_charge, 'time', mode='previous')
 
-    Parameters
-    ----------
-    data:
-        Corrected detector data to be normalized.
-    proton_charge:
-        Proton charge data for normalization.
-    """
-
-    charge_per_frame = _compute_proton_charge_per_exposure(
-        data, proton_charge, exposure_time
+    return FluxNormalizedDetector[RunType](
+        data / proton_charge_lookup[data.coords['time']]
     )
 
-    # Here we preserve the time dimension of the sample data and the proton charge.
-    return FluxNormalizedDetector[SampleRun](data / charge_per_frame)
 
-
-providers = (
-    load_exposure_time,
-    normalize_by_proton_charge_orca,
-    normalize_by_proton_charge_orca_sample,
-)
+orca_providers = (normalize_by_proton_charge_orca,)
 
 
 def default_parameters() -> dict:
@@ -156,7 +80,7 @@ def OrcaNormalizedImagesWorkflow(**kwargs) -> sl.Pipeline:
     for provider in (
         *imaging.normalization.providers,
         *imaging.masking.providers,
-        *providers,
+        *orca_providers,
     ):
         wf.insert(provider)
     for key, param in default_parameters().items():

packages/essimaging/tests/tbl/orca_test.py

@@ -17,6 +17,7 @@
     Filename,
     FluxNormalizedDetector,
     MaskingRules,
+    MeanDarkFrame,
     NeXusDetectorName,
     NormalizedImage,
     OpenBeamRun,
@@ -73,28 +74,26 @@ def test_workflow_applies_masks(workflow, run):
     assert da.sum().value > masked_da.sum().value
 
 
-@pytest.mark.parametrize("run", [OpenBeamRun, DarkBackgroundRun])
+@pytest.mark.parametrize("run", [OpenBeamRun, DarkBackgroundRun, SampleRun])
 def test_workflow_normalizes_by_proton_charge(workflow, run):
     da = workflow.compute(FluxNormalizedDetector[run])
-    # Dark and open beam runs have been averaged over time dimension
-    assert da.ndim == 2
-    assert "time" not in da.dims
-    # TODO: should it be "counts / uC"?
-    assert da.unit == "1 / uC"
-
-
-def test_workflow_normalizes_sample_by_proton_charge(workflow):
-    da = workflow.compute(FluxNormalizedDetector[SampleRun])
-    # Sample run retains time dimension
     assert da.ndim == 3
     assert "time" in da.dims
     # TODO: should it be "counts / uC"?
     assert da.unit == "1 / uC"
 
 
+def test_workflow_computes_mean_dark_frame(workflow):
+    dark_frames = workflow.compute(FluxNormalizedDetector[DarkBackgroundRun])
+    mean_dark_frame = workflow.compute(MeanDarkFrame)
+    assert mean_dark_frame.ndim == 2
+    assert "time" not in mean_dark_frame.dims
+    assert mean_dark_frame.unit == dark_frames.unit
+
+
 @pytest.mark.parametrize("run", [OpenBeamRun, SampleRun])
 def test_workflow_subtracts_dark_background(workflow, run):
-    background = workflow.compute(FluxNormalizedDetector[DarkBackgroundRun])
+    background = workflow.compute(MeanDarkFrame)
     before = workflow.compute(FluxNormalizedDetector[run])
     after = workflow.compute(BackgroundSubtractedDetector[run])
 
@@ -106,5 +105,7 @@ def test_workflow_computes_normalized_image(workflow):
     open_beam = workflow.compute(BackgroundSubtractedDetector[OpenBeamRun])
     normalized = workflow.compute(NormalizedImage)
 
-    assert_identical(sc.values(normalized), sc.values(sample) / sc.values(open_beam))
+    assert_identical(
+        sc.values(normalized), sc.values(sample) / sc.values(open_beam.mean('time'))
+    )
     assert normalized.unit == sc.units.dimensionless