diff --git a/docs/user-guide/dream/dream-powder-reduction.ipynb b/docs/user-guide/dream/dream-powder-reduction.ipynb index cea8b537..8bc3e98f 100644 --- a/docs/user-guide/dream/dream-powder-reduction.ipynb +++ b/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", + ")" ] }, { diff --git a/src/ess/dream/data.py b/src/ess/dream/data.py index 82c73ce0..0cee8a67 100644 --- a/src/ess/dream/data.py +++ b/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: diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index ece8c317..3316923f 100644 --- a/src/ess/dream/io/geant4.py +++ b/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) diff --git a/src/ess/dream/workflow.py b/src/ess/dream/workflow.py index 7525f1fb..6504e821 100644 --- a/src/ess/dream/workflow.py +++ b/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,9 +115,25 @@ 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): @@ -120,6 +141,8 @@ def DreamGeant4Workflow(*, run_norm: RunNormalization) -> sciline.Pipeline: 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() ): diff --git a/src/ess/powder/correction.py b/src/ess/powder/correction.py index f9ca00e1..3fb53e1b 100644 --- a/src/ess/powder/correction.py +++ b/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, diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index cf127809..ae277061 100644 --- a/src/ess/powder/types.py +++ b/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.""" diff --git a/tests/dream/geant4_reduction_test.py b/tests/dream/geant4_reduction_test.py index 51603cd2..1a252f61 100644 --- a/tests/dream/geant4_reduction_test.py +++ b/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()