fit sample
[1]:
import nexus as nx
import numpy as np
import matplotlib.pyplot as plt
load data
[2]:
data = np.loadtxt('sample_spectrum.dat')
processed_data = data
velocities = np.loadtxt("velocity_calibration.dat")
[3]:
plt.plot(data)
plt.xlabel("data point")
plt.ylabel("counts")
plt.show()
print('number of data points: {}'.format(len(data)))
number of data points: 1024
data folding
[4]:
# use same settings as for velocity calibration
# in verison 1
#processed_data, lag = nx.data.AutoFold(processed_data,
# flip = 'right', # determies which part of the spectrum is fliped
# factor = 100, # interpolation on a 100 times finer grid
# method = "cubic") # or cubic, PChip, Akima
# since version 2
processed_data, lag, additional_data = nx.data.AutoFold(processed_data,
flip = 'right', # determies which part of the spectrum is fliped
factor = 100, # interpolation on a 100 times finer grid
method = "cubic") # or cubic, PChip, Akima
print('lag: {}'.format(lag))
lag: -1.06
normalization
[5]:
processed_data, norm_factor = nx.data.Normalize(data = processed_data,
method = "baseline", # 'max', 'min', 'value'
# value = 0, # only used when method is value
# left_point=20,
# right_point=250,
# poly_order = 1 # use 1 in case you want to correct a tilted baseline
)
[6]:
plt.plot(velocities, processed_data)
plt.xlabel("velocitiy (mm/s)")
plt.ylabel("data point")
plt.show()
setup experiment
hyperfine sites
[7]:
# we assume to have two sites and an isotropic distribution (e.g. powder sample)
site1 = nx.Hyperfine(id = "site 1",
weight = nx.Var(1, min = 0, max = 1, fit = True, id = "weight 1"),
isomer = nx.Var(0, min = -1, max = 0, fit = True, id = "isomer 1"),
magnetic_field = nx.Var(15, min = 10, max = 25, fit = True, id = "magnetic field 1"),
quadrupole = nx.Var(0, min = -2, max = 2, fit = True, id = "quadrupole 1"),
# texture = 1, # only used if isotropic = False
# broadening = 1 # model free broadening parameter
isotropic = True)
site2 = nx.Hyperfine(id = "site 2",
weight = nx.Var(1, min = 0, max = 1, fit = True, id = "weight 2"),
isomer = nx.Var(0, min = -1, max = 1, fit = True, id = "isomer 2"),
magnetic_field = 0,
quadrupole = nx.Var(0, min = -1, max = 1, fit = True, id = "quadrupole 2"),
isotropic = True)
sample and experiment
always fit the sample thickness for a Moessabuer experiment
it serves as the scaling factor for the absorption
[8]:
# always fit the sample thickness, it serves as the scaling factor for the absorption
sample = nx.SimpleSample(id = "sample",
thickness = nx.Var(3000, min=2000, max=4000, fit=True, id="t"),
composition = [["Fe", 1]],
density = 7.874,
isotope = nx.lib.moessbauer.Fe57,
abundance = 0.02119,
lamb_moessbauer = 0.796,
hyperfine_sites = [site1, site2]
)
beam = nx.Beam(id = "upolarized beam",
polarization = 0,
# mixing_angle = 0, # ony for polarized beams
# canting_angle = 0 # only for polarized beams
# profile = "g" # only for grazing incidence
# fwhm = 0 # only for grazing incidence
)
#beam.Unpolarized() # not needed as polarization has been set to zero
exp = nx.Experiment(id = "experiment",
beam = beam,
objects = [sample],
isotope = nx.lib.moessbauer.Fe57)
spectrum = nx.MoessbauerSpectrum(id= "Moessbauer spectrum",
experiment = exp,
velocity = velocities,
intensity_data = processed_data,
# scaling = "auto" # use a Var in csae you want to specify this value
# background = "auto" # use a Var in csae you want to specify this value
# residual = nx.lib.residual.Sqrt() # in case you want to change to Gaussian statistics use nx.lib.residual.StdDev()
# kernel_type = "Lorentz" # is standard for Moessbauer spectrum
# resolution = 2.197) # value from calibration file or fit value as Var object # only in version 1
instrument = nx.lib.instrument.Gaussian(2.197) # since version 2
)
intensity = spectrum.Calculate()
spectrum.Plot(velocity = True
# name=None # provide a filename string to save under this filename
# data=True,
# residuals=True,
# datacolor='black'
# ... many more style options, see API
# legend=True,
# errors=False,
# errorcap=2
)
setup fit
[9]:
# create a fit object with a list of measurements to be fit in parallel.
fit = nx.Fit(id = "fit sample",
measurements = [spectrum],
# external_fit_variables = [] # for advanced fitting
# inequalities = None # for advanced fitting
)
fit.options.method = "PagmoDiffEvol" # or ""LevMar", "Annealing", "NaturalEvol", ...
# many fit.options area vailable, see API
# run the fit
fit.Evaluate()
# the fit creates an output file "fit_id_fit sample.txt" with the results
Run Fit instance with id: fit sample
Starting fit with 1 measurement data set(s) and 10 fit parameter(s):
no. | id | initial value | min | max
0 | ES scaling | 1.00599 | 1e-16 | 100.599
1 | ES backgr | 0.0992719 | 0 | 9.92719
2 | t | 3000 | 2000 | 4000
3 | weight 1 | 0.5 | 0 | 1
4 | isomer 1 | 0 | -1 | 0
5 | magnetic field 1 | 15 | 10 | 25
6 | quadrupole 1 | 0 | -2 | 2
7 | weight 2 | 0.5 | 0 | 1
8 | isomer 2 | 0 | -1 | 1
9 | quadrupole 2 | 0 | -1 | 1
Using 0 equality constraint(s) on parameter(s):
Using 0 inequality constraint(s).
Calling Pagmo solver with fit method PagmoDiffEvol
population: 110
iterations: 100
cost = 2.520352e-03
Calling ceres solver with fit method LevMar
Ceres Solver Report: Iterations: 7, Initial cost: 2.520352e-03, Final cost: 7.551236e-04, Termination: CONVERGENCE
Gradient error analysis.
Fit performed with algorithm:
PagmoDiffEvol
Local algorithm:
LevMar
Error analysis:
Gradient
Using 10 fit parameter(s):
no. | id | fit value | +/- std dev | initial value | min | max
0 | ES scaling | 0.753773 | 0.00734606 | 1.00599 | 1e-16 | 100.599
1 | ES backgr | 0.380518 | 0.00596679 | 0.0992719 | 0 | 9.92719
2 | t | 3993.38 | 1.84941e-06 | 3000 | 2000 | 4000
3 | weight 1 | 0.487548 | 0.00510111 | 0.5 | 0 | 1
4 | isomer 1 | -0.497012 | 0.00304509 | 0 | -1 | 0
5 | magnetic field 1 | 20.0165 | 0.0216377 | 15 | 10 | 25
6 | quadrupole 1 | 1.39912 | 0.00898236 | 0 | -2 | 2
7 | weight 2 | 0.512452 | 0.00485321 | 0.5 | 0 | 1
8 | isomer 2 | 0.725023 | 0.00206514 | 0 | -1 | 1
9 | quadrupole 2 | 0.792529 | 0.00396664 | 0 | -1 | 1
Values at boundaries:
Correlation matrix:
no. | 0 1 2 3 4 5 6 7 8 9
-----|------------------------------------------------------------------------------------------
0 | 1.000 -1.000 1.000 0.303 0.005 0.014 -0.005 -0.303 -0.026 0.068
1 | -1.000 1.000 -1.000 -0.302 -0.005 -0.014 0.005 0.302 0.026 -0.068
2 | 1.000 -1.000 1.000 0.315 0.004 0.012 -0.003 -0.315 -0.027 0.068
3 | 0.303 -0.302 0.315 1.000 -0.066 -0.040 0.072 -1.000 -0.046 0.075
4 | 0.005 -0.005 0.004 -0.066 1.000 0.314 -0.374 0.066 0.205 -0.084
5 | 0.014 -0.014 0.012 -0.040 0.314 1.000 -0.243 0.040 0.155 -0.062
6 | -0.005 0.005 -0.003 0.072 -0.374 -0.243 1.000 -0.072 -0.234 0.096
7 | -0.303 0.302 -0.315 -1.000 0.066 0.040 -0.072 1.000 0.046 -0.075
8 | -0.026 0.026 -0.027 -0.046 0.205 0.155 -0.234 0.046 1.000 -0.029
9 | 0.068 -0.068 0.068 0.075 -0.084 -0.062 0.096 -0.075 -0.029 1.000
Using 0 equality constraint(s) on parameter(s):
and 0 inequality constraint(s).
Total cost = 2.950e-06
cost for each FitMeasurement is:
no. | id | cost | %
0 | Moessbauer spectrum | 2.950e-06 | 100.000
Fit instance with id "fit sample" finished.
[10]:
spectrum.Plot(velocity=True)
# in case the fit does not look good:
# restart when a global mehtod has been used (as here), maybe change parameter boundaries
# when a local method has been used, change start parameter (maybe also parameter boundaries)
plot individual contributions
[11]:
plt.plot(velocities, processed_data, '.', color='black')
#get intensities from individual components
intensity_sites = nx.tools.GetSiteSpectra(spectrum, sample)
# intensity_sites is a 3D array with first dimension being the layer, the second the site of the layer, and the 3rd the intensity
print(np.array(intensity_sites).shape) # 1 layer, 2 sites, 512 data points
# go through intensity and layers of the sample in parallel to get intensity and id of the site
for intens_lay, lay in zip(intensity_sites, sample.layers):
for intens_site, site in zip(intens_lay, lay.material.hyperfine_sites):
plt.plot(velocities, intens_site, label = site.id)
plt.fill_between(velocities, intens_site, intens_site[0], alpha = 0.3)
plt.plot(velocities, spectrum.result, color='red', label = "total fit")
plt.xlabel("norm. intensity")
plt.ylabel("velocitiy (mm/s)")
plt.legend()
plt.show()
(1, 2, 512)
area ratios of intensities
[12]:
area_sites = nx.tools.AreaSites(spectrum,
sample = sample,
norm = True)
print(np.array(area_sites).shape) # 1 layer, 2 sites
print("")
print("Relative areas of the sites")
print(area_sites)
(1, 2)
Relative areas of the sites
[[0.49816666 0.50183334]]
[ ]: