[essreduce] Compute wavelength LUT from chopper cascade

packages/essdiffraction/docs/user-guide/dream/dream-make-wavelength-lookup-table.ipynb

@@ -63,7 +63,6 @@
     "wf = unwrap.LookupTableWorkflow()\n",
     "\n",
     "wf[unwrap.LtotalRange] = sc.scalar(5.0, unit=\"m\"), sc.scalar(80.0, unit=\"m\")\n",
-    "wf[unwrap.NumberOfSimulatedNeutrons] = 200_000  # Increase this number for more reliable results\n",
     "wf[unwrap.SourcePosition] = sc.vector([0, 0, 0], unit='m')\n",
     "wf[unwrap.DiskChoppers[AnyRun]] = disk_choppers\n",
     "wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit=\"m\")\n",

packages/essdiffraction/tests/dream/geant4_reduction_test.py

@@ -123,8 +123,6 @@ def dream_lookup_table():
         dream.beamline.InstrumentConfiguration.high_flux_BC215
     )
     lut_wf[unwrap.SourcePosition] = sc.vector(value=[0, 0, -76.55], unit="m")
-    lut_wf[unwrap.NumberOfSimulatedNeutrons] = 500_000
-    lut_wf[unwrap.SimulationSeed] = 555
     lut_wf[unwrap.PulseStride] = 1
     lut_wf[unwrap.LtotalRange] = (
         sc.scalar(60.0, unit="m"),

packages/essimaging/docs/odin/odin-make-wavelength-lookup-table.ipynb

@@ -42,8 +42,6 @@
     "wf = unwrap.LookupTableWorkflow()\n",
     "wf[unwrap.DiskChoppers[AnyRun]] = disk_choppers\n",
     "wf[unwrap.SourcePosition] = source_position\n",
-    "wf[unwrap.NumberOfSimulatedNeutrons] = 200_000  # Increase this number for more reliable results\n",
-    "wf[unwrap.SimulationSeed] = 1234\n",
     "wf[unwrap.PulseStride] = 2\n",
     "wf[unwrap.LtotalRange] = sc.scalar(5.0, unit=\"m\"), sc.scalar(65.0, unit=\"m\")\n",
     "wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit=\"m\")\n",

packages/essimaging/docs/tbl/tbl-make-wavelength-lookup-table.ipynb

@@ -40,8 +40,6 @@
     "wf = unwrap.LookupTableWorkflow()\n",
     "wf[unwrap.DiskChoppers[AnyRun]] = {}\n",
     "wf[unwrap.SourcePosition] = source_position\n",
-    "wf[unwrap.NumberOfSimulatedNeutrons] = 200_000  # Increase this number for more reliable results\n",
-    "wf[unwrap.SimulationSeed] = 1234\n",
     "wf[unwrap.PulseStride] = 1\n",
     "wf[unwrap.LtotalRange] = sc.scalar(5.0, unit=\"m\"), sc.scalar(35.0, unit=\"m\")\n",
     "wf[unwrap.DistanceResolution] = sc.scalar(0.1, unit=\"m\")\n",

packages/essimaging/tests/tbl/data_reduction_test.py

@@ -29,8 +29,6 @@ def wavelength_lookup_table() -> sl.Pipeline:
     lut_wf = unwrap.LookupTableWorkflow()
     lut_wf[unwrap.DiskChoppers[AnyRun]] = {}
     lut_wf[unwrap.SourcePosition] = sc.vector([0, 0, 0], unit="m")
-    lut_wf[unwrap.NumberOfSimulatedNeutrons] = 200_000
-    lut_wf[unwrap.SimulationSeed] = 333
     lut_wf[unwrap.PulseStride] = 1
     lut_wf[unwrap.LtotalRange] = (sc.scalar(25.0, unit="m"), sc.scalar(35.0, unit="m"))
     return lut_wf.compute(LookupTable)

packages/essnmx/src/ess/nmx/workflows.py

@@ -5,10 +5,11 @@
 import sciline
 import scipp as sc
 import scippnexus as snx
-import tof
 from scippneutron.conversion.tof import tof_from_wavelength
 
 from ess.reduce.nexus.types import (
+    AnyRun,
+    DiskChoppers,
     EmptyDetector,
     Filename,
     NeXusComponent,
@@ -18,15 +19,13 @@
     SampleRun,
 )
 from ess.reduce.unwrap import (
-    BeamlineComponentReading,
     GenericUnwrapWorkflow,
     LookupTableFilename,
     LookupTableRelativeErrorThreshold,
     LookupTableWorkflow,
     LtotalRange,
-    NumberOfSimulatedNeutrons,
-    SimulationResults,
-    SimulationSeed,
+    SourcePosition,
+    SourcePulse,
     WavelengthDetector,
 )
 from ess.reduce.workflow import register_workflow
@@ -45,60 +44,11 @@
     TofSimulationMaxWavelength: sc.scalar(3.6, unit='angstrom'),
     TofSimulationMinWavelength: sc.scalar(1.8, unit='angstrom'),
     LookupTableRelativeErrorThreshold: {f'detector_panel_{i}': 0.1 for i in range(5)},
+    # TODO: This should become DiskChoppers[RunType] once we add choppers
+    DiskChoppers[AnyRun]: {},
 }
 
 
-def _simulate_fixed_wavelength_tof(
-    wmin: TofSimulationMinWavelength,
-    wmax: TofSimulationMaxWavelength,
-    neutrons: NumberOfSimulatedNeutrons,
-    seed: SimulationSeed,
-) -> SimulationResults:
-    """
-    Simulate a pulse of neutrons propagating through the instrument using the
-    ``tof`` package (https://scipp.github.io/tof/).
-    This runs a simulation assuming there are no choppers in the instrument.
-
-    Parameters
-    ----------
-    wmin:
-        Minimum wavelength of the simulated neutrons.
-    wmax:
-        Maximum wavelength of the simulated neutrons.
-    neutrons:
-        Number of neutrons to simulate.
-    seed:
-        Random seed for the simulation.
-    """
-    source = tof.Source(
-        facility="ess",
-        neutrons=neutrons,
-        pulses=1,
-        seed=seed,
-        wmax=wmax,
-        wmin=wmin,
-    )
-    events = source.data.squeeze().flatten(to="event")
-
-    return SimulationResults(
-        readings={
-            "source": BeamlineComponentReading(
-                time_of_arrival=events.coords["birth_time"],
-                wavelength=events.coords["wavelength"],
-                weight=events.data,
-                distance=source.distance,
-            )
-        },
-        choppers=None,
-    )
-
-
-def _merge_panels(*da: sc.DataArray) -> sc.DataArray:
-    """Merge multiple DataArrays representing different panels into one."""
-    merged = sc.concat(da, dim='panel')
-    return merged
-
-
 def select_detector_names(*, detector_ids: Iterable[int] = (0, 1, 2)):
     import os
 
@@ -259,15 +209,32 @@ def compute_detector_tof(da: WavelengthDetector[RunType]) -> TofDetector[RunType
     )
 
 
+def _source_position_to_SourcePosition(
+    source_position: Position[snx.NXsource, SampleRun],
+) -> SourcePosition:
+    """
+    This is a temporary provider to convert the source position from the Nexus file to
+    the SourcePosition type used in the unwrapping workflow.
+    In the next iteration, we will directly use the source position from the Nexus file
+    in the unwrapping workflow and remove this provider.
+    """
+    return SourcePosition(source_position)
+
+
 @register_workflow
 def NMXWorkflow() -> sciline.Pipeline:
     generic_wf = GenericUnwrapWorkflow(run_types=[SampleRun], monitor_types=[])
 
-    generic_wf.insert(_retrieve_crystal_rotation)
-    generic_wf.insert(assemble_sample_metadata)
-    generic_wf.insert(assemble_source_metadata)
-    generic_wf.insert(assemble_detector_metadata)
-    generic_wf.insert(compute_detector_tof)
+    for provider in (
+        _retrieve_crystal_rotation,
+        assemble_sample_metadata,
+        assemble_source_metadata,
+        assemble_detector_metadata,
+        compute_detector_tof,
+        _source_position_to_SourcePosition,
+    ):
+        generic_wf.insert(provider)
+
     for key, value in default_parameters.items():
         generic_wf[key] = value
 
@@ -317,12 +284,12 @@ def initialize_nmx_workflow(*, config: WorkflowConfig) -> sciline.Pipeline:
         wf[LookupTableFilename] = config.lookup_table_file_path
     else:
         wf = _merge_workflows(base_wf=wf, merged_wf=LookupTableWorkflow())
-        wf.insert(_simulate_fixed_wavelength_tof)
         wmax = sc.scalar(config.tof_simulation_max_wavelength, unit='angstrom')
         wmin = sc.scalar(config.tof_simulation_min_wavelength, unit='angstrom')
-        wf[TofSimulationMaxWavelength] = wmax
-        wf[TofSimulationMinWavelength] = wmin
-        wf[SimulationSeed] = config.tof_simulation_seed
+        wf[SourcePulse] = SourcePulse(
+            time=(sc.scalar(0.0, unit='ms'), sc.scalar(5.0, unit='ms')),
+            wavelength=(wmin, wmax),
+        )
         ltotal_min = sc.scalar(value=config.tof_simulation_min_ltotal, unit='m')
         ltotal_max = sc.scalar(value=config.tof_simulation_max_ltotal, unit='m')
         wf[LtotalRange] = LtotalRange((ltotal_min, ltotal_max))

packages/essnmx/tests/executable_test.py

@@ -25,6 +25,7 @@
 from ess.nmx.configurations import TimeBinCoordinate, TimeBinUnit, to_command_arguments
 from ess.nmx.executables import reduction
 from ess.nmx.types import Compression, NMXLauetof
+from ess.reduce.unwrap.types import SourcePosition
 
 
 def _build_arg_list_from_pydantic_instance(*instances: pydantic.BaseModel) -> list[str]:
@@ -410,6 +411,7 @@ def lut_file_path(tmp_path: pathlib.Path):
 
     # Simply use the default workflow for testing.
     workflow = initialize_nmx_workflow(config=WorkflowConfig())
+    workflow[SourcePosition] = sc.vector([0, 0, 0], unit='m')
     lut: LookupTable = workflow.compute(LookupTable)
 
     # Change the tof range a bit for testing.