The fluorescence/scattered intensity $I_f$ due to irradiation with flux $I_s$
$$ \begin{equation} I_f=I_s(Hz)\Delta t(s)\Omega(sr)\sum_{i,j}\epsilon_{i,j} c_{i,j}(sr^{-1}) \end{equation} $$where $I$ the detected intensity (sum of selected fluorescence and or scattering lines), $\Delta t$ the exposure time, $\Omega$ the solid angle of the detector, $\epsilon_{i,j}$ a product of filter transmission and detector absorbance and $c_{i,j}$ the rate of line $j$ with energy $E_j$ due to source line $i$ with energy $E_i$ (depending on sample composition).
As an example, the fluorescence rate of a flat-multilayer sample can be written as (only primary interactions)
$$ \begin{equation} \begin{split} c_{i,j}(sr^{-1})=&\frac{d\mu_{i,j}}{d\Omega}\sum_k w_{j,k}\rho_k t_k^\prime(E_i,E_j)\\ \frac{d\mu_{i,j}^{fluo}}{d\Omega} =& \frac{\mu_j(E_i)}{4\pi}\\ \frac{d\mu_{i,j}^R}{d\Omega} =& r_e^2 K_R(\phi,\theta) \frac{N_A}{M_j}f_j^2(E_i,\theta)\\ \frac{d\mu_{i,j}^C}{d\Omega} =& r_e^2 K_C(\phi,\theta) \frac{N_A}{M_j}S_j(E_i,\theta) \end{split} \end{equation} $$where $k$ loops over the layers. Note that $j$ refers to a particular interaction type (fluorescence of element $Z$, elastic or inelastic scattering). See polarization for the definition of the differential scattering cross-sections (SI units $cm^2/g/sr$ and $M_j$ the molar mass of the atom ($g/mol$), $N_A$ the Avogadro constant ($1/mol$), $f$ the atomic form factor and $S$ the incoherent scattering function of the atom).
The corrected layer thickness $t_k^\prime$ takes attenuation of primary X-rays and fluorescence/scattering into account. For a single layer in reflection geometry it can be written as
$$ \begin{equation} \begin{split} t^\prime(E_i,E_j) =& \frac{e^{\chi(E_i,E_j) t}-1}{\chi(E_i,E_j)\cos\alpha_{in}}\\ \chi(E_i,E_j) =& \rho\left(\frac{\mu(E_j)}{\cos\alpha_{out}}-\frac{\mu(E_i)}{\cos\alpha_{in}}\right) \end{split} \end{equation} $$where $\alpha$ the angle between the sample surface normal (pointing away from the source) and the incident(in) or fluorescence/scattering(out) direction ($\alpha_{out}>90^\circ$ in reflection geometry). Note that $\lim_{\chi\to\infty}t^\prime=\frac{t}{\cos\alpha_{in}}$.
See the notebook on diodes on how $I_s$ is measured. We will assume the detector has a centric-cone geometry with solid angle
$$ \begin{equation} \Omega=2\pi\left(1-\frac{x+d_0}{\sqrt{\frac{A}{\pi}+\left(x+d_0\right)^2}}\right) \end{equation} $$where $A(mm^2)$ the active area of the detector, $x(mm)$ the position of the detector and $d_0(mm)$ the distance to the sample for $x=0$. To determine $A$ and $d_0$ we can measure the fluorescence of any sample as function of $x$:
$$ \begin{equation} I_f(x,c,d_0,A)=c\Omega(x,d_0,A) \end{equation} $$As an illustration we will define a detector geometry and multilayer sample. A thin-film standard is used here but any other material can be considered:
In [60]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
from spectrocrunch.materials import xrfstandards
from spectrocrunch.detectors import xrf as xrfdetectors
from spectrocrunch.geometries import xrf as xrfgeometries
from spectrocrunch.sources import xray as xraysources
source = xraysources.factory("synchrotron")
detector = xrfdetectors.factory("leia")
geometry = xrfgeometries.factory("sxm120",detectorposition=-10,positionunits="mm",\
detector=detector,source=source)
addnoise = True # add noise to simulations done below
method = "fisx" # XRF simulation method
realsample = xrfstandards.factory("RF7-200-S2371-03",geometry=geometry,\
filmthickness=10e-7) # 10 nm
print(realsample)
Simulate a detector scan $I_f(x,c,d_0,A)$:
In [61]:
from spectrocrunch.common import units
from spectrocrunch.math import noisepropagation
from spectrocrunch.materials import pymca
# Geometry at which the data is collected
geometry.zerodistance = units.Quantity(5.,"cm")
detector.activearea = units.Quantity(70.,"mm^2")
print("\nTheoretical geometry:")
print(" Zero-distance: {:~}".format(geometry.zerodistance.to("cm")))
print(" Active area: {:~}".format(detector.activearea.to("mm^2")))
# Simulate measurement at current distance
energy = 7.3
flux = 1e9
time = 5
pymcahandle = pymca.PymcaHandle(sample=realsample,energy=energy,flux=flux,time=time,\
linear=True,escape=False,continuum=False,scatter=False)
mcaref = pymcahandle.mca(histogram=True,scattering=False,method=method)
# Simulate detector scan
n = 100
x = units.Quantity(np.linspace(-20,60,n),"mm")
I0 = np.full(n,flux*time)
solidangle = geometry.detector.solidangle_calc(activearea=detector.activearea,distance=x+geometry.zerodistance)
fluo = mcaref.sum()/geometry.solidangle*solidangle
if addnoise:
I0 = np.random.poisson(np.round(I0).astype(int))
fluo = np.random.poisson(np.round(fluo).astype(int))
fig,axs = plt.subplots(1,2,figsize=(12,5))
u = x.units
plt.sca(axs[0])
plt.plot(x,fluo/I0.astype(float))
xref = geometry.detectorposition.to(u).magnitude
iref = mcaref.sum()/(flux*time)
lines = plt.plot([xref,xref,x[0].magnitude],[0,iref,iref])
color = lines[0].get_color()
plt.ylabel("Normalized fluorescence")
plt.xlabel("Motor position ({:~})".format(u))
plt.sca(axs[1])
plt.plot(mcaref,color=color)
plt.gca().set_yscale('log', basey=10)
plt.xlim([0,len(mcaref)-1])
plt.ylim([1,np.max(mcaref)*1.1])
plt.ylabel("ph/channel")
plt.xlabel("MCA channels")
plt.title("\nSpectrum at x ={:~}:".format(geometry.detectorposition.to("mm")))
plt.show()
Calibrate the geometry (starting from different values as the ones used to simulate the data):
In [62]:
# Calibration resources
intensities = noisepropagation.poisson(fluo)/noisepropagation.poisson(I0)
calibrc = {"signal":noisepropagation.E(intensities),\
"var":noisepropagation.VAR(intensities),\
"detectorposition":x.magnitude,\
"positionunits":x.units}
# Calibrate the geometry (starting from wrong values)
geometry.zerodistance += units.Quantity(-5.,"cm")
detector.activearea += units.Quantity(10.,"mm^2")
print("\nInitial geometry:")
print("Zero-distance: {:~}".format(geometry.zerodistance.to("cm")))
print("Active area: {:f~}".format(detector.activearea.to("mm^2")))
geometry.calibrate(calibrc=calibrc,plot=True,fit=True,fixedactivearea=False)
plt.show()
print("Calibrated geometry:")
print("Zero-distance: {:~}".format(geometry.zerodistance_rv.to("cm")))
print("Active area: {:f~}".format(detector.activearea_rv.to("mm^2")))
The correlation between $c$, $A$ and $d_0$ is too high to provide a usable result (also when fixing the active area).
To solve the correlation issue, we determine $\Omega_{ref}$ at a particular motor position $x=x_{ref}$ by fitting the fluorescence spectrum of a standard measured with the detector at this position:
$$ \begin{equation} \Omega_{ref}=\frac{I_{f}}{I_s(Hz)\Delta t(s)\sum_{i,j}\epsilon_{i,j} c_{i,j}(sr^{-1})} \end{equation} $$where $I_s$, $\Delta t$, $\epsilon_{i,j}$ and $c_{i,j}$ are assumed to be known (flux measured by calibrated diodes, known sample and filter composition).
This provides a fixed relationship between $A$ and $d_0$ which can be substituted in the expression used for calibrating the geometry
$$ \begin{equation} \begin{split} I_f(x,c,d_0)=&c\Omega(x,d_0,A(d_0))\\ A(d_0)=&\pi\left(\frac{\left(x_{ref}+d_0\right)^2}{\left(1-\frac{\Omega_{ref}}{2\pi}\right)^2}-\left(x_{ref}+d_0\right)^2\right) \end{split} \end{equation} $$When using a thin-film standard, the thickness and density of the film are unknown but the areal densities of the elements in the film are known. For elements only present in the film and assuming absorption and other secondary effects are negligible, we can write
$$ \begin{equation} \begin{split} c_{i,j}^{film}=&\frac{d\mu_{i,j}}{d\Omega} w_{j}^{film}\rho_{film} t_{film}^{\prime}(E_i,E_j)\\ \approx&\frac{1}{\cos\alpha_{in}}\frac{d\mu_{i,j}}{d\Omega} w_{Z}^{film}\rho_{film} t_{film}\\ =&\frac{1}{\cos\alpha_{in}}\frac{d\mu_{i,j}}{d\Omega}\rho_{Z,A}^{film} \end{split} \end{equation} $$Hence for elements in the thin-film, it is enough to know their areal densities. In practice however we use mass fractions calculated from the areal densities using the density and the thickness of the substrate. The mass fractions obtained are physically meaningless but valid for the purpose of calculating $\Omega_{ref}$.
For elements in the substrate, density and thickness need to be known if self-absorption is not non-negligible:
$$ \begin{equation} t_{subs}^{\prime}(E_i,E_j)\neq \frac{t_{subs}}{\cos\alpha_{in}} \end{equation} $$Simulate and fit an XRF spectrum of a thin-film standard (simulation and fit are done with a different $d_0$ and $A$; scattering, escape and sum peaks are omitted):
In [63]:
# Thin film standards have an unknown film thickness and density,
# only the areal densities of the different elements and the
# composition and thickness of the substrate are known.
thinfilmapprox = True
# Geometry at which the data is collected
geometry.zerodistance = units.Quantity(5.,"cm")
detector.activearea = units.Quantity(70.,"mm^2")
print("\nTheoretical geometry:")
print(" Zero-distance: {:~}".format(geometry.zerodistance_rv.to("cm")))
print(" Active area: {:~}".format(detector.activearea_rv.to("mm^2")))
# Simulate measurement (use the sample with known film thickness)
energy = 7.3
flux = 1e9
time = 1000 # the sum spectrum of a 2D map
pymcahandle = pymca.PymcaHandle(sample=realsample,energy=energy,flux=flux,time=time,\
linear=True,escape=False,continuum=False,scatter=False)
mca = pymcahandle.mca(histogram=True,scattering=False,method=method)
if method=="fisx":
mca1 = mca
mca2 = pymcahandle.mca(histogram=True,scattering=False,method="analytical")
if addnoise:
mca = np.random.poisson(np.round(mca).astype(int))
# Calibrate with unknown film thickness
if thinfilmapprox:
thinfilmsample = xrfstandards.factory("RF7-200-S2371-03",geometry=geometry)
pymcahandle.sample = thinfilmsample
# Initialize fit with the wrong geometry
pymcahandle.setdata(mca)
geometry.zerodistance += units.Quantity(5.,"cm")
detector.activearea += units.Quantity(-10.,"mm^2")
pymcahandle.addtopymca(fresh=True)
# Adapt config manually if needed:
#config = pymcahandle.mcafit.getConfiguration()
#config["fit"]["stripflag"] = 0
#...
#pymcahandle.mcafit.configure(config)
# Perform fit
fitresult = pymcahandle.fit()
# Print errors
def strwerror(e,wfrac,exwfrac):
error = (wfrac-exwfrac)/exwfrac
return " {}: {:6.02f} wt% (expected: {:6.02f} wt%, error: {:.02f}%)".\
format(e,wfrac*100,exwfrac*100,error*100)
def straderror(e,ad,exad):
error = (ad-exad)/exad
return " {}: {:6.02f} ng/mm^2 (expected: {:6.02f} ng/mm^2, error: {:.02f}%)".\
format(e,ad*1e7,exad*1e7,error*100)
def printerrors(fitresult,sample):
out = {}
if thinfilmapprox:
exarealdensities = sample.arealdensity()
rho = sample[0].density
t = sample[0].thickness
for k,wfrac in fitresult["massfractions"].items():
element = k.element
ad = wfrac*rho*t
if element in exarealdensities:
exad = exarealdensities[element]
exwfrac = exad/(rho*t)
out[element] = {"ad":straderror(element,ad,exad),\
"wfrac":strwerror(element,wfrac,exwfrac)}
else:
exarealdensities = sample.arealdensity()
arealdensities = {}
massfractions = {}
exmassfractions = {}
exarealdensities = {}
for layer,wfracs in zip(sample,fitresult["lmassfractions"]):
rho = layer.density
t = layer.thickness
exwfracs = layer.elemental_massfractions()
exad = layer.arealdensity()
for k,wfrac in wfracs.items():
if wfrac!=0:
element = k.element
arealdensities[k] = wfrac*rho*t
massfractions[k] = wfrac
exmassfractions[k] = exwfracs[element]
exarealdensities[k] = exad[element]
for k,wfrac in massfractions.items():
if k in exmassfractions:
element = k.element
exwfrac = exmassfractions[k]
exad = exarealdensities[k]
ad = arealdensities[k]
out[element] = {"ad":straderror(element,ad,exad),\
"wfrac":strwerror(element,wfrac,exwfrac)}
print(" Mass fractions and areal densities (within one layer):")
for k in out:
print out[k]["wfrac"]
print out[k]["ad"]
print("\nFitted vs. theory (before geometry calibration):")
printerrors(fitresult,pymcahandle.sample)
# Plot fit
def plotfit(fitresult):
plt.plot(fitresult["energy"],fitresult["y"],label='data')
plt.plot(fitresult["energy"],fitresult["yfit"],label='pymca fit')
backfunc = fitresult["interpol_energy"](fitresult["yback"])
plt.plot(fitresult["energy"],backfunc(fitresult["energy"]),label='background')
plt.gca().set_yscale('log', basey=10)
plt.ylim([1,np.max(fitresult["y"])*1.1])
plt.ylabel("ph/channel")
plt.xlabel("Energy (keV)")
plotfit(fitresult)
plt.show()
if method=="fisx":
plt.plot(mca1,label="fisx")
plt.plot(mca2,label="xraylib")
plt.gca().set_yscale('log', basey=10)
plt.ylim([1,np.max(mca)*1.1])
plt.ylabel("ph/ch")
plt.xlabel("Channels")
plt.legend()
plt.show()
Determine $\Omega_{ref}$ by comparing the fitted and theoretical fluorescence intensities:
In [64]:
caliblines = ["Ca"]
useline = lambda k: any(str(k).startswith(e) for e in caliblines)
# rate = Ifluo/I0 with I0 = flux * time
Rfit = {k:v for k,v in fitresult["fitrates"].items() if useline(k)}
Rinit = {k:v for k,v in fitresult["rates"].items() if useline(k)}
if thinfilmapprox:
# for an element within the film:
# - pymca mass fraction = 1
# - substrate density and thicknes
rho = pymcahandle.sample[0].density
t = pymcahandle.sample[0].thickness
arealdensities = pymcahandle.sample.arealdensity()
substrate = pymcahandle.sample[0].elements
for k in Rinit:
el = k.element
if el not in substrate:
Rinit[k] *= arealdensities[el]/(rho*t)
solidangleref = geometry.solidangle * sum(Rfit.values())/sum(Rinit.values())
Calibrate the geometry ($d_0$ and $A$) with a known $[\Omega_{ref},x_{ref}]$ pair as constraint:
In [65]:
geometry.calibrate(calibrc=calibrc,solidanglecalib=solidangleref,\
plot=True,fit=True,fixedactivearea=False)
# Force to real values for testing:
#geometry.zerodistance = units.Quantity(5,"cm")
#detector.activearea = units.Quantity(70,"mm^2")
print("\nCalibrate geometry using {}:".format(caliblines))
print(" Zero-distance: {:~}".format(geometry.zerodistance_rv.to("cm")))
print(" Active area: {:~}".format(detector.activearea_rv.to("mm^2")))
print("\nCurrent distance:")
print " Motor position = {:~}".format(geometry.detectorposition)
print(" Distance: {:~}".format(geometry.distance_rv.to("cm")))
The correlation between the two unknowns $c$ and $d_0$ is low enough to provide estimates of $d_0$ and $A$ with acceptable uncertainty. Known and fitted areal densities should be the same:
In [66]:
pymcahandle.addtopymca(fresh=False)
fitresult = pymcahandle.fit()
print("\nFitted vs. theory (after geometry calibration):")
printerrors(fitresult,pymcahandle.sample)
plt.figure(figsize=(12,5))
plotfit(fitresult)
spectrum = realsample.xrayspectrum(energy,emin=1,emax = energy+0.5,scattering=False,method=method)
matplotlib.rcParams.update({'font.size': 15})
spectrum.plot(histogram=True,decompose=True,fluxtime=pymcahandle.I0,\
legend=False,forcelines=True)
matplotlib.rcParams.update({'font.size': 14})
plt.show()
Errors originate from least-squares fitting, the thin-film approximation (if enabled) and discrepancies between the cross-sections used to similate fluorescence and the ones used for fitting (if method not "fisx").