Fit - multiple mesurements
[1]:
import nexus as nx
import numpy as np
import matplotlib.pyplot as plt
# load data
data_A = np.loadtxt('spectrum_A.txt')
velocities_A = data_A[:,0]
intensity_A = data_A[:,1]
data_B = np.loadtxt('spectrum_B.txt')
velocities_B = data_B[:,0]
intensity_B = data_B[:,1]
data_C = np.loadtxt('spectrum_C.txt')
velocities_C = data_C[:,0]
intensity_C = data_C[:,1]
# setup common objects
beam = nx.Beam()
beam.Unpolarized()
# define a measurement for each temperature
# there are common fit parameters as the isomer shift and the quadrupole splitting
# they should be the same for every measurement
# so we define them independently of the measurement
isomer_for_all = nx.Var(1.1, min = 0, max = 1.2, fit = True, id = "isomer")
quadrupole_for_all = nx.Var(1.0, min = 0.5, max = 1.2, fit = True, id = "quadrupole")
# the magnetic field and the Lamb-Moessbauer factor are sample specific and have to be defined for each measurement
# ---- iron experiment at temperature A ----
site_A = nx.Hyperfine(weight = 1,
isomer = isomer_for_all,
magnetic_field = nx.Var(27, min = 15, max = 40, fit = True, id = "magnetic field A"),
quadrupole = quadrupole_for_all,
isotropic = True)
mat_iron_A = nx.Material(composition = [["Fe", 1]],
density = 7.8,
isotope = nx.lib.moessbauer.Fe57,
abundance = 0.02119,
lamb_moessbauer = nx.Var(0.81, min = 0.5, max = 1, fit = True, id = "Lamb Moessbauer A"),
hyperfine_sites = [site_A])
layer_iron_A = nx.Layer(thickness = 1000,
material = mat_iron_A)
sample_A = nx.Sample(layers = [layer_iron_A],
geometry = "f")
experiment_A = nx.Experiment(beam = beam,
objects = [sample_A],
isotope = nx.lib.moessbauer.Fe57)
moessbauer_spec_A = nx.MoessbauerSpectrum(experiment = experiment_A,
velocity = velocities_A,
intensity_data = intensity_A)
# ---- iron experiment at temperature B ----
site_B = nx.Hyperfine(weight = 1,
isomer = isomer_for_all,
magnetic_field = nx.Var(23, min = 15, max = 40, fit = True, id = "magnetic field B"),
quadrupole = quadrupole_for_all,
isotropic = True)
mat_iron_B = nx.Material(composition = [["Fe", 1]],
density = 7.8,
isotope = nx.lib.moessbauer.Fe57,
abundance = 0.02119,
lamb_moessbauer = nx.Var(0.63, min = 0.5, max = 1, fit = True, id = "Lamb Moessbauer B"),
hyperfine_sites = [site_B])
layer_iron_B = nx.Layer(thickness = 1000,
material = mat_iron_B)
sample_B = nx.Sample(layers = [layer_iron_B],
geometry = "f")
experiment_B = nx.Experiment(beam = beam,
objects = [sample_B],
isotope = nx.lib.moessbauer.Fe57)
moessbauer_spec_B = nx.MoessbauerSpectrum(experiment = experiment_B,
velocity = velocities_B,
intensity_data = intensity_B)
# ---- iron experiment at temperature C ----
site_C = nx.Hyperfine(weight = 1,
isomer = isomer_for_all,
magnetic_field = nx.Var(8, min = 0, max = 25, fit = True, id = "magnetic field C"),
quadrupole = quadrupole_for_all,
isotropic = True)
mat_iron_C = nx.Material(composition = [["Fe", 1]],
density = 7.8,
isotope = nx.lib.moessbauer.Fe57,
abundance = 0.02119,
lamb_moessbauer = nx.Var(0.6, min = 0.5, max = 1, fit = True, id = "Lamb Moessbauer C"),
hyperfine_sites = [site_C])
layer_iron_C = nx.Layer(thickness = 1000,
material = mat_iron_C)
sample_C = nx.Sample(layers = [layer_iron_C],
geometry = "f")
experiment_C = nx.Experiment(beam = beam,
objects = [sample_C],
isotope = nx.lib.moessbauer.Fe57)
moessbauer_spec_C = nx.MoessbauerSpectrum(experiment = experiment_C,
velocity = velocities_C,
intensity_data = intensity_C)
[2]:
plt.plot(velocities_A, intensity_A)
plt.plot(velocities_A, moessbauer_spec_A())
plt.xlabel('velocity (mm/s)')
plt.ylabel('intensity')
plt.show()
plt.plot(velocities_B, intensity_B)
plt.plot(velocities_B, moessbauer_spec_B())
plt.xlabel('velocity (mm/s)')
plt.ylabel('intensity')
plt.show()
plt.plot(velocities_C, intensity_C)
plt.plot(velocities_C, moessbauer_spec_C())
plt.xlabel('velocity (mm/s)')
plt.ylabel('intensity')
plt.show()
[3]:
# set up the fit
fit = nx.Fit(measurements = [moessbauer_spec_A, moessbauer_spec_B, moessbauer_spec_C])
fit.options.method = "Annealing"
fit()
Run Fit instance with id:
Starting fit with 3 measurement data set(s) and 14 fit parameter(s):
no. | id | initial value | min | max
0 | ES scaling | 955405 | 1e-16 | 9.55405e+07
1 | ES backgr | 94961.7 | 0 | 9.49617e+06
2 | Lamb Moessbauer A | 0.81 | 0.5 | 1
3 | isomer | 1.1 | 0 | 1.2
4 | magnetic field A | 27 | 15 | 40
5 | quadrupole | 1 | 0.5 | 1.2
6 | ES scaling | 954576 | 1e-16 | 9.54576e+07
7 | ES backgr | 94991.8 | 0 | 9.49918e+06
8 | Lamb Moessbauer B | 0.63 | 0.5 | 1
9 | magnetic field B | 23 | 15 | 40
10 | ES scaling | 954905 | 1e-16 | 9.54905e+07
11 | ES backgr | 95025.2 | 0 | 9.50252e+06
12 | Lamb Moessbauer C | 0.6 | 0.5 | 1
13 | magnetic field C | 8 | 0 | 25
Using 0 equality constraint(s) on parameter(s):
Using 0 inequality constraint(s).
Calling Pagmo solver with fit method Annealing
population: 150
iterations: 100
cost = 9.990074e+02
Calling ceres solver with fit method LevMar
Ceres Solver Report: Iterations: 1, Initial cost: 9.990074e+02, Final cost: 9.990074e+02, Termination: CONVERGENCE
Gradient error analysis.
Fit performed with algorithm:
Annealing
Local algorithm:
LevMar
Error analysis:
Gradient
Using 14 fit parameter(s):
no. | id | fit value | +/- std dev | initial value | min | max
0 | ES scaling | 741179 | 53.5402 | 955405 | 1e-16 | 9.55405e+07
1 | ES backgr | 246242 | 56.2298 | 94961.7 | 0 | 9.49617e+06
2 | Lamb Moessbauer A | 0.999818 | 0.0183316 | 0.81 | 0.5 | 1
3 | isomer | 0.999672 | 0.00184813 | 1.1 | 0 | 1.2
4 | magnetic field A | 33.0281 | 0.0165562 | 27 | 15 | 40
5 | quadrupole | 0.800043 | 0.00402797 | 1 | 0.5 | 1.2
6 | ES scaling | 318999 | 53.465 | 954576 | 1e-16 | 9.54576e+07
7 | ES backgr | 647259 | 56.1509 | 94991.8 | 0 | 9.49918e+06
8 | Lamb Moessbauer B | 0.999995 | 0.042474 | 0.63 | 0.5 | 1
9 | magnetic field B | 24.983 | 0.0387327 | 23 | 15 | 40
10 | ES scaling | 459683 | 52.1788 | 954905 | 1e-16 | 9.54905e+07
11 | ES backgr | 514012 | 54.7989 | 95025.2 | 0 | 9.50252e+06
12 | Lamb Moessbauer C | 0.999985 | 0.0245795 | 0.6 | 0.5 | 1
13 | magnetic field C | 12.0035 | 0.0240921 | 8 | 0 | 25
Values at boundaries:
Lamb Moessbauer A at UPPER boundary: 1
Lamb Moessbauer B at UPPER boundary: 1
Lamb Moessbauer C at UPPER boundary: 1
Correlation matrix:
no. | 0 1 2 3 4 5 6 7 8 9 10 11 12 13
-----|------------------------------------------------------------------------------------------------------------------------------
0 | 1.000 1.000 0.387 0.001 -0.000 0.002 0.000 0.000 0.000 0.000 -0.000 -0.000 -0.000 0.000
1 | 1.000 1.000 0.387 0.001 -0.000 0.002 0.000 0.000 0.000 0.000 -0.000 -0.000 -0.000 0.000
2 | 0.387 0.387 1.000 0.009 -0.002 0.006 0.000 0.000 0.000 0.000 -0.000 -0.000 -0.000 0.001
3 | 0.001 0.001 0.009 1.000 -0.007 -0.284 0.001 0.001 0.006 0.001 0.002 0.002 0.004 0.152
4 | -0.000 -0.000 -0.002 -0.007 1.000 0.003 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.000 -0.001
5 | 0.002 0.002 0.006 -0.284 0.003 1.000 -0.001 -0.001 -0.002 0.005 -0.008 -0.008 -0.023 -0.091
6 | 0.000 0.000 0.000 0.001 -0.000 -0.001 1.000 1.000 0.384 0.001 0.000 0.000 0.000 0.000
7 | 0.000 0.000 0.000 0.001 -0.000 -0.001 1.000 1.000 0.384 0.001 0.000 0.000 0.000 0.000
8 | 0.000 0.000 0.000 0.006 -0.000 -0.002 0.384 0.384 1.000 0.002 0.000 0.000 0.000 0.001
9 | 0.000 0.000 0.000 0.001 -0.000 0.005 0.001 0.001 0.002 1.000 -0.000 -0.000 -0.000 -0.000
10 | -0.000 -0.000 -0.000 0.002 -0.000 -0.008 0.000 0.000 0.000 -0.000 1.000 1.000 0.323 0.005
11 | -0.000 -0.000 -0.000 0.002 -0.000 -0.008 0.000 0.000 0.000 -0.000 1.000 1.000 0.323 0.005
12 | -0.000 -0.000 -0.000 0.004 -0.000 -0.023 0.000 0.000 0.000 -0.000 0.323 0.323 1.000 0.010
13 | 0.000 0.000 0.001 0.152 -0.001 -0.091 0.000 0.000 0.001 -0.000 0.005 0.005 0.010 1.000
Using 0 equality constraint(s) on parameter(s):
and 0 inequality constraint(s).
Total cost = 3.902e+00
cost for each FitMeasurement is:
no. | id | cost | %
0 | | 3.119e-01 | 7.992
1 | | 2.819e+00 | 72.242
2 | | 7.713e-01 | 19.766
Fit instance with id "" finished.
[4]:
plt.plot(velocities_A, intensity_A)
plt.plot(velocities_A, moessbauer_spec_A.result)
plt.xlabel('velocity (mm/s)')
plt.ylabel('intensity')
plt.show()
plt.plot(velocities_B, intensity_B)
plt.plot(velocities_B, moessbauer_spec_B.result)
plt.xlabel('velocity (mm/s)')
plt.ylabel('intensity')
plt.show()
plt.plot(velocities_C, intensity_C)
plt.plot(velocities_C, moessbauer_spec_C.result)
plt.xlabel('velocity (mm/s)')
plt.ylabel('intensity')
plt.show()
