Add inelastic sample component that changes neutron wavelengths

[nvaytet]

This PR adds a new type of component: InelasticSample.
It modifies the neutron energies according to a user-defined function for energy shift.

It is used as follows:

# Sample 1: double-peak at min and max
def double_peak(e_i):
    # Either -0.2 or 0.2 meV
    de = sc.array(dims=e_i.dims, values=rng.choice([-0.2, 0.2], size=e_i.shape), unit='meV')
    return e_i.to(unit='meV', copy=False) - de

sample1 = tof.InelasticSample(
    distance=28.0 * meter,
    name="sample",
    delta_e=double_peak,
)

# Sample 2: normal distribution
def normal_deltae(e_i):
    de = sc.array(dims=e_i.dims, values=rng.normal(scale=0.05, size=e_i.shape), unit='meV')
    return e_i.to(unit='meV', copy=False) - de

sample2 = tof.InelasticSample(
    distance=28.0 * meter,
    name="sample",
    delta_e=normal_deltae,
)

tof.Model(source=source, components=choppers + detectors + [sample1])

Internal changes

Quite a bit of the internals had to be refactored, because of the way the chopper and detector logic was hard-coded inside the Model (especially the run function).
Instead of taking in choppers and detectors, the Model now accepts a flat list of components.
The effects of each component are applied onto the neutrons in succession, and the logic for those effects has been moved inside the components themselves.

I think it makes for a much cleaner structure in the model. Logic for plotting was also moved into the respective components.

Fixes #117

[nvaytet]

@bingli621 please take a look (and also try it out!)

[nvaytet]

This is where we need to also update the dashboard to support samples. This will be done in a follow up PR.

[nvaytet]

Comments from in-person discussion:

• Use a callable instead of a data array with probability distribution to apply the sample effects on the neutrons (more flexible)
• Prevent neutrons from becoming unphysical (large DeltaE for large wavelengths could mean the final energy is less than zero)

[nvaytet]

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?

[nvaytet]

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

[bingli621]

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

[nvaytet]

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.

[bingli621]

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...

[nvaytet]

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?

[nvaytet]

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

[bingli621]

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.

[nvaytet]

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.

[nvaytet]

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.

[nvaytet]

@bingli621 if you have time to take another look, it would be great :-)

[bingli621]

I will try it out, I have a great use case in mind, but it needs a bit more work on the code I'm working on, will update asap.