Analysing non-spinflip polarised NR data

Currently (5/Aug/2021) refnx does not have the ability to fully analyse polarised neutron reflectometry data. However, it can be used to analyse the non-spinflip channels of a polarised neutron reflectometry measurement. Here we demonstrate how to do this using auxiliary Parameter. The datasets of interest have the structure:

Si | SiO2 | Permalloy | Au | 2-mercaptoethanol | D2O

# some necessary imports
import os.path
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import refnx
from refnx.analysis import Parameter, Objective, CurveFitter, GlobalObjective
from refnx.reflect import SLD, Slab, Structure, ReflectModel
from refnx.dataset import Data1D
from refnx._lib import flatten

# create datasets from the NSF PNR data
pth = os.path.dirname(refnx.__file__)
dd = 'c_PLP0007882.dat'
uu = 'c_PLP0007885.dat'

file_path_uu = os.path.join(pth, 'reflect', 'test', uu)
file_path_dd = os.path.join(pth, 'reflect', 'test', dd)

data_uu = Data1D(file_path_uu)
data_dd = Data1D(file_path_dd)

# create SLD (Scattering Length Density) objects for each of the materials
si = SLD(2.07, name="Si")
sio2 = SLD(3.47, name="SiO2")
au = SLD(4.66, name="Au")
mercapto = SLD(3.49, name="2-mercaptoethanol")
d2o = SLD(6.35, name="d2o")

# to describe the Permalloy layer we're going to create parameters to describe the nuclear
# and magnetic parts of the SLD. Instead of using the magnetic moment and angle we
# could just use a magnetic SLD

nuclear_py = Parameter(9.0, name="Py nuclear part")
mag_moment_py = Parameter(600, name="Py emu/cc")
angle = Parameter(0, name="angle", bounds=(0, 90.))

# Now create two SLD objects for the Permalloy layer, one for the UU channel, one for the DD channel.
# don't worry that the SLD is set to zero to start with

py_dd = SLD(0.0, name="Py DD SLD")
py_uu = SLD(0.0, name="Py UU SLD")

# Now we make constraints for the SLD objects. Each SLD object has two parameters,
# SLD.real and SLD.imag. The conversion factor of 2.85e-3 converts the magnetic moment
# from emu/cc to 10**-6 Å**-2
py_dd.real.constraint = nuclear_py - mag_moment_py * 2.85e-3 * np.cos(angle*np.pi/180)
py_uu.real.constraint = nuclear_py + mag_moment_py * 2.85e-3 * np.cos(angle*np.pi/180)

# let's check on the SLDs. Observe that the SLDs obey the constraints

complex(py_dd), complex(py_uu)

((7.29+0j), (10.71+0j))

# let's try altering the magnetic part and see if the values are updated in the SLDs.
# `magnetic_Fe` is a Parameter, and Parameter values are modified like this:

mag_moment_py.value = 700

# note how both the SLD objects are updated.

complex(py_dd), complex(py_uu)

((7.005+0j), (10.995000000000001+0j))

# Now make Slabs that describe each layer. These can either be made from SLD objects,
# or by using the `Slab` constructor directly.

# sio2 slab has a thickness of 20 and roughness of 4 with the Si fronting medium
sio2_l = sio2(20, 4)

au_l = au(215, 4)
mercapto_l = mercapto(8, 4)
d2o_l = d2o(0, 4)

# now make the Fe layers for each of the spin channels. Note that we create
# Parameter for the thickness and roughness which will be shared over both spin channels.

py_thickness = Parameter(50, name="Py thickness")
py_roughness = Parameter(5, name="Py roughness")

py_dd_l = Slab(py_thickness, py_dd, py_roughness, name="Py dd slab")
py_uu_l = Slab(py_thickness, py_uu, py_roughness, name="Py uu slab")

# now we make structures for each of the spin channels
# note that we use the same `sio2_l` for each of the structures. This
# will share the same sio2 thickness and roughness in the structures
# and will reduce parameter numbers in a fit

s_dd = si | sio2_l | py_dd_l | au_l | mercapto_l | d2o_l
s_uu = si | sio2_l | py_uu_l | au_l | mercapto_l | d2o_l

# The total number of parameters is reduced by sharing Parameter/SLD/Slab objects
# over the two structures.

# what are the number of parameters in each of the structures?
print(len(list(flatten(s_dd.parameters))), len(list(flatten(s_uu.parameters))))

# now what are the number of unique parameters in both parameter sets?
combined_set = set(flatten(s_dd.parameters)).union(set(flatten(s_uu.parameters)))
print(len(combined_set))

# this shows that the unique number of parameters over both datasets is 36, reduced from 60.
# i.e. there are parameters that are joint over both datasets

30 30
36

# now place the Structures into a ReflectModel. ReflectModel applies resolution smearing, etc.

model_dd = ReflectModel(s_dd)
model_uu = ReflectModel(s_uu)

objective_dd = Objective(model_dd, data_dd,
                         auxiliary_params=(nuclear_py, mag_moment_py, angle))
objective_uu = Objective(model_uu, data_uu,
                         auxiliary_params=(nuclear_py, mag_moment_py, angle))

global_objective = GlobalObjective([objective_dd, objective_uu])

fig, ax = plt.subplots()
plt.scatter(data_dd.x, data_dd.y, label="dd", s=4)
plt.plot(data_dd.x, objective_dd.generative())

plt.scatter(data_uu.x, data_uu.y, label="uu", s=4)
plt.plot(data_uu.x, objective_uu.generative())

plt.ylabel("R")
plt.xlabel("Q / $\AA^{-1}$")
plt.yscale('log')
plt.legend();

plt.plot(*s_dd.sld_profile(), label='dd')
plt.plot(*s_uu.sld_profile(), label="uu")
plt.ylabel("SLD")
plt.xlabel("z / $\AA$")
plt.legend();

# select the parameters to be fitted and their bounds
model_uu.scale.setp(vary=True, bounds=(0.9, 1.1))
model_uu.bkg.setp(vary=True, bounds=(1e-7, 5e-6))
model_dd.scale.setp(vary=True, bounds=(0.9, 1.1))
model_dd.bkg.setp(vary=True, bounds=(1e-7, 5e-6))

sio2_l.thick.setp(vary=True, bounds=(10, 25))
sio2_l.rough.setp(vary=True, bounds=(1, 8))

py_thickness.setp(vary=True, bounds=(38, 55))
py_roughness.setp(vary=True, bounds=(1, 8))
nuclear_py.setp(vary=True, bounds=(8.5, 9.5))
mag_moment_py.setp(vary=True, bounds=(500, 800))

au_l.thick.setp(vary=True, bounds=(200, 240))
au_l.rough.setp(vary=True, bounds=(1, 8))
au.real.setp(vary=True, bounds=(4.5, 4.66))

mercapto_l.thick.setp(vary=True, bounds=(5, 15))
mercapto_l.rough.setp(vary=True, bounds=(1, 8))
mercapto.real.setp(vary=True, bounds=(3, 4))

d2o_l.rough.setp(vary=True, bounds=(1, 8))
d2o.real.setp(vary=True, bounds=(6.2, 6.36))

fitter = CurveFitter(global_objective)

fitter.fit('differential_evolution', seed=1);

77it [00:51,  1.51it/s]

fig, ax = plt.subplots()
plt.scatter(data_dd.x, data_dd.y, label="dd", s=4)
plt.plot(data_dd.x, objective_dd.generative())

plt.scatter(data_uu.x, data_uu.y, label="uu", s=4)
plt.plot(data_uu.x, objective_uu.generative())

plt.ylabel("R")
plt.xlabel("Q / $\AA^{-1}$")
plt.yscale('log')
plt.legend();

plt.plot(*s_dd.sld_profile(), label='dd')
plt.plot(*s_uu.sld_profile(), label="uu")
plt.ylabel("SLD")
plt.xlabel("z / $\AA$")
plt.legend();

By printing out the objectives we can see what the parameters are. Here we see that the Permalloy magnetic moment is 615 emu/cc, with a nuclear SLD of \(9.15\times10^{-6}Å^{-2}\)

print(objective_dd)

________________________________________________________________________________
Objective - 140698445453824
Dataset = c_PLP0007882
datapoints = 94
chi2 = 301.3613165434196
Weighted = True
Transform = None
________________________________________________________________________________
Parameters:      None
<Parameter:'Py nuclear part', value=9.14537 +/- 0.0239, bounds=[8.5, 9.5]>
<Parameter:  'Py emu/cc'  , value=614.615 +/- 4.88 , bounds=[500.0, 800.0]>
<Parameter:    'angle'    , value=0 (fixed)  , bounds=[0.0, 90.0]>
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    , value=0.995443 +/- 0.00293, bounds=[0.9, 1.1]>
<Parameter:     'bkg'     , value=1.47825e-06 +/- 1.9e-07, bounds=[1e-07, 5e-06]>
<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'q_offset'   , value=0 (fixed)  , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'Structure - '
________________________________________________________________________________
Parameters:      'Si'
<Parameter: 'Si - thick'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - sld'   , value=2.07 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'Si - rough'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'Si - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' , value=17.1657 +/- 0.723, bounds=[10.0, 25.0]>
<Parameter: 'SiO2 - sld'  , value=3.47 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'SiO2 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - rough' , value=5.79942 +/- 1.54 , bounds=[1.0, 8.0]>
<Parameter:'SiO2 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:  'Py dd slab'
<Parameter:'Py thickness' , value=43.7966 +/- 0.138, bounds=[38.0, 55.0]>
<Parameter:'Py DD SLD - sld', value=7.39371          , bounds=[-inf, inf], constraint=<Parameter:     None      , value=7.39371 (fixed)  , bounds=None>>
<Parameter:'Py DD SLD - isld', value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'Py roughness' , value=3.81817 +/- 0.542, bounds=[1.0, 8.0]>
<Parameter:'Py dd slab - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:      'Au'
<Parameter: 'Au - thick'  , value=209.656 +/- 34.5 , bounds=[200.0, 240.0]>
<Parameter:  'Au - sld'   , value=4.5 +/- 0.00966, bounds=[4.5, 4.66]>
<Parameter:  'Au - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'Au - rough'  , value=4.88107 +/- 0.29 , bounds=[1.0, 8.0]>
<Parameter:'Au - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters: '2-mercaptoethanol'
<Parameter:'2-mercaptoethanol - thick', value=13.8574 +/- 54.2 , bounds=[5.0, 15.0]>
<Parameter:'2-mercaptoethanol - sld', value=3.26745 +/- 7.36 , bounds=[3.0, 4.0]>
<Parameter:'2-mercaptoethanol - isld', value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'2-mercaptoethanol - rough', value=5.62675 +/- 14.8 , bounds=[1.0, 8.0]>
<Parameter:'2-mercaptoethanol - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:      'd2o'
<Parameter: 'd2o - thick' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'd2o - sld'  , value=6.24009 +/- 0.000567, bounds=[6.2, 6.36]>
<Parameter: 'd2o - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'd2o - rough' , value=6.75234 +/- 8.35 , bounds=[1.0, 8.0]>
<Parameter:'d2o - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>