%matplotlib widget
Fitting SANS data¶
Some small angle neutron scattering (SANS) 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
def fetch_data(name: str) -> str:
"""
Fetch pre-prepared data from a remote source and return the path to the file.
"""
import pooch
return pooch.retrieve(
url=f'https://public.esss.dk/groups/scipp/dmsc-summer-school/2025/{name}',
known_hash=None,
)
filename = fetch_data('4-reduction/sans_iofq.dat')
Downloading data from 'https://public.esss.dk/groups/scipp/dmsc-summer-school/2025/4-reduction/sans_iofq.dat' to file '/home/runner/.cache/pooch/4b17da1f86c526e73f9c150933e50ad9-sans_iofq.dat'.
/home/runner/work/core/core/.pixi/envs/default/lib/python3.13/site-packages/requests/__init__.py:113: RequestsDependencyWarning: urllib3 (2.6.3) or chardet (7.0.1)/charset_normalizer (3.4.4) doesn't match a supported version! warnings.warn(
SHA256 hash of downloaded file: fb29bc63f793cac8cb2b9832afa159469b8525cd39a1e134efef6357be4be67b Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
def load(filename: str):
"""
Load data from a file. Filter out any NaN values.
"""
x, y, e = np.loadtxt(filename, unpack=True)
sel = np.isfinite(y)
return x[sel], y[sel], e[sel]
q, i, di = load(filename)
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(yscale='log', xlabel='$q$/Å^-1', ylabel='I(q)')
plt.show()
We now want to consider the mathematical model to be used in the analysis. There are SANS models a range of different systems, see for instance the models in SasView. However, initially, we will assume that our data has arisen from a spherical system.
The mathematical model for a sphere is
$$ I(q) = \frac{\text{scale}}{V} \bigg(\frac{3 V \Delta \rho [\sin{(qr)} - qr \cos{(qr)}]}{(qr)^3}\bigg)^2 + \text{bkg}, $$ (sphere)
where $\text{scale}$ is a scale factor, $V$ is the volume of the sphere, $\Delta \rho$ is the difference between the solvent and particle scattering length density, $r$ is the radius of the sphere, $\text{bkg}$ is a uniform background, and $q$ is the q-vector that the intensity is being calculated for.
Exercise 1: simplify the expression¶
The mathematical model described in Eqn. {eq}sphere has five parameters.
What simple mathematical simplification can be performed to reduce this to four?
Solution:
The volume of a sphere is related to the radius of the sphere as
$$ V = \frac{4}{3} \pi r^3. $$ (volume-sphere)
Therefore, the parameter $V$ can be replaced with Eqn. {eq}volume-sphere.
Exercise 2: write a function that computes for the form factor of a sphere¶
Four parameters is a suitable number for modelling.
Therefore, we should write a function that implements your reduced dimensionality version of Eqn. {eq}sphere.
def sphere(q):
"""
The function for the form factor of a sphere.
Parameters
----------
q:
q-vectors to calculate for.
Returns
-------
:
The modelled intensity.
"""
qr = q * r.value
V = 4 / 3 * np.pi * r.value**3
return (
scale.value
/ V
* (3 * V * delta_rho.value * (np.sin(qr) - qr * np.cos(qr)) / ((qr) ** 3)) ** 2
+ bkg.value
)
Exercise 3: create fitting parameters¶
Create four Parameter objects, for $\text{scale}$, $\Delta \rho$,
$r$, and $\text{bkg}$.
Each should have an initial value and a uniform prior distribution
based on the values given in {numref}sans-parameters, except for the
$\text{scale}$ which should be fixed to a value of 1.4 × 10-7.
| Parameter | Initial Value | Min | Max |
|---|---|---|---|
| $\text{scale}$ | 1.4 × 10-7 | N/A | N/A |
| $\Delta \rho$ | 3 | 0 | 10 |
| $r$ | 80 | 10 | 1000 |
| $\text{bkg}$ | 3.0 × 10-3 | 1.0 × 10-3 | 1.0 × 10-2 |
Solution:
from easyscience import Parameter
scale = Parameter(name='scale', value=1.4e-7, fixed=True)
delta_rho = Parameter(name='delta_rho', value=3, fixed=False, min=0, max=10)
r = Parameter(name='r', value=80, fixed=False, min=0, max=1000)
bkg = Parameter(name='bkg', value=0.01, fixed=False, min=0.001, 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, sphere(q), '-k', zorder=10)
ax.set(yscale='log', xlabel='$q$/Å^-1', ylabel='I(q)')
plt.show()
Exercise 4: find maximum likelihood estimates¶
Using easyscience, obtain maximum likelihood estimates for the four parameters of the model from comparison with the data.
Solution:
from easyscience import Fitter
from easyscience import ObjBase
params = ObjBase(name='params', delta_rho=delta_rho, r=r, bkg=bkg)
f = Fitter(params, sphere)
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, sphere(q), '-k', zorder=10)
ax.set(yscale='log', xlabel='$q$/Å^-1', ylabel='I(q)')
plt.show()
delta_rho, r, bkg
(<Parameter 'delta_rho': 3.5213 ± 0.0244, bounds=[0.0:10.0]>, <Parameter 'r': 90.6172 ± 0.2628, bounds=[0.0:1000.0]>, <Parameter 'bkg': 0.0186 ± 0.0004, bounds=[0.001:0.1]>)