import numpy as np
import escape as esc
from escape.utils import show
import escape.utils.mdb as mdb
In this notebook we are going to fit the data published in the following article [https://doi.org/10.1103/PhysRevB.97.144511]. The given sample is a Nb(25nm)/Gd(7nm)/Nb(25nm) trilayer on a Al2O3 substrate and covered with Ta/Cu thin films. In the previous notebook about specular reflectivity we have demonstarted how the module works, where as an experimental data we took the curve generated with Parratt32 program. The idea of this notebook is to demonstrate how to work with real experimental data in terms of ESCAPE. The result obtained here could be slightly different from what is published, because we do not have all details about published model.
We first define our parameters, which we are going to fit. The names of these parameters are self-explanatory.
wl = 2.29;
mdTa=esc.par("Ta Density", 16.65, units="g/cm^3", userlim=[2, 16.8])
mdCu=esc.par("Cu Density", 8.96, units="g/cm^3", userlim=[7, 9])
mdNb1=esc.par("Nb1 Density", 8.57, units="g/cm^3", userlim=[7.5, 8.8])
mdGd=esc.par("Gd Density", 7.90, units="g/cm^3", userlim=[6.0, 8.0])
mdNb2=esc.par("Nb2 Density", 8.57, units="g/cm^3", userlim=[7, 8.8])
mdAl2O3=esc.par("Al2O3 Density", 3.95, units="g/cm^3", userlim=[3.5, 4.0], fixed=True)
thknTa = esc.par("Ta Thkn", 3.0, units="nm", userlim=[1, 5])
roughTa = esc.par("Ta Roughness", 0.1, units="nm", userlim=[0, 3])
thknCu = esc.par("Cu Thkn", 4.0, units="nm", userlim=[3, 7])
roughCu = esc.par("Cu Roughness", 0.5, units="nm", userlim=[0, 4])
thknGd = esc.par("Gd Thkn", 7.6, units="nm", userlim=[0, 10])
roughGd = esc.par("Gd Roughness", 0.5, units="nm", userlim=[0, 5])
thknNb1 = esc.par("Nb1 Thkn", 25.0, units="nm", userlim=[20, 30])
roughNb1 = esc.par("Nb1 Roughness", 1.0, units="nm", userlim=[0, 5])
thknNb2 = esc.par("Nb2 Thkn", 25.0, units="nm", userlim=[20, 30])
roughNb2 = esc.par("Nb2 Roughness", 1.0, units="nm", userlim=[0, 5])
roughAl2O3 = esc.par("Al2O3 Roughness", 1.0, units="nm", userlim=[0, 1])
Next we create materials, layers and finally the sample objects.
ta = mdb.amorphous("Ta", sldscale=1e-4, density=mdTa, source=mdb.Source.XRAYS, wavelength=wl)
cu = mdb.amorphous("Cu", sldscale=1e-4, density=mdCu, source=mdb.Source.XRAYS, wavelength=wl)
gd = mdb.amorphous("Gd", sldscale=1e-4, density=mdGd, source=mdb.Source.XRAYS, wavelength=wl)
nb1 = mdb.amorphous("Nb1", sldscale=1e-4, density=mdNb1, source=mdb.Source.XRAYS, wavelength=wl)
nb2 = mdb.amorphous("Nb2", sldscale=1e-4, density=mdNb2, source=mdb.Source.XRAYS, wavelength=wl)
al2o3=mdb.amorphous("Al2O3", sldscale=1e-4, density=mdAl2O3, source=mdb.Source.XRAYS, wavelength=wl)
ta_layer = esc.layer("Layer: Ta", ta, thknTa, roughTa)
cu_layer = esc.layer("Layer: Cu", cu, thknCu, roughCu)
gd_layer = esc.layer("Layer: Gd", gd, thknGd, roughGd)
nb1_layer = esc.layer("Layer: Nb1", nb1, thknNb1, roughNb1)
nb2_layer = esc.layer("Layer: Nb2", nb2, thknNb2, roughNb2)
sub = esc.substrate("Substrate: Al2O3", al2o3, roughAl2O3)
sample = esc.multilayer("Ta/Cu/Nb1/Gd/Nb2/Al2O3", esc.air("Air"), sub)
sample.add(ta_layer)
sample.add(cu_layer)
sample.add(nb1_layer)
sample.add(gd_layer)
sample.add(nb2_layer)
show(sample)
The contribution of resolution function is negligible, but as a demonstration we will take it into account. The beam FWHM can be also defined as a parameter. This can help, if resolution function cannot be estimated properly. Normally this paramater is fixed.
Qz=esc.var("qz")
Qz0=esc.var("qz0")
fwhm=esc.par("FWHM", 2.0e-2, userlim=[0.001, 0.02], fixed=True)
R = esc.specrefl("Specrefl", Qz, sample, "matrix")
#we integrate over the variable Qz, thus intensity is a function of I(Qz0)
R = esc.average_normal(R, fwhm, Qz, Qz0)
Below we take into account sample size. At low angles the X-rays or neutrons beam doesn't cover the whole sample leading to the reduction of intensity. Usually this correction is done for the data, but can be also added to the model as a fit parameter in the case if some characteristics are not known.
I0=esc.par("I0", 1, scale=1e7, userlim=[0.1, 10], units="Cnt")
B=esc.par("Bgr", 10, userlim=[0, 30], units="Cnt")
h=0.05 #size of a beam
L=5 # size of a sample
Qmax=esc.par("Qmax", 4*np.pi/wl*h/L*10, fixed=True, userlim=[0, 1], units="1/nm")
I0f=esc.conditional(Qz0>Qmax, esc.func("", Qz0, I0), I0/(Qmax/Qz0))
I=I0f*R+B
#Opening the experimental data and creating the data object
theta, y=np.loadtxt("data/GdNb/TriLayer/XRR/xrr2.dat", unpack=True)
qz=4*np.pi/wl*np.sin((theta/2)*np.pi/180.0)*10 #*10 for units conversion wl is in Angstroms
qz_r=qz[(qz>=0.15) & (qz<=3)]
y_r=y[(qz>=0.15) & (qz<=3)]
dobj = esc.data("Trilayer.dat", qz_r, y_r, copy=True)
#Creating the model
mobj = esc.model("Model", I, dobj, residuals_scale="q2", weight_type="data")
#show(mobj, ylog=True, xlog=False, xlabel="Q[1/nm]", ylabel="R")
#Now we can perform optimization with DE algorithm
opt = esc.diffevol("DiffEvol", mobj, popsize=7, maxiter=150,
mutation=0.5, crossover=0.5, minconv=1e-3, mincost=1e-10,
nupdate=1,
polish_final_maxiter=150, polish_candidate_maxiter=2)
opt()
show(opt, ylog=True, xlog=False, xlabel="Q[1/nm]", ylabel="R")
show(sample, xlabel="Z[nm]", ylabel="SLDRe[1/nm]", yaxis="sld0re")
opt
Name: DiffEvol Parameters number: 15 Parameter Value +- Error Units Fixed Qmax 0.54875 +- 0 1/nm 1 I0 1.0855 +- 9.114e-05 x1e+07 Cnt 0 FWHM 0.02 +- 0 arbitr. 1 Ta Thkn 3.1649 +- 0.0015669 nm 0 Ta Roughness 1.0753 +- 0.00033265 nm 0 Cu Thkn 5.1026 +- 0.0072938 nm 0 Cu Roughness 2.3498 +- 0.002391 nm 0 Nb1 Thkn 24.484 +- 0.0074052 nm 0 Nb1 Roughness 2.205 +- 0.008042 nm 0 Gd Thkn 7.6388 +- 0.0027389 nm 0 Gd Roughness 0.54613 +- 0.0024425 nm 0 Nb2 Thkn 26.771 +- 0.0027079 nm 0 Nb2 Roughness 0.58904 +- 0.0029044 nm 0 Al2O3 Roughness 0.51177 +- 0.0019026 nm 0 Bgr 1.4901e-08 +- 0 Cnt 0