diff --git a/docs/ess/dream.ipynb b/docs/ess/dream.ipynb new file mode 100644 index 0000000..1f44217 --- /dev/null +++ b/docs/ess/dream.ipynb @@ -0,0 +1,360 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "2b54cd42-dfc2-4ad7-b9af-b165d5709888", + "metadata": {}, + "source": [ + "# DREAM in WFM mode\n", + "\n", + "This is a simulation of the DREAM chopper cascade in WFM mode.\n", + "We also show how one can convert the neutron arrival times at the detector to wavelength." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69e77a24-0c1b-49dd-855b-0c5b5bbb936e", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import scipp as sc\n", + "import plopp as pp\n", + "import tof\n", + "\n", + "Hz = sc.Unit('Hz')\n", + "deg = sc.Unit('deg')\n", + "meter = sc.Unit('m')\n", + "AA = sc.Unit('angstrom')" + ] + }, + { + "cell_type": "markdown", + "id": "72a1ba68-dae3-4d61-8313-f8801bb69bb6", + "metadata": {}, + "source": [ + "## Create a source\n", + "\n", + "We first create an ESS source with 2 pulses containing 500,000 neutrons each." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "659a8cff-06ab-43f2-9566-c1679e693915", + "metadata": {}, + "outputs": [], + "source": [ + "source = tof.Source(facility='ess', neutrons=500_000, pulses=2)\n", + "source.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "305d1424-637c-4d4a-9d1a-8c4e28102691", + "metadata": {}, + "source": [ + "## Component set-up\n", + "\n", + "## Choppers\n", + "\n", + "The DREAM chopper cascade consists of:\n", + "\n", + "- Two counter-rotating pulse-shaping choppers (PSC) that are very close to each other, located ~6m from the source\n", + "- An overlap chopper placed right after the two PSCs\n", + "- A band control chopper\n", + "- A T0 chopper" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b01dbf4-9f05-4eea-82b4-89c310dbd340", + "metadata": {}, + "outputs": [], + "source": [ + "choppers = [\n", + " tof.Chopper(\n", + " frequency=14 * Hz,\n", + " direction=tof.AntiClockwise,\n", + " centers=sc.array(\n", + " dims=['cutout'],\n", + " values=[0, 72, 86.4, 115.2, 172.8, 273.6, 288.0, 302.4],\n", + " unit='deg',\n", + " ),\n", + " widths=sc.array(\n", + " dims=['cutout'],\n", + " values=[2.46, 3.02, 3.27, 3.27, 5.02, 3.93, 3.93, 2.46],\n", + " unit='deg',\n", + " ),\n", + " phase=(286 - 180) * deg,\n", + " distance=6.145 * meter,\n", + " name=\"PSC1\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=14 * Hz,\n", + " direction=tof.Clockwise,\n", + " centers=sc.array(\n", + " dims=['cutout'],\n", + " values=[0, 28.8, 57.6, 144, 158.4, 216, 259.2, 316.8],\n", + " unit='deg',\n", + " ),\n", + " widths=sc.array(\n", + " dims=['cutout'],\n", + " values=[2.46, 3.60, 3.60, 3.23, 3.27, 3.77, 3.94, 2.62],\n", + " unit='deg',\n", + " ),\n", + " phase=236 * deg,\n", + " distance=6.155 * meter,\n", + " name=\"PSC2\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=14 * Hz,\n", + " direction=tof.AntiClockwise,\n", + " centers=sc.array(\n", + " dims=['cutout'],\n", + " values=[0.0],\n", + " unit='deg',\n", + " ),\n", + " widths=sc.array(\n", + " dims=['cutout'],\n", + " values=[27.6],\n", + " unit='deg',\n", + " ),\n", + " phase=(297 - 180 - 90) * deg,\n", + " distance=6.174 * meter,\n", + " name=\"OC\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=112 * Hz,\n", + " direction=tof.AntiClockwise,\n", + " centers=sc.array(\n", + " dims=['cutout'],\n", + " values=[0.0, 180.0],\n", + " unit='deg',\n", + " ),\n", + " widths=sc.array(\n", + " dims=['cutout'],\n", + " values=[73.75, 73.75],\n", + " unit='deg',\n", + " ),\n", + " phase=(240 - 180) * deg,\n", + " distance=9.78 * meter,\n", + " name=\"BC\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=28 * Hz,\n", + " direction=tof.AntiClockwise,\n", + " centers=sc.array(\n", + " dims=['cutout'],\n", + " values=[0.0],\n", + " unit='deg',\n", + " ),\n", + " widths=sc.array(\n", + " dims=['cutout'],\n", + " values=[314.9],\n", + " unit='deg',\n", + " ),\n", + " phase=(280 - 180) * deg,\n", + " distance=13.05 * meter,\n", + " name=\"T0\",\n", + " ),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "ab93843a-0da5-414d-a215-efe3fa78d569", + "metadata": {}, + "source": [ + "### Detector banks and monitors\n", + "\n", + "DREAM has 5 detector banks: the Mantle, two End-caps, a High-resolution detector and a SANS detector.\n", + "\n", + "For each detector bank, we use a single mean distance (in practice, one could have a different distance for each pixel)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a720e32-5fef-41ce-9afe-7415a539a779", + "metadata": {}, + "outputs": [], + "source": [ + "sample_position = 76.55 * meter\n", + "\n", + "detectors = [\n", + " tof.Detector(distance=sample_position + 1.125 * meter, name='mantle'),\n", + " tof.Detector(distance=sample_position + 1.125 * meter, name='end-cap'),\n", + " tof.Detector(distance=sample_position + 2.5 * meter, name='high-resolution'),\n", + " tof.Detector(distance=sample_position + 2.5 * meter, name='sans'),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "741cd7ae-075b-4a04-b71a-f9c977fab9a8", + "metadata": {}, + "source": [ + "## Run the simulation\n", + "\n", + "We propagate our pulse of neutrons through the chopper cascade and inspect the results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4de7426f-90be-45da-89ac-5de1a2d3bb6f", + "metadata": {}, + "outputs": [], + "source": [ + "model = tof.Model(source=source, choppers=choppers, detectors=detectors)\n", + "results = model.run()\n", + "results.plot(blocked_rays=5000)" + ] + }, + { + "cell_type": "markdown", + "id": "ad7f2d08-a50f-4a19-b65f-5010173bc3b5", + "metadata": {}, + "source": [ + "## Wavelength as a function of time-of-arrival\n", + "\n", + "### Plotting wavelength vs time-of-arrival\n", + "\n", + "Since we know the true wavelength of our neutrons,\n", + "as well as the time at which the neutrons arrive at the detector\n", + "(coordinate named `toa` in the detector reading),\n", + "we can plot an image of the wavelengths as a function of time-of-arrival:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6fdf820-fb80-43e8-995f-e43f049d503b", + "metadata": {}, + "outputs": [], + "source": [ + "events = sc.DataGroup()\n", + "for key, da in results.detectors.items():\n", + " bank = da.data.flatten(to='event')\n", + " events[key] = bank[~bank.masks['blocked_by_others']]\n", + "\n", + "# Histogram and plot\n", + "events['mantle'].hist(wavelength=500, toa=500).plot(norm='log', grid=True)" + ] + }, + { + "cell_type": "markdown", + "id": "2ffa0a00-cdc7-492d-a0b2-191325980f43", + "metadata": {}, + "source": [ + "### Defining a conversion from `toa` to `wavelength`\n", + "\n", + "The image above shows that there is a pretty tight correlation between time-of-arrival and wavelength.\n", + "\n", + "We compute the mean wavelength inside a given `toa` bin to define a relation between `toa` and `wavelength`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "afc57e9d-9766-450e-aa60-d2a86fb6c9da", + "metadata": {}, + "outputs": [], + "source": [ + "binned = events.bin(tof=500)\n", + "\n", + "# Weighted mean of wavelength inside each bin\n", + "mu = sc.DataGroup(\n", + " {\n", + " key: (da.bins.data * da.bins.coords['wavelength']).bins.sum() / da.bins.sum()\n", + " for key, da in binned.items()\n", + " }\n", + ")\n", + "\n", + "mu.plot(grid=True)" + ] + }, + { + "cell_type": "markdown", + "id": "97112112-efd3-47b6-bd7d-2c045ef660fe", + "metadata": {}, + "source": [ + "## Computing wavelengths\n", + "\n", + "We set up an interpolator that will compute wavelengths given an array of `toas`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd01c024-dff6-4c4a-b733-7ee26acacdf9", + "metadata": {}, + "outputs": [], + "source": [ + "from scipp.scipy.interpolate import interp1d\n", + "\n", + "wavelengths = sc.DataGroup()\n", + "\n", + "for key in mu:\n", + " # Set up interpolator\n", + " y = mu[key].copy()\n", + " y.coords['tof'] = sc.midpoints(y.coords['tof'])\n", + " f = interp1d(y, 'tof', bounds_error=False)\n", + "\n", + " # Compute wavelengths\n", + " wavs = f(events[key].coords['tof'].rename_dims(event='tof'))\n", + " wavelengths[key] = sc.DataArray(\n", + " data=sc.ones(sizes=wavs.sizes, unit='counts'), coords={'wavelength': wavs.data}\n", + " ).rename_dims(tof='event')\n", + "\n", + "wavelengths" + ] + }, + { + "cell_type": "markdown", + "id": "35aff0f0-a1a1-4247-9f40-e75f0a73ba8d", + "metadata": {}, + "source": [ + "We can now compare our computed wavelengths to the true wavelengths of the neutrons." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5256b1e7-b0b2-464c-889c-0cb14da400b5", + "metadata": {}, + "outputs": [], + "source": [ + "pp.plot(\n", + " {\n", + " 'wfm': wavelengths['mantle'].hist(wavelength=300),\n", + " 'original': events['mantle'].hist(wavelength=300),\n", + " }\n", + ")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/ess/index.rst b/docs/ess/index.rst new file mode 100644 index 0000000..6cd0de3 --- /dev/null +++ b/docs/ess/index.rst @@ -0,0 +1,7 @@ +ESS instruments +*************** + +.. toctree:: + :maxdepth: 2 + + dream diff --git a/docs/index.rst b/docs/index.rst index 522c2fa..b2fbb44 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -42,5 +42,6 @@ You can install from ``conda`` using components multiple-pulses wfm + ess/index api Release notes diff --git a/src/tof/chopper.py b/src/tof/chopper.py index 7e36844..395c22b 100644 --- a/src/tof/chopper.py +++ b/src/tof/chopper.py @@ -60,10 +60,10 @@ def __init__( frequency: sc.Variable, distance: sc.Variable, phase: sc.Variable, - open: sc.Variable | None = None, - close: sc.Variable | None = None, - centers: sc.Variable | None = None, - widths: sc.Variable | None = None, + open: Optional[sc.Variable] = None, + close: Optional[sc.Variable] = None, + centers: Optional[sc.Variable] = None, + widths: Optional[sc.Variable] = None, direction: Direction = Clockwise, name: str = "", ):