From 97e43b2dda3fe76f2d27250f3467aa1650214c51 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Wed, 4 Jun 2025 15:34:52 +0200 Subject: [PATCH 01/10] Begin adding Zoom sample workflow --- docs/user-guide/zoom.ipynb | 157 ++++++++++++++++++++++++----------- src/ess/polarization/zoom.py | 75 +++++++++-------- 2 files changed, 149 insertions(+), 83 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index a45232c..9f93032 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -23,6 +23,7 @@ "from ess import polarization as pol\n", "from ess import sans\n", "from ess import isissans as isis\n", + "import ess.isissans.data # noqa: F401\n", "from ess.sans.types import *" ] }, @@ -33,8 +34,14 @@ "metadata": {}, "outputs": [], "source": [ - "sans_workflow = isis.zoom.ZoomWorkflow()\n", - "sans_workflow = sans.with_pixel_mask_filenames(sans_workflow, []) # no masks" + "from pathlib import Path\n", + "\n", + "data_folder = Path('zoom_polarized_data')\n", + "# Runs with analyzer at 4 different times\n", + "cell_runs = [str(data_folder / f'ZOOM00022{run}.nxs') for run in [710, 712, 714, 716]]\n", + "empty_run = data_folder / 'ZOOM00034787.nxs'\n", + "depolarized_run = data_folder / 'ZOOM00022718.nxs'\n", + "cell_runs = [*cell_runs, depolarized_run]" ] }, { @@ -44,14 +51,9 @@ "metadata": {}, "outputs": [], "source": [ - "from pathlib import Path\n", - "\n", - "data_folder = Path('zoom_polarized_data')\n", - "# Runs with analyzer at 4 different times\n", - "cell_runs = [str(data_folder / f'ZOOM00022{run}.nxs') for run in [710, 712, 714, 716]]\n", - "empty_run = data_folder / 'ZOOM00034787.nxs'\n", - "depolarized_run = data_folder / 'ZOOM00022718.nxs'\n", - "cell_runs = [*cell_runs, depolarized_run]" + "data_folder = Path('/home/simon/instruments/polarization/2024-08')\n", + "sample_runs = [str(data_folder / f'ZOOM00038{run}.nxs') for run in [513]]\n", + "empty_run = data_folder / 'ZOOM00038336.nxs'" ] }, { @@ -69,20 +71,79 @@ "metadata": {}, "outputs": [], "source": [ - "from ess.polarization.zoom import ZoomTransmissionFractionWorkflow\n", + "params = {\n", + " Filename[EmptyBeamRun]: str(empty_run),\n", + " WavelengthBins: sc.geomspace(\n", + " 'wavelength', start=1.75, stop=16.5, num=141, unit='Å'\n", + " ),\n", + " isis.mantidio.Period: 0,\n", + " isis.mantidio.CalibrationWorkspace: None,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "from ess.polarization.zoom import (\n", + " ZoomTransmissionFractionWorkflow,\n", + " ZoomPolarizedWorkflow,\n", + ")\n", "\n", "sans_workflow = ZoomTransmissionFractionWorkflow(cell_runs)\n", - "sans_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", - "sans_workflow[WavelengthBins] = sc.geomspace(\n", - " 'wavelength', start=1.75, stop=16.5, num=141, unit='Å'\n", + "for key, value in params.items():\n", + " sans_workflow[key] = value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "sample_workflow = ZoomPolarizedWorkflow()\n", + "for key, value in params.items():\n", + " sample_workflow[key] = value\n", + "masks = isis.data.zoom_tutorial_mask_filenames()\n", + "sample_workflow = sans.with_pixel_mask_filenames(sample_workflow, masks)\n", + "sample_workflow.insert(isis.io.transmission_from_background_run)\n", + "sample_workflow.insert(isis.io.transmission_from_sample_run)\n", + "sample_workflow[QBins] = sc.geomspace(\n", + " dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom'\n", + ")\n", + "sample_workflow[NonBackgroundWavelengthRange] = sc.array(\n", + " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", ")\n", - "sans_workflow[isis.mantidio.Period] = 0\n", - "sans_workflow[isis.mantidio.CalibrationWorkspace] = None" + "sample_workflow[CorrectForGravity] = False\n", + "sample_workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop\n", + "sample_workflow[ReturnEvents] = False\n", + "\n", + "# TODO\n", + "# sample_workflow[DirectBeamFilename]\n", + "sample_workflow[DirectBeam] = None\n", + "\n", + "sample_workflow[Filename[SampleRun]] = sample_runs[0]\n", + "# sample_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", + "# sample_workflow[isis.SampleOffset] = sc.vector([0.0, 0.0, 0.11], unit='m')\n", + "# sample_workflow[isis.DetectorBankOffset] = sc.vector([0.0, 0.0, 0.5], unit='m')\n", + "sample_workflow[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m')\n", + "sample_workflow.visualize(IofQ[SampleRun], graph_attr={'rankdir': 'LR'})\n", + "# sample_workflow.compute(NeXusComponent[Incident, SampleRun])\n", + "# sample_workflow.compute(NeXusComponent[Transmission, SampleRun])\n", + "# sample_workflow.compute(DetectorData[SampleRun])\n", + "# sample_workflow.compute(IofQ[SampleRun])\n", + "# ws = sample_workflow.compute(isis.mantidio.DataWorkspace[SampleRun])\n", + "dg = sample_workflow.compute(isis.mantidio.LoadedFileContents[SampleRun])\n", + "dg" ] }, { "cell_type": "markdown", - "id": "6", + "id": "8", "metadata": {}, "source": [ "## Inspect data for one of the runs with analyzer cell" @@ -91,7 +152,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7", + "id": "9", "metadata": {}, "outputs": [], "source": [ @@ -106,7 +167,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -115,7 +176,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "We can load the combined time-dependent incident and transmission monitors.\n", @@ -125,7 +186,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -145,7 +206,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "The task graph for computing the transmission fraction is:" @@ -154,7 +215,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -163,7 +224,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "## Compute transmission fractions\n", @@ -175,7 +236,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -184,7 +245,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "We can plot the computed transmission fractions:" @@ -193,7 +254,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -209,7 +270,7 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "19", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +291,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -244,7 +305,7 @@ }, { "cell_type": "markdown", - "id": "19", + "id": "21", "metadata": {}, "source": [ "We can now setup the polarization analysis workflow.\n", @@ -254,7 +315,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -284,7 +345,7 @@ }, { "cell_type": "markdown", - "id": "21", + "id": "23", "metadata": {}, "source": [ "The workflow can compute the transmission function:" @@ -293,7 +354,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -302,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "23", + "id": "25", "metadata": {}, "source": [ "We can evaluate this transmission function at desired time and wavelength points:" @@ -311,7 +372,7 @@ { "cell_type": "code", "execution_count": null, - "id": "24", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -326,7 +387,7 @@ { "cell_type": "code", "execution_count": null, - "id": "25", + "id": "27", "metadata": {}, "outputs": [], "source": [ @@ -339,7 +400,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -352,7 +413,7 @@ { "cell_type": "code", "execution_count": null, - "id": "27", + "id": "29", "metadata": {}, "outputs": [], "source": [ @@ -362,7 +423,7 @@ { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -372,7 +433,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -381,7 +442,7 @@ }, { "cell_type": "markdown", - "id": "30", + "id": "32", "metadata": {}, "source": [ "## Correction workflow\n", @@ -394,7 +455,7 @@ { "cell_type": "code", "execution_count": null, - "id": "31", + "id": "33", "metadata": {}, "outputs": [], "source": [ @@ -404,7 +465,7 @@ }, { "cell_type": "markdown", - "id": "32", + "id": "34", "metadata": {}, "source": [ "We will use a second-order polynomial supermirror efficiency function:" @@ -413,7 +474,7 @@ { "cell_type": "code", "execution_count": null, - "id": "33", + "id": "35", "metadata": {}, "outputs": [], "source": [ @@ -430,7 +491,7 @@ { "cell_type": "code", "execution_count": null, - "id": "34", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -442,7 +503,7 @@ }, { "cell_type": "markdown", - "id": "35", + "id": "37", "metadata": {}, "source": [ "For a single channel, the complete workflow looks as follows:" @@ -451,7 +512,7 @@ { "cell_type": "code", "execution_count": null, - "id": "36", + "id": "38", "metadata": {}, "outputs": [], "source": [ @@ -463,7 +524,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "mantid310", "language": "python", "name": "python3" }, diff --git a/src/ess/polarization/zoom.py b/src/ess/polarization/zoom.py index 4fd3933..2f565fb 100644 --- a/src/ess/polarization/zoom.py +++ b/src/ess/polarization/zoom.py @@ -16,7 +16,6 @@ from ess.isissans.mantidio import DataWorkspace, Period from ess.reduce.nexus.types import Position from ess.sans.types import ( - EmptyBeamRun, Filename, Incident, MonitorType, @@ -106,7 +105,9 @@ class MonitorSpectrumNumber(Generic[MonitorType]): def get_monitor_data( - dg: LoadedFileContents[RunType], nexus_name: NeXusMonitorName[MonitorType] + dg: LoadedFileContents[RunType], + nexus_name: NeXusMonitorName[MonitorType], + spectrum_number: MonitorSpectrumNumber[MonitorType], ) -> NeXusComponent[MonitorType, RunType]: """ Same as :py:func:`ess.isissans.get_monitor_data` but dropping variances. @@ -114,29 +115,22 @@ def get_monitor_data( Dropping variances is a workaround required since ESSsans does not handle variance broadcasting when combining monitors. In our case some of the monitors are time-dependent, so this is required for now. + + If the raw files is histogram data, Mantid stores this as a Workspace2D, where some + or all spectra corresponds to monitors. """ - # See https://github.com/scipp/sciline/issues/52 why copy needed - mon = dg['monitors'][nexus_name]['data'].copy() + if 'monitors' in dg: + # From EventWorkspace + # See https://github.com/scipp/sciline/issues/52 why copy needed + mon = dg['monitors'][nexus_name]['data'].copy() + else: + # From Workspace2D + mon = sc.values(dg["data"]["spectrum", sc.index(spectrum_number.value)]).copy() return NeXusComponent[MonitorType, RunType]( sc.DataGroup(data=sc.values(mon), position=mon.coords['position']) ) -def get_monitor_data_from_empty_beam_run( - dg: LoadedFileContents[EmptyBeamRun], - spectrum_number: MonitorSpectrumNumber[MonitorType], -) -> NeXusComponent[MonitorType, EmptyBeamRun]: - """ - Extract incident or transmission monitor from ZOOM empty beam run - - The files in this case do not contain detector data, only monitor data. Mantid - stores this as a Workspace2D, where each spectrum corresponds to a monitor. - """ - # Note we index with a scipp.Variable, i.e., by the spectrum number used at ISIS - monitor = sc.values(dg["data"]["spectrum", sc.index(spectrum_number.value)]).copy() - return sc.DataGroup(data=monitor, position=monitor.coords['position']) - - def get_monitor_data_from_transmission_run( dg: LoadedFileContents[TransmissionRun[RunType]], spectrum_number: MonitorSpectrumNumber[MonitorType], @@ -153,6 +147,32 @@ def get_monitor_data_from_transmission_run( return sc.DataGroup(data=monitor, position=monitor.coords['position']) +def ZoomPolarizedWorkflow() -> sl.Pipeline: + """ + Zoom workflow for polarized data. + + This workflow is based on `ess.sans.isis.zoom.ZoomWorkflow` but is adapted to work + with histogrammed data. + """ + workflow = isis.zoom.ZoomWorkflow() + workflow.insert(get_monitor_data) + workflow.insert(load_histogrammed_run) + + # We are dealing with two different types of files, and monitors are identified + # differently in each case, so there is some duplication here. + workflow[MonitorSpectrumNumber[Incident]] = MonitorSpectrumNumber[Incident](3) + workflow[MonitorSpectrumNumber[Transmission]] = MonitorSpectrumNumber[Transmission]( + 4 + ) + workflow[NeXusMonitorName[Incident]] = NeXusMonitorName[Incident]("monitor3") + workflow[NeXusMonitorName[Transmission]] = NeXusMonitorName[Transmission]( + "monitor4" + ) + workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound + + return workflow + + def ZoomTransmissionFractionWorkflow(runs: Sequence[str]) -> sl.Pipeline: """ Workflow computing time-dependent SANS transmission fraction from ZOOM data. @@ -171,11 +191,8 @@ def ZoomTransmissionFractionWorkflow(runs: Sequence[str]) -> sl.Pipeline: runs: List of filenames of the runs to use for the transmission fraction. """ - workflow = isis.zoom.ZoomWorkflow() - workflow.insert(get_monitor_data) - workflow.insert(get_monitor_data_from_empty_beam_run) + workflow = ZoomPolarizedWorkflow() workflow.insert(get_monitor_data_from_transmission_run) - workflow.insert(load_histogrammed_run) mapped = workflow.map({Filename[TransmissionRun[SampleRun]]: runs}) for mon_type in (Incident, Transmission): @@ -186,16 +203,4 @@ def ZoomTransmissionFractionWorkflow(runs: Sequence[str]) -> sl.Pipeline: Position[NXsource, TransmissionRun[SampleRun]] ].reduce(func=_get_unique_position) - # We are dealing with two different types of files, and monitors are identified - # differently in each case, so there is some duplication here. - workflow[MonitorSpectrumNumber[Incident]] = MonitorSpectrumNumber[Incident](3) - workflow[MonitorSpectrumNumber[Transmission]] = MonitorSpectrumNumber[Transmission]( - 4 - ) - workflow[NeXusMonitorName[Incident]] = NeXusMonitorName[Incident]("monitor3") - workflow[NeXusMonitorName[Transmission]] = NeXusMonitorName[Transmission]( - "monitor4" - ) - workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound - return workflow From 9785bc158e9302a7115fe5646f111e2099eabc6b Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Wed, 4 Jun 2025 15:38:08 +0200 Subject: [PATCH 02/10] Comment --- docs/user-guide/zoom.ipynb | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 9f93032..67bf4e1 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -105,6 +105,8 @@ "metadata": {}, "outputs": [], "source": [ + "# TODO Figure out the filetypes we have as inputs. Can we not use the workflow\n", + "# from ESSsans? If the problem is histogrammed files, can we add support there?\n", "sample_workflow = ZoomPolarizedWorkflow()\n", "for key, value in params.items():\n", " sample_workflow[key] = value\n", From 4ed1e85436ea2fd7d45b1bbcf0ae0fb401120840 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Fri, 6 Jun 2025 07:11:44 +0200 Subject: [PATCH 03/10] Move ISIS hist file handling to ESSsans, add sample workflow --- docs/user-guide/zoom.ipynb | 282 +++++++++++++++++++++++------------ src/ess/polarization/zoom.py | 102 +------------ 2 files changed, 190 insertions(+), 194 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 67bf4e1..9922835 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -18,6 +18,7 @@ "outputs": [], "source": [ "%matplotlib widget\n", + "import plopp\n", "import sciline\n", "import scipp as sc\n", "from ess import polarization as pol\n", @@ -44,21 +45,9 @@ "cell_runs = [*cell_runs, depolarized_run]" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "3", - "metadata": {}, - "outputs": [], - "source": [ - "data_folder = Path('/home/simon/instruments/polarization/2024-08')\n", - "sample_runs = [str(data_folder / f'ZOOM00038{run}.nxs') for run in [513]]\n", - "empty_run = data_folder / 'ZOOM00038336.nxs'" - ] - }, { "cell_type": "markdown", - "id": "4", + "id": "3", "metadata": {}, "source": [ "## Setup SANS workflow for computing transmission fraction" @@ -67,85 +56,46 @@ { "cell_type": "code", "execution_count": null, - "id": "5", + "id": "4", "metadata": {}, "outputs": [], "source": [ "params = {\n", - " Filename[EmptyBeamRun]: str(empty_run),\n", " WavelengthBins: sc.geomspace(\n", " 'wavelength', start=1.75, stop=16.5, num=141, unit='Å'\n", " ),\n", - " isis.mantidio.Period: 0,\n", " isis.mantidio.CalibrationWorkspace: None,\n", + " # Spectrum numbers used for histogram files (Workspace2D)\n", + " isis.MonitorSpectrumNumber[Incident]: isis.MonitorSpectrumNumber[Incident](3),\n", + " isis.MonitorSpectrumNumber[Transmission]: isis.MonitorSpectrumNumber[Transmission](\n", + " 4\n", + " ),\n", + " # Monitor names used for event files (EventWorkspace)\n", + " NeXusMonitorName[Incident]: NeXusMonitorName[Incident](\"monitor3\"),\n", + " NeXusMonitorName[Transmission]: NeXusMonitorName[Transmission](\"monitor4\"),\n", + " UncertaintyBroadcastMode: UncertaintyBroadcastMode.upper_bound,\n", "}" ] }, { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "5", "metadata": {}, "outputs": [], "source": [ - "from ess.polarization.zoom import (\n", - " ZoomTransmissionFractionWorkflow,\n", - " ZoomPolarizedWorkflow,\n", - ")\n", + "from ess.polarization.zoom import ZoomTransmissionFractionWorkflow\n", "\n", "sans_workflow = ZoomTransmissionFractionWorkflow(cell_runs)\n", "for key, value in params.items():\n", - " sans_workflow[key] = value" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7", - "metadata": {}, - "outputs": [], - "source": [ - "# TODO Figure out the filetypes we have as inputs. Can we not use the workflow\n", - "# from ESSsans? If the problem is histogrammed files, can we add support there?\n", - "sample_workflow = ZoomPolarizedWorkflow()\n", - "for key, value in params.items():\n", - " sample_workflow[key] = value\n", - "masks = isis.data.zoom_tutorial_mask_filenames()\n", - "sample_workflow = sans.with_pixel_mask_filenames(sample_workflow, masks)\n", - "sample_workflow.insert(isis.io.transmission_from_background_run)\n", - "sample_workflow.insert(isis.io.transmission_from_sample_run)\n", - "sample_workflow[QBins] = sc.geomspace(\n", - " dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom'\n", - ")\n", - "sample_workflow[NonBackgroundWavelengthRange] = sc.array(\n", - " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", - ")\n", - "sample_workflow[CorrectForGravity] = False\n", - "sample_workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop\n", - "sample_workflow[ReturnEvents] = False\n", - "\n", - "# TODO\n", - "# sample_workflow[DirectBeamFilename]\n", - "sample_workflow[DirectBeam] = None\n", - "\n", - "sample_workflow[Filename[SampleRun]] = sample_runs[0]\n", - "# sample_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", - "# sample_workflow[isis.SampleOffset] = sc.vector([0.0, 0.0, 0.11], unit='m')\n", - "# sample_workflow[isis.DetectorBankOffset] = sc.vector([0.0, 0.0, 0.5], unit='m')\n", - "sample_workflow[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m')\n", - "sample_workflow.visualize(IofQ[SampleRun], graph_attr={'rankdir': 'LR'})\n", - "# sample_workflow.compute(NeXusComponent[Incident, SampleRun])\n", - "# sample_workflow.compute(NeXusComponent[Transmission, SampleRun])\n", - "# sample_workflow.compute(DetectorData[SampleRun])\n", - "# sample_workflow.compute(IofQ[SampleRun])\n", - "# ws = sample_workflow.compute(isis.mantidio.DataWorkspace[SampleRun])\n", - "dg = sample_workflow.compute(isis.mantidio.LoadedFileContents[SampleRun])\n", - "dg" + " sans_workflow[key] = value\n", + "sans_workflow[isis.mantidio.Period] = 0\n", + "sans_workflow[Filename[EmptyBeamRun]] = str(empty_run)" ] }, { "cell_type": "markdown", - "id": "8", + "id": "6", "metadata": {}, "source": [ "## Inspect data for one of the runs with analyzer cell" @@ -154,7 +104,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9", + "id": "7", "metadata": {}, "outputs": [], "source": [ @@ -169,7 +119,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +128,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "9", "metadata": {}, "source": [ "We can load the combined time-dependent incident and transmission monitors.\n", @@ -188,7 +138,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -208,7 +158,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "11", "metadata": {}, "source": [ "The task graph for computing the transmission fraction is:" @@ -217,7 +167,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -226,7 +176,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "13", "metadata": {}, "source": [ "## Compute transmission fractions\n", @@ -238,7 +188,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -247,7 +197,7 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "15", "metadata": {}, "source": [ "We can plot the computed transmission fractions:" @@ -256,7 +206,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -272,7 +222,7 @@ { "cell_type": "code", "execution_count": null, - "id": "19", + "id": "17", "metadata": {}, "outputs": [], "source": [ @@ -293,7 +243,7 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -307,7 +257,7 @@ }, { "cell_type": "markdown", - "id": "21", + "id": "19", "metadata": {}, "source": [ "We can now setup the polarization analysis workflow.\n", @@ -317,7 +267,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -347,7 +297,7 @@ }, { "cell_type": "markdown", - "id": "23", + "id": "21", "metadata": {}, "source": [ "The workflow can compute the transmission function:" @@ -356,7 +306,7 @@ { "cell_type": "code", "execution_count": null, - "id": "24", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -365,7 +315,7 @@ }, { "cell_type": "markdown", - "id": "25", + "id": "23", "metadata": {}, "source": [ "We can evaluate this transmission function at desired time and wavelength points:" @@ -374,7 +324,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -389,7 +339,7 @@ { "cell_type": "code", "execution_count": null, - "id": "27", + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -402,7 +352,7 @@ { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -415,7 +365,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "27", "metadata": {}, "outputs": [], "source": [ @@ -425,7 +375,7 @@ { "cell_type": "code", "execution_count": null, - "id": "30", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -435,7 +385,7 @@ { "cell_type": "code", "execution_count": null, - "id": "31", + "id": "29", "metadata": {}, "outputs": [], "source": [ @@ -444,7 +394,7 @@ }, { "cell_type": "markdown", - "id": "32", + "id": "30", "metadata": {}, "source": [ "## Correction workflow\n", @@ -457,7 +407,7 @@ { "cell_type": "code", "execution_count": null, - "id": "33", + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -467,7 +417,7 @@ }, { "cell_type": "markdown", - "id": "34", + "id": "32", "metadata": {}, "source": [ "We will use a second-order polynomial supermirror efficiency function:" @@ -476,7 +426,7 @@ { "cell_type": "code", "execution_count": null, - "id": "35", + "id": "33", "metadata": {}, "outputs": [], "source": [ @@ -493,7 +443,7 @@ { "cell_type": "code", "execution_count": null, - "id": "36", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -505,7 +455,7 @@ }, { "cell_type": "markdown", - "id": "37", + "id": "35", "metadata": {}, "source": [ "For a single channel, the complete workflow looks as follows:" @@ -514,7 +464,7 @@ { "cell_type": "code", "execution_count": null, - "id": "38", + "id": "36", "metadata": {}, "outputs": [], "source": [ @@ -522,6 +472,142 @@ " pol.PolarizationCorrectedData[pol.Up, pol.Up], graph_attr={'rankdir': 'LR'}\n", ")" ] + }, + { + "cell_type": "markdown", + "id": "37", + "metadata": {}, + "source": [ + "## Sample data\n", + "\n", + "We now proceed to process the sample data.\n", + "We use the Zoom workflow from ESSsans, configured as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38", + "metadata": {}, + "outputs": [], + "source": [ + "def make_sample_workflow() -> sciline.Pipeline:\n", + " wf = isis.zoom.ZoomWorkflow()\n", + " wf.insert(isis.io.transmission_from_background_run)\n", + " wf.insert(isis.io.transmission_from_sample_run)\n", + "\n", + " # TODO Are these masks correct? Provide correct XML files!\n", + " masks = isis.data.zoom_tutorial_mask_filenames()\n", + " wf = sans.with_pixel_mask_filenames(wf, masks)\n", + "\n", + " for key, value in params.items():\n", + " wf[key] = value\n", + " wf[QBins] = sc.geomspace(dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom')\n", + " wf[NonBackgroundWavelengthRange] = sc.array(\n", + " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", + " )\n", + " wf[CorrectForGravity] = True\n", + " wf[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", + " wf[ReturnEvents] = True\n", + "\n", + " # TODO Configure offsets and beam center correctly!\n", + " # sample_workflow[isis.SampleOffset] = sc.vector([0.0, 0.0, 0.11], unit='m')\n", + " # sample_workflow[isis.DetectorBankOffset] = sc.vector([0.0, 0.0, 0.5], unit='m')\n", + " wf[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m')\n", + "\n", + " # TODO Use direct beam!\n", + " # sample_workflow[DirectBeamFilename]\n", + " wf[DirectBeam] = None\n", + "\n", + " return wf" + ] + }, + { + "cell_type": "markdown", + "id": "39", + "metadata": {}, + "source": [ + "The sample runs are event data, so we can simply merge all runs to obtain a single event-mode $I(Q)$.\n", + "We use `sans.with_sample_runs` to obtain a workflow that merges the runs and computes the $I(Q)$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40", + "metadata": {}, + "outputs": [], + "source": [ + "water_folder = Path('zoom_water_data')\n", + "sample_runs = [str(water_folder / f'ZOOM00038{run}.nxs') for run in range(439, 466, 2)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41", + "metadata": {}, + "outputs": [], + "source": [ + "sample_workflow = make_sample_workflow()\n", + "# TODO Is this the correct empty beam run?\n", + "sample_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", + "multi_run_sample_workflow = sans.with_sample_runs(sample_workflow, runs=sample_runs)" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": {}, + "source": [ + "The sample workflow looks like this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43", + "metadata": {}, + "outputs": [], + "source": [ + "multi_run_sample_workflow.visualize(\n", + " IofQ[SampleRun], graph_attr={'rankdir': 'LR'}, compact=True\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44", + "metadata": {}, + "outputs": [], + "source": [ + "periods = []\n", + "for period in range(4):\n", + " multi_run_sample_workflow[isis.mantidio.Period] = period\n", + " periods.append(multi_run_sample_workflow.compute(IofQ[SampleRun]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "45", + "metadata": {}, + "outputs": [], + "source": [ + "dg = sc.DataGroup({f'period {period}': periods[period] for period in range(4)})\n", + "plopp.slicer(sc.values(dg.hist(pulse_time=10)), keep='Q')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46", + "metadata": {}, + "outputs": [], + "source": [ + "periods[0].hist(pulse_time=100).plot(norm='log', scale={'Q': 'log'})" + ] } ], "metadata": { diff --git a/src/ess/polarization/zoom.py b/src/ess/polarization/zoom.py index 2f565fb..6398e68 100644 --- a/src/ess/polarization/zoom.py +++ b/src/ess/polarization/zoom.py @@ -1,19 +1,13 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) -import threading from collections.abc import Sequence -from dataclasses import dataclass -from typing import Generic -import mantid.api as _mantid_api import sciline as sl import scipp as sc -from mantid import simpleapi as _mantid_simpleapi from scippnexus import NXsource import ess.isissans as isis from ess.isissans.io import LoadedFileContents -from ess.isissans.mantidio import DataWorkspace, Period from ess.reduce.nexus.types import Position from ess.sans.types import ( Filename, @@ -25,7 +19,6 @@ SampleRun, Transmission, TransmissionRun, - UncertaintyBroadcastMode, ) # In this case the "sample" is the analyzer cell, of which we want to measure @@ -33,45 +26,6 @@ sample_run_type = RunType -def load_histogrammed_run( - filename: Filename[sample_run_type], period: Period -) -> DataWorkspace[sample_run_type]: - """Load a non-event-data ISIS file""" - # Loading many small files with Mantid is, for some reason, very slow when using - # the default number of threads in the Dask threaded scheduler (1 thread worked - # best, 2 is a bit slower but still fast). We can either limit that thread count, - # or add a lock here, which is more specific. - with load_histogrammed_run.lock: - loaded = _mantid_simpleapi.Load(Filename=str(filename), StoreInADS=False) - if isinstance(loaded, _mantid_api.Workspace): - # A single workspace - data_ws = loaded - if isinstance(data_ws, _mantid_api.WorkspaceGroup): - if period is None: - raise ValueError( - f'Needs {Period} to be set to know what ' - 'section of the event data to load' - ) - data_ws = data_ws.getItem(period) - else: - # Separate data and monitor workspaces - data_ws = loaded.OutputWorkspace - if isinstance(data_ws, _mantid_api.WorkspaceGroup): - if period is None: - raise ValueError( - f'Needs {Period} to be set to know what ' - 'section of the event data to load' - ) - data_ws = data_ws.getItem(period) - data_ws.setMonitorWorkspace(loaded.MonitorWorkspace.getItem(period)) - else: - data_ws.setMonitorWorkspace(loaded.MonitorWorkspace) - return DataWorkspace[sample_run_type](data_ws) - - -load_histogrammed_run.lock = threading.Lock() - - def _get_time(dg: sc.DataGroup) -> sc.Variable: start = sc.datetime(dg['run_start'].value) end = sc.datetime(dg['run_end'].value) @@ -99,41 +53,23 @@ def _get_unique_position(*positions: sc.DataArray) -> sc.DataArray: return unique -@dataclass -class MonitorSpectrumNumber(Generic[MonitorType]): - value: int - - def get_monitor_data( dg: LoadedFileContents[RunType], nexus_name: NeXusMonitorName[MonitorType], - spectrum_number: MonitorSpectrumNumber[MonitorType], + spectrum_number: isis.MonitorSpectrumNumber[MonitorType], ) -> NeXusComponent[MonitorType, RunType]: """ Same as :py:func:`ess.isissans.get_monitor_data` but dropping variances. - - Dropping variances is a workaround required since ESSsans does not handle - variance broadcasting when combining monitors. In our case some of the monitors - are time-dependent, so this is required for now. - - If the raw files is histogram data, Mantid stores this as a Workspace2D, where some - or all spectra corresponds to monitors. """ - if 'monitors' in dg: - # From EventWorkspace - # See https://github.com/scipp/sciline/issues/52 why copy needed - mon = dg['monitors'][nexus_name]['data'].copy() - else: - # From Workspace2D - mon = sc.values(dg["data"]["spectrum", sc.index(spectrum_number.value)]).copy() - return NeXusComponent[MonitorType, RunType]( - sc.DataGroup(data=sc.values(mon), position=mon.coords['position']) + monitor = isis.general.get_monitor_data( + dg, nexus_name=nexus_name, spectrum_number=spectrum_number ) + return NeXusComponent[MonitorType, RunType](sc.values(monitor)) def get_monitor_data_from_transmission_run( dg: LoadedFileContents[TransmissionRun[RunType]], - spectrum_number: MonitorSpectrumNumber[MonitorType], + spectrum_number: isis.MonitorSpectrumNumber[MonitorType], ) -> NeXusComponent[MonitorType, TransmissionRun[RunType]]: """ Extract incident or transmission monitor from ZOOM direct-beam run @@ -147,32 +83,6 @@ def get_monitor_data_from_transmission_run( return sc.DataGroup(data=monitor, position=monitor.coords['position']) -def ZoomPolarizedWorkflow() -> sl.Pipeline: - """ - Zoom workflow for polarized data. - - This workflow is based on `ess.sans.isis.zoom.ZoomWorkflow` but is adapted to work - with histogrammed data. - """ - workflow = isis.zoom.ZoomWorkflow() - workflow.insert(get_monitor_data) - workflow.insert(load_histogrammed_run) - - # We are dealing with two different types of files, and monitors are identified - # differently in each case, so there is some duplication here. - workflow[MonitorSpectrumNumber[Incident]] = MonitorSpectrumNumber[Incident](3) - workflow[MonitorSpectrumNumber[Transmission]] = MonitorSpectrumNumber[Transmission]( - 4 - ) - workflow[NeXusMonitorName[Incident]] = NeXusMonitorName[Incident]("monitor3") - workflow[NeXusMonitorName[Transmission]] = NeXusMonitorName[Transmission]( - "monitor4" - ) - workflow[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound - - return workflow - - def ZoomTransmissionFractionWorkflow(runs: Sequence[str]) -> sl.Pipeline: """ Workflow computing time-dependent SANS transmission fraction from ZOOM data. @@ -191,7 +101,7 @@ def ZoomTransmissionFractionWorkflow(runs: Sequence[str]) -> sl.Pipeline: runs: List of filenames of the runs to use for the transmission fraction. """ - workflow = ZoomPolarizedWorkflow() + workflow = isis.zoom.ZoomWorkflow() workflow.insert(get_monitor_data_from_transmission_run) mapped = workflow.map({Filename[TransmissionRun[SampleRun]]: runs}) From ad62c580c3507834d37f54986a9c25a6ad940e5b Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Fri, 6 Jun 2025 10:27:33 +0200 Subject: [PATCH 04/10] Return variable instead of data array from trans functions --- src/ess/polarization/correction.py | 3 ++ src/ess/polarization/he3.py | 9 ++--- src/ess/polarization/supermirror.py | 12 +++---- src/ess/polarization/types.py | 2 +- src/ess/polarization/zoom.py | 6 ++-- tests/correction_test.py | 24 +++++++++++++ tests/he3/opacity_test.py | 6 ++-- tests/he3/polarization_test.py | 54 ++++++++++++++++++++++------- tests/supermirror_test.py | 5 +-- 9 files changed, 87 insertions(+), 34 deletions(-) diff --git a/src/ess/polarization/correction.py b/src/ess/polarization/correction.py index 83528f0..3c1c5fd 100644 --- a/src/ess/polarization/correction.py +++ b/src/ess/polarization/correction.py @@ -103,6 +103,9 @@ def compute_polarizing_element_correction( : Correction matrix coefficients. """ + if isinstance(channel, sc.Variable | sc.DataArray) and channel.bins is not None: + channel = channel.bins + t_plus = transmission.apply(channel, 'plus') t_minus = transmission.apply(channel, 'minus') t_minus *= -1 diff --git a/src/ess/polarization/he3.py b/src/ess/polarization/he3.py index 47153be..42483e1 100644 --- a/src/ess/polarization/he3.py +++ b/src/ess/polarization/he3.py @@ -99,10 +99,7 @@ def __call__(self, wavelength: sc.Variable) -> sc.Variable: scale = broadcast_with_upper_bound_variances( self.opacity0, prototype=wavelength ) - return sc.DataArray( - (scale * wavelength).to(unit='', copy=False), - coords={'wavelength': wavelength}, - ) + return (scale * wavelength).to(unit='', copy=False) def he3_opacity_from_cell_params( @@ -192,7 +189,7 @@ def T1(self) -> sc.Variable: return self._T1 def __call__(self, time: sc.Variable) -> sc.Variable: - return sc.DataArray(self.C * sc.exp(-time / self.T1), coords={'time': time}) + return self.C * sc.exp(-time / self.T1) @dataclass @@ -214,7 +211,7 @@ def __call__( polarization *= -plus_minus return self.transmission_empty_glass * sc.exp(-opacity * (1.0 + polarization)) - def apply(self, data: sc.DataArray, plus_minus: PlusMinus) -> sc.DataArray: + def apply(self, data: sc.DataArray, plus_minus: PlusMinus) -> sc.Variable: return self( time=data.coords['time'], wavelength=data.coords['wavelength'], diff --git a/src/ess/polarization/supermirror.py b/src/ess/polarization/supermirror.py index d983ecd..1276117 100644 --- a/src/ess/polarization/supermirror.py +++ b/src/ess/polarization/supermirror.py @@ -17,7 +17,7 @@ class SupermirrorEfficiencyFunction(Generic[PolarizingElement], ABC): """Base class for supermirror efficiency functions""" @abstractmethod - def __call__(self, *, wavelength: sc.Variable) -> sc.DataArray: + def __call__(self, *, wavelength: sc.Variable) -> sc.Variable: """Return the efficiency of a supermirror for a given wavelength""" @@ -44,7 +44,7 @@ class SecondDegreePolynomialEfficiency( b: sc.Variable c: sc.Variable - def __call__(self, *, wavelength: sc.Variable) -> sc.DataArray: + def __call__(self, *, wavelength: sc.Variable) -> sc.Variable: """Return the efficiency of a supermirror for a given wavelength""" return ( (self.a * wavelength**2).to(unit='', copy=False) @@ -71,9 +71,9 @@ def __post_init__(self): table = self.table if self.table.variances is None else sc.values(self.table) self._lut = sc.lookup(table, 'wavelength') - def __call__(self, *, wavelength: sc.Variable) -> sc.DataArray: + def __call__(self, *, wavelength: sc.Variable) -> sc.Variable: """Return the efficiency of a supermirror for a given wavelength""" - return sc.DataArray(self._lut(wavelength), coords={'wavelength': wavelength}) + return self._lut(wavelength) @classmethod def from_file( @@ -107,7 +107,7 @@ class SupermirrorTransmissionFunction(TransmissionFunction[PolarizingElement]): def __call__( self, *, wavelength: sc.Variable, plus_minus: PlusMinus - ) -> sc.DataArray: + ) -> sc.Variable: """Return the transmission fraction for a given wavelength""" efficiency = self.efficiency_function(wavelength=wavelength) if plus_minus == 'plus': @@ -115,7 +115,7 @@ def __call__( else: return 0.5 * (1 - efficiency) - def apply(self, data: sc.DataArray, plus_minus: PlusMinus) -> sc.DataArray: + def apply(self, data: sc.DataArray, plus_minus: PlusMinus) -> sc.Variable: """Apply the transmission function to a data array""" return self(wavelength=data.coords['wavelength'], plus_minus=plus_minus) diff --git a/src/ess/polarization/types.py b/src/ess/polarization/types.py index 16a28c8..2ab50ae 100644 --- a/src/ess/polarization/types.py +++ b/src/ess/polarization/types.py @@ -32,7 +32,7 @@ class TransmissionFunction(Generic[PolarizingElement], ABC): """Wavelength- and time-dependent transmission for a given cell.""" @abstractmethod - def apply(self, data: sc.DataArray, plus_minus: PlusMinus) -> sc.DataArray: ... + def apply(self, data: sc.DataArray, plus_minus: PlusMinus) -> sc.Variable: ... @dataclass diff --git a/src/ess/polarization/zoom.py b/src/ess/polarization/zoom.py index 6398e68..384fc5f 100644 --- a/src/ess/polarization/zoom.py +++ b/src/ess/polarization/zoom.py @@ -53,7 +53,7 @@ def _get_unique_position(*positions: sc.DataArray) -> sc.DataArray: return unique -def get_monitor_data( +def get_monitor_data_no_variances( dg: LoadedFileContents[RunType], nexus_name: NeXusMonitorName[MonitorType], spectrum_number: isis.MonitorSpectrumNumber[MonitorType], @@ -64,7 +64,8 @@ def get_monitor_data( monitor = isis.general.get_monitor_data( dg, nexus_name=nexus_name, spectrum_number=spectrum_number ) - return NeXusComponent[MonitorType, RunType](sc.values(monitor)) + monitor['data'] = sc.values(monitor['data']) + return NeXusComponent[MonitorType, RunType](monitor) def get_monitor_data_from_transmission_run( @@ -102,6 +103,7 @@ def ZoomTransmissionFractionWorkflow(runs: Sequence[str]) -> sl.Pipeline: List of filenames of the runs to use for the transmission fraction. """ workflow = isis.zoom.ZoomWorkflow() + workflow.insert(get_monitor_data_no_variances) workflow.insert(get_monitor_data_from_transmission_run) mapped = workflow.map({Filename[TransmissionRun[SampleRun]]: runs}) diff --git a/tests/correction_test.py b/tests/correction_test.py index c7c9562..fe356d3 100644 --- a/tests/correction_test.py +++ b/tests/correction_test.py @@ -68,6 +68,30 @@ def test_compute_polarizing_element_correction() -> None: assert_allclose(off_diag, -transmission_minus / denom) +def test_compute_polarizing_element_correction_binned_data() -> None: + time = sc.linspace('event', 1, 10, 10, unit='') + wavelength = sc.linspace('event', 0.1, 1, 10, unit='') + events = sc.DataArray( + sc.arange('event', 10), + coords={'time': time, 'wavelength': wavelength}, + ) + binned = sc.bins( + data=events, dim='event', begin=sc.array(dims=['Q'], values=[0, 3], unit=None) + ) + transmission = SimpleTransmissionFunction() + + result = compute_polarizing_element_correction( + channel=binned, transmission=transmission + ) + diag = result.diag + off_diag = result.off_diag + transmission_plus = transmission(time, wavelength, 'plus') + transmission_minus = transmission(time, wavelength, 'minus') + denom = transmission_plus**2 - transmission_minus**2 + assert_allclose(diag.bins.concat().value, transmission_plus / denom) + assert_allclose(off_diag.bins.concat().value, -transmission_minus / denom) + + class FakeTransmissionFunction: def __init__(self, coeffs: np.ndarray) -> None: self.coeffs = coeffs diff --git a/tests/he3/opacity_test.py b/tests/he3/opacity_test.py index 8ddaaab..3f4af08 100644 --- a/tests/he3/opacity_test.py +++ b/tests/he3/opacity_test.py @@ -18,7 +18,7 @@ def test_opacity_from_cell_params() -> None: pressure=pressure, length=length, temperature=temperature ) opacity_function = he3.he3_opacity_function_from_cell_opacity(opacity0) - opacity = opacity_function(wavelength).data + opacity = opacity_function(wavelength) assert_identical(2 * opacity['pressure', 0], opacity['pressure', 1]) assert_identical(2 * opacity['cell_length', 0], opacity['cell_length', 1]) assert_identical(2 * opacity['wavelength', 0], opacity['wavelength', 1]) @@ -38,7 +38,7 @@ def test_opacity_from_cell_params_reproduces_literature_value() -> None: pressure=pressure, length=length, temperature=temperature ) opacity_function = he3.he3_opacity_function_from_cell_opacity(opacity0) - opacity = opacity_function(wavelength).data + opacity = opacity_function(wavelength) assert sc.isclose(opacity, sc.scalar(0.0733, unit=''), rtol=sc.scalar(1e-3)) @@ -64,7 +64,7 @@ def test_opacity_from_beam_data() -> None: transmission_fraction=transmission, opacity0_initial_guess=opacity0 * 1.23, # starting guess imperfect ) - opacity = opacity_function(wavelength).data + opacity = opacity_function(wavelength) assert sc.isclose( opacity_function.opacity0, opacity0.to(unit=opacity_function.opacity0.unit) ) diff --git a/tests/he3/polarization_test.py b/tests/he3/polarization_test.py index f19aaba..49a9312 100644 --- a/tests/he3/polarization_test.py +++ b/tests/he3/polarization_test.py @@ -18,10 +18,16 @@ def test_incoming_unpolarized_reproduces_input_params_within_errors() -> None: transmission_empty_glass = sc.scalar(0.9) opacity = opacity_function(wavelength) polarization = polarization_function(time) - transmission = he3.transmission_incoming_unpolarized( - transmission_empty_glass=transmission_empty_glass, - opacity=opacity, - polarization=polarization, + transmission = sc.DataArray( + he3.transmission_incoming_unpolarized( + transmission_empty_glass=transmission_empty_glass, + opacity=opacity, + polarization=polarization, + ), + coords={ + 'time': time, + 'wavelength': wavelength, + }, ) result = he3.get_he3_transmission_incoming_unpolarized_from_fit_to_direct_beam( @@ -67,11 +73,23 @@ def test_incoming_polarized_reproduces_input_params_within_errors() -> None: polarization_function=polarization_function, ) # State switch at 456th time point, cut further below into 4 channels total. - plus = transmission_function( - time=time[:456], wavelength=wavelength, plus_minus='plus' + plus = sc.DataArray( + transmission_function( + time=time[:456], wavelength=wavelength, plus_minus='plus' + ), + coords={ + 'time': time[:456], + 'wavelength': wavelength, + }, ) - minus = transmission_function( - time=time[456:], wavelength=wavelength, plus_minus='minus' + minus = sc.DataArray( + transmission_function( + time=time[456:], wavelength=wavelength, plus_minus='minus' + ), + coords={ + 'time': time[456:], + 'wavelength': wavelength, + }, ) result = he3.get_he3_transmission_incoming_polarized_from_fit_to_direct_beam( @@ -121,11 +139,23 @@ def test_incoming_polarized_raises_if_plus_minus_coord_is_bad() -> None: polarization_function=polarization_function, ) # State switch at 456th time point, cut further below into 4 channels total. - plus = transmission_function( - time=time[:456], wavelength=wavelength, plus_minus='plus' + plus = sc.DataArray( + transmission_function( + time=time[:456], wavelength=wavelength, plus_minus='plus' + ), + coords={ + 'time': time[:456], + 'wavelength': wavelength, + }, ) - minus = transmission_function( - time=time[456:], wavelength=wavelength, plus_minus='minus' + minus = sc.DataArray( + transmission_function( + time=time[456:], wavelength=wavelength, plus_minus='minus' + ), + coords={ + 'time': time[456:], + 'wavelength': wavelength, + }, ) with pytest.raises(ValueError, match='plus-minus coordinate of plus channel'): diff --git a/tests/supermirror_test.py b/tests/supermirror_test.py index a9afd44..2ff4de3 100644 --- a/tests/supermirror_test.py +++ b/tests/supermirror_test.py @@ -78,10 +78,7 @@ def test_EfficiencyLookupTable_returns_expected_result(): ) ) x = sc.midpoints(sc.linspace('wavelength', 0, 1, 10, unit='angstrom')) - assert_allclose( - tab(wavelength=x), - sc.DataArray(sc.values(tab.table.data), coords={'wavelength': x}), - ) + assert_allclose(tab(wavelength=x), sc.values(tab.table.data)) def test_EfficiencyLookupTable_load_from_file(): From a2c3ad955e6bbf4249d93775f1124a0ea6733240 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Fri, 6 Jun 2025 10:40:17 +0200 Subject: [PATCH 05/10] "Finalize" Zoom notebook --- docs/user-guide/zoom.ipynb | 234 ++++++++++++++++++++++++------------- 1 file changed, 152 insertions(+), 82 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 9922835..a15170c 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -40,7 +40,7 @@ "data_folder = Path('zoom_polarized_data')\n", "# Runs with analyzer at 4 different times\n", "cell_runs = [str(data_folder / f'ZOOM00022{run}.nxs') for run in [710, 712, 714, 716]]\n", - "empty_run = data_folder / 'ZOOM00034787.nxs'\n", + "empty_run = data_folder / 'ZOOM00038238.nxs'\n", "depolarized_run = data_folder / 'ZOOM00022718.nxs'\n", "cell_runs = [*cell_runs, depolarized_run]" ] @@ -171,7 +171,9 @@ "metadata": {}, "outputs": [], "source": [ - "sans_workflow.visualize(TransmissionFraction[SampleRun], graph_attr={'rankdir': 'LR'})" + "sans_workflow.visualize(\n", + " TransmissionFraction[SampleRun], graph_attr={'rankdir': 'LR'}, compact=True\n", + ")" ] }, { @@ -343,7 +345,13 @@ "metadata": {}, "outputs": [], "source": [ - "trans = func(wavelength=wavelength, time=time, plus_minus='plus')\n", + "trans = sc.DataArray(\n", + " func(wavelength=wavelength, time=time, plus_minus='plus'),\n", + " coords={\n", + " 'wavelength': wavelength,\n", + " 'time': time,\n", + " },\n", + ")\n", "sc.DataGroup(\n", " {f\"{time:c}\": trans['time', time] for time in trans.coords['time'][::20]}\n", ").plot(norm='linear', linestyle='solid', marker=None)" @@ -356,7 +364,13 @@ "metadata": {}, "outputs": [], "source": [ - "trans = func(wavelength=wavelength, time=time, plus_minus='plus')\n", + "trans = sc.DataArray(\n", + " func(wavelength=wavelength, time=time, plus_minus='plus'),\n", + " coords={\n", + " 'wavelength': wavelength,\n", + " 'time': time,\n", + " },\n", + ")\n", "sc.DataGroup(\n", " {f\"{wav:c}\": trans['wavelength', wav] for wav in trans.coords['wavelength'][::20]}\n", ").plot(norm='log', linestyle='solid', marker=None)" @@ -397,11 +411,10 @@ "id": "30", "metadata": {}, "source": [ - "## Correction workflow\n", + "## Sample data\n", "\n", - "In the previous section we have setup the workflow for the analyzer.\n", - "We also computed the transmission function there, but in production this will be done implicitly by running the entire workflow we will setup here.\n", - "We can combine this with the workflow for the polarizer to obtain the full correction workflow:" + "We now proceed to process the sample data.\n", + "We use the Zoom workflow from ESSsans, configured as follows:" ] }, { @@ -411,8 +424,35 @@ "metadata": {}, "outputs": [], "source": [ - "supermirror_workflow = pol.SupermirrorWorkflow()\n", - "supermirror_workflow.visualize(pol.TransmissionFunction[pol.Polarizer])" + "def make_sample_workflow() -> sciline.Pipeline:\n", + " wf = isis.zoom.ZoomWorkflow()\n", + " wf.insert(isis.io.transmission_from_background_run)\n", + " wf.insert(isis.io.transmission_from_sample_run)\n", + "\n", + " # TODO Are these masks correct? Provide correct XML files!\n", + " masks = isis.data.zoom_tutorial_mask_filenames()\n", + " wf = sans.with_pixel_mask_filenames(wf, masks)\n", + "\n", + " for key, value in params.items():\n", + " wf[key] = value\n", + " wf[QBins] = sc.geomspace(dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom')\n", + " wf[NonBackgroundWavelengthRange] = sc.array(\n", + " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", + " )\n", + " wf[CorrectForGravity] = True\n", + " wf[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", + " wf[ReturnEvents] = True\n", + "\n", + " # TODO Configure offsets and beam center correctly!\n", + " # sample_workflow[isis.SampleOffset] = sc.vector([0.0, 0.0, 0.11], unit='m')\n", + " # sample_workflow[isis.DetectorBankOffset] = sc.vector([0.0, 0.0, 0.5], unit='m')\n", + " wf[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m')\n", + "\n", + " # TODO Use direct beam!\n", + " # sample_workflow[DirectBeamFilename]\n", + " wf[DirectBeam] = None\n", + "\n", + " return wf" ] }, { @@ -420,7 +460,8 @@ "id": "32", "metadata": {}, "source": [ - "We will use a second-order polynomial supermirror efficiency function:" + "The sample runs are event data, so we can simply merge all runs to obtain a single event-mode $I(Q)$.\n", + "We use `sans.with_sample_runs` to obtain a workflow that merges the runs and computes the $I(Q)$:" ] }, { @@ -430,14 +471,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Note that these coefficients are meaningless, please fill in correct values!\n", - "supermirror_workflow[pol.SupermirrorEfficiencyFunction[pol.Polarizer]] = (\n", - " pol.SecondDegreePolynomialEfficiency(\n", - " a=0.5 * sc.Unit('1/angstrom**2'),\n", - " b=0.4 * sc.Unit('1/angstrom'),\n", - " c=0.3 * sc.Unit('dimensionless'),\n", - " )\n", - ")" + "water_folder = Path('zoom_water_data')\n", + "sample_runs = [str(water_folder / f'ZOOM00038{run}.nxs') for run in range(439, 466, 2)]" ] }, { @@ -447,10 +482,10 @@ "metadata": {}, "outputs": [], "source": [ - "workflow = pol.PolarizationAnalysisWorkflow(\n", - " polarizer_workflow=supermirror_workflow,\n", - " analyzer_workflow=he3_workflow,\n", - ")" + "sample_workflow = make_sample_workflow()\n", + "# TODO Is this the correct empty beam run?\n", + "sample_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", + "multi_run_sample_workflow = sans.with_sample_runs(sample_workflow, runs=sample_runs)" ] }, { @@ -458,7 +493,7 @@ "id": "35", "metadata": {}, "source": [ - "For a single channel, the complete workflow looks as follows:" + "The sample workflow looks like this:" ] }, { @@ -468,8 +503,8 @@ "metadata": {}, "outputs": [], "source": [ - "workflow.visualize(\n", - " pol.PolarizationCorrectedData[pol.Up, pol.Up], graph_attr={'rankdir': 'LR'}\n", + "multi_run_sample_workflow.visualize(\n", + " IofQ[SampleRun], graph_attr={'rankdir': 'LR'}, compact=True\n", ")" ] }, @@ -478,10 +513,9 @@ "id": "37", "metadata": {}, "source": [ - "## Sample data\n", - "\n", - "We now proceed to process the sample data.\n", - "We use the Zoom workflow from ESSsans, configured as follows:" + "ESSsans does not support processing multiple periods (and there are no plans for this, as this is ISIS-specific).\n", + "We therefore run the workflow once per period.\n", + "Note that is is inefficient as it loads every file four times." ] }, { @@ -491,44 +525,31 @@ "metadata": {}, "outputs": [], "source": [ - "def make_sample_workflow() -> sciline.Pipeline:\n", - " wf = isis.zoom.ZoomWorkflow()\n", - " wf.insert(isis.io.transmission_from_background_run)\n", - " wf.insert(isis.io.transmission_from_sample_run)\n", - "\n", - " # TODO Are these masks correct? Provide correct XML files!\n", - " masks = isis.data.zoom_tutorial_mask_filenames()\n", - " wf = sans.with_pixel_mask_filenames(wf, masks)\n", - "\n", - " for key, value in params.items():\n", - " wf[key] = value\n", - " wf[QBins] = sc.geomspace(dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom')\n", - " wf[NonBackgroundWavelengthRange] = sc.array(\n", - " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", - " )\n", - " wf[CorrectForGravity] = True\n", - " wf[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n", - " wf[ReturnEvents] = True\n", - "\n", - " # TODO Configure offsets and beam center correctly!\n", - " # sample_workflow[isis.SampleOffset] = sc.vector([0.0, 0.0, 0.11], unit='m')\n", - " # sample_workflow[isis.DetectorBankOffset] = sc.vector([0.0, 0.0, 0.5], unit='m')\n", - " wf[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m')\n", + "# WARNING: Running this cell takes a couple of minutes when all sample runs are used.\n", "\n", - " # TODO Use direct beam!\n", - " # sample_workflow[DirectBeamFilename]\n", - " wf[DirectBeam] = None\n", + "# TODO Is this correct? We get very large times. Are we using the correct transmission\n", + "# runs above?\n", + "time_base = transmission.coords['datetime'].min()\n", "\n", - " return wf" + "periods = []\n", + "for period in range(4):\n", + " multi_run_sample_workflow[isis.mantidio.Period] = period\n", + " iofq = multi_run_sample_workflow.compute(IofQ[SampleRun])\n", + " iofq = iofq.bins.assign_coords(\n", + " time=iofq.bins.coords['pulse_time'].to(unit=time_base.unit) - time_base\n", + " ).drop_coords('wavelength') # Wavelength bounds get in the way\n", + " periods.append(iofq)" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "id": "39", "metadata": {}, + "outputs": [], "source": [ - "The sample runs are event data, so we can simply merge all runs to obtain a single event-mode $I(Q)$.\n", - "We use `sans.with_sample_runs` to obtain a workflow that merges the runs and computes the $I(Q)$:" + "dg = sc.DataGroup({f'period {period}': periods[period] for period in range(4)})\n", + "plopp.slicer(sc.values(dg.hist(pulse_time=10)), keep='Q')" ] }, { @@ -538,75 +559,124 @@ "metadata": {}, "outputs": [], "source": [ - "water_folder = Path('zoom_water_data')\n", - "sample_runs = [str(water_folder / f'ZOOM00038{run}.nxs') for run in range(439, 466, 2)]" + "periods[0].hist(pulse_time=100).plot(norm='log', scale={'Q': 'log'})" + ] + }, + { + "cell_type": "markdown", + "id": "41", + "metadata": {}, + "source": [ + "## Correction workflow\n", + "\n", + "In the previous section we have setup the workflow for the analyzer.\n", + "We also computed the transmission function there, but in production this will be done implicitly by running the entire workflow we will setup here.\n", + "We can combine this with the workflow for the polarizer to obtain the full correction workflow:" ] }, { "cell_type": "code", "execution_count": null, - "id": "41", + "id": "42", "metadata": {}, "outputs": [], "source": [ - "sample_workflow = make_sample_workflow()\n", - "# TODO Is this the correct empty beam run?\n", - "sample_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", - "multi_run_sample_workflow = sans.with_sample_runs(sample_workflow, runs=sample_runs)" + "supermirror_workflow = pol.SupermirrorWorkflow()\n", + "supermirror_workflow.visualize(pol.TransmissionFunction[pol.Polarizer])" ] }, { "cell_type": "markdown", - "id": "42", + "id": "43", "metadata": {}, "source": [ - "The sample workflow looks like this:" + "We will use a second-order polynomial supermirror efficiency function:" ] }, { "cell_type": "code", "execution_count": null, - "id": "43", + "id": "44", "metadata": {}, "outputs": [], "source": [ - "multi_run_sample_workflow.visualize(\n", - " IofQ[SampleRun], graph_attr={'rankdir': 'LR'}, compact=True\n", + "from ess.polarization.data import example_polarization_efficiency_table\n", + "\n", + "supermirror_workflow[pol.SupermirrorEfficiencyFunction[pol.Polarizer]] = (\n", + " pol.EfficiencyLookupTable.from_file(\n", + " example_polarization_efficiency_table(),\n", + " wavelength_colname='# X ',\n", + " efficiency_colname=' Y ',\n", + " )\n", ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "44", + "id": "45", "metadata": {}, "outputs": [], "source": [ - "periods = []\n", - "for period in range(4):\n", - " multi_run_sample_workflow[isis.mantidio.Period] = period\n", - " periods.append(multi_run_sample_workflow.compute(IofQ[SampleRun]))" + "workflow = pol.PolarizationAnalysisWorkflow(\n", + " polarizer_workflow=supermirror_workflow,\n", + " analyzer_workflow=he3_workflow,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "46", + "metadata": {}, + "source": [ + "For a single channel, the complete workflow looks as follows:" ] }, { "cell_type": "code", "execution_count": null, - "id": "45", + "id": "47", "metadata": {}, "outputs": [], "source": [ - "dg = sc.DataGroup({f'period {period}': periods[period] for period in range(4)})\n", - "plopp.slicer(sc.values(dg.hist(pulse_time=10)), keep='Q')" + "workflow.visualize(\n", + " pol.PolarizationCorrectedData[pol.Up, pol.Up], graph_attr={'rankdir': 'LR'}\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "48", + "metadata": {}, + "source": [ + "We can now set the reduced sample data for each measured channel:" ] }, { "cell_type": "code", "execution_count": null, - "id": "46", + "id": "49", "metadata": {}, "outputs": [], "source": [ - "periods[0].hist(pulse_time=100).plot(norm='log', scale={'Q': 'log'})" + "# TODO Which periods correspond to which spin channels?\n", + "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Up]] = periods[0]\n", + "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Down]] = periods[1]\n", + "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Up]] = periods[2]\n", + "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Down]] = periods[3]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO This needs debugging, currently we get all NaN values. Please check all TODOs\n", + "# above!\n", + "channels = workflow.compute(pol.PolarizationCorrectedData[pol.Up, pol.Up])\n", + "sc.DataGroup(vars(channels)).nanhist()" ] } ], From 4b410df9aeeabf30b79e973e06a38ba4cc09d3d4 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 21 Aug 2025 11:16:32 +0200 Subject: [PATCH 06/10] Update to use correct Zoom data files --- docs/user-guide/zoom.ipynb | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index a15170c..214c36b 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -39,9 +39,12 @@ "\n", "data_folder = Path('zoom_polarized_data')\n", "# Runs with analyzer at 4 different times\n", - "cell_runs = [str(data_folder / f'ZOOM00022{run}.nxs') for run in [710, 712, 714, 716]]\n", - "empty_run = data_folder / 'ZOOM00038238.nxs'\n", - "depolarized_run = data_folder / 'ZOOM00022718.nxs'\n", + "cell_runs = [\n", + " str(data_folder / f'ZOOM00038{run}.nxs')\n", + " for run in [423, 425, 427, 429, 431, 433, 435, 437]\n", + "]\n", + "empty_run = data_folder / 'ZOOM00038336.nxs'\n", + "depolarized_run = data_folder / 'ZOOM00038470.nxs'\n", "cell_runs = [*cell_runs, depolarized_run]" ] }, From 57a4a45f3bac3470869953b4917f4623dee50a02 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 21 Aug 2025 11:26:49 +0200 Subject: [PATCH 07/10] Use correct period->channel mapping --- docs/user-guide/zoom.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 214c36b..0ceb589 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -664,8 +664,8 @@ "source": [ "# TODO Which periods correspond to which spin channels?\n", "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Up]] = periods[0]\n", - "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Down]] = periods[1]\n", - "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Up]] = periods[2]\n", + "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Up]] = periods[1]\n", + "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Down]] = periods[2]\n", "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Down]] = periods[3]" ] }, From 178611f1f43d8bce0a664fd28b0af8603381e6f9 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 21 Aug 2025 12:59:25 +0200 Subject: [PATCH 08/10] Fix paths, compute Qxy, plots --- docs/user-guide/zoom.ipynb | 32 +++++++++++++++++++------------- src/ess/polarization/types.py | 6 ++++++ 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 0ceb589..0a1974f 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -37,14 +37,14 @@ "source": [ "from pathlib import Path\n", "\n", - "data_folder = Path('zoom_polarized_data')\n", + "data_folder = Path('zoom-data-2025-07')\n", "# Runs with analyzer at 4 different times\n", "cell_runs = [\n", - " str(data_folder / f'ZOOM00038{run}.nxs')\n", + " str(data_folder / '3He_unpolarized_in' / f'ZOOM00038{run}.nxs')\n", " for run in [423, 425, 427, 429, 431, 433, 435, 437]\n", "]\n", - "empty_run = data_folder / 'ZOOM00038336.nxs'\n", - "depolarized_run = data_folder / 'ZOOM00038470.nxs'\n", + "empty_run = data_folder / 'Background_files' / 'ZOOM00038336.nxs'\n", + "depolarized_run = data_folder / 'Background_files' / 'ZOOM00038470.nxs'\n", "cell_runs = [*cell_runs, depolarized_run]" ] }, @@ -439,6 +439,8 @@ " for key, value in params.items():\n", " wf[key] = value\n", " wf[QBins] = sc.geomspace(dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom')\n", + " wf[QxBins] = sc.linspace('Qx', start=-0.4, stop=0.4, num=101, unit='1/angstrom')\n", + " wf[QyBins] = sc.linspace('Qy', start=-0.4, stop=0.4, num=101, unit='1/angstrom')\n", " wf[NonBackgroundWavelengthRange] = sc.array(\n", " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", " )\n", @@ -474,7 +476,7 @@ "metadata": {}, "outputs": [], "source": [ - "water_folder = Path('zoom_water_data')\n", + "water_folder = data_folder / 'Sample_water_cell3'\n", "sample_runs = [str(water_folder / f'ZOOM00038{run}.nxs') for run in range(439, 466, 2)]" ] }, @@ -537,7 +539,7 @@ "periods = []\n", "for period in range(4):\n", " multi_run_sample_workflow[isis.mantidio.Period] = period\n", - " iofq = multi_run_sample_workflow.compute(IofQ[SampleRun])\n", + " iofq = multi_run_sample_workflow.compute(IofQxy[SampleRun])\n", " iofq = iofq.bins.assign_coords(\n", " time=iofq.bins.coords['pulse_time'].to(unit=time_base.unit) - time_base\n", " ).drop_coords('wavelength') # Wavelength bounds get in the way\n", @@ -551,8 +553,7 @@ "metadata": {}, "outputs": [], "source": [ - "dg = sc.DataGroup({f'period {period}': periods[period] for period in range(4)})\n", - "plopp.slicer(sc.values(dg.hist(pulse_time=10)), keep='Q')" + "periods[0]" ] }, { @@ -562,7 +563,10 @@ "metadata": {}, "outputs": [], "source": [ - "periods[0].hist(pulse_time=100).plot(norm='log', scale={'Q': 'log'})" + "dg = sc.DataGroup({f'period {period}': periods[period] for period in range(1)})\n", + "plopp.slicer(\n", + " sc.values(sc.concat(periods, 'period').hist(pulse_time=10)), keep=('Qx', 'Qy')\n", + ")" ] }, { @@ -662,7 +666,6 @@ "metadata": {}, "outputs": [], "source": [ - "# TODO Which periods correspond to which spin channels?\n", "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Up]] = periods[0]\n", "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Up]] = periods[1]\n", "workflow[pol.ReducedSampleDataBySpinChannel[pol.Up, pol.Down]] = periods[2]\n", @@ -676,10 +679,13 @@ "metadata": {}, "outputs": [], "source": [ - "# TODO This needs debugging, currently we get all NaN values. Please check all TODOs\n", - "# above!\n", "channels = workflow.compute(pol.PolarizationCorrectedData[pol.Up, pol.Up])\n", - "sc.DataGroup(vars(channels)).nanhist()" + "tiled = plopp.tiled(2, 2)\n", + "tiled[0, 0] = channels.upup.nanhist().plot()\n", + "tiled[0, 1] = channels.updown.nanhist().plot()\n", + "tiled[1, 0] = channels.downup.nanhist().plot()\n", + "tiled[1, 1] = channels.downdown.nanhist().plot()\n", + "tiled" ] } ], diff --git a/src/ess/polarization/types.py b/src/ess/polarization/types.py index 2ab50ae..55ec9a0 100644 --- a/src/ess/polarization/types.py +++ b/src/ess/polarization/types.py @@ -85,6 +85,12 @@ class PolarizationCorrectedData(Generic[PolarizerSpin, AnalyzerSpin]): downup: sc.DataArray downdown: sc.DataArray + def __post_init__(self): + self.upup.name = '(up, up)' + self.updown.name = '(up, down)' + self.downup.name = '(down, up)' + self.downdown.name = '(down, down)' + """The sum of polarization corrected data from all flipper state channels.""" TotalPolarizationCorrectedData = NewType( From b4ca601fdb97d7b108ea3f1a6db56315c1ae901c Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 21 Aug 2025 13:40:01 +0200 Subject: [PATCH 09/10] Add total corrected result plot --- docs/user-guide/zoom.ipynb | 70 ++++++++++++++++++++++++-------------- 1 file changed, 45 insertions(+), 25 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 0a1974f..2cdcd17 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -439,8 +439,8 @@ " for key, value in params.items():\n", " wf[key] = value\n", " wf[QBins] = sc.geomspace(dim='Q', start=0.004, stop=0.8, num=141, unit='1/angstrom')\n", - " wf[QxBins] = sc.linspace('Qx', start=-0.4, stop=0.4, num=101, unit='1/angstrom')\n", - " wf[QyBins] = sc.linspace('Qy', start=-0.4, stop=0.4, num=101, unit='1/angstrom')\n", + " wf[QxBins] = sc.linspace('Qx', start=-0.3, stop=0.3, num=101, unit='1/angstrom')\n", + " wf[QyBins] = sc.linspace('Qy', start=-0.3, stop=0.3, num=101, unit='1/angstrom')\n", " wf[NonBackgroundWavelengthRange] = sc.array(\n", " dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n", " )\n", @@ -488,7 +488,6 @@ "outputs": [], "source": [ "sample_workflow = make_sample_workflow()\n", - "# TODO Is this the correct empty beam run?\n", "sample_workflow[Filename[EmptyBeamRun]] = str(empty_run)\n", "multi_run_sample_workflow = sans.with_sample_runs(sample_workflow, runs=sample_runs)" ] @@ -553,17 +552,6 @@ "metadata": {}, "outputs": [], "source": [ - "periods[0]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "40", - "metadata": {}, - "outputs": [], - "source": [ - "dg = sc.DataGroup({f'period {period}': periods[period] for period in range(1)})\n", "plopp.slicer(\n", " sc.values(sc.concat(periods, 'period').hist(pulse_time=10)), keep=('Qx', 'Qy')\n", ")" @@ -571,7 +559,7 @@ }, { "cell_type": "markdown", - "id": "41", + "id": "40", "metadata": {}, "source": [ "## Correction workflow\n", @@ -584,7 +572,7 @@ { "cell_type": "code", "execution_count": null, - "id": "42", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -594,7 +582,7 @@ }, { "cell_type": "markdown", - "id": "43", + "id": "42", "metadata": {}, "source": [ "We will use a second-order polynomial supermirror efficiency function:" @@ -603,7 +591,7 @@ { "cell_type": "code", "execution_count": null, - "id": "44", + "id": "43", "metadata": {}, "outputs": [], "source": [ @@ -621,7 +609,7 @@ { "cell_type": "code", "execution_count": null, - "id": "45", + "id": "44", "metadata": {}, "outputs": [], "source": [ @@ -633,7 +621,7 @@ }, { "cell_type": "markdown", - "id": "46", + "id": "45", "metadata": {}, "source": [ "For a single channel, the complete workflow looks as follows:" @@ -642,7 +630,7 @@ { "cell_type": "code", "execution_count": null, - "id": "47", + "id": "46", "metadata": {}, "outputs": [], "source": [ @@ -653,7 +641,7 @@ }, { "cell_type": "markdown", - "id": "48", + "id": "47", "metadata": {}, "source": [ "We can now set the reduced sample data for each measured channel:" @@ -662,7 +650,7 @@ { "cell_type": "code", "execution_count": null, - "id": "49", + "id": "48", "metadata": {}, "outputs": [], "source": [ @@ -672,6 +660,14 @@ "workflow[pol.ReducedSampleDataBySpinChannel[pol.Down, pol.Down]] = periods[3]" ] }, + { + "cell_type": "markdown", + "id": "49", + "metadata": {}, + "source": [ + "The result for a single input channel (up-up) is:" + ] + }, { "cell_type": "code", "execution_count": null, @@ -687,11 +683,35 @@ "tiled[1, 1] = channels.downdown.nanhist().plot()\n", "tiled" ] + }, + { + "cell_type": "markdown", + "id": "51", + "metadata": {}, + "source": [ + "The combination of all four channels is:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52", + "metadata": {}, + "outputs": [], + "source": [ + "channels = workflow.compute(pol.TotalPolarizationCorrectedData)\n", + "tiled = plopp.tiled(2, 2)\n", + "tiled[0, 0] = channels.upup.nanhist().plot()\n", + "tiled[0, 1] = channels.updown.nanhist().plot()\n", + "tiled[1, 0] = channels.downup.nanhist().plot()\n", + "tiled[1, 1] = channels.downdown.nanhist().plot()\n", + "tiled" + ] } ], "metadata": { "kernelspec": { - "display_name": "mantid310", + "display_name": "dev311", "language": "python", "name": "python3" }, @@ -705,7 +725,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.11.13" }, "nbsphinx": { "execute": "never" From dfaf18b29df81b129a4f63fc75af6aa8b809fb68 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 21 Aug 2025 13:42:52 +0200 Subject: [PATCH 10/10] Update wavelength band --- docs/user-guide/zoom.ipynb | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/docs/user-guide/zoom.ipynb b/docs/user-guide/zoom.ipynb index 2cdcd17..7a89e42 100644 --- a/docs/user-guide/zoom.ipynb +++ b/docs/user-guide/zoom.ipynb @@ -252,9 +252,8 @@ "metadata": {}, "outputs": [], "source": [ - "# TODO Which wavelength bounds should be used?\n", - "wav_min = 2.2 * sc.Unit('angstrom')\n", - "wav_max = 2.8 * sc.Unit('angstrom')\n", + "wav_min = 1.75 * sc.Unit('angstrom')\n", + "wav_max = 16.5 * sc.Unit('angstrom')\n", "transmission_truncated = raw_transmission['wavelength', wav_min:wav_max]\n", "transmission_depolarized = transmission_truncated['time', -1].copy()\n", "transmission = transmission_truncated['time', :-1].copy()"