Fitting QENS data

Fitting QENS data#

Previously, some quasi-elastic neutron scattering (QENS) data has been simulated and reduced and can now be analysed with easyscience. Before the analysis can begin, it is necessary to load the experimental data and check that it looks reasonable. The data can be loaded with np.loadtxt as the data has been stored in a simple space-separated column file.

import numpy as np
from load import load

q, i, di = load('../4-reduction/qens_energy_transfer_unknown_quasi_elastic_3_pulse.dat')

With the data read in, we can produce a quick plot.

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.errorbar(q, i, di, fmt='.')
ax.set(xlabel='$\omega$/meV', ylabel='$I(\omega)$')
plt.show()
../_images/8c932da12e914f17fb23d4f406afcb8e98fa2a07f16dbefc7866b252b289666f.png

The QENS data covers a slightly broader range than we want to investigate. Therefore, we will crop the data so that we are only using points between -0.06 and 0.06 meV.

sel = (q > -0.06) & (q < 0.06)
q, i, di = q[sel], i[sel], di[sel]

We now want to consider the mathematical model to be used in the analysis. QENS data is typically analysed by fitting a Lorentzian function to the data, which has the functional form

(5)#\[ I(\omega) = \frac{A\gamma}{\pi\big[(\omega - \omega_0) ^ 2 + \gamma ^ 2\big]}, \]

where \(A\) is a scale factor, \(\gamma\) is the half-width at half-maximum, and \(\omega_0\) is the centre offset.

Exercise 1: implement a Lorentzian function#

Write a function that implements Eqn. (5), called lorentzian.

Solution:

Hide code cell content 
def lorentzian(omega: np.ndarray) -> np.ndarray:
    """
    A Lorentzian function.
    
    :omega: the energy transfer values to calculate over.

    :return: intensity values.
    """
    return A.value / np.pi * gamma.value / ((omega - omega_0.value) ** 2 + gamma.value ** 2)

The instrument has a finite resolution, and to take this into account, the Lorentzian function must be convolved with this resolution. In our case, the resolution is a Gaussian distribution that is centred at zero with width, \(\sigma\). In real experiments there are ways to measure \(\sigma\), but for this project, we will simply model it.

from scipy.stats import norm

def intensity(omega):
    """
    The convolution of a Gaussian and the Lorentzian.
    
    :omega: the energy transfer values to calculate over.

    :return: intensity values.
    """
    gauss = norm(0, sigma.value).pdf(omega)
    gauss /= gauss.sum()
    return np.convolve(lorentzian(omega), gauss, 'same')

This means that there are a total of four parameters in our model.

Exercise 2: create fitting parameters#

Create four Parameter objects, for \(A\), \(\gamma\), \(\omega_0\), and \(\sigma\).

Each should have an initial value and a uniform prior distribution, based on the values given in Table 2.

Table 2 Parameter values.#

Parameter

Initial Value

Min

Max

\(A\)

10

1

100

\(\gamma\)

8.0 × 10-3

1.0 × 10-4

1.0 × 10-2

\(\omega_0\)

1.0 × 10-3

0

2.0 × 10-3

\(\sigma\)

1.0 × 10-3

1.0 × 10-5

1.0 × 10-1

Solution:

Hide code cell content 
from easyscience.Objects.new_variable import Parameter

A = Parameter(name='A', value=10, fixed=False, min=0.01, max=100)
gamma = Parameter(name='gamma', value=8e-3, fixed=False, min=1e-4, max=10e-3)
omega_0 = Parameter(name='omega_0', value=0.001, fixed=False, min=0, max=0.002)
sigma = Parameter(name='sigma', value=0.001, fixed=False, min=1e-5, max=0.1)

It is now possible to compare our model at the initial estimates to the simulated data.

fig, ax = plt.subplots()
ax.errorbar(q, i, di, fmt='.')
ax.plot(q, intensity(q), '-k', zorder=10)
ax.set(xlabel='$\omega$/meV', ylabel='$I(\omega)$')
plt.show()
../_images/b169397df84a9a8b4711bb85fc1a1155578229c6301289c3c2fc0c9e1b794262.png

Exercise 3: find maximum likelihood estimates#

Using easyscience, obtain maximum likelihood estimates for the four parameters of the model from comparison with the data.

Solution:

Hide code cell content 
from easyscience.Objects.ObjectClasses import BaseObj
from easyscience.fitting import Fitter

params = BaseObj(name='params', A=A, gamma=gamma, omega_0=omega_0, sigma=sigma)
f = Fitter(params, intensity)

res = f.fit(x=q, y=i, weights=1/di)

We can then plot the model and the data together as before and print the values of the parameters along with their uncertainties.

fig, ax = plt.subplots()
ax.errorbar(q, i, di, fmt='.')
ax.plot(q, intensity(q), '-k', zorder=10)
ax.set(xlabel='$\omega$/meV', ylabel='$I(\omega)$')
plt.show()
../_images/222441c8627ac169fc6b6c5e9c57c57504366d0cdc51e862d1e832802310d333.png 
A, gamma, omega_0, sigma

(<Parameter 'A': 15.5713, bounds=[0.01:100.0]>,
 <Parameter 'gamma': 0.0071, bounds=[0.0001:0.01]>,
 <Parameter 'omega_0': 0.0008, bounds=[0.0:0.002]>,
 <Parameter 'sigma': 0.0000, bounds=[1e-05:0.1]>)
Contents