Reduction of TOFTOF data
Introduction
This notebook is based on the paper by D. Noferini et al., “Role of the doping level in localized proton motions in acceptor-doped barium zirconate proton conductors”, PCCP 20, 13697 (2018). https://doi.org/10.1039/C7CP07340B
We kindly thank the authors for allowing us to use their data.
In this work, the authors have used Quasi Elastic Neutron Scattering (QENS), among other techniques, to investigate the microscopic details of the proton motions in two proton conducting, acceptor-doped, perovskites, namely BaZr\(_{0.9}\)In\(_{0.1}\)O\(_{2.95}\) and BaZr\(_{0.8}\)In\(_{0.2}\)O\(_{2.9}\). Both samples were measured at the time-of-flight spectrometer TOFTOF at MLZ using an incident wavelength of 2.5 \(\unicode{x212B}\), giving a Gaussian-shaped resolution function with a full width at half maximum (FWHM) of approximately 450 \(\mu\)eV and a \(Q\)-range between 0.36 and 4.7 \(\unicode{x212B}^{-1}\). However, as some of the detectors were contaminated by Bragg reflections, only a set of 9 Q-cuts were retained in the analysis: \(Q = 0.6, 0.8, 1.0, 1.2, 1.7, 1.9, 2.3, 3.5, 3.9\) \(\unicode{x212B}^{-1}\).
This experimental setup provides information about the fast (in the order of picoseconds) local dynamics of the protons in the sample, i.e., how the proton jumps between adjacent oxygens or how the O-H group reorients. This is done by fitting the measured incoherent dynamic structure factor, \(S(Q, \omega)\), as a sum of an elastic (represented by a Dirac delta function) and a quasielastic (represented by a Lorentzian function) contributions, from which the elastic incoherent structure factor (EISF) is extracted. The latter can then be compared to different theoretical models (e.g., jumps between two equivalent sites vs jumps over four sites) in order to determine the microscopic details of the initial proton diffusion mechanism at the origin of the long-range diffusion (observable on a larger time scale of the order of a few nanoseconds, and therefore requiring sub-\(\mu\)eV resolution to measure it) and hence influencing the protonic conductivity of these systems. See the article for more details.
The two samples have been measured at several temperatures (200, 350, 450 and 540 K) in order also to determine the activation energy of this process.
In this notebook we will analyse only the data of one sample (\(BaZr_{0.8}In_{0.2}O_{3}H_{0.2}\)) at 350, 450 and 540 K, starting from the fully corrected \(S(Q,\omega)\) data derived from the raw counts after applying the following reduction steps: - normalization of the detectors using a vanadium standard - suppression of noisy or blind detectors - subtraction of the scattering from an empty sample container - correction of neutron energy-dependent efficiency of the neutron counters
During the analysis, the fit was done in a limited exchanged energy range: [-5, 0.7] meV, in order to avoid the elastic multiple scattering artefacts coming from scattering of the sample environment in the Stokes part of the spectra.
Explanations to preparation of fitting
Reading input files
Each sample file corresponds to a temperature. They contain reduced data at several \(q\)-values. There are four files (named In20_2p5_###K_tt_cr.txt) containing \(S(Q,\hbar\omega)\) at T = 200, 350, 450, and 540 K, and a fifth one containing the signal of a vanadium reference sample serving as a measurement of the resolution function.
The files are written in the INX format (F. Rieutord, INX - Program for time-of-flight data reduction, ILL Technical Report 90RI17T, 1990). For each \(Q\) we have a block consisting of a header of 4 lines, followed by 3 columns of data corresponding to \(\hbar\omega\) (in meV), \(S(Q,\hbar\omega)\) for that particular \(Q\) value, and the associated error.
The header is organized as follows:
- Line 1: 8 integers.The first integer indicates the number \(N\) of following lines that correspond to the present block, and the last integer is the number of points in the spectrum, i.e., \(N-3\).
Line 2: String containing the experiment title
- Line 3: 5 floats and 1 integer.The first one corresponds to the block index (or alternatively it can contain the value of \(Q\) or \(2\theta\) in the case the file contains \(S(2\theta, \hbar\omega)\) instead of \(S(Q,\hbar\omega)\)).The second and third values correspond to the incident energy (in meV) and wavevector (in \(\unicode{x212B}^{-1}\)), respectively.The fourth float number is the temperature.Here, the remaining float and integer values do not contain relevant information.
- Line 4: 3 floats.Not used here, although the second number can contain the bin width (normally in microseconds) of the time channels used during the time-of-flight acquisition.
Expression of fitting model
The spectrum of a vanadium standard was used as a resolution function \(R(Q, \hbar \omega)\).
\[S_{meas}(Q, \hbar \omega) = S(Q, \hbar \omega) \times [1 + n(\hbar \omega)] \hbar \omega \beta \circledast R(Q, \hbar \omega ),\]
where $ S(Q, \hbar `:nbsphinx-math:omega`) = a_D(Q) \delta`(:nbsphinx-math:hbar :nbsphinx-math:omega`) + a_L(Q) L(Q, \hbar `:nbsphinx-math:omega`) + bkg(Q)$, \(n(\hbar \omega) = \frac{1}{\exp(\frac{\hbar \omega}{k_B T}) - 1}\) is the Bose factor, and \(\beta = \frac{1}{k_B T}\) the reciprocal of the thermodynamic temperature. \(\delta\) and \(L\) are a Dirac peak and a Lorentzian, respectively. The background is constant, only Q-dependent. \(\circledast\) is the convolution symbol.
To avoid elastic multiple scattering artefacts from the environment present in the Stokes part of the spectra, the fit range was limited to \(-5\leq \hbar \omega \leq 0.7\) meV. We further focus on the peak and use a fitting range of \(-1\leq \hbar \omega \leq 0.7\) meV at temperatures of 350, 450 amd 540 K.
In order to investigate the geometry of the observed proton dynamics, we analyze the Q-dependence of the elastic and quasielastic scattering intensities. Comparing \(a_D(Q)\) and \(a_L(Q)\) normalized for their sum, with the expressions expected within the framework of a model for a jump diffusion over 2 or 4 equivalent sites located on a circle with radius \(r\): jump length, mean residence time (Noferini et al., “Localized proton motions in acceptor-doped Barium Zirconates” J. Phys. Chem. C 121, 7088 (2017))
Calculating variances for \(\frac{a_D}{a_D+a_L}\) and \(\frac{a_L}{a_D+a_L}\)
Combinations of refined parameters will be plotted after the fit in order to compare with the above publication: \(\frac{a_D}{a_D + a_L}\) and \(\frac{a_L}{a_D + a_L}\). Therefore, we have to calculate the related uncertainty. Straightforward calculations of error propagation, assuming uncorrelated parameters, lead to the following expressions:
\[\sigma^2_{\frac{a_D}{a_D + a_L}} = \frac{a_L^2}{(a_D + a_L)^4}\sigma^2_{a_D} + \frac{a_D^2}{(a_D + a_L)^4}\sigma^2_{a_L}\]
\[\sigma^2_{\frac{a_L}{a_D + a_L}} = \frac{a_L^2}{(a_D + a_L)^4}\sigma^2_{a_D} + \frac{a_D^2}{(a_D + a_L)^4}\sigma^2_{a_L}\]
Importing libraries
[1]:
import numpy as np
import matplotlib.pyplot as plt
import QENSmodels
import scipy
import lmfit
from lmfit import Model, CompositeModel
from scipy.interpolate import interp1d
import pandas as pd
%matplotlib widget
q-values and range of \(\hbar \omega\)
[2]:
q_tof = [0.6, 0.8, 1.0, 1.2, 1.7, 1.9, 2.3, 3.5, 3.9] # in 1/Angstrom
# Range of energy in meV
hbar_omega_min = -1.
hbar_omega_max = 0.7
Reading reduced TOFTOF files
[3]:
def read_toftof_file(filename):
"""
Read reduced files from TOFTOF
The file contains hbar omega, I, I_err for different values of q
Some metadata are stored as the header of the array in each case
Parameters
----------
filename: str
file containing TOFTOF reduced data
Returns
-------
numpy.array: 4 columns: hbar omega, I, I_err and index corresponding to different values of q
"""
with open(filename, 'r') as fin:
file = fin.readlines()
for i, item in enumerate(file):
# get number of recording per spectrum
if i == 0:
REC_PER_SPEC = int(item.split()[0]) + 1
if (i % REC_PER_SPEC)==0:
number_records = int(item.split()[0])
elif (i % REC_PER_SPEC)==2:
q_value = item.strip()[0]
elif (i % REC_PER_SPEC) == 1 or (i % REC_PER_SPEC) == 3:
pass
elif i == 4:
list_2_store = [float(elt) for elt in item.strip().split()]
list_2_store.append(float(q_value))
data = np.array(list_2_store)
elif (i % REC_PER_SPEC) == 4 and i > 4:
pass
else:
list_2_store = [float(elt) for elt in item.strip().split()]
list_2_store.append(float(q_value))
new_row = np.array(list_2_store)
data = np.vstack([data, new_row])
return data
[4]:
data_350K = read_toftof_file('./data/toftof/In20_2p5_350K__tt_cr.txt')
data_450K = read_toftof_file('./data/toftof/In20_2p5_450K__tt_cr.txt')
data_540K = read_toftof_file('./data/toftof/In20_2p5_540K__tt_cr.txt')
vana = read_toftof_file('./data/toftof/Vana_2p5_cool__tt_cr.txt')
We define a single dictionary to contain all experimental data.
The keys correspond to either a temperature for the sample or to the vanadium run.
The values are also dictionaries storing \(\hbar \omega\), \(S\) and its uncertainty for each value of \(q\).
For example,
dict_toftof['350K']['q_0.6']['x'] stores the values of \(\hbar \omega\) at \(q=0.6\) 1/Angstrom and \(T=350K\).The ‘y’ and ‘e’ keys of
dict_toftof['350K']['q_0.6'] correspond to the values of \(S\) and the associated uncertainty, respectively.[5]:
dict_toftof = {}
dict_data_350K = {}
dict_data_450K = {}
dict_data_540K = {}
dict_vana = {}
for i in range(1, 10):
key = f"q_{q_tof[i-1]}"
# 350K
indx = (data_350K[:, 3]==i)
dict_data_350K[key] = {'x': data_350K[indx][:, 0],
'y': data_350K[indx][:, 1],
'e': data_350K[indx][:, 2]}
# 450K
indx = (data_450K[:, 3]==i)
dict_data_450K[key] = {'x': data_450K[indx][:, 0],
'y': data_450K[indx][:, 1],
'e': data_450K[indx][:, 2]}
# 540K
indx = (data_540K[:, 3]==i)
dict_data_540K[key] = {'x': data_540K[indx][:, 0],
'y': data_540K[indx][:, 1],
'e': data_540K[indx][:, 2]}
# vanadium
indx = (vana[:, 3]==i)
dict_vana[key] = {'x': vana[indx][:, 0],
'y': vana[indx][:, 1],
'e': vana[indx][:, 2]}
dict_toftof['350K'] = dict_data_350K
dict_toftof['450K'] = dict_data_450K
dict_toftof['540K'] = dict_data_540K
dict_toftof['vana'] = dict_vana
[6]:
fig, ax = plt.subplots(nrows=2,
ncols=2,
figsize=(10, 10),
sharex=True,
sharey=True)
[ax[0, 0].semilogy(dict_toftof['350K'][q_value]['x'],
dict_toftof['350K'][q_value]['y'],
label=f"{q_value}") for q_value in dict_toftof['350K']]
[ax[0, 1].semilogy(dict_toftof['450K'][q_value]['x'],
dict_toftof['450K'][q_value]['y'],
label=f"{q_value}") for q_value in dict_toftof['450K']]
[ax[1, 0].semilogy(dict_toftof['540K'][q_value]['x'],
dict_toftof['540K'][q_value]['y'],
label=f"{q_value}") for q_value in dict_toftof['540K']]
[ax[1, 1].semilogy(dict_toftof['vana'][q_value]['x'],
dict_toftof['vana'][q_value]['y'],
label=f"{q_value}") for q_value in dict_toftof['vana']]
for i in range(2):
for j in range(2):
ax[i, j].grid()
ax[i, j].legend()
fig.suptitle("All TOFTOF sample data")
ax[0, 0].set_title('350K')
ax[0, 1].set_title('450K')
ax[1, 0].set_title('540K')
ax[1, 1].set_title('vanadium')
ax[1, 0].set_xlabel(r"$\hbar \omega$ (meV)")
ax[1, 1].set_xlabel(r"$\hbar \omega$ (meV)")
ax[0, 0].set_ylabel(r"S(Q, $\hbar \omega$) (arb. units)")
ax[1, 0].set_ylabel(r"S(Q, $\hbar \omega$) (arb. units)");
Fitting
Creating the Model
Dirac and Lorentzian functions
[7]:
def toftof_delta_lorentz(w, aD, aL, cst_bgd, center, gamma):
""" Define S(Q, omega) using the QENS models
Parameters
----------
w: float or list or :class:`~numpy:numpy.ndarray`
Energy (meV)
aD: float
Scaling coefficient of Delta function
aL: float
Scaling coefficient of Lorentzian function
cst_bgd: float
Constant coefficient of the background
center: float
Center of the Delta and Lorentzian functions
gamma: float
Half Width at Half Maximum of Lorentzian function
Returns
-------
numpy.array: 1 column corresponding to Dirac + Lorentzian + constant_background
"""
return QENSmodels.delta(w,
scale=aD,
center=center) + \
QENSmodels.lorentzian(w,
scale=aL,
center=center,
hwhm=gamma) + \
QENSmodels.background_polynomials(w,
[cst_bgd])
def S_meas_no_convol(w, T, aD, aL, cst_bgd, center, gamma):
""" Define S_meas = S(Q, omega) multiplied by Bose factor
Parameters
----------
w: float or list or :class:`~numpy:numpy.ndarray`
Energy (meV)
T: float
Temperature (K)
aD: float
Scaling coefficient of Delta function
aL: float
Scaling coefficient of Lorentzian function
cst_bgd: float
Constant coefficient of the background
center: float
Center of the Delta and Lorentzian functions
gamma: float
Half Width at Half Maximum of Lorentzian function
Returns
-------
numpy.array: 1 column corresponding to (Dirac + Lorentzian + constant_background) * Bose factor
"""
beta = scipy.constants.physical_constants['electron volt-joule relationship'][0] / (1000 * scipy.constants.Boltzmann * T)
return (w * beta * toftof_delta_lorentz(w, aD, aL, cst_bgd, center, gamma)
* (1. + 1. / (np.exp(beta * w) - 1)))
Convolution
From https://github.com/jmborr/qef
[8]:
def convolve(model, resolution):
r"""
Convolution of resolution with model data
.. math:: (model \otimes resolution)[n] = dx * \sum_m model[m] * resolution[m-n]
"""
c = np.convolve(model, resolution, mode='valid')
if len(model) % len(resolution) == 0:
c = c[:-1]
return c
class Convolve(CompositeModel):
"""
Convolution between model and resolution.
It is assumed that the resolution FWHM is energy independent.
Non-symmetric energy ranges are allowed (when the range of negative values
is different than that of positive values).
The convolution requires multiplication by the X-spacing to preserve
normalization
"""
def __init__(self, resolution, model, **kws):
super(Convolve, self).__init__(resolution, model, convolve, **kws)
self.resolution = resolution
self.model = model
def eval(self, params=None, **kwargs):
res_data = self.resolution.eval(params=params, **kwargs)
# evaluate model on an extended energy range to avoid boundary effects
independent_var = self.resolution.independent_vars[0]
e = kwargs[independent_var] # energy values
neg_e = min(e) - np.flip(e[np.where(e > 0)], axis=0)
pos_e = max(e) - np.flip(e[np.where(e < 0)], axis=0)
e = np.concatenate((neg_e, e, pos_e))
kwargs.update({independent_var: e})
model_data = self.model.eval(params=params, **kwargs)
# Multiply by the X-spacing to preserve normalization
de = (e[-1] - e[0])/(len(e) - 1) # energy spacing
return de * convolve(model_data, res_data)
Complete fitting model and procedure
[9]:
def build_and_fit_S_meas_single_q(q_value, ini_vals, plot=True):
"""
Define the complete fitting model S_meas
Parameters
----------
q_value: float
momentum transfer (non-fitting, in 1/Angstrom)
ini_vals: dictionary of initial guesses for the parameters.
For example, {'T': 350 , 'aD': 2, 'aL': 10, 'cst_bgd': 0.1, 'center': 0, 'gamma': 0.4}
Note that the temperature, T, is not refined
plot: boolean
If True, plots initial model, experimental data and fitting result
Returns
-------
lmfit.model.ModelResult: structure containing the output of a fit using lmfit for a single value of q
"""
temperature_key = str(ini_vals['T']) + 'K'
def irf_gate(w):
"""
Function defined from the interpolation of instrument resolution data
Used to define fitting model and plot
"""
f_interp = interp1d(dict_toftof['vana'][q_value]['x'],
dict_toftof['vana'][q_value]['y']/np.sum(dict_toftof['vana'][q_value]['y']),
kind='cubic', bounds_error=False, fill_value='extrapolate')
return f_interp(w)
gmodel = Convolve(Model(irf_gate), Model(S_meas_no_convol))
pars = lmfit.Parameters()
pars.add('T', value=ini_vals['T'], vary=False)
pars.add('aD', value=ini_vals['aD'], min=0, vary=True)
pars.add('aL', value=ini_vals['aL'], min=0.01, vary=True)
pars.add('cst_bgd', value=ini_vals['cst_bgd'], min=0, vary=True)
pars.add('center', value=ini_vals['center'], min=-1, max=1, vary=True)
pars.add('gamma', value=ini_vals['gamma'], min=0.01, vary=True)
xx = dict_toftof[temperature_key][q_value]['x']
mask = np.intersect1d(np.where(xx>=hbar_omega_min), np.where(xx<=hbar_omega_max))
# remove NaNs due to singularity on model when xx=0
mask = np.intersect1d(mask, np.where(abs(xx)>=1e-12))
# remove negative S(Q, hbar omega) values
mask = np.intersect1d(mask, np.where(dict_toftof[temperature_key][q_value]['y']>0))
out = gmodel.fit(dict_toftof[temperature_key][q_value]['y'][mask],
pars,
method='leastsq',
w=xx[mask],
weights=1/dict_toftof[temperature_key][q_value]['e'][mask],
fit_kws={'xtol': 1e-9, 'ftol': 1e-9},
)
if plot:
# plot initial model, fit results and experimental data
fig1, ax1 = plt.subplots()
ax1.semilogy(xx[mask],
dict_toftof[temperature_key][q_value]['y'][mask],
'.-',
label='exp data')
ax1.semilogy(xx[mask], out.init_fit, '--', label='initial fit')
ax1.semilogy(xx[mask], out.best_fit, 'r-', label='best fit')
ax1.set_title(f'{q_value}, reduced chi_sq={out.redchi:.2f}')
ax1.set_xlabel(r"$\hbar \omega$ (meV)")
ax1.set_ylabel(r"S(Q, $\hbar \omega$) (arb. units)")
ax1.legend(); ax1.grid()
plt.show()
return out
def fit_S_meas_all_q(ini_vals, plot=True):
"""
Fit S_meas (convolution with resolution function) using lmfit
and plot experimental data, initial model and fit for all q-values
(One value of q per plot)
Parameters
----------
ini_vals: dictionary of initial guesses for the parameters.
For example, {'T': 350 , 'aD': 2, 'aL': 10, 'cst_bgd': 0.1, 'center': 0, 'gamma': 0.4}
Note that the temperature, T, is not refined
plot: boolean
If True, plots initial model, experimental data and fitting result
Returns
-------
dictionary of lmfit.model.ModelResult
Each key corresponds to a value of q. The value is the output of a fit using lmfit for a single value of q
"""
result_fit = {} # store fits for all spectra
temperature_key = str(ini_vals['T']) + 'K'
# Loop over q-values
for i, q_value in enumerate(dict_toftof[temperature_key]):
result_fit[str(q_value)] = build_and_fit_S_meas_single_q(q_value, ini_vals, plot=plot)
if plot:
fig3, ax3 = plt.subplots(nrows=3, ncols=2, sharex=True, figsize=(10, 10))
ax3[0, 0].plot([key for key in result_fit.keys()],
[result_fit[key].params['aD'].value for key in result_fit.keys()],
marker='o')
ax3[0, 0].set_ylabel('aD')
ax3[0, 1].plot([key for key in result_fit.keys()],
[result_fit[key].params['aL'].value for key in result_fit.keys()],
marker='o')
ax3[0, 1].set_ylabel('aL')
ax3[1, 0].plot([key for key in result_fit.keys()],
[result_fit[key].params['center'].value for key in result_fit.keys()],
marker='o')
ax3[1, 0].set_ylabel('center')
ax3[1, 1].plot([key for key in result_fit.keys()],
[result_fit[key].params['gamma'].value for key in result_fit.keys()],
marker='o')
ax3[1, 1].set_ylabel('gamma')
ax3[2, 0].plot([key for key in result_fit.keys()],
[result_fit[key].params['cst_bgd'].value for key in result_fit.keys()],
marker='o')
ax3[2, 0].set_ylabel('cst_bgd'); ax3[2, 0].set_xlabel('q')
ax3[2, 1].set_ylabel('reduced chi squared'); ax3[2, 1].set_xlabel('q values')
ax3[2, 1].bar(range(len(q_tof)),
[result_fit[key].redchi for key in result_fit.keys()],
align='center',
tick_label=list(result_fit.keys()))
plt.tight_layout()
for i in range(3):
for j in range(2):
ax3[i, j].grid()
return result_fit
[10]:
def dict_ref_params_for_a_temperature(result_fit):
"""
Creates a pandas.DataFrame to store the values and errors of the refined parameters at a given temperature.
Each row corresponds to a value of q.
Note that we replace NaN by zero.
Parameters
----------
result_fit : dictionary of lmfit.model.ModelResult for a given temperature
Each key corresponds to a value of q.
The value is the output of a fit using lmfit for the corresponding q
Returns
-------
pandas.DataFrame: stores the value and error of the refined parameters
"""
q_keys = result_fit.keys()
data_frame_one_temp = pd.DataFrame(
data= {"aD": [result_fit[q].params['aD'].value for q in q_keys],
"aD_err": [result_fit[q].params['aD'].stderr for q in q_keys],
"aL": [result_fit[q].params['aL'].value for q in q_keys],
"aL_err": [result_fit[q].params['aL'].stderr for q in q_keys],
"cst_bgd": [result_fit[q].params['cst_bgd'].value for q in q_keys],
"cst_bgd_err": [result_fit[q].params['cst_bgd'].stderr for q in q_keys],
"center": [result_fit[q].params['center'].value for q in q_keys],
"center_err": [result_fit[q].params['center'].stderr for q in q_keys],
"gamma": [result_fit[q].params['gamma'].value for q in q_keys],
"gamma_err": [result_fit[q].params['gamma'].stderr for q in q_keys],
"red_chi_squared": [result_fit[q].redchi for q in q_keys],
},
index=q_keys).fillna(0)
return data_frame_one_temp
Creating and defining parameters and functions to collect and display fitting results
[11]:
# This dictionary will store the value and error of the refined parameters for each temperature as a pandas dataframe.
dict_ref_params = {}
[12]:
def plot_ratio_coefficients(dict_to_plot):
"""
Plot ratios of aD and aL to compare with publication
Parameters
----------
dict_to_plot: select key of dictionary to plot: for example, dict_ref_params['350K']
Returns
-------
None
"""
std_err_ratio = np.sqrt((dict_to_plot['aD'].values**2 * dict_to_plot['aL_err'].values**2
+ dict_to_plot['aL'].values**2 * dict_to_plot['aD_err'].values**2)
/ (dict_to_plot['aD'].values + dict_to_plot['aL'].values)**4)
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax[0].errorbar(dict_to_plot.index,
dict_to_plot['aD'].values / (dict_to_plot['aD'].values + dict_to_plot['aL'].values),
std_err_ratio,
marker='o',
linestyle='--')
ax[0].set_ylabel('$a_D/(a_D+a_L)$')
ax[0].set_ylim((0, 1))
ax[1].errorbar(dict_to_plot.index,
dict_to_plot['aL'].values / (dict_to_plot['aD'].values + dict_to_plot['aL'].values),
std_err_ratio,
marker='o',
linestyle='--')
ax[1].set_ylabel('$a_L/(a_D+a_L$)')
ax[1].set_ylim((0, 1))#0.5))
for i in range(2):
ax[i].grid()
ax[i].set_xlabel(r'$Q (1/\AA)$')
for tick in ax[i].get_xticklabels():
tick.set_rotation(45)
plt.tight_layout()
plt.show()
return
Fitting at T = 350K
[13]:
ini_vals_350K = {'T': 350,
'aD': 1,
'aL': 40,
'cst_bgd': 0.04,
'center': -1e-7,
'gamma': 0.1}
fit350 = fit_S_meas_all_q(ini_vals_350K)
[14]:
for item in fit350.keys():
print(f'----------\nResult of fit for {item}:\n',
fit350[item].fit_report())
----------
Result of fit for q_0.6:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 227
# data points = 34
# variables = 5
chi-square = 51.6875996
reduced chi-square = 1.78233102
Akaike info crit = 24.2411508
Bayesian info crit = 31.8729534
R-squared = 0.94612305
## Warning: uncertainties could not be estimated:
aL: at boundary
[[Variables]]
T: 350 (fixed)
aD: 162.021246 (init = 1)
aL: 0.01000001 (init = 40)
cst_bgd: 0.32552634 (init = 0.04)
center: -0.05029700 (init = -1e-07)
gamma: 0.01002385 (init = 0.1)
----------
Result of fit for q_0.8:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 2544
# data points = 34
# variables = 5
chi-square = 7.78665625
reduced chi-square = 0.26850539
Akaike info crit = -40.1142658
Bayesian info crit = -32.4824632
R-squared = 0.99149584
## Warning: uncertainties could not be estimated:
cst_bgd: at boundary
gamma: at boundary
[[Variables]]
T: 350 (fixed)
aD: 140.600595 (init = 1)
aL: 94.2895752 (init = 40)
cst_bgd: 1.2922e-09 (init = 0.04)
center: -5.3771e-10 (init = -1e-07)
gamma: 0.01000000 (init = 0.1)
----------
Result of fit for q_1.0:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 437
# data points = 34
# variables = 5
chi-square = 41.0859028
reduced chi-square = 1.41675527
Akaike info crit = 16.4363543
Bayesian info crit = 24.0681570
R-squared = 0.95599024
## Warning: uncertainties could not be estimated:
gamma: at boundary
[[Variables]]
T: 350 (fixed)
aD: 140.215403 (init = 1)
aL: 104.529441 (init = 40)
cst_bgd: 0.05944927 (init = 0.04)
center: -2.0172e-09 (init = -1e-07)
gamma: 0.01000000 (init = 0.1)
----------
Result of fit for q_1.2:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 1301
# data points = 34
# variables = 5
chi-square = 15.9822260
reduced chi-square = 0.55111124
Akaike info crit = -15.6660321
Bayesian info crit = -8.03422946
R-squared = 0.98277193
## Warning: uncertainties could not be estimated:
cst_bgd: at boundary
gamma: at boundary
[[Variables]]
T: 350 (fixed)
aD: 136.150152 (init = 1)
aL: 123.231659 (init = 40)
cst_bgd: 5.3832e-09 (init = 0.04)
center: -1.8017e-09 (init = -1e-07)
gamma: 0.01000000 (init = 0.1)
----------
Result of fit for q_1.7:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 461
# data points = 34
# variables = 5
chi-square = 11.9454578
reduced chi-square = 0.41191234
Akaike info crit = -25.5643201
Bayesian info crit = -17.9325175
R-squared = 0.98616150
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 350 (fixed)
aD: 134.012085 (init = 1)
aL: 40.4025544 (init = 40)
cst_bgd: 0.58316481 (init = 0.04)
center: -2.6425e-11 (init = -1e-07)
gamma: 0.03111282 (init = 0.1)
----------
Result of fit for q_1.9:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 247
# data points = 34
# variables = 5
chi-square = 5.77290980
reduced chi-square = 0.19906586
Akaike info crit = -50.2882653
Bayesian info crit = -42.6564627
R-squared = 0.99323979
[[Variables]]
T: 350 (fixed)
aD: 126.304360 +/- 4.05392517 (3.21%) (init = 1)
aL: 139.484223 +/- 274.572499 (196.85%) (init = 40)
cst_bgd: 0.52247309 +/- 0.13280326 (25.42%) (init = 0.04)
center: -0.00271709 +/- 0.00391600 (144.12%) (init = -1e-07)
gamma: 0.01000003 +/- 0.01838278 (183.83%) (init = 0.1)
[[Correlations]] (unreported correlations are < 0.100)
C(aL, gamma) = -0.999
C(aD, aL) = -0.701
C(aD, center) = 0.700
C(aD, gamma) = 0.662
C(aD, cst_bgd) = 0.476
C(cst_bgd, center) = 0.400
C(aL, center) = -0.399
C(center, gamma) = 0.371
C(cst_bgd, gamma) = -0.196
C(aL, cst_bgd) = 0.151
----------
Result of fit for q_2.3:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 422
# data points = 34
# variables = 5
chi-square = 12.6151137
reduced chi-square = 0.43500392
Akaike info crit = -23.7098075
Bayesian info crit = -16.0780049
R-squared = 0.98477413
## Warning: uncertainties could not be estimated:
gamma: at boundary
[[Variables]]
T: 350 (fixed)
aD: 127.214601 (init = 1)
aL: 139.725264 (init = 40)
cst_bgd: 0.73536362 (init = 0.04)
center: -5.5085e-09 (init = -1e-07)
gamma: 0.01000000 (init = 0.1)
----------
Result of fit for q_3.5:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 413
# data points = 34
# variables = 5
chi-square = 10.8341633
reduced chi-square = 0.37359184
Akaike info crit = -28.8843078
Bayesian info crit = -21.2525052
R-squared = 0.98448776
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 350 (fixed)
aD: 117.247075 (init = 1)
aL: 36.9115245 (init = 40)
cst_bgd: 0.43955199 (init = 0.04)
center: -5.7454e-09 (init = -1e-07)
gamma: 0.05710292 (init = 0.1)
----------
Result of fit for q_3.9:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 359
# data points = 34
# variables = 5
chi-square = 14.7833730
reduced chi-square = 0.50977148
Akaike info crit = -18.3171523
Bayesian info crit = -10.6853496
R-squared = 0.97741417
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 350 (fixed)
aD: 107.387591 (init = 1)
aL: 60.7418892 (init = 40)
cst_bgd: 0.71141647 (init = 0.04)
center: -3.4166e-09 (init = -1e-07)
gamma: 0.03121186 (init = 0.1)
Storing fitted parameters and quality factors
[15]:
dict_ref_params['350K'] = dict_ref_params_for_a_temperature(fit350)
dict_ref_params['350K']
[15]:
| aD | aD_err | aL | aL_err | cst_bgd | cst_bgd_err | center | center_err | gamma | gamma_err | red_chi_squared | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| q_0.6 | 162.021246 | 0.000000 | 0.010000 | 0.000000 | 3.255263e-01 | 0.000000 | -5.029700e-02 | 0.000000 | 0.010024 | 0.000000 | 1.782331 |
| q_0.8 | 140.600595 | 0.000000 | 94.289575 | 0.000000 | 1.292243e-09 | 0.000000 | -5.377085e-10 | 0.000000 | 0.010000 | 0.000000 | 0.268505 |
| q_1.0 | 140.215403 | 0.000000 | 104.529441 | 0.000000 | 5.944927e-02 | 0.000000 | -2.017207e-09 | 0.000000 | 0.010000 | 0.000000 | 1.416755 |
| q_1.2 | 136.150152 | 0.000000 | 123.231659 | 0.000000 | 5.383250e-09 | 0.000000 | -1.801730e-09 | 0.000000 | 0.010000 | 0.000000 | 0.551111 |
| q_1.7 | 134.012085 | 0.000000 | 40.402554 | 0.000000 | 5.831648e-01 | 0.000000 | -2.642497e-11 | 0.000000 | 0.031113 | 0.000000 | 0.411912 |
| q_1.9 | 126.304360 | 4.053925 | 139.484223 | 274.572499 | 5.224731e-01 | 0.132803 | -2.717089e-03 | 0.003916 | 0.010000 | 0.018383 | 0.199066 |
| q_2.3 | 127.214601 | 0.000000 | 139.725264 | 0.000000 | 7.353636e-01 | 0.000000 | -5.508479e-09 | 0.000000 | 0.010000 | 0.000000 | 0.435004 |
| q_3.5 | 117.247075 | 0.000000 | 36.911525 | 0.000000 | 4.395520e-01 | 0.000000 | -5.745414e-09 | 0.000000 | 0.057103 | 0.000000 | 0.373592 |
| q_3.9 | 107.387591 | 0.000000 | 60.741889 | 0.000000 | 7.114165e-01 | 0.000000 | -3.416569e-09 | 0.000000 | 0.031212 | 0.000000 | 0.509771 |
Plotting \(\frac{a_D}{a_D+a_L}\) and \(\frac{a_L}{a_D+a_L}\) at 350K
[16]:
plot_ratio_coefficients(dict_ref_params['350K'])
Tuning fits for some selected q-values (optional)
If the above fit is not satisfactory for only a few values of \(q\), the fit at these particular \(q\) can be rerun using the following method.
[17]:
dict_ref_params['350K']
[17]:
| aD | aD_err | aL | aL_err | cst_bgd | cst_bgd_err | center | center_err | gamma | gamma_err | red_chi_squared | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| q_0.6 | 162.021246 | 0.000000 | 0.010000 | 0.000000 | 3.255263e-01 | 0.000000 | -5.029700e-02 | 0.000000 | 0.010024 | 0.000000 | 1.782331 |
| q_0.8 | 140.600595 | 0.000000 | 94.289575 | 0.000000 | 1.292243e-09 | 0.000000 | -5.377085e-10 | 0.000000 | 0.010000 | 0.000000 | 0.268505 |
| q_1.0 | 140.215403 | 0.000000 | 104.529441 | 0.000000 | 5.944927e-02 | 0.000000 | -2.017207e-09 | 0.000000 | 0.010000 | 0.000000 | 1.416755 |
| q_1.2 | 136.150152 | 0.000000 | 123.231659 | 0.000000 | 5.383250e-09 | 0.000000 | -1.801730e-09 | 0.000000 | 0.010000 | 0.000000 | 0.551111 |
| q_1.7 | 134.012085 | 0.000000 | 40.402554 | 0.000000 | 5.831648e-01 | 0.000000 | -2.642497e-11 | 0.000000 | 0.031113 | 0.000000 | 0.411912 |
| q_1.9 | 126.304360 | 4.053925 | 139.484223 | 274.572499 | 5.224731e-01 | 0.132803 | -2.717089e-03 | 0.003916 | 0.010000 | 0.018383 | 0.199066 |
| q_2.3 | 127.214601 | 0.000000 | 139.725264 | 0.000000 | 7.353636e-01 | 0.000000 | -5.508479e-09 | 0.000000 | 0.010000 | 0.000000 | 0.435004 |
| q_3.5 | 117.247075 | 0.000000 | 36.911525 | 0.000000 | 4.395520e-01 | 0.000000 | -5.745414e-09 | 0.000000 | 0.057103 | 0.000000 | 0.373592 |
| q_3.9 | 107.387591 | 0.000000 | 60.741889 | 0.000000 | 7.114165e-01 | 0.000000 | -3.416569e-09 | 0.000000 | 0.031212 | 0.000000 | 0.509771 |
[18]:
# selected q-value for which to fit again
q_indx = [0, 1]
redo_q_keys = [list(fit350.keys())[q] for q in q_indx]
new_ini_vals = ini_vals_350K.copy()
stored_params = dict_ref_params['350K'].columns
data_frame_selected_qs = pd.DataFrame()
for indx, item in enumerate(redo_q_keys):
# store previous fitting result
old_fit_result = dict_ref_params['350K'].loc[item]
# define new initial parameters' values
new_ini_vals['aD'] = 2
new_ini_vals['aL'] = 20
new_ini_vals['gamma'] = 0.2
new_ini_vals['cst_bgd'] = 1
# rerun fitting
new_out = build_and_fit_S_meas_single_q(item, new_ini_vals, plot=True)
new_row = pd.DataFrame(
data = {"aD": [new_out.params['aD'].value],
"aD_err": [new_out.params['aD'].stderr],
"aL": [new_out.params['aL'].value],
"aL_err": [new_out.params['aL'].stderr],
"cst_bgd": [new_out.params['cst_bgd'].value],
"cst_bgd_err": [new_out.params['cst_bgd'].stderr],
"center": [new_out.params['center'].value],
"center_err": [new_out.params['center'].stderr],
"gamma": [new_out.params['gamma'].value],
"gamma_err": [new_out.params['gamma'].stderr],
"red_chi_squared": [new_out.redchi],
},
index=[item]).fillna(0)
data_frame_selected_qs = pd.concat([data_frame_selected_qs, new_row])
print((f'\nPrevious initial values: {ini_vals_350K}\n'),
(f'\nPrevious fit result at {item}:\n{old_fit_result}\n'),
('\n==================================================\n'),
(f'\nNew initial values: {new_ini_vals}\n'),
(f'\nNew fit result at {item}:\n{data_frame_selected_qs.loc[item]}')
)
Previous initial values: {'T': 350, 'aD': 1, 'aL': 40, 'cst_bgd': 0.04, 'center': -1e-07, 'gamma': 0.1}
Previous fit result at q_0.6:
aD 162.021246
aD_err 0.000000
aL 0.010000
aL_err 0.000000
cst_bgd 0.325526
cst_bgd_err 0.000000
center -0.050297
center_err 0.000000
gamma 0.010024
gamma_err 0.000000
red_chi_squared 1.782331
Name: q_0.6, dtype: float64
==================================================
New initial values: {'T': 350, 'aD': 2, 'aL': 20, 'cst_bgd': 1, 'center': -1e-07, 'gamma': 0.2}
New fit result at q_0.6:
aD 162.022131
aD_err 2.187066
aL 0.010001
aL_err 2.227643
cst_bgd 0.324623
cst_bgd_err 0.555146
center -0.003586
center_err 43.514908
gamma 0.152750
gamma_err 125.515161
red_chi_squared 1.782579
Name: q_0.6, dtype: float64
Previous initial values: {'T': 350, 'aD': 1, 'aL': 40, 'cst_bgd': 0.04, 'center': -1e-07, 'gamma': 0.1}
Previous fit result at q_0.8:
aD 1.406006e+02
aD_err 0.000000e+00
aL 9.428958e+01
aL_err 0.000000e+00
cst_bgd 1.292243e-09
cst_bgd_err 0.000000e+00
center -5.377085e-10
center_err 0.000000e+00
gamma 1.000000e-02
gamma_err 0.000000e+00
red_chi_squared 2.685054e-01
Name: q_0.8, dtype: float64
==================================================
New initial values: {'T': 350, 'aD': 2, 'aL': 20, 'cst_bgd': 1, 'center': -1e-07, 'gamma': 0.2}
New fit result at q_0.8:
aD 1.406006e+02
aD_err 0.000000e+00
aL 9.428956e+01
aL_err 0.000000e+00
cst_bgd 6.193712e-11
cst_bgd_err 0.000000e+00
center -6.281897e-10
center_err 0.000000e+00
gamma 1.000000e-02
gamma_err 0.000000e+00
red_chi_squared 2.685054e-01
Name: q_0.8, dtype: float64
If the new fitting results are better than the old ones, the dictionary containing the fitted parameters can be changed using the following method.
[19]:
for item in redo_q_keys:
dict_ref_params['350K'].loc[item,:] = data_frame_selected_qs.loc[item,:]
Fitting at T = 450K
[20]:
ini_vals_450K = {'T': 450,
'aD': 2,
'aL': 20,
'cst_bgd': 0.01,
'center': -1e-7,
'gamma': 0.5}
fit450 = fit_S_meas_all_q(ini_vals_450K)
/tmp/ipykernel_1950/3601744967.py:65: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
fig1, ax1 = plt.subplots()
[21]:
for key in fit450.keys():
print(fit450[key].params['gamma'].value)
0.2993277049065741
0.010000000067557524
0.010000000162895706
0.02255181713742438
0.20277187557712328
0.05329829058595026
0.04997755395595749
0.09754213986639626
0.10001211628251294
[22]:
for item in fit450.keys():
print(f'----------\nResult of fit for {item}:\n',
fit450[item].fit_report())
----------
Result of fit for q_0.6:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 187
# data points = 34
# variables = 5
chi-square = 43.7914165
reduced chi-square = 1.51004884
Akaike info crit = 18.6046283
Bayesian info crit = 26.2364309
R-squared = 0.95312794
[[Variables]]
T: 450 (fixed)
aD: 158.677519 +/- 4.73805081 (2.99%) (init = 2)
aL: 3.70476658 +/- 3.04405310 (82.17%) (init = 20)
cst_bgd: 8.6099e-05 +/- 1.27204820 (1477431.19%) (init = 0.01)
center: -1.0256e-08 +/- 0.06987671 (681320801.02%) (init = -1e-07)
gamma: 0.29932770 +/- 1.10032650 (367.60%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
C(aD, gamma) = 0.918
C(cst_bgd, gamma) = -0.917
C(aD, cst_bgd) = -0.746
C(aD, aL) = -0.405
C(aL, cst_bgd) = -0.237
C(aL, center) = 0.185
C(aL, gamma) = -0.140
C(aD, center) = -0.113
----------
Result of fit for q_0.8:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 546
# data points = 34
# variables = 5
chi-square = 15.1673664
reduced chi-square = 0.52301263
Akaike info crit = -17.4452881
Bayesian info crit = -9.81348546
R-squared = 0.98294157
## Warning: uncertainties could not be estimated:
gamma: at boundary
[[Variables]]
T: 450 (fixed)
aD: 134.276380 (init = 2)
aL: 112.589524 (init = 20)
cst_bgd: 0.24954321 (init = 0.01)
center: -6.9180e-09 (init = -1e-07)
gamma: 0.01000000 (init = 0.5)
----------
Result of fit for q_1.0:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 3009
# data points = 34
# variables = 5
chi-square = 34.3359012
reduced chi-square = 1.18399659
Akaike info crit = 10.3342528
Bayesian info crit = 17.9660555
R-squared = 0.96183671
## Warning: uncertainties could not be estimated:
cst_bgd: at boundary
gamma: at boundary
[[Variables]]
T: 450 (fixed)
aD: 121.460581 (init = 2)
aL: 197.409880 (init = 20)
cst_bgd: 1.0383e-10 (init = 0.01)
center: -3.9375e-09 (init = -1e-07)
gamma: 0.01000000 (init = 0.5)
----------
Result of fit for q_1.2:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 434
# data points = 34
# variables = 5
chi-square = 14.7052511
reduced chi-square = 0.50707763
Akaike info crit = -18.4972997
Bayesian info crit = -10.8654971
R-squared = 0.98280484
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 450 (fixed)
aD: 114.323599 (init = 2)
aL: 106.909694 (init = 20)
cst_bgd: 0.07078031 (init = 0.01)
center: -2.8608e-09 (init = -1e-07)
gamma: 0.02255182 (init = 0.5)
----------
Result of fit for q_1.7:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 129
# data points = 34
# variables = 5
chi-square = 17.6298751
reduced chi-square = 0.60792673
Akaike info crit = -12.3300307
Bayesian info crit = -4.69822812
R-squared = 0.97704487
[[Variables]]
T: 450 (fixed)
aD: 135.993697 +/- 5.13273875 (3.77%) (init = 2)
aL: 18.5453410 +/- 4.93325581 (26.60%) (init = 20)
cst_bgd: 0.74049209 +/- 0.59224950 (79.98%) (init = 0.01)
center: -1.1185e-08 +/- 0.01776411 (158816260.16%) (init = -1e-07)
gamma: 0.20277188 +/- 0.12636601 (62.32%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
C(aD, gamma) = 0.957
C(aD, aL) = -0.954
C(aL, gamma) = -0.880
C(cst_bgd, gamma) = -0.850
C(aD, cst_bgd) = -0.707
C(aL, center) = -0.614
C(aD, center) = 0.549
C(aL, cst_bgd) = 0.519
C(center, gamma) = 0.506
C(cst_bgd, center) = -0.251
----------
Result of fit for q_1.9:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 472
# data points = 34
# variables = 5
chi-square = 11.9184640
reduced chi-square = 0.41098152
Akaike info crit = -25.6412388
Bayesian info crit = -18.0094362
R-squared = 0.98411107
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 450 (fixed)
aD: 110.809424 (init = 2)
aL: 55.5261452 (init = 20)
cst_bgd: 1.05991021 (init = 0.01)
center: -5.2051e-09 (init = -1e-07)
gamma: 0.05329829 (init = 0.5)
----------
Result of fit for q_2.3:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 384
# data points = 34
# variables = 5
chi-square = 14.4682817
reduced chi-square = 0.49890627
Akaike info crit = -19.0496592
Bayesian info crit = -11.4178566
R-squared = 0.97918965
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 450 (fixed)
aD: 99.8529875 (init = 2)
aL: 72.1356688 (init = 20)
cst_bgd: 0.99903951 (init = 0.01)
center: -1.4892e-09 (init = -1e-07)
gamma: 0.04997755 (init = 0.5)
----------
Result of fit for q_3.5:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 485
# data points = 34
# variables = 5
chi-square = 13.2723296
reduced chi-square = 0.45766654
Akaike info crit = -21.9830907
Bayesian info crit = -14.3512880
R-squared = 0.97573346
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 450 (fixed)
aD: 99.3803260 (init = 2)
aL: 40.2196736 (init = 20)
cst_bgd: 1.56052431 (init = 0.01)
center: -3.6556e-09 (init = -1e-07)
gamma: 0.09754214 (init = 0.5)
----------
Result of fit for q_3.9:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 82
# data points = 34
# variables = 5
chi-square = 8.87675892
reduced chi-square = 0.30609514
Akaike info crit = -35.6594167
Bayesian info crit = -28.0276141
R-squared = 0.98261423
[[Variables]]
T: 450 (fixed)
aD: 98.0932671 +/- 7.27546336 (7.42%) (init = 2)
aL: 36.0073800 +/- 10.5143041 (29.20%) (init = 20)
cst_bgd: 1.44590678 +/- 0.24645269 (17.04%) (init = 0.01)
center: -0.00544586 +/- 0.00690306 (126.76%) (init = -1e-07)
gamma: 0.10001212 +/- 0.04095525 (40.95%) (init = 0.5)
[[Correlations]] (unreported correlations are < 0.100)
C(aD, aL) = -0.998
C(aD, gamma) = 0.984
C(aL, gamma) = -0.983
C(aD, center) = 0.828
C(aL, center) = -0.827
C(center, gamma) = 0.811
C(cst_bgd, gamma) = -0.718
C(aD, cst_bgd) = -0.611
C(aL, cst_bgd) = 0.597
C(cst_bgd, center) = -0.452
Storing fitted parameters
[23]:
dict_ref_params['450K'] = dict_ref_params_for_a_temperature(fit450)
dict_ref_params['450K']
[23]:
| aD | aD_err | aL | aL_err | cst_bgd | cst_bgd_err | center | center_err | gamma | gamma_err | red_chi_squared | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| q_0.6 | 158.677519 | 4.738051 | 3.704767 | 3.044053 | 8.609864e-05 | 1.272048 | -1.025607e-08 | 0.069877 | 0.299328 | 1.100326 | 1.510049 |
| q_0.8 | 134.276380 | 0.000000 | 112.589524 | 0.000000 | 2.495432e-01 | 0.000000 | -6.918045e-09 | 0.000000 | 0.010000 | 0.000000 | 0.523013 |
| q_1.0 | 121.460581 | 0.000000 | 197.409880 | 0.000000 | 1.038298e-10 | 0.000000 | -3.937475e-09 | 0.000000 | 0.010000 | 0.000000 | 1.183997 |
| q_1.2 | 114.323599 | 0.000000 | 106.909694 | 0.000000 | 7.078031e-02 | 0.000000 | -2.860767e-09 | 0.000000 | 0.022552 | 0.000000 | 0.507078 |
| q_1.7 | 135.993697 | 5.132739 | 18.545341 | 4.933256 | 7.404921e-01 | 0.592250 | -1.118532e-08 | 0.017764 | 0.202772 | 0.126366 | 0.607927 |
| q_1.9 | 110.809424 | 0.000000 | 55.526145 | 0.000000 | 1.059910e+00 | 0.000000 | -5.205118e-09 | 0.000000 | 0.053298 | 0.000000 | 0.410982 |
| q_2.3 | 99.852988 | 0.000000 | 72.135669 | 0.000000 | 9.990395e-01 | 0.000000 | -1.489175e-09 | 0.000000 | 0.049978 | 0.000000 | 0.498906 |
| q_3.5 | 99.380326 | 0.000000 | 40.219674 | 0.000000 | 1.560524e+00 | 0.000000 | -3.655621e-09 | 0.000000 | 0.097542 | 0.000000 | 0.457667 |
| q_3.9 | 98.093267 | 7.275463 | 36.007380 | 10.514304 | 1.445907e+00 | 0.246453 | -5.445864e-03 | 0.006903 | 0.100012 | 0.040955 | 0.306095 |
Plotting \(\frac{a_D}{a_D+a_L}\) and \(\frac{a_L}{a_D+a_L}\) at 450K
[24]:
plot_ratio_coefficients(dict_ref_params['450K'])
Fitting at T = 540K
[25]:
ini_vals_540K = {'T': 540,
'aD': 2,
'aL': 20, #10
'cst_bgd': 0.01,
'center': -1e-7,
'gamma': 0.4}
fit540 = fit_S_meas_all_q(ini_vals_540K)
[26]:
for item in fit540.keys():
print(f'----------\nResult of fit for {item}:\n',
fit540[item].fit_report())
----------
Result of fit for q_0.6:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 207
# data points = 34
# variables = 5
chi-square = 40.0521297
reduced chi-square = 1.38110792
Akaike info crit = 15.5699250
Bayesian info crit = 23.2017277
R-squared = 0.95629801
[[Variables]]
T: 540 (fixed)
aD: 158.547645 +/- 3.31754022 (2.09%) (init = 2)
aL: 1.13586531 +/- 6.28228601 (553.08%) (init = 20)
cst_bgd: 0.90628707 +/- 2.09175202 (230.80%) (init = 0.01)
center: -6.8205e-08 +/- 0.14557907 (213444758.97%) (init = -1e-07)
gamma: 0.42871318 +/- 5.40344561 (1260.39%) (init = 0.4)
[[Correlations]] (unreported correlations are < 0.100)
C(aL, cst_bgd) = -0.969
C(cst_bgd, gamma) = -0.955
C(aD, gamma) = 0.871
C(aL, gamma) = 0.869
C(aD, cst_bgd) = -0.756
C(aD, aL) = 0.610
C(center, gamma) = 0.163
C(cst_bgd, center) = -0.107
----------
Result of fit for q_0.8:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 260
# data points = 34
# variables = 5
chi-square = 19.4101642
reduced chi-square = 0.66931601
Akaike info crit = -9.05916462
Bayesian info crit = -1.42736200
R-squared = 0.97748690
[[Variables]]
T: 540 (fixed)
aD: 152.158615 +/- 3.96264939 (2.60%) (init = 2)
aL: 6.91754542 +/- 8.19579624 (118.48%) (init = 20)
cst_bgd: 0.33482953 +/- 2.72626626 (814.23%) (init = 0.01)
center: -7.3475e-09 +/- 0.04693234 (638749373.36%) (init = -1e-07)
gamma: 0.42587342 +/- 1.11870240 (262.68%) (init = 0.4)
[[Correlations]] (unreported correlations are < 0.100)
C(aL, cst_bgd) = -0.974
C(cst_bgd, gamma) = -0.952
C(aD, gamma) = 0.902
C(aL, gamma) = 0.873
C(aD, cst_bgd) = -0.777
C(aD, aL) = 0.643
C(aD, center) = 0.434
C(center, gamma) = 0.421
C(cst_bgd, center) = -0.294
C(aL, center) = 0.225
----------
Result of fit for q_1.0:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 490
# data points = 34
# variables = 5
chi-square = 20.8968702
reduced chi-square = 0.72058173
Akaike info crit = -6.54987834
Bayesian info crit = 1.08192428
R-squared = 0.97494468
## Warning: uncertainties could not be estimated:
gamma: at boundary
[[Variables]]
T: 540 (fixed)
aD: 114.190594 (init = 2)
aL: 205.967553 (init = 20)
cst_bgd: 0.92020046 (init = 0.01)
center: -1.3638e-09 (init = -1e-07)
gamma: 0.01000000 (init = 0.4)
----------
Result of fit for q_1.2:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 448
# data points = 34
# variables = 5
chi-square = 15.3307995
reduced chi-square = 0.52864826
Akaike info crit = -17.0808872
Bayesian info crit = -9.44908461
R-squared = 0.98123665
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 540 (fixed)
aD: 103.745406 (init = 2)
aL: 123.363272 (init = 20)
cst_bgd: 0.89614028 (init = 0.01)
center: -2.2194e-09 (init = -1e-07)
gamma: 0.02339889 (init = 0.4)
----------
Result of fit for q_1.7:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 151
# data points = 34
# variables = 5
chi-square = 11.9950215
reduced chi-square = 0.41362143
Akaike info crit = -25.4235403
Bayesian info crit = -17.7917377
R-squared = 0.98229696
[[Variables]]
T: 540 (fixed)
aD: 121.896291 +/- 4.26929006 (3.50%) (init = 2)
aL: 31.8858848 +/- 3.58005066 (11.23%) (init = 20)
cst_bgd: 0.46858781 +/- 0.60929435 (130.03%) (init = 0.01)
center: -0.00552346 +/- 0.00793701 (143.70%) (init = -1e-07)
gamma: 0.21309504 +/- 0.06871882 (32.25%) (init = 0.4)
[[Correlations]] (unreported correlations are < 0.100)
C(aD, gamma) = 0.961
C(aD, aL) = -0.936
C(cst_bgd, gamma) = -0.876
C(aL, gamma) = -0.851
C(aD, cst_bgd) = -0.751
C(aL, cst_bgd) = 0.519
C(aD, center) = 0.446
C(aL, center) = -0.440
C(center, gamma) = 0.433
C(cst_bgd, center) = -0.253
----------
Result of fit for q_1.9:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 392
# data points = 34
# variables = 5
chi-square = 17.2774862
reduced chi-square = 0.59577539
Akaike info crit = -13.0165123
Bayesian info crit = -5.38470971
R-squared = 0.97324611
## Warning: uncertainties could not be estimated:
[[Variables]]
T: 540 (fixed)
aD: 84.2325105 (init = 2)
aL: 94.6015439 (init = 20)
cst_bgd: 2.31440952 (init = 0.01)
center: -3.3465e-09 (init = -1e-07)
gamma: 0.04122589 (init = 0.4)
----------
Result of fit for q_2.3:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 145
# data points = 34
# variables = 5
chi-square = 12.0183987
reduced chi-square = 0.41442754
Akaike info crit = -25.3573419
Bayesian info crit = -17.7255393
R-squared = 0.97930081
[[Variables]]
T: 540 (fixed)
aD: 63.0237050 +/- 32.8521185 (52.13%) (init = 2)
aL: 129.646254 +/- 117.148004 (90.36%) (init = 20)
cst_bgd: 2.26345187 +/- 0.27127280 (11.98%) (init = 0.01)
center: -0.00924788 +/- 0.00569512 (61.58%) (init = -1e-07)
gamma: 0.03691462 +/- 0.03594815 (97.38%) (init = 0.4)
[[Correlations]] (unreported correlations are < 0.100)
C(aL, gamma) = -0.999
C(aD, aL) = -0.998
C(aD, gamma) = 0.994
C(aD, center) = 0.955
C(aL, center) = -0.946
C(center, gamma) = 0.941
C(cst_bgd, gamma) = -0.570
C(aL, cst_bgd) = 0.534
C(aD, cst_bgd) = -0.496
C(cst_bgd, center) = -0.439
----------
Result of fit for q_3.5:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 78
# data points = 34
# variables = 5
chi-square = 13.4329791
reduced chi-square = 0.46320617
Akaike info crit = -21.5740224
Bayesian info crit = -13.9422197
R-squared = 0.96915486
[[Variables]]
T: 540 (fixed)
aD: 85.5216563 +/- 7.18528936 (8.40%) (init = 2)
aL: 44.0400079 +/- 9.27356864 (21.06%) (init = 20)
cst_bgd: 2.40985757 +/- 0.36477943 (15.14%) (init = 0.01)
center: -8.4728e-04 +/- 0.00745222 (879.55%) (init = -1e-07)
gamma: 0.12317551 +/- 0.04126871 (33.50%) (init = 0.4)
[[Correlations]] (unreported correlations are < 0.100)
C(aD, aL) = -0.995
C(aD, gamma) = 0.978
C(aL, gamma) = -0.969
C(aD, center) = 0.817
C(aL, center) = -0.817
C(center, gamma) = 0.795
C(cst_bgd, gamma) = -0.743
C(aD, cst_bgd) = -0.620
C(aL, cst_bgd) = 0.576
C(cst_bgd, center) = -0.444
----------
Result of fit for q_3.9:
[[Model]]
(Model(irf_gate) <function convolve at 0x7fe9b0304670> Model(S_meas_no_convol))
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 108
# data points = 34
# variables = 5
chi-square = 11.8142365
reduced chi-square = 0.40738746
Akaike info crit = -25.9398782
Bayesian info crit = -18.3080756
R-squared = 0.97144293
[[Variables]]
T: 540 (fixed)
aD: 84.9276971 +/- 8.52714225 (10.04%) (init = 2)
aL: 40.8223414 +/- 10.4848622 (25.68%) (init = 20)
cst_bgd: 1.91641835 +/- 0.45973725 (23.99%) (init = 0.01)
center: -0.01962833 +/- 0.00441652 (22.50%) (init = -1e-07)
gamma: 0.12484801 +/- 0.05359118 (42.93%) (init = 0.4)
[[Correlations]] (unreported correlations are < 0.100)
C(aD, aL) = -0.997
C(aD, gamma) = 0.987
C(aL, gamma) = -0.981
C(cst_bgd, gamma) = -0.864
C(aD, cst_bgd) = -0.794
C(aL, cst_bgd) = 0.763
C(aD, center) = 0.359
C(aL, center) = -0.352
C(center, gamma) = 0.347
C(cst_bgd, center) = -0.224
Storing fitted parameters
[27]:
dict_ref_params['540K'] = dict_ref_params_for_a_temperature(fit540)
dict_ref_params['540K']
[27]:
| aD | aD_err | aL | aL_err | cst_bgd | cst_bgd_err | center | center_err | gamma | gamma_err | red_chi_squared | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| q_0.6 | 158.547645 | 3.317540 | 1.135865 | 6.282286 | 0.906287 | 2.091752 | -6.820457e-08 | 0.145579 | 0.428713 | 5.403446 | 1.381108 |
| q_0.8 | 152.158615 | 3.962649 | 6.917545 | 8.195796 | 0.334830 | 2.726266 | -7.347535e-09 | 0.046932 | 0.425873 | 1.118702 | 0.669316 |
| q_1.0 | 114.190594 | 0.000000 | 205.967553 | 0.000000 | 0.920200 | 0.000000 | -1.363848e-09 | 0.000000 | 0.010000 | 0.000000 | 0.720582 |
| q_1.2 | 103.745406 | 0.000000 | 123.363272 | 0.000000 | 0.896140 | 0.000000 | -2.219353e-09 | 0.000000 | 0.023399 | 0.000000 | 0.528648 |
| q_1.7 | 121.896291 | 4.269290 | 31.885885 | 3.580051 | 0.468588 | 0.609294 | -5.523458e-03 | 0.007937 | 0.213095 | 0.068719 | 0.413621 |
| q_1.9 | 84.232511 | 0.000000 | 94.601544 | 0.000000 | 2.314410 | 0.000000 | -3.346523e-09 | 0.000000 | 0.041226 | 0.000000 | 0.595775 |
| q_2.3 | 63.023705 | 32.852118 | 129.646254 | 117.148004 | 2.263452 | 0.271273 | -9.247881e-03 | 0.005695 | 0.036915 | 0.035948 | 0.414428 |
| q_3.5 | 85.521656 | 7.185289 | 44.040008 | 9.273569 | 2.409858 | 0.364779 | -8.472765e-04 | 0.007452 | 0.123176 | 0.041269 | 0.463206 |
| q_3.9 | 84.927697 | 8.527142 | 40.822341 | 10.484862 | 1.916418 | 0.459737 | -1.962833e-02 | 0.004417 | 0.124848 | 0.053591 | 0.407387 |
Plotting \(\frac{a_D}{a_D+a_L}\) and \(\frac{a_L}{a_D+a_L}\) at 540K
[28]:
plot_ratio_coefficients(dict_ref_params['540K'])
Plotting ratio of coefficients for all temperatures
[29]:
# Define errors for each ratio
std_err_ratio_350K = np.sqrt((dict_ref_params['350K']['aD'].values**2 * dict_ref_params['350K']['aL_err'].values**2
+ dict_ref_params['350K']['aL'].values**2 * dict_ref_params['350K']['aD_err'].values**2)
/ (dict_ref_params['350K']['aD'].values + dict_ref_params['350K']['aL'].values)**4)
std_err_ratio_450K = np.sqrt((dict_ref_params['450K']['aD'].values**2 * dict_ref_params['450K']['aL_err'].values**2
+ dict_ref_params['450K']['aL'].values**2 * dict_ref_params['450K']['aD_err'].values**2)
/ (dict_ref_params['450K']['aD'].values + dict_ref_params['450K']['aL'].values)**4)
std_err_ratio_540K = np.sqrt((dict_ref_params['540K']['aD'].values**2 * dict_ref_params['540K']['aL_err'].values**2
+ dict_ref_params['540K']['aL'].values**2 * dict_ref_params['540K']['aD_err'].values**2)
/ (dict_ref_params['540K']['aD'].values + dict_ref_params['540K']['aL'].values)**4)
[30]:
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 4))
ax[0].errorbar(dict_ref_params['350K'].index,
dict_ref_params['350K']['aD'].values / (dict_ref_params['350K']['aD'].values + dict_ref_params['350K']['aL'].values),
std_err_ratio_350K,
label='350K',
marker='o',
linestyle='--')
ax[0].errorbar(dict_ref_params['450K'].index,
dict_ref_params['450K']['aD'].values / (dict_ref_params['450K']['aD'].values + dict_ref_params['450K']['aL'].values),
std_err_ratio_450K,
label='450K',
marker='x',
linestyle='--')
ax[0].errorbar(dict_ref_params['540K'].index,
dict_ref_params['540K']['aD'].values / (dict_ref_params['540K']['aD'].values + dict_ref_params['540K']['aL'].values),
std_err_ratio_540K,
label='540K',
marker='+',
linestyle='--')
ax[0].set_ylabel('$a_D/(a_D+a_L)$')
ax[0].set_ylim((0, 1))
ax[1].errorbar(dict_ref_params['350K'].index,
dict_ref_params['350K']['aL'].values / (dict_ref_params['350K']['aD'].values + dict_ref_params['350K']['aL'].values),
std_err_ratio_350K,
label='350K',
marker='o',
linestyle='--')
ax[1].errorbar(dict_ref_params['450K'].index,
dict_ref_params['450K']['aL'].values / (dict_ref_params['450K']['aD'].values + dict_ref_params['450K']['aL'].values),
std_err_ratio_450K,
label='450K',
marker='x',
linestyle='--')
ax[1].errorbar(dict_ref_params['540K'].index,
dict_ref_params['540K']['aL'].values / (dict_ref_params['540K']['aD'].values + dict_ref_params['540K']['aL'].values),
std_err_ratio_540K,
label='540K',
marker='+',
linestyle='--')
ax[1].set_ylabel('$a_L/(a_D+a_L$)')
ax[1].set_ylim((0, 0.5))
for i in range(2):
ax[i].grid()
ax[i].legend()
ax[i].set_xlabel(r'$Q (1/\AA)$')
for tick in ax[i].get_xticklabels():
tick.set_rotation(45)
fig.tight_layout()
plt.show()
Checking Q-dependence of \(\Gamma\) (width of Lorentzian)
We plot \(\Gamma\), i.e., the HWHM of the Lorentzian for all temperatures (x-axis) and all values of \(q\) (y-axis).
The colorscale is associated with the value of \(\Gamma\), which is also displayed next to the bullet point with a 2-digit precision.
[31]:
temperaturesK = [350, 450, 540]
fig = plt.figure()
ax = fig.add_subplot()
# T = 350K
scat0 = ax.scatter(
temperaturesK[0] * np.ones(len(q_tof)),
q_tof,
label='350K',
c=dict_ref_params['350K']['gamma'].values)
for i, txt in enumerate(dict_ref_params['350K']['gamma'].values):
ax.annotate(f"$\Gamma=${txt:.2f}",
(temperaturesK[0] * 1.01,
q_tof[i]))
# T = 450K
ax.scatter(
temperaturesK[1] * np.ones(len(q_tof)),
q_tof,
c=dict_ref_params['450K']['gamma'].values,
label='450K')
for i, txt in enumerate(dict_ref_params['450K']['gamma'].values):
ax.annotate(f"$\Gamma=${txt:.2f}",
(temperaturesK[1] * 1.01,
q_tof[i]))
# T = 540K
ax.scatter(
temperaturesK[2] * np.ones(len(q_tof)),
q_tof,
c=dict_ref_params['540K']['gamma'].values,
label='540K')
for i, txt in enumerate(dict_ref_params['540K']['gamma'].values):
ax.annotate(f"$\Gamma=${txt:.2f}",
(temperaturesK[2] * 1.01,
q_tof[i]))
ax.set_xlim((340, 580))
ax.grid()
ax.set_xlabel('T (K)')
ax.set_ylabel('q ($\AA^{-1}$)')
fig.colorbar(scat0, label=r'$\Gamma$')
plt.show()
\(\Gamma\) shows some \(q\)-dependences.
[ ]: