Skip to content

Add inelastic sample component that changes neutron wavelengths #124

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

Open
nvaytet wants to merge 31 commits into
base: main
Choose a base branch
from
Open

Add inelastic sample component that changes neutron wavelengths #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
9f527aa
start modifying the way model runs
nvaytet Feb 20, 2026
6f49224
update model to work with components
nvaytet Feb 20, 2026
b43a3e3
move more logic into components
nvaytet Feb 21, 2026
e23e753
fix imports
nvaytet Feb 21, 2026
50da6f8
model runs and results are plotting
nvaytet Feb 22, 2026
58f0b25
fix plotting of rays
nvaytet Feb 22, 2026
ddb0eb8
remove old reading file
nvaytet Feb 22, 2026
7a8101d
start adding inelastic sample
nvaytet Feb 22, 2026
1bf824c
fix deltae
nvaytet Feb 23, 2026
0933abb
plot wavelength color for each section between components
nvaytet Feb 23, 2026
4029e49
fix ray colors
nvaytet Feb 24, 2026
b2c6bed
add inelastic sample section to components notebooks
nvaytet Feb 24, 2026
047d685
cleanup
nvaytet Feb 24, 2026
87fefc4
fix plotting blocked rays
nvaytet Feb 24, 2026
a199784
start fixing tests
nvaytet Feb 25, 2026
acb156b
fix remaining tests
nvaytet Feb 25, 2026
d10610f
remove commented test
nvaytet Feb 25, 2026
0e1e99d
fix dashboard
nvaytet Feb 25, 2026
0d59007
add tests for inelastic sample
nvaytet Feb 25, 2026
cfde0c4
finish sample tests
nvaytet Feb 25, 2026
13ab32c
static analysis
nvaytet Feb 25, 2026
7201e0c
small style change
nvaytet Feb 25, 2026
c1aade2
use callable to compute deltaE for more flexibility, and prevent fina…
nvaytet Mar 4, 2026
0c141a2
add test for final energy positive
nvaytet Mar 4, 2026
a5d209c
add rays in one go so that they share the same colorbar
nvaytet Mar 4, 2026
e1c4fc1
fix tests
nvaytet Mar 4, 2026
504ad6e
Merge pull request #130 from scipp/fix-cbar-zoom
nvaytet Mar 5, 2026
2435559
change function to return final energy
nvaytet Mar 5, 2026
83e770b
update tests
nvaytet Mar 5, 2026
c585f81
format noteookb
nvaytet Mar 5, 2026
b54f955
update text in notebook
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
5 changes: 4 additions & 1 deletion docs/api-reference/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,17 @@

Chopper
ChopperReading
Component
ComponentReading
Dashboard
Detector
DetectorReading
InelasticSample
Model
ReadingField
Result
Source
SourceParameters
SourceReading
```

## Top-level functions
Expand All @@ -39,5 +41,6 @@
:template: module-template.rst
:recursive:

facilities
utils
```
220 changes: 202 additions & 18 deletions docs/components.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,10 @@
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import scipp as sc\n",
"import plopp as pp\n",
"import tof\n",
"\n",
"meter = sc.Unit('m')\n",
Expand Down Expand Up @@ -412,7 +415,193 @@
"id": "33",
"metadata": {},
"source": [
"## Loading from a JSON file\n",
"## Inelastic sample\n",
"\n",
"Placing an `InelasticSample` in the instrument will change the energy of the incoming neutrons by a $\\Delta E$ defined in a function.\n",
"That function takes in the incident energy and returns the final energy after the energy transfer took place.\n",
"\n",
"To give an equal chance for a random $\\Delta E$ between -0.2 and 0.2 meV, we can use Numpy's uniform random sampling:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "34",
"metadata": {},
"outputs": [],
"source": [
"rng = np.random.default_rng(seed=83)\n",
"\n",
"\n",
"def uniform_deltae(e_i):\n",
Copy link
Member Author

@nvaytet nvaytet Mar 5, 2026

Choose a reason for hiding this comment

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

I"m wondering if it would make more sense to have a function take in the initial energy and return the final energy instead of the energy transfer?

It could look something like:

def apply_energy_transfer(e_i):
    de = sc.array(dims=e_i.dims, values=rng.uniform(-0.2, 0.2, size=e_i.shape), unit='meV')
    return e_i.to(unit='meV') - de

?

I initially went with returning DeltaE so that the user writing the function would not have to worry about handling unit conversions and unphysical final energies. But the unit conversions are not too much code to add, and we can still fix final energies inside the package code after the final energy is returned by the function.

This mechanism here would also allow the user to fix the negative energies in the way they want: either throw them away or give them a zero energy.

@bingli621 any opinions?

Copy link
Member Author

@nvaytet nvaytet Mar 5, 2026

Choose a reason for hiding this comment

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

In the end, I went with the function that returns final energy, I think it makes more sense.

Copy link
Contributor

@bingli621 bingli621 Mar 5, 2026

Choose a reason for hiding this comment

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

from the simulation's perspective, having Ef will make the calculation much easier. but the user should be able to change the energy transfer as well.

uniform_sampler = lambda size: rng.uniform(-0.2, 0.2, size=size)

def apply_energy_transfer(e_i, energy_transfer_sampler= uniform_sampler):
de = sc.array(dims=e_i.dims, values=sampler(size=e_i.shape), unit='meV')
return e_i.to(unit='meV') - de

Copy link
Member Author

@nvaytet nvaytet Mar 5, 2026

Choose a reason for hiding this comment

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

Not sure I understand.
The user supplied the whole function, so they can put in there what they want.

The question was: should the user supply a function that outputs energy transfer, or should it output final_energy?

In the first case, computing the final_energy from the initial_energy and the energy_transfer would be done internally in the package.

I think the second case make more sense; a function that would output the energy_transfer feels unusual/strange.

Copy link
Contributor

@bingli621 bingli621 Mar 5, 2026

Choose a reason for hiding this comment

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

At T-REX and CSPEC, one sample can be measured with different Ei. In this sense, I would imagine user should supply a function describing the energy transfer, which is the property of the sample, regardless of what Ei is used to measure it. And the Sample component should stay the same when we change the condition how it is measured. I hope this what you are asking...

" # Uniform sampling between -0.2 and 0.2 meV\n",
" de = sc.array(\n",
" dims=e_i.dims, values=rng.uniform(-0.2, 0.2, size=e_i.shape), unit='meV'\n",
" )\n",
" return e_i.to(unit='meV', copy=False) - de\n",
"\n",
"\n",
"sample = tof.InelasticSample(\n",
" distance=28.0 * meter, name=\"sample\", delta_e=uniform_deltae\n",
")\n",
"sample"
]
},
{
"cell_type": "markdown",
"id": "35",
"metadata": {},
"source": [
"We then make a single fast-rotating chopper with one small opening,\n",
"to select a narrow wavelength range at every rotation.\n",
"\n",
"We also add a monitor before the sample, and a detector after the sample so we can follow the changes in energies/wavelengths."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "36",
"metadata": {},
"outputs": [],
"source": [
"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=[1.0], unit='deg'),\n",
" phase=0.0 * deg,\n",
" distance=20.0 * meter,\n",
" name=\"fastchopper\",\n",
" ),\n",
"]\n",
"\n",
"detectors = [\n",
" tof.Detector(distance=26.0 * meter, name='monitor'),\n",
" tof.Detector(distance=32.0 * meter, name='detector'),\n",
"]\n",
"\n",
"source = tof.Source(facility='ess', neutrons=5_000_000)\n",
"\n",
"model = tof.Model(source=source, components=choppers + detectors + [sample])\n",
"res = model.run()\n",
"\n",
"\n",
"fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))\n",
"\n",
"dw = sc.scalar(0.1, unit='angstrom')\n",
"pp.plot(\n",
" {\n",
" 'monitor': res['monitor'].data.hist(wavelength=dw),\n",
" 'detector': res['detector'].data.hist(wavelength=dw),\n",
" },\n",
" title=\"With inelastic sample\",\n",
" xmin=4,\n",
" xmax=20,\n",
" ymin=-20,\n",
" ymax=400,\n",
" ax=ax[1],\n",
")\n",
"\n",
"res.plot(visible_rays=10000, ax=ax[0])"
]
},
{
"cell_type": "markdown",
"id": "37",
"metadata": {},
"source": [
"### Non-uniform energy-transfer distributions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "38",
"metadata": {},
"outputs": [],
"source": [
"# Sample 1: double-peak at min and max\n",
"def double_peak(e_i):\n",
" # Either -0.2 or 0.2 meV\n",
" de = sc.array(\n",
" dims=e_i.dims, values=rng.choice([-0.2, 0.2], size=e_i.shape), unit='meV'\n",
" )\n",
" return e_i.to(unit='meV', copy=False) - de\n",
"\n",
"\n",
"sample1 = tof.InelasticSample(\n",
" distance=28.0 * meter,\n",
" name=\"sample\",\n",
" delta_e=double_peak,\n",
")\n",
"\n",
"\n",
"# Sample 2: normal distribution\n",
"def normal_deltae(e_i):\n",
" de = sc.array(\n",
" dims=e_i.dims, values=rng.normal(scale=0.05, size=e_i.shape), unit='meV'\n",
" )\n",
" return e_i.to(unit='meV', copy=False) - de\n",
"\n",
"\n",
"sample2 = tof.InelasticSample(\n",
" distance=28.0 * meter,\n",
" name=\"sample\",\n",
" delta_e=normal_deltae,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "39",
"metadata": {},
"outputs": [],
"source": [
"model1 = tof.Model(source=source, components=choppers + detectors + [sample1])\n",
"model2 = tof.Model(source=source, components=choppers + detectors + [sample2])\n",
"\n",
"res1 = model1.run()\n",
"res2 = model2.run()\n",
"\n",
"fig, ax = plt.subplots(2, 2, figsize=(12, 9))\n",
"\n",
"res1.plot(visible_rays=10000, ax=ax[0, 0], title=\"Sample 1\")\n",
"pp.plot(\n",
" {\n",
" 'monitor': res1['monitor'].data.hist(wavelength=dw),\n",
" 'detector': res1['detector'].data.hist(wavelength=dw),\n",
" },\n",
" title=\"Sample 1\",\n",
" xmin=4,\n",
" xmax=20,\n",
" ymin=-20,\n",
" ymax=400,\n",
" ax=ax[0, 1],\n",
")\n",
"\n",
"res2.plot(visible_rays=10000, ax=ax[1, 0], title=\"Sample 2\")\n",
"_ = pp.plot(\n",
" {\n",
" 'monitor': res2['monitor'].data.hist(wavelength=dw),\n",
" 'detector': res2['detector'].data.hist(wavelength=dw),\n",
" },\n",
" title=\"Sample 2\",\n",
" xmin=4,\n",
" xmax=20,\n",
" ymin=-20,\n",
" ymax=400,\n",
" ax=ax[1, 1],\n",
")"
]
},
{
"cell_type": "markdown",
"id": "40",
"metadata": {},
"source": [
"## Loading components from a JSON file\n",
"\n",
"It is also possible to load components from a JSON file,\n",
"which can be very useful to quickly load a pre-configured instrument.\n",
Expand All @@ -423,19 +612,14 @@
{
"cell_type": "code",
"execution_count": null,
"id": "34",
"id": "41",
"metadata": {},
"outputs": [],
"source": [
"import json\n",
"\n",
"params = {\n",
" \"source\": {\n",
" \"type\": \"source\",\n",
" \"facility\": \"ess\",\n",
" \"neutrons\": 1e6,\n",
" \"pulses\": 1\n",
" },\n",
" \"source\": {\"type\": \"source\", \"facility\": \"ess\", \"neutrons\": 1e6, \"pulses\": 1},\n",
" \"chopper1\": {\n",
" \"type\": \"chopper\",\n",
" \"frequency\": {\"value\": 56.0, \"unit\": \"Hz\"},\n",
Expand Down Expand Up @@ -470,7 +654,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "35",
"id": "42",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -479,7 +663,7 @@
},
{
"cell_type": "markdown",
"id": "36",
"id": "43",
"metadata": {},
"source": [
"We now use the `tof.Model.from_json()` method to load our instrument:"
Expand All @@ -488,7 +672,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "37",
"id": "44",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -498,7 +682,7 @@
},
{
"cell_type": "markdown",
"id": "38",
"id": "45",
"metadata": {},
"source": [
"We can see that all components have been read in correctly.\n",
Expand All @@ -508,7 +692,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "39",
"id": "46",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -517,7 +701,7 @@
},
{
"cell_type": "markdown",
"id": "40",
"id": "47",
"metadata": {},
"source": [
"### Modifying the source\n",
Expand All @@ -531,7 +715,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "41",
"id": "48",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -542,7 +726,7 @@
},
{
"cell_type": "markdown",
"id": "42",
"id": "49",
"metadata": {},
"source": [
"## Saving to JSON\n",
Expand All @@ -553,7 +737,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "43",
"id": "50",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -577,7 +761,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.7"
"version": "3.12.12"
}
},
"nbformat": 4,
Expand Down
5 changes: 3 additions & 2 deletions src/tof/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,13 @@
submodules=['facilities'],
submod_attrs={
'chopper': ['AntiClockwise', 'Chopper', 'ChopperReading', 'Clockwise'],
'component': ['Component', 'ComponentReading', 'ReadingField'],
'dashboard': ['Dashboard'],
'detector': ['Detector', 'DetectorReading'],
'inelastic': ['InelasticSample', 'InelasticSampleReading'],
'model': ['Model'],
'reading': ['ComponentReading', 'ReadingField'],
'result': ['Result'],
'source': ['Source', 'SourceParameters'],
'source': ['Source', 'SourceReading'],
},
)

Expand Down
Loading
Loading