Basic fitting in EasyDynamics¶

We here show how to use EasyDynamics to fit a more advanced synthetic data set with multiple signals and a background, as well as an energy offset that could from from slight instrument misalignment.

The general procedure is to create an Experiment object to hold the data, a SampleModel to describe the model, and Analysis to fit the model to the data. The Analysis object now also contains an InstrumentModel, which describes the background and energy offset.

In [1]:

# Imports
import pooch

import easydynamics as edyn
import easydynamics.sample_model as sm

# Make the plots interactive
%matplotlib widget

# Imports import pooch import easydynamics as edyn import easydynamics.sample_model as sm # Make the plots interactive %matplotlib widget

As before, we createa an Experiment to hold the data, and rebin it.

In [2]:

# Load the data
experiment = edyn.Experiment(display_name='Tutorial')

file_path = pooch.retrieve(
    url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial0/docs/docs/tutorials/data/fake_advanced_data.hdf5',
    known_hash='ee7310249df71a312ebc219f3e16b8da4e9aa37d29df919bbcaa541a38e1c39f',
)

experiment.load_hdf5(filename=file_path)
experiment.rebin({'Q': 16, 'energy': 256})
experiment.plot_data(slicer=True)

# Load the data experiment = edyn.Experiment(display_name='Tutorial') file_path = pooch.retrieve( url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial0/docs/docs/tutorials/data/fake_advanced_data.hdf5', known_hash='ee7310249df71a312ebc219f3e16b8da4e9aa37d29df919bbcaa541a38e1c39f', ) experiment.load_hdf5(filename=file_path) experiment.rebin({'Q': 16, 'energy': 256}) experiment.plot_data(slicer=True)

Downloading data from 'https://github.com/easyscience/dynamics-lib/raw/refs/heads/tutorial0/docs/docs/tutorials/data/fake_advanced_data.hdf5' to file '/home/runner/.cache/pooch/b435b6d060ac2c8248a84531ddd80995-fake_advanced_data.hdf5'.

Out[2]:

This time, we see a sharp Gaussian near zero energy transfer, but with broad tails. There is also a pair of symmetrically placed peaks near +-1.5 meV. The background is not zero, and looks like it has a slope. We first focus on the three contributions to the signal: the Gaussian, its tails, which can be described by a Lorentzian, and the peaks near +- 1.5 meV, which can be described by a damped harmonic oscillator.

We create and place these in a ComponentCollection like this:

In [3]:

gaussian = sm.Gaussian(display_name='Gaussian', area=3, width=0.05)
lorentzian = sm.Lorentzian(display_name='Lorentzian', area=2, width=0.3)
dho = sm.DampedHarmonicOscillator(display_name='DHO', area=1.5, width=0.2, center=1.5)

collection = sm.ComponentCollection()
collection.append_component(gaussian)
collection.append_component(lorentzian)
collection.append_component(dho)

gaussian = sm.Gaussian(display_name='Gaussian', area=3, width=0.05) lorentzian = sm.Lorentzian(display_name='Lorentzian', area=2, width=0.3) dho = sm.DampedHarmonicOscillator(display_name='DHO', area=1.5, width=0.2, center=1.5) collection = sm.ComponentCollection() collection.append_component(gaussian) collection.append_component(lorentzian) collection.append_component(dho)

Note

Note that we do not give the Gaussian and Lorentzian a `center` value. In this case, the center is at zero (with a potential offset, see below), and fixed. This is typically what you want when describing elastic (the Gaussian) and quasielastic (the Lorentzian) signals.

We pass this `ComponentCollection` to our `SampleModel`, just like before when we passed only a single component.

In [4]:

model = sm.SampleModel(components=collection)

model = sm.SampleModel(components=collection)

Next, we want to handle the two instrument effects. The first is the small shift in energy: if you zoom in on the peak at the center, you'll see that it is in fact not centered at zero - it is offset. The offset applies to elastic, inelastic and quasielastic scattering. This is very typical in real experiments, and can be caused by slight misalignment or imperfections in the instrument. To handle this, we create an InstrumentModel, and give it an energy_offset.

Note

As noted in the previous tutorial, the energy_offset is always present, even when you do not explicitly set an InstrumentModel. This may give problems in two cases: 1) If all peaks are given a center value and allowed to fit it. In the simplest case of a single peak, this is problematic because there are now two parameters determining the center of the peak. The peak center/energy offset are equivalent in this case, and cannot be uniquely determined. 2) If there is no central peak, e.g. if only background scattering is present. In either case, the energy_offset should be fixed. You can do this in the InstrumentModel using InstrumentModel.fix_energy_offset() or Analysis using Analysis.fix_energy_offset().

In [5]:

instrument = sm.InstrumentModel(energy_offset=0.05)

instrument = sm.InstrumentModel(energy_offset=0.05)

Second, we want to handle the background. We do this with a BackgroundModel, which works the same as the SampleModel. We give it a Polynomial component:

💡 Tip

SampleModel, BackgroundModel and ResolutionModel (introduced later) all take components in several different ways: you can add a single component like in the previous tutorial, append components using the `append_component` method, or give a ComponentCollection like we do here.

In [6]:

background = sm.BackgroundModel(components=sm.Polynomial(coefficients=[1.2, 0.05]))

background = sm.BackgroundModel(components=sm.Polynomial(coefficients=[1.2, 0.05]))

We add the background to the instrument model:

💡 Tip

We could also have created the background model first, and then added it to the InstrumentModel when creating it, i.e. InstrumentModel(energy_offset = 0.05, background_model = background)

In [7]:

instrument.background_model = background

instrument.background_model = background

Now we create out Analysis object and give it the instrument model:

In [8]:

analysis = edyn.Analysis(
    experiment=experiment,
    sample_model=model,
    instrument_model=instrument,
)

analysis = edyn.Analysis( experiment=experiment, sample_model=model, instrument_model=instrument, )

It is always good to check if the model somewhat matches the data.

In [9]:

analysis.plot_data_and_model()

analysis.plot_data_and_model()

Out[9]:

As before, let us first fit a single $Q$ index and plot the data and model to see how it looks. We choose an arbitrary $Q$ and plot only that one

In [10]:

fit_result_independent_single_Q = analysis.fit(Q_index=5)
analysis.plot_data_and_model(Q_index=5)

fit_result_independent_single_Q = analysis.fit(Q_index=5) analysis.plot_data_and_model(Q_index=5)

Out[10]:

The fit looks very good. We can again get a list of the parameters for this fit by accesing the corresponding Analysis1d object. Note ethat the Gaussian and Lorentzian centers are both zero, but that the energy_offset is non-zero.

In [11]:

analysis.analysis_list[5].get_all_parameters()

analysis.analysis_list[5].get_all_parameters()

Out[11]:

[<Parameter 'Gaussian area': 3.5594 ± 0.0229 meV, bounds=[0.0:inf]>,
 <Parameter 'Gaussian center': 0.0000 meV (fixed), bounds=[-inf:inf]>,
 <Parameter 'Gaussian width': 0.0503 ± 0.0003 meV, bounds=[1e-10:inf]>,
 <Parameter 'Lorentzian area': 10.3617 ± 0.0724 meV, bounds=[0.0:inf]>,
 <Parameter 'Lorentzian center': 0.0000 meV (fixed), bounds=[-inf:inf]>,
 <Parameter 'Lorentzian width': 0.3947 ± 0.0049 meV, bounds=[1e-10:inf]>,
 <Parameter 'DHO area': 2.5554 ± 0.0369 meV, bounds=[0.0:inf]>,
 <Parameter 'DHO center': 1.4832 ± 0.0013 meV, bounds=[1e-10:inf]>,
 <Parameter 'DHO width': 0.0978 ± 0.0019 meV, bounds=[1e-10:inf]>,
 <Parameter 'energy_offset': 0.1007 ± 0.0002 meV, bounds=[-inf:inf]>,
 <Parameter 'Polynomial_c0': 1.2031 ± 0.0150, bounds=[-inf:inf]>,
 <Parameter 'Polynomial_c1': 0.1958 ± 0.0047, bounds=[-inf:inf]>]

Since the fit looked good, we can now fit all $Q$. We also plot the result, again using the slicer.

In [12]:

fit_result_all_Q = analysis.fit()
analysis.plot_data_and_model()

fit_result_all_Q = analysis.fit() analysis.plot_data_and_model()

Out[12]:

Let us plot a few parameters as a function of Q using the plot_parameters method.

In [13]:

analysis.plot_parameters(names=['Gaussian area', 'Lorentzian area', 'DHO area'])

analysis.plot_parameters(names=['Gaussian area', 'Lorentzian area', 'DHO area'])

Out[13]: