diff --git a/packages/essimaging/src/ess/imaging/normalization.py b/packages/essimaging/src/ess/imaging/normalization.py index 7dc953e5c..a730ed1fa 100644 --- a/packages/essimaging/src/ess/imaging/normalization.py +++ b/packages/essimaging/src/ess/imaging/normalization.py @@ -10,6 +10,7 @@ BackgroundSubtractedDetector, DarkBackgroundRun, FluxNormalizedDetector, + MeanDarkFrame, NormalizedImage, OpenBeamRun, SampleRun, @@ -17,9 +18,23 @@ ) +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, diff --git a/packages/essimaging/src/ess/imaging/types.py b/packages/essimaging/src/ess/imaging/types.py index e17217a3b..f101648f0 100644 --- a/packages/essimaging/src/ess/imaging/types.py +++ b/packages/essimaging/src/ess/imaging/types.py @@ -87,13 +87,13 @@ 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.""" -class ExposureTime(sciline.Scope[RunType, sc.DataArray], sc.DataArray): - """Exposure time of each frame recorded by the camera detector.""" - - del sc, sciline, NewType diff --git a/packages/essimaging/src/ess/tbl/data.py b/packages/essimaging/src/ess/tbl/data.py index aa5278a97..3bec21786 100644 --- a/packages/essimaging/src/ess/tbl/data.py +++ b/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", }, ) diff --git a/packages/essimaging/src/ess/tbl/orca.py b/packages/essimaging/src/ess/tbl/orca.py index fb74592b0..41eec9975 100644 --- a/packages/essimaging/src/ess/tbl/orca.py +++ b/packages/essimaging/src/ess/tbl/orca.py @@ -7,19 +7,13 @@ 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 from .. import imaging from ..imaging.types import ( CorrectedDetector, DarkBackgroundRun, - ExposureTime, FluxNormalizedDetector, OpenBeamRun, ProtonCharge, @@ -28,58 +22,14 @@ ) -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 ---------- @@ -87,58 +37,30 @@ def normalize_by_proton_charge_orca( 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: """Return the default NeXus names and detector name for the ORCA workflow.""" return { NeXusDetectorName: 'orca_detector', - NeXusName[ExposureTime]: '/entry/instrument/orca_detector/camera_exposure', } @@ -156,7 +78,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(): diff --git a/packages/essimaging/tests/tbl/orca_test.py b/packages/essimaging/tests/tbl/orca_test.py index 86069202e..f9bd8a033 100644 --- a/packages/essimaging/tests/tbl/orca_test.py +++ b/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 diff --git a/packages/essimaging/tools/make-tbl-images-from-ymir.ipynb b/packages/essimaging/tools/make-tbl-images-from-ymir.ipynb index d30aeb4d3..487defbaf 100644 --- a/packages/essimaging/tools/make-tbl-images-from-ymir.ipynb +++ b/packages/essimaging/tools/make-tbl-images-from-ymir.ipynb @@ -96,13 +96,15 @@ "exposure_time = int(5e8) # 0.5s (made up)\n", "period = 1e9/14\n", "\n", - "dark_charge_time = np.arange(float(dark_time[0]), float(dark_time[-1] + exposure_time), period).astype('uint64')\n", + "# Need to start a bit before the start of the time axis of the data\n", + "# so that sc.lookup in 'previous' mode can find the pulse time\n", + "dark_charge_time = np.arange(float(dark_time[0]) - exposure_time, float(dark_time[-1] + exposure_time), period).astype('uint64')\n", "dark_charge_value = np.random.uniform(0.99, 1.01, size=len(dark_charge_time)).astype('float64')\n", "\n", - "ob_charge_time = np.arange(float(ob_time[0]), float(ob_time[-1] + exposure_time), period).astype('uint64')\n", + "ob_charge_time = np.arange(float(ob_time[0]) - exposure_time, float(ob_time[-1] + exposure_time), period).astype('uint64')\n", "ob_charge_value = np.random.uniform(0.99, 1.01, size=len(ob_charge_time)).astype('float64')\n", "\n", - "sample_charge_time = np.arange(float(sample_time[0]), float(sample_time[-1] + exposure_time), period).astype('uint64')\n", + "sample_charge_time = np.arange(float(sample_time[0]) - exposure_time, float(sample_time[-1] + exposure_time), period).astype('uint64')\n", "sample_charge_value = np.random.uniform(0.99, 1.01, size=len(sample_charge_time)).astype('float64')" ] }, @@ -149,8 +151,10 @@ " import shutil\n", " shutil.copyfile(template_nexus_file, outfile)\n", "\n", - " # detector_path = 'entry/instrument/histogram_mode_detectors/orca' # ODIN\n", - " detector_path = 'entry/instrument/orca_detector/' # TBL\n", + " instr_path = 'entry/instrument/'\n", + "\n", + " # detector_path = instr_path + 'histogram_mode_detectors/orca' # ODIN\n", + " detector_path = instr_path + 'orca_detector/' # TBL\n", " proton_charge_path = 'entry/neutron_prod_info/pulse_charge'\n", "\n", " with h5.File(outfile, \"r+\") as f:\n", @@ -158,33 +162,36 @@ " detector_data = f[detector_path+\"data\"]\n", " replace_dataset(detector_data, name=\"value\", values=data[\"value\"])\n", " replace_dataset(detector_data, name=\"time\", values=data[\"time\"])\n", + " detector_data.attrs['axes'] = ['time', 'y_pixel_offset', 'x_pixel_offset']\n", "\n", " detector = f[detector_path]\n", - " detector.copy('intensifier', detector, 'camera_exposure')\n", - " exp = f[detector_path+\"camera_exposure\"]\n", - " replace_dataset(exp, name=\"value\", values=np.array([exposure_time]))\n", - " replace_dataset(exp, name=\"average_value\", values=exposure_time)\n", - " replace_dataset(exp, name=\"minimum_value\", values=exposure_time)\n", - " replace_dataset(exp, name=\"maximum_value\", values=exposure_time)\n", - " f[detector_path+\"camera_exposure/value\"].attrs['units'] = 'ns'\n", - " f[detector_path+\"camera_exposure/average_value\"].attrs['units'] = 'ns'\n", - " f[detector_path+\"camera_exposure/minimum_value\"].attrs['units'] = 'ns'\n", - " f[detector_path+\"camera_exposure/maximum_value\"].attrs['units'] = 'ns'\n", - "\n", + " del detector['x_pixel_offset']\n", " xoff = detector.create_dataset('x_pixel_offset', data=x_pixel_offset.values)\n", " xoff.attrs['units'] = str(x_pixel_offset.unit)\n", " xoff.attrs['axis'] = 1\n", + " del detector['y_pixel_offset']\n", " yoff = detector.create_dataset('y_pixel_offset', data=y_pixel_offset.values)\n", " yoff.attrs['units'] = str(y_pixel_offset.unit)\n", " yoff.attrs['axis'] = 2\n", - " detector.attrs['axes'] = ['dim_2', 'dim_1']\n", "\n", " pcharge = f[proton_charge_path]\n", " replace_dataset(pcharge, name=\"value\", values=proton_charge[\"value\"])\n", " replace_dataset(pcharge, name=\"time\", values=proton_charge[\"time\"])\n", " replace_dataset(pcharge, name=\"average_value\", values=proton_charge[\"value\"].mean())\n", " replace_dataset(pcharge, name=\"minimum_value\", values=proton_charge[\"value\"].min())\n", - " replace_dataset(pcharge, name=\"maximum_value\", values=proton_charge[\"value\"].max())" + " replace_dataset(pcharge, name=\"maximum_value\", values=proton_charge[\"value\"].max())\n", + "\n", + " for det in ('he3_detector_bank0', 'he3_detector_bank1', 'ngem_detector', 'multiblade_detector', 'timepix3_detector'):\n", + " del f[instr_path + det]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "template_nexus_file = \"coda_tbl_999999_00010814.hdf\"" ] }, { @@ -193,8 +200,6 @@ "metadata": {}, "outputs": [], "source": [ - "template_nexus_file = \"857127_00000212.hdf\"\n", - "\n", "# Dark run\n", "make_new_file(\n", " template_nexus_file=template_nexus_file,\n", @@ -251,7 +256,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.7" + "version": "3.12.12" } }, "nbformat": 4,