diff --git a/.readthedocs.yaml b/.readthedocs.yaml index e884b0d..18ca1dd 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -3,7 +3,7 @@ version: 2 build: os: "ubuntu-22.04" tools: - python: "3.9" + python: "3.10" sphinx: configuration: docs/conf.py diff --git a/docs/index.rst b/docs/index.rst index 5c1f414..a69bf0c 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -41,5 +41,6 @@ You can install from ``conda`` using sources components multiple-pulses + wfm-stitching api Release notes diff --git a/docs/wfm-stitching.ipynb b/docs/wfm-stitching.ipynb new file mode 100644 index 0000000..d887a6e --- /dev/null +++ b/docs/wfm-stitching.ipynb @@ -0,0 +1,618 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6b645603-12c0-4b3d-ab3c-529ad4f373e3", + "metadata": {}, + "source": [ + "# Stitching WFM data\n", + "\n", + "Wavelength-frame-multiplication (WFM) is a technique commonly used at long-pulse facilities to improve the resolution of the results measured at the neutron detectors.\n", + "See for example the article by [Schmakat et al. (2020)](https://www.sciencedirect.com/science/article/pii/S0168900220308640) for a description of how WFM works.\n", + "\n", + "In this notebook, we show how `tof` can be used to find the boundaries of the WFM frames, and apply a time correction to each frame,\n", + "in order to obtain more accurate wavelengths." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "899a245a-0ae7-45dd-b684-3b9d1c034cc6", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import scipp as sc\n", + "from scippneutron.conversion.graph.beamline import beamline\n", + "from scippneutron.conversion.graph.tof import elastic\n", + "import plopp as pp\n", + "import tof\n", + "\n", + "Hz = sc.Unit(\"Hz\")\n", + "deg = sc.Unit(\"deg\")\n", + "meter = sc.Unit(\"m\")" + ] + }, + { + "cell_type": "markdown", + "id": "3890d24c-93ab-4faa-b653-c9208b5eba23", + "metadata": {}, + "source": [ + "## Create a source pulse\n", + "\n", + "We first create a source with one pulse containing 1 million neutrons whose distribution follows the ESS time and wavelength profiles (both thermal and cold neutrons are included)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "872800a5-4cf7-424e-8d2d-0db81292456c", + "metadata": {}, + "outputs": [], + "source": [ + "source = tof.Source(facility=\"ess\", neutrons=1_000_000)\n", + "source.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "b5ffdab6-10cd-4719-8ff2-86e5e1d10571", + "metadata": {}, + "source": [ + "## Chopper set-up\n", + "\n", + "We create a list of choppers that will be included in our beamline.\n", + "In our case, we make two WFM choppers, and two frame-overlap choppers.\n", + "All choppers have 6 openings.\n", + "\n", + "Finally, we also add a pulse-overlap chopper with a single opening.\n", + "These choppers are copied after the [V20 ESS beamline at HZB](https://www.sciencedirect.com/science/article/pii/S0168900216309597)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34314bcc-978c-468c-b5a1-d84ce33e53d5", + "metadata": {}, + "outputs": [], + "source": [ + "choppers = [\n", + " tof.Chopper(\n", + " frequency=70.0 * Hz,\n", + " open=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[98.71, 155.49, 208.26, 257.32, 302.91, 345.3],\n", + " unit=\"deg\",\n", + " ),\n", + " close=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[109.7, 170.79, 227.56, 280.33, 329.37, 375.0],\n", + " unit=\"deg\",\n", + " ),\n", + " phase=47.10 * deg,\n", + " distance=6.6 * meter,\n", + " name=\"WFM1\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=70 * Hz,\n", + " open=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[80.04, 141.1, 197.88, 250.67, 299.73, 345.0],\n", + " unit=\"deg\",\n", + " ),\n", + " close=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[91.03, 156.4, 217.18, 269.97, 322.74, 375.0],\n", + " unit=\"deg\",\n", + " ),\n", + " phase=76.76 * deg,\n", + " distance=7.1 * meter,\n", + " name=\"WFM2\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=56 * Hz,\n", + " open=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[74.6, 139.6, 194.3, 245.3, 294.8, 347.2],\n", + " unit=\"deg\",\n", + " ),\n", + " close=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[95.2, 162.8, 216.1, 263.1, 310.5, 371.6],\n", + " unit=\"deg\",\n", + " ),\n", + " phase=62.40 * deg,\n", + " distance=8.8 * meter,\n", + " name=\"Frame-overlap 1\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=28 * Hz,\n", + " open=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[98.0, 154.0, 206.8, 254.0, 299.0, 344.65],\n", + " unit=\"deg\",\n", + " ),\n", + " close=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[134.6, 190.06, 237.01, 280.88, 323.56, 373.76],\n", + " unit=\"deg\",\n", + " ),\n", + " phase=12.27 * deg,\n", + " distance=15.9 * meter,\n", + " name=\"Frame-overlap 2\",\n", + " ),\n", + " tof.Chopper(\n", + " frequency=7 * Hz,\n", + " open=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[30.0],\n", + " unit=\"deg\",\n", + " ),\n", + " close=sc.array(\n", + " dims=[\"cutout\"],\n", + " values=[140.0],\n", + " unit=\"deg\",\n", + " ),\n", + " phase=0 * deg,\n", + " distance=22 * meter,\n", + " name=\"Pulse-overlap\",\n", + " ),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "bd6e1399-7f05-4d5c-83e3-7d249d9e0a61", + "metadata": {}, + "source": [ + "## Detector set-up\n", + "\n", + "We add a detector 32 meters from the source." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a377e75c-51a5-4994-af9c-92139d27bee3", + "metadata": {}, + "outputs": [], + "source": [ + "detectors = [\n", + " tof.Detector(distance=32.0 * meter, name=\"detector\"),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "fcb52d39-56d6-4cfb-8bcd-be9f45f4cc70", + "metadata": {}, + "source": [ + "## Find WFM frame edges\n", + "\n", + "To compute the frame edges, we use one of the openings in the WFM choppers at a time,\n", + "run the `tof` simulation, and find the min and max tof of the neutrons that make it through the chopper cascade." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20597dbd-5204-49fc-8867-9a9ba73fd8de", + "metadata": {}, + "outputs": [], + "source": [ + "nframes = len(choppers[0].open)\n", + "results = []\n", + "\n", + "print(\"Computing frames\")\n", + "\n", + "# Run the model with a single frame at a time\n", + "for i in range(nframes):\n", + " wfm_choppers = [\n", + " tof.Chopper(\n", + " frequency=choppers[j].frequency,\n", + " open=choppers[j].open[i : i + 1],\n", + " close=choppers[j].close[i : i + 1],\n", + " phase=choppers[j].phase,\n", + " distance=choppers[j].distance,\n", + " name=choppers[j].name,\n", + " )\n", + " for j in (0, 1)\n", + " ]\n", + " new_choppers = wfm_choppers + choppers[2:]\n", + " model = tof.Model(source=source, choppers=new_choppers, detectors=detectors)\n", + " res = model.run()\n", + " results.append(res)\n", + "\n", + "invalid_frames = True\n", + "fact = 0.01\n", + "while invalid_frames:\n", + " print(f\"Searching for bounds: threshold={fact:.2f}\")\n", + " frame_bounds = []\n", + " for res in results:\n", + " tofs = res.detectors[\"detector\"].tofs.data[\"visible\"][\"pulse:0\"].hist(tof=500)\n", + " tofs.coords[\"tof\"] = sc.midpoints(tofs.coords[\"tof\"])\n", + " # We need to filter out the outliers because some stray rays from other frames make it through\n", + " filtered = tofs[tofs.data > fact * tofs.data.max()].coords[\"tof\"]\n", + " frame_bounds.append((filtered.min(), filtered.max()))\n", + " if all(\n", + " frame_bounds[k][1] < frame_bounds[k + 1][0]\n", + " for k in range(len(frame_bounds) - 1)\n", + " ):\n", + " invalid_frames = False\n", + " else:\n", + " fact += 0.01" + ] + }, + { + "cell_type": "markdown", + "id": "866a02f3-6314-40b5-ac1c-cd8ddacee366", + "metadata": {}, + "source": [ + "The edges of the frames are the following" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed050525-7f88-4bb2-af2f-827d1d7afe92", + "metadata": {}, + "outputs": [], + "source": [ + "frame_bounds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5283e4c-699d-4669-a1ff-3f86a751cd68", + "metadata": {}, + "outputs": [], + "source": [ + "frames = sc.DataGroup(\n", + " {\n", + " \"time_min\": sc.concat([b[0] for b in frame_bounds], dim=\"frame\"),\n", + " \"time_max\": sc.concat([b[1] for b in frame_bounds], dim=\"frame\"),\n", + " }\n", + ")\n", + "frames" + ] + }, + { + "cell_type": "markdown", + "id": "3891b351-4c17-42c2-9c88-68facd82b2aa", + "metadata": {}, + "source": [ + "### Inspecting the frames\n", + "\n", + "As a consistency check, we can run the model with all of the chopper openings, and overlay the frame bounds,\n", + "to verify that there is no overlap between the frames." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d425e97a-247f-47ac-ad85-10b6cb1fea0e", + "metadata": {}, + "outputs": [], + "source": [ + "model = tof.Model(source=source, choppers=choppers, detectors=detectors)\n", + "res = model.run()\n", + "f = res.detectors[\"detector\"].tofs.plot(legend=False)\n", + "for i, bound in enumerate(frame_bounds):\n", + " col = f\"C{i + 1}\"\n", + " f.ax.axvspan(bound[0].value, bound[1].value, alpha=0.1, color=col)\n", + " f.ax.axvline(bound[0].value, color=col)\n", + " f.ax.axvline(bound[1].value, color=col)\n", + "f" + ] + }, + { + "cell_type": "markdown", + "id": "5fd905ab-f9b2-437f-9c3a-1bad599447d5", + "metadata": {}, + "source": [ + "### Time-distance diagram\n", + "\n", + "Another way of verifying the frames that were computed is to find the fastest and slowest neutron in each frame,\n", + "and propagate those through the choppers and show the time-distance diagram." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35f59c87-7435-4248-b8df-c05a1eda0c5e", + "metadata": {}, + "outputs": [], + "source": [ + "frame_min = []\n", + "frame_max = []\n", + "\n", + "for i in range(nframes):\n", + " print(\"frame\", i)\n", + " da = results[i][\"detector\"].data.squeeze()\n", + " da = da[~da.masks[\"blocked_by_others\"]]\n", + " ts = da.coords[\"tof\"]\n", + " sel = (ts > frames[\"time_min\"][i]) & (ts < frames[\"time_max\"][i])\n", + " filtered = da[sel]\n", + " frame_min.append(filtered[np.argmin(filtered.coords[\"tof\"])])\n", + " frame_max.append(filtered[np.argmax(filtered.coords[\"tof\"])])\n", + "\n", + "# Create a source by manually setting neutron birth times and wavelengths\n", + "birth_times = sc.concat(\n", + " [f.coords[\"time\"] for f in frame_min] + [f.coords[\"time\"] for f in frame_max],\n", + " dim=\"event\",\n", + ")\n", + "wavelengths = sc.concat(\n", + " [f.coords[\"wavelength\"] for f in frame_min]\n", + " + [f.coords[\"wavelength\"] for f in frame_max],\n", + " dim=\"event\",\n", + ")\n", + "source_min_max = tof.Source.from_neutrons(birth_times=birth_times, wavelengths=wavelengths)\n", + "source_min_max" + ] + }, + { + "cell_type": "markdown", + "id": "3ce82c32-2b0d-4e19-a17e-e280caa734c4", + "metadata": {}, + "source": [ + "We can see that the source has 12 neutrons, which is 2 per frame.\n", + "Re-running the model with those yields" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a797542b-130d-4571-8d3f-6b81fb4781df", + "metadata": {}, + "outputs": [], + "source": [ + "model = tof.Model(source=source_min_max, choppers=choppers, detectors=detectors)\n", + "model.run().plot()" + ] + }, + { + "cell_type": "markdown", + "id": "54b85305-c89b-4077-97c2-10c447793725", + "metadata": {}, + "source": [ + "What is interesting, and also a contract of WFM, is that the wavelength of the slowest neutron in one frame is very close to the wavelength of the fastest neutron in the next frame." + ] + }, + { + "cell_type": "markdown", + "id": "c808287b-e569-48c6-a964-ee19d06380ef", + "metadata": {}, + "source": [ + "## Stitching the data: computing a new time-of-flight\n", + "\n", + "Using WFM choppers allows us to re-define the burst time of the neutrons, and compute a more accurate wavelength.\n", + "\n", + "In the following, we use the boundaries of the frames to select neutrons in each frame,\n", + "and apply a correction to the time-of-flight of those neutrons which corresponds to the time when the WFM choppers are open.\n", + "\n", + "### Computing wavelengths from the naive time-of-flight\n", + "\n", + "We first begin by computing the neutron wavelengths as if there were no WFM choppers.\n", + "We take the distance from the source to the detector, and use the neutron arrival time at the detector to compute a speed and hence a wavelength." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "979a7fd7-3597-4522-ab41-9ddf6098c241", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract the tof data of the events that make it through to the detector\n", + "tofs = res[\"detector\"].tofs.data[\"visible\"][\"pulse:0\"].copy()\n", + "tofs.coords[\"source_position\"] = sc.vector([0.0, 0.0, 0.0], unit=\"m\")\n", + "tofs.coords[\"position\"] = sc.vector(\n", + " [0.0, 0.0, detectors[0].distance.value], unit=detectors[0].distance.unit\n", + ")\n", + "tofs" + ] + }, + { + "cell_type": "markdown", + "id": "1584ab73-9a08-43e5-9047-9f8ec15e56e4", + "metadata": {}, + "source": [ + "Converting the time-of-flight to wavelength is done using Scipp's `transform_coords`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dbb79512-efbc-4f9c-92e4-e1125f18b07f", + "metadata": {}, + "outputs": [], + "source": [ + "# Make a coordinate transformation graph to compute wavelength from tof\n", + "graph = {**beamline(scatter=False), **elastic(\"tof\")}\n", + "wav_naive = tofs.transform_coords(\"wavelength\", graph=graph)\n", + "wav_naive.hist(wavelength=300).plot()" + ] + }, + { + "cell_type": "markdown", + "id": "f72c4e94-2052-488d-90d3-113af1334287", + "metadata": {}, + "source": [ + "### Computing time-of-flight from WFM choppers to detector\n", + "\n", + "Instead of using the source as the departure point of the neutrons, we use the WFM choppers.\n", + "This means that the distance used for the flight is from the WFM choppers to the detector,\n", + "and the flight time is from when the choppers open to the arrival time at detector.\n", + "\n", + "We first get the times when the choppers are open (mid-point between open and close times for the 2 WFM choppers):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d792986-194a-4714-9e75-f7a038102f43", + "metadata": {}, + "outputs": [], + "source": [ + "times_wfm1 = choppers[0].open_close_times()\n", + "times_wfm2 = choppers[1].open_close_times()\n", + "\n", + "corrections = [\n", + " sc.concat(\n", + " [\n", + " times_wfm1[0][i + nframes], # open wfm1\n", + " times_wfm1[1][i + nframes], # close wfm1\n", + " times_wfm2[0][i + nframes], # open wfm2\n", + " times_wfm2[1][i + nframes], # close wfm2\n", + " ],\n", + " dim=\"x\",\n", + " ).mean()\n", + " for i in range(nframes)\n", + "]\n", + "\n", + "frames[\"time_correction\"] = sc.concat(corrections, dim=\"frame\")\n", + "frames" + ] + }, + { + "cell_type": "markdown", + "id": "8c7daac3-b38f-41e4-a90d-b9daf2d18a1b", + "metadata": {}, + "source": [ + "We apply the correction which effectively 'stitches' the data back together:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7da162f-bb0b-4bf9-812e-ba14b12ad1fd", + "metadata": {}, + "outputs": [], + "source": [ + "def stitch(\n", + " data: sc.DataArray,\n", + " frames: sc.DataGroup,\n", + " dim: str,\n", + ") -> sc.DataArray:\n", + " edges = sc.flatten(\n", + " sc.transpose(\n", + " sc.concat([frames[\"time_min\"], frames[\"time_max\"]], \"dummy\"),\n", + " dims=[\"frame\", \"dummy\"],\n", + " ),\n", + " to=dim,\n", + " )\n", + "\n", + " binned = data.bin({dim: edges})\n", + "\n", + " for i in range(frames.sizes[\"frame\"]):\n", + " binned[dim, i * 2].bins.coords[dim] -= frames[\"time_correction\"][\"frame\", i]\n", + "\n", + " binned.masks[\"frame_gaps\"] = (\n", + " sc.arange(dim, 2 * frames.sizes[\"frame\"] - 1) % 2\n", + " ).astype(bool)\n", + " binned.masks[\"frame_gaps\"].unit = None\n", + " return binned.bins.concat()\n", + "\n", + "\n", + "wfm_tofs = stitch(data=tofs, frames=frames, dim=\"tof\")" + ] + }, + { + "cell_type": "markdown", + "id": "8895fba5-c777-4190-ab02-9d7be7acbe8e", + "metadata": {}, + "source": [ + "Finally, we change the `source_position` to now be the mid-point between the WFM choppers:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "160412bd-6253-4ac3-87ff-058e9db2132f", + "metadata": {}, + "outputs": [], + "source": [ + "wfm_tofs.coords[\"source_position\"] = sc.vector(\n", + " [0.0, 0.0, 0.5 * (choppers[0].distance.value + choppers[1].distance.value)],\n", + " unit=choppers[0].distance.unit,\n", + ")\n", + "wfm_tofs" + ] + }, + { + "cell_type": "markdown", + "id": "c2bdf70a-a2c9-4d93-a17b-5191d3f8c784", + "metadata": {}, + "source": [ + "We can now compute wavelengths using the `transform_coords`, as before:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea7a25f0-4ce8-4ace-9bec-2cece7e4ee08", + "metadata": {}, + "outputs": [], + "source": [ + "wav_wfm = wfm_tofs.transform_coords(\"wavelength\", graph=graph)" + ] + }, + { + "cell_type": "markdown", + "id": "ea677910-7670-4588-991b-baa3fc54f0eb", + "metadata": {}, + "source": [ + "### Comparison between naive and WFM computations\n", + "\n", + "We compare the wavelengths computed using the naive approach, the WFM approach, and also with the true wavelengths of the neutrons:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "893a2d62-82d0-4bd2-acf2-688cd0d5fec9", + "metadata": {}, + "outputs": [], + "source": [ + "pp.plot(\n", + " {\n", + " \"naive\": wav_naive.hist(wavelength=300),\n", + " \"wfm\": wav_wfm.hist(wavelength=300),\n", + " \"truth\": res[\"detector\"]\n", + " .wavelengths.data[\"visible\"][\"pulse:0\"]\n", + " .hist(wavelength=300),\n", + " }\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "da6a161b-2051-4afa-a49a-7eb045356d01", + "metadata": {}, + "source": [ + "As we can see, the WFM approach vastly outperforms the naive approach." + ] + } + ], + "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/requirements/docs.in b/requirements/docs.in index cb2ee66..729f9d8 100644 --- a/requirements/docs.in +++ b/requirements/docs.in @@ -2,6 +2,7 @@ gitpython ipykernel ipywidgets nbsphinx +scippneutron sphinx sphinx-autodoc-typehints sphinx-book-theme diff --git a/requirements/docs.txt b/requirements/docs.txt index 4cbc32b..c1ee3fe 100644 --- a/requirements/docs.txt +++ b/requirements/docs.txt @@ -1,4 +1,4 @@ -# SHA1:132cbd83d2a0310697d4b6af67651e1fc07f84c9 +# SHA1:6fd4d9201f4bf5006252fb9dc1d498ef452c7015 # # This file is autogenerated by pip-compile-multi # To update, run: @@ -33,6 +33,10 @@ comm==0.2.2 # via # ipykernel # ipywidgets +contourpy==1.2.1 + # via matplotlib +cycler==0.12.1 + # via matplotlib debugpy==1.8.2 # via ipykernel decorator==5.1.1 @@ -50,10 +54,16 @@ executing==2.0.1 # via stack-data fastjsonschema==2.20.0 # via nbformat +fonttools==4.53.0 + # via matplotlib gitdb==4.0.11 # via gitpython gitpython==3.1.43 # via -r docs.in +h5py==3.11.0 + # via + # scippneutron + # scippnexus idna==3.7 # via requests imagesize==1.4.1 @@ -92,16 +102,24 @@ jupyterlab-pygments==0.3.0 # via nbconvert jupyterlab-widgets==3.0.11 # via ipywidgets +kiwisolver==1.4.5 + # via matplotlib markupsafe==2.1.5 # via # jinja2 # nbconvert +matplotlib==3.9.0 + # via + # mpltoolbox + # plopp matplotlib-inline==0.1.7 # via # ipykernel # ipython mistune==3.0.2 # via nbconvert +mpltoolbox==24.5.1 + # via scippneutron nbclient==0.10.0 # via nbconvert nbconvert==7.16.4 @@ -115,10 +133,21 @@ nbsphinx==0.9.4 # via -r docs.in nest-asyncio==1.6.0 # via ipykernel +numpy==2.0.0 + # via + # contourpy + # h5py + # matplotlib + # mpltoolbox + # scipp + # scippneutron + # scipy packaging==24.1 # via # ipykernel + # matplotlib # nbconvert + # pooch # pydata-sphinx-theme # sphinx pandocfilters==1.5.1 @@ -127,8 +156,16 @@ parso==0.8.4 # via jedi pexpect==4.9.0 # via ipython +pillow==10.3.0 + # via matplotlib platformdirs==4.2.2 - # via jupyter-core + # via + # jupyter-core + # pooch +plopp==24.6.0 + # via scippneutron +pooch==1.8.2 + # via scippneutron prompt-toolkit==3.0.47 # via ipython psutil==6.0.0 @@ -146,8 +183,13 @@ pygments==2.18.0 # nbconvert # pydata-sphinx-theme # sphinx +pyparsing==3.1.2 + # via matplotlib python-dateutil==2.9.0.post0 - # via jupyter-client + # via + # jupyter-client + # matplotlib + # scippnexus pyzmq==26.0.3 # via # ipykernel @@ -157,11 +199,25 @@ referencing==0.35.1 # jsonschema # jsonschema-specifications requests==2.32.3 - # via sphinx + # via + # pooch + # sphinx rpds-py==0.18.1 # via # jsonschema # referencing +scipp==24.6.0 + # via + # scippneutron + # scippnexus +scippneutron==24.6.0 + # via -r docs.in +scippnexus==24.6.0 + # via scippneutron +scipy==1.14.0 + # via + # scippneutron + # scippnexus six==1.16.0 # via # asttokens diff --git a/src/tof/source.py b/src/tof/source.py index a46fd0d..38896de 100644 --- a/src/tof/source.py +++ b/src/tof/source.py @@ -13,17 +13,17 @@ from . import facilities from .utils import Plot, wavelength_to_speed -TIME_UNIT = 'us' -WAV_UNIT = 'angstrom' +TIME_UNIT = "us" +WAV_UNIT = "angstrom" def _default_frequency(frequency: Union[None, sc.Variable], pulses: int) -> sc.Variable: if frequency is None: if pulses > 1: raise ValueError( - 'If pulses is greater than one, a frequency must be supplied.' + "If pulses is greater than one, a frequency must be supplied." ) - frequency = 1.0 * sc.Unit('Hz') + frequency = 1.0 * sc.Unit("Hz") return frequency @@ -71,8 +71,8 @@ def _make_pulses( wmax: Maximum neutron wavelength. """ - t_dim = 'time' - w_dim = 'wavelength' + t_dim = "time" + w_dim = "wavelength" p_time = _convert_coord(p_time, unit=TIME_UNIT, coord=t_dim) p_wav = _convert_coord(p_wav, unit=WAV_UNIT, coord=w_dim) @@ -85,8 +85,8 @@ def _make_pulses( if wmax is None: wmax = p_wav.coords[w_dim].max() - time_interpolator = interp1d(p_time, dim=t_dim, fill_value='extrapolate') - wav_interpolator = interp1d(p_wav, dim=w_dim, fill_value='extrapolate') + time_interpolator = interp1d(p_time, dim=t_dim, fill_value="extrapolate") + wav_interpolator = interp1d(p_wav, dim=w_dim, fill_value="extrapolate") x_time = sc.linspace( dim=t_dim, start=tmin.value, @@ -145,13 +145,13 @@ def _make_pulses( wavs.append(w[mask]) n += mask.sum() - dim = 'event' + dim = "event" birth_time = sc.array( dims=[dim], values=np.concatenate(times), unit=TIME_UNIT, - ).fold(dim=dim, sizes={'pulse': pulses, dim: neutrons}) + ( - sc.arange('pulse', pulses) / frequency + ).fold(dim=dim, sizes={"pulse": pulses, dim: neutrons}) + ( + sc.arange("pulse", pulses) / frequency ).to( unit=TIME_UNIT, copy=False ) @@ -160,12 +160,12 @@ def _make_pulses( dims=[dim], values=np.concatenate(wavs), unit=WAV_UNIT, - ).fold(dim=dim, sizes={'pulse': pulses, dim: neutrons}) + ).fold(dim=dim, sizes={"pulse": pulses, dim: neutrons}) speed = wavelength_to_speed(wavelength) return { - 'time': birth_time, - 'wavelength': wavelength, - 'speed': speed, + "time": birth_time, + "wavelength": wavelength, + "speed": speed, } @@ -221,11 +221,11 @@ def __init__( wmax=wmax, ) self.data = sc.DataArray( - data=sc.ones(sizes=pulse_params['time'].sizes, unit='counts'), + data=sc.ones(sizes=pulse_params["time"].sizes, unit="counts"), coords={ - 'time': pulse_params['time'], - 'wavelength': pulse_params['wavelength'], - 'speed': pulse_params['speed'], + "time": pulse_params["time"], + "wavelength": pulse_params["wavelength"], + "speed": pulse_params["speed"], }, ) @@ -257,7 +257,7 @@ def from_neutrons( source = cls(facility=None, neutrons=len(birth_times), pulses=pulses) source.frequency = _default_frequency(frequency, pulses) - birth_times = (sc.arange('pulse', pulses) / source.frequency).to( + birth_times = (sc.arange("pulse", pulses) / source.frequency).to( unit=TIME_UNIT, copy=False ) + birth_times.to(unit=TIME_UNIT, copy=False) wavelengths = sc.broadcast( @@ -265,11 +265,11 @@ def from_neutrons( ) source.data = sc.DataArray( - data=sc.ones(sizes=birth_times.sizes, unit='counts'), + data=sc.ones(sizes=birth_times.sizes, unit="counts"), coords={ - 'time': birth_times, - 'wavelength': wavelengths, - 'speed': wavelength_to_speed(wavelengths).to(unit='m/s', copy=False), + "time": birth_times, + "wavelength": wavelengths, + "speed": wavelength_to_speed(wavelengths).to(unit="m/s", copy=False), }, ) @@ -321,17 +321,17 @@ def from_distribution( sampling=sampling, ) source.data = sc.DataArray( - data=sc.ones(sizes=pulse_params['time'].sizes, unit='counts'), + data=sc.ones(sizes=pulse_params["time"].sizes, unit="counts"), coords={ - 'time': pulse_params['time'], - 'wavelength': pulse_params['wavelength'], - 'speed': pulse_params['speed'], + "time": pulse_params["time"], + "wavelength": pulse_params["wavelength"], + "speed": pulse_params["speed"], }, ) return source def __len__(self) -> int: - return self.data.sizes['pulse'] + return self.data.sizes["pulse"] def plot(self, bins: int = 300) -> tuple: """ @@ -343,7 +343,7 @@ def plot(self, bins: int = 300) -> tuple: Number of bins to use for histogramming the neutrons. """ fig, ax = plt.subplots(1, 2) - dim = (set(self.data.dims) - {'pulse'}).pop() + dim = (set(self.data.dims) - {"pulse"}).pop() collapsed = sc.collapse(self.data, keep=dim) pp.plot( {k: da.hist(time=bins) for k, da in collapsed.items()}, @@ -353,8 +353,7 @@ def plot(self, bins: int = 300) -> tuple: {k: da.hist(wavelength=bins) for k, da in collapsed.items()}, ax=ax[1], ) - size = fig.get_size_inches() - fig.set_size_inches(size[0] * 2, size[1]) + fig.set_size_inches(10, 4) fig.tight_layout() return Plot(fig=fig, ax=ax)