Skip to content

Commit

Permalink
Merge pull request #61 from scipp/notebook-for-dream
Browse files Browse the repository at this point in the history 
Notebook for dream
  • Loading branch information
@nvaytet
nvaytet authored Nov 11, 2024
2 parents 7d1b716 + 671dbeb commit f93c00d
Show file tree
Hide file tree
Showing 4 changed files with 372 additions and 4 deletions.
360 changes: 360 additions & 0 deletions docs/ess/dream.ipynb
View file
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
View file
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
View file
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
View file
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,
direction: Direction = Clockwise,
name: str = "",
):
Expand Down

0 comments on commit f93c00d

Please sign in to comment.