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

Merged
nvaytet merged 56 commits into from
Mar 11, 2026
Merged

Add inelastic sample component that changes neutron wavelengths #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
56 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
cfb8093
add section on dealing with unphysical DeltaE
nvaytet Mar 6, 2026
7011a02
prettier printing
nvaytet Mar 6, 2026
c85fcb3
NaN negative final energies instead of thresholding
nvaytet Mar 6, 2026
00248de
delta_e -> func
nvaytet Mar 9, 2026
3f17914
get the time_limit from neutron toas instead of passing it to all com…
nvaytet Mar 10, 2026
f6620c7
remove sample plot method
nvaytet Mar 10, 2026
da7261c
missing space
nvaytet Mar 10, 2026
65425dd
remove tof property in component reading
nvaytet Mar 10, 2026
29769fa
fix return type
nvaytet Mar 10, 2026
e2d2de4
replace _ with space
nvaytet Mar 10, 2026
857c747
replace _ with space again
nvaytet Mar 10, 2026
046f6c2
fix test name
nvaytet Mar 10, 2026
542abe6
use assign instead of manually making shallow copies
nvaytet Mar 10, 2026
6537cce
fix docstring for Model.add() method
nvaytet Mar 10, 2026
bca6503
fix Result.data to also add samples to the data group
nvaytet Mar 10, 2026
54187fc
skip inelastic sample in json (de)serialization
nvaytet Mar 10, 2026
bb66427
component base class fixes
nvaytet Mar 10, 2026
10af7f1
make kind an abstract property of components
nvaytet Mar 10, 2026
3e49e2f
require exact match to extract components from model and result
nvaytet Mar 10, 2026
7d460df
add sample __eq__
nvaytet Mar 10, 2026
9f00518
add sample eq test
nvaytet Mar 10, 2026
0d2290e
use assign methods instead of manually shallow-copying
nvaytet Mar 10, 2026
11749a8
add edge cases tests with multiple samples and sample between choppers
nvaytet Mar 10, 2026
84512e3
Apply automatic formatting
pre-commit-ci-lite[bot] Mar 10, 2026
ce486d0
remove unused imports
nvaytet Mar 10, 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
```
216 changes: 199 additions & 17 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,191 @@
"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...

Copy link
Member Author

@nvaytet nvaytet Mar 6, 2026

Choose a reason for hiding this comment

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

I think it is a good way to deal with the unphysical energy transfer.

Note that I am still applying a fix to the final energies here, but the function gives the opportunity to the user to fix it the way they want.

Is it safer to remove that check to avoid doing something silently that the user did not expect?

Copy link
Member Author

@nvaytet nvaytet Mar 6, 2026

Choose a reason for hiding this comment

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

I now added a section in the docs that shows how to manually deal with unphysical neutrons.

Copy link
Contributor

@bingli621 bingli621 Mar 6, 2026

Choose a reason for hiding this comment

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

The documentation helps! I have two comments. The first is, tof.InelasticSample takes a keyword argument "delta_e", but this actually is a function that takes Ei as input and outputs Ef. This is what confused me. The logic is correct here, just the argument name is a bit misleading. The second comment is, when encountering the case with unphysical energy transfer, by replacing the negative Ef with a small ( 10^-30) positive value, we are basically disregarding them, because neutrons with such small Ef will pile up at a TOA value that is out of the time range of interests. In this sense, I suppose both ways to handle the unphysical Ef should work.

Copy link
Member Author

@nvaytet nvaytet Mar 6, 2026

Choose a reason for hiding this comment

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

neutrons with such small Ef will pile up at a TOA value that is out of the time range of interests

I think the 1e-30 leads to an issue: the wavelengths are so large that it ends up being problematic to plot results (histograms and time-distance diagrams).
So I think by default I will NaN the numbers instead of thresholding at 1e-30.

Copy link
Member Author

@nvaytet nvaytet Mar 6, 2026

Choose a reason for hiding this comment

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

but this actually is a function that takes Ei as input and outputs Ef. This is what confused me

Ok let's find a better name for the argument.

" # 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(distance=28.0 * meter, name=\"sample\", func=uniform_deltae)\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",
" func=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",
" func=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 +610,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 +652,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "35",
"id": "42",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -479,7 +661,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 +670,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "37",
"id": "44",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -498,7 +680,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 +690,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "39",
"id": "46",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -517,7 +699,7 @@
},
{
"cell_type": "markdown",
"id": "40",
"id": "47",
"metadata": {},
"source": [
"### Modifying the source\n",
Expand All @@ -531,7 +713,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "41",
"id": "48",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -542,7 +724,7 @@
},
{
"cell_type": "markdown",
"id": "42",
"id": "49",
"metadata": {},
"source": [
"## Saving to JSON\n",
Expand All @@ -553,7 +735,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "43",
"id": "50",
"metadata": {},
"outputs": [],
"source": [
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