Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Notebook for dream #61

Merged
merged 3 commits into from
Nov 11, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
360 changes: 360 additions & 0 deletions docs/ess/dream.ipynb
Original file line number Diff line number Diff line change
@@ -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
}
7 changes: 7 additions & 0 deletions docs/ess/index.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
ESS instruments
***************

.. toctree::
:maxdepth: 2

dream
1 change: 1 addition & 0 deletions docs/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -42,5 +42,6 @@ You can install from ``conda`` using
components
multiple-pulses
wfm
ess/index
api
Release notes <https://github.com/scipp/tof/releases>
8 changes: 4 additions & 4 deletions src/tof/chopper.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Comment on lines +63 to +66
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why did you change this part...? Doesn't ruff complain about this...?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Because we use optional in the rest of the codebase and I didn't want to make it incompatible with older versions of python, as people outside of DMSC are also using it ;-)

direction: Direction = Clockwise,
name: str = "",
):
Expand Down