Skip to content

Document how to use a component reading to create a new source #131

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.

Sign up for GitHub

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

Draft
nvaytet wants to merge 3 commits into
base: main
Choose a base branch
from
Draft

Document how to use a component reading to create a new source #131

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
193363e
add section docs to show how to use a component reading to create a m…
nvaytet Mar 5, 2026
d074703
fix typos
nvaytet Mar 5, 2026
89aa0fa
fix function taking in reading
nvaytet Mar 5, 2026
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
276 changes: 269 additions & 7 deletions docs/sources.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,12 @@
"source": [
"import numpy as np\n",
"import scipp as sc\n",
"import plopp as pp\n",
"import tof\n",
"\n",
"Hz = sc.Unit('Hz')"
"meter = sc.Unit('m')\n",
"Hz = sc.Unit('Hz')\n",
"deg = sc.Unit('deg')"
]
},
{
Expand All @@ -41,7 +44,12 @@
"\n",
"The default way of creating a source is to choose a facility name and a number of neutrons.\n",
"Each facility defines time and wavelength probability distributions for the neutrons, as well as a pulse frequency.\n",
"By default, a single pulse is generated."
"By default, a single pulse is generated.\n",
"\n",
"Here is a list of the currently supported sources:\n",
"\n",
"- `'ess'`: the standard ESS source profile, applicable for all ESS instruments\n",
"- `'ess-odin'`: a source specific to the ESS Odin instrument, sampled at the location where cold and thermal neutrons are combined (~2.35m away from the surface of the moderator)"
]
},
{
Expand Down Expand Up @@ -342,12 +350,266 @@
"id": "23",
"metadata": {},
"source": [
"## List of available sources/facilities\n",
"## Using a component reading to optimize source efficiency\n",
"\n",
"Here is a list of the currently supported sources, that can be selected via the `facility` argument:\n",
"In many cases, only a small fraction of neutrons make it through the entire chopper cascade,\n",
"while most are blocked by choppers.\n",
"\n",
"- `'ess'`: the standard ESS source profile, applicable for all ESS instruments\n",
"- `'ess-odin'`: a source specific to the ESS Odin instrument, sampled at the location where cold and thermal neutrons are combined (~2.35m away from the surface of the moderator)"
"Increasing the signal on a detector that lies after a chopper with very small openings can get very wasteful/expensive.\n",
"\n",
"Here we show a trick that can be used to only use the parts of the source pulse that make it to the detector,\n",
"and sample neutrons from that instead of the full birth time and wavelength distribution.\n",
"\n",
"**Note:** This section will only work for a model where the source only has **1 pulse!**\n",
"\n",
"We set up a model with a chopper that blocks most neutrons coming from the source:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "24",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"choppers = [\n",
" tof.Chopper(\n",
" frequency=70.0 * Hz,\n",
" open=sc.array(dims=['cutout'], values=[0.0], unit='deg'),\n",
" close=sc.array(dims=['cutout'], values=[150.0], unit='deg'),\n",
" phase=240.0 * deg,\n",
" distance=20.0 * meter,\n",
" name=\"chopchop\",\n",
" ),\n",
"]\n",
"\n",
"detectors = [\n",
" tof.Detector(distance=32.0 * meter, name='detector'),\n",
"]\n",
"\n",
"source = tof.Source(facility='ess', neutrons=500_000)\n",
"\n",
"model = tof.Model(source=source, components=choppers + detectors)\n",
"res = model.run()\n",
"\n",
"fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
"\n",
"res.plot(blocked_rays=5000, ax=ax[0])\n",
"_ = res['detector'].toa.plot(ax=ax[1])"
]
},
{
"cell_type": "markdown",
"id": "25",
"metadata": {},
"source": [
"As shown in the figure above, most neutrons from the source get blocked by the chopper.\n",
"\n",
"We can count the fraction of neutrons that make it through:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "26",
"metadata": {},
"outputs": [],
"source": [
"float(res['detector'].data.sum().value / source.data.sum().value) * 100"
]
},
{
"cell_type": "markdown",
"id": "27",
"metadata": {},
"source": [
"Only ~40% of the neutrons make it through.\n",
"\n",
"We now histogram the detector reading (this also works if we take the reading from any other component, such as the chopper)\n",
"in two dimension `(wavelength, birth_time)`, showing us where from the source the neutrons that arrive at the detector came from.\n",
"\n",
"We also show in the first panel a histogram with `toa` instead of `birth_time` to show what the data looks like at the detector."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28",
"metadata": {},
"outputs": [],
"source": [
"da = res['detector'].data\n",
"\n",
"p1 = da.hist(\n",
" wavelength=sc.scalar(0.025, unit='angstrom'), toa=sc.scalar(100, unit='us')\n",
")\n",
"p2 = da.hist(\n",
" wavelength=sc.scalar(0.025, unit='angstrom'), birth_time=sc.scalar(20, unit='us')\n",
")\n",
"\n",
"fig, ax = plt.subplots(1, 3, figsize=(16, 4))\n",
"\n",
"p1.plot(logc=True, title=\"Neutrons at the detector\", ax=ax[0])\n",
"p2.plot(logc=True, title=\"Neutrons at the source\", ax=ax[1])\n",
"source.data.hist(\n",
" wavelength=sc.scalar(0.025, unit='angstrom'), birth_time=sc.scalar(20, unit='us')\n",
").plot(\n",
" cmap='gray',\n",
" logc=True,\n",
" cbar=False,\n",
" title=\"Detected neutrons overlaid on original source\",\n",
" ax=ax[2],\n",
")\n",
"_ = p2.plot(logc=True, ax=ax[2])"
]
},
{
"cell_type": "markdown",
"id": "29",
"metadata": {},
"source": [
"Now if we simply sampled from the histogram of neutrons at the source,\n",
"we would get counts at the detector which are just a noisy as our original run above.\n",
"\n",
"So before we do this, we use a Gaussian kernel to smooth out the distribution:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "30",
"metadata": {},
"outputs": [],
"source": [
"from scipp.scipy.ndimage import gaussian_filter\n",
"\n",
"p = gaussian_filter(\n",
" p2,\n",
" sigma={\n",
" 'birth_time': sc.scalar(100, unit='us'),\n",
" 'wavelength': sc.scalar(0.05, unit='angstrom'),\n",
" },\n",
")\n",
"p.plot(logc=True, title=\"New source\")"
]
},
{
"cell_type": "markdown",
"id": "31",
"metadata": {},
"source": [
"We can now sample from this new distribution and create a new source with only 200K neutrons (instead of 500K above)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "32",
"metadata": {},
"outputs": [],
"source": [
"s2 = tof.Source.from_distribution(neutrons=200_000, p=p)\n",
"s2.data.hist(wavelength=300, birth_time=300).plot(logc=True)"
]
},
{
"cell_type": "markdown",
"id": "33",
"metadata": {},
"source": [
"Finally, we re-run the model using the new source:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "34",
"metadata": {},
"outputs": [],
"source": [
"m2 = tof.Model(source=s2, components=choppers + detectors)\n",
"r2 = m2.run()\n",
"\n",
"fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
"\n",
"r2.plot(blocked_rays=5000, ax=ax[0])\n",
"dt = sc.scalar(200, unit='us')\n",
"_ = pp.plot(\n",
" {\"original\": da.hist(toa=dt), \"new\": r2['detector'].data.hist(toa=dt)}, ax=ax[1]\n",
")"
]
},
{
"cell_type": "markdown",
"id": "35",
"metadata": {},
"source": [
"Now we can see (in the right panel) that we have obtained very similar results to the original model (same level of noise in the statistics) but using less than half the neutrons.\n",
"\n",
"We can also see in the left panel that very few neutrons get blocked by the chopper (grey lines).\n",
"\n",
"For convenice, we can package all this in a single function:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "36",
"metadata": {},
"outputs": [],
"source": [
"def source_from_reading(\n",
" reading,\n",
" neutrons: int,\n",
" time_binning: sc.Variable,\n",
" wavelength_binning: sc.Variable,\n",
" time_smooth_kernel: sc.Variable | None = None,\n",
" wavelength_smooth_kernel: sc.Variable | None = None,\n",
"):\n",
" from scipp.scipy.ndimage import gaussian_filter\n",
"\n",
" p = reading.data.hist(wavelength=wavelength_binning, birth_time=time_binning)\n",
" if time_smooth_kernel is None:\n",
" time_smooth_kernel = 2 * time_binning\n",
" if wavelength_smooth_kernel is None:\n",
" wavelength_smooth_kernel = 2 * wavelength_binning\n",
" p = gaussian_filter(\n",
" p,\n",
" sigma={\n",
" 'birth_time': time_smooth_kernel,\n",
" 'wavelength': wavelength_smooth_kernel,\n",
" },\n",
" )\n",
" return tof.Source.from_distribution(neutrons=neutrons, p=p)"
]
},
{
"cell_type": "markdown",
"id": "37",
"metadata": {},
"source": [
"Which we can then use as follows:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "38",
"metadata": {},
"outputs": [],
"source": [
"s3 = source_from_reading(\n",
" res['chopchop'],\n",
" neutrons=500_000,\n",
" time_binning=sc.scalar(40, unit='us'),\n",
" wavelength_binning=sc.scalar(0.025, unit='angstrom'),\n",
")\n",
"\n",
"m3 = tof.Model(source=s3, components=choppers + detectors)\n",
"r3 = m3.run()\n",
"r3['detector'].plot(bins=800)"
]
}
],
Expand All @@ -367,7 +629,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.7"
"version": "3.12.12"
}
},
"nbformat": 4,
Expand Down
7 changes: 3 additions & 4 deletions src/tof/component.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,15 @@

import plopp as pp
import scipp as sc

from .utils import Plot
from plopp.core.typing import FigureLike


@dataclass(frozen=True)
class ReadingField:
data: sc.DataArray
dim: str

def plot(self, bins: int = 300, **kwargs):
def plot(self, bins: int = 300, **kwargs) -> FigureLike:
by_pulse = sc.collapse(self.data, keep="event")
to_plot = {}
color = {}
Expand Down Expand Up @@ -120,7 +119,7 @@ def speed(self) -> ReadingField:
"""
return _make_reading_field(self.data, dim="speed")

def plot(self, bins: int = 300) -> Plot:
def plot(self, bins: int = 300) -> FigureLike:
"""
Plot both the toa and wavelength data side by side.

Expand Down
3 changes: 2 additions & 1 deletion src/tof/source.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
import numpy as np
import plopp as pp
import scipp as sc
from plopp.core.typing import FigureLike

from .component import ComponentReading
from .utils import wavelength_to_speed
Expand Down Expand Up @@ -454,7 +455,7 @@ def from_distribution(
def __len__(self) -> int:
return self.data.sizes["pulse"]

def plot(self, bins: int = 300) -> tuple:
def plot(self, bins: int = 300) -> FigureLike:
"""
Plot the pulses of the source.

Expand Down
Loading