ADR Suggestion ResolutionHandler

[henrikjacobsenfys]

General

The instrument resolution must be included when fitting QENS data, and often also for INS. We shall here only consider fitting 1d cuts of the data, with intensity as function of energy. The fit function is then the convolution of the model with the resolution, plus a background. The model is most commonly a sum of one or more Lorentzian functions, but can also include Gaussians, damped harmonic oscillators (DHO), and many other functions.

The resolution function is determined using measurements of either vanadium or the sample itself at low temperature. Either way, there are two common ways to determine it, and both must be supported. The first way (a) is to simply use the measurement directly, after cleaning it up a bit (detailed below). The convolution is then calculated numerically using, e.g., scipy.signal.fftconvolve. This leads to some potential errors when the model function is sharply peaked or if the data is not linearly spaced.

The second way is to fit the resolution to an analytical model. Unless the resolution is particularly ugly (which we do not expect it to be), this can be done using a sum of Gaussians and Lorentzians. In general, the convolution with the model can then be calculated numerically as above, with the same potential errors. However, the convolution of a Gaussian with a Gaussian is another Gaussian, the convolution of a Lorentzian with another Lorentzian is another Lorentzian, and the convolution of a Gaussian with a Lorentzian is a Voigt profile, which is also available, e.g., in scipy.special.voigt_profile. This can be used to improve accuracy and make the calculations faster.

Current implementation

Nothing is implemented yet

Proposed implementation

I propose a class, ResolutionHandler to handle this. It takes a SampleModel and ResolutionModel and produces the convolution of the two. The Analysis class can use this for the fitting, together with a BackgroundModel (to be defined elsewhere). Let's talk about the SampleModel and ResolutionModel

SampleModel

This holds the model, which consists of components, which can be Gaussian, Lorentzian, DHO, etc., or user defined. (Care must be taken to also allow 2d models often used in QENS, such as JumpDiffusion. This will probably be handled by a separate class, which can then generate a SampleModel for each probed Q).
It has methods to create and remove components.

ResolutionModel

In case (a), this holds the measured resolution. It will have methods to clean the data, such as smoothing and background subtraction.
In case (b), this will be a SampleModel that holds the fitted resolution.
It has methods to create and remove components.

ResolutionHandler

This will loop over the components in SampleModel and calculate their convolution with the ResolutionModel, using the analytical method where possible, and numerical methods otherwise. (If the analytical method is not possible, the calculation can be sped up by summing the SampleModel components first. Let's not worry about this unless the fitting is slow).

It will take x as input, i.e. the points where the convoluted model will be evaluated, as well as a SampleModel and ResolutionModel.
Optionally, it can take
upsampling: how many more points to evaluate on compared to the input (default: probably around 4) (perhaps also a num_points to let the user choose how many points instead)

It will need some methods in addition to the convolute method to combat numerical instabilities:

interpolate: If the data is not evenly spaced, the fft convolution method produces incorrect results. The model must therefore be evaluated at evenly spaced points and interpolated onto the input.

It will also need to check if the input is symmetric, and handle it properly.

[henrikjacobsenfys]

Two problems need to be carefully addressed:

1. Detailed Balancing: It is easier to create excitations than to absorb them at low temperature, and as a result, the intensities are different on the energy gain and energy loss side. This means that analytical convolutions give incorrect results when the peaks are broad.
2. Sharp peaks give numerical instabilities when using numerical convolution.

For now, I will stick with numerical convolutions, and be mindful of the instabilities. They can to some extent be overcome with the upsampling and interpolate methods

[henrikjacobsenfys]

I finally did some tests. Numerical convolution is definitely required at low temperature. At some point I'll look at this in terms of real units.

[henrikjacobsenfys]

Typical backscatterig (MIRACLESS):

Typical ToF (CSPEC):

Conclusion:
For backscattering, we can use analytical convolution with no problems, unless we use a dilution fridge.
For ToF, we have to be careful at low T. Sensible, overly cautious, checks can be made automatically.

[damskii9992]

I would propose some slight changes to this structure. Instead of a ResolutionModel and BackgroundModel class, I would implement a single InstrumentModel class with 2 attributes: resolution and background.

Now for the resolution attribute, instead of having a single ResolutionModel class, I would propose to let this resolution attribute either be an instance of a SampleModel or an instance of a Measurement class.
Since the Measurement class is designed to hold mesurement data, this class should be ideal for case a)
Then when you perform the convolution, you can simply query whether the resolution attribute is either a Measurement class instance, or a SampleModel class instance to determine whether the convolution should be numerical or analytic.

[henrikjacobsenfys]

All very sensible, thanks

[henrikjacobsenfys]

Update on the analytical vs numerical convolution: it's essentially impossible to define a region where it's OK to use the analytical convolution with DBF. It can therefore only be reliably used without DBF. Only deep blue in this color map is safe for convolution of Gaussian resolution with Lorentzian detailed balance