Pilatus on a goniometer at ID28¶
Nguyen Thanh Tra who was post-doc at ESRF-ID28 enquired about a potential bug in pyFAI in October 2016: he calibrated 3 images taken with a Pilatus-1M detector at various detector angles: 0, 17 and 45 degrees. While everything looked correct, in first approximation, one peak did not overlap properly with itself depending on the detector angle. This peak correspond to the peak in the angle of the detector, at 23.6° …
This notebook will guide you through the calibration of the goniometer setup.
Let’s first retrieve the images and initialize the environment:
[1]:
%pylab nbagg
Populating the interactive namespace from numpy and matplotlib
[2]:
import os, sys
print(sys.executable)
/users/kieffer/.venv/py37/bin/python3
[3]:
import fabio, pyFAI, os
print("Using pyFAI version:", pyFAI.version)
from os.path import basename
from pyFAI.gui import jupyter
from pyFAI.calibrant import get_calibrant
from silx.resources import ExternalResources
downloader = ExternalResources("thick", "http://www.silx.org/pub/pyFAI/testimages")
all_files = downloader.getdir("gonio_ID28.tar.bz2")
for afile in all_files:
print(basename(afile))
Using pyFAI version: 0.19.0-beta0
gonio_ID28
det130_g45_0001p.npt
det130_g0_0001p.cbf
det130_g17_0001p.npt
det130_g0_0001p.npt
det130_g45_0001p.cbf
det130_g17_0001p.cbf
There are 3 images stored as CBF files and the associated control points as npt files.
[4]:
images = [i for i in all_files if i.endswith("cbf")]
images.sort()
mask = None
fig, ax = subplots(1,3, figsize=(9,3))
for i, cbf in enumerate(images):
fimg = fabio.open(cbf)
jupyter.display(fimg.data, label=basename(cbf), ax=ax[i])
if mask is None:
mask = fimg.data<0
else:
mask |= fimg.data<0
numpy.save("mask.npy", mask)
To be able to calibrate the detector position, the calibrant used is LaB6 and the wavelength was 0.69681e-10m
[5]:
wavelength=0.6968e-10
calibrant = get_calibrant("LaB6")
calibrant.wavelength = wavelength
print(calibrant)
detector = pyFAI.detector_factory("Pilatus1M")
LaB6 Calibrant with 120 reflections at wavelength 6.968e-11
[6]:
# Define the function that extracts the angle from the filename:
def get_angle(basename):
"""Takes the basename (like det130_g45_0001.cbf ) and returns the angle of the detector"""
return float(os.path.basename((basename.split("_")[-2][1:])))
for afile in images:
print('filename', afile, "angle:",get_angle(afile))
filename /tmp/thick_testdata_kieffer/gonio_ID28.tar.bz2__content/gonio_ID28/det130_g0_0001p.cbf angle: 0.0
filename /tmp/thick_testdata_kieffer/gonio_ID28.tar.bz2__content/gonio_ID28/det130_g17_0001p.cbf angle: 17.0
filename /tmp/thick_testdata_kieffer/gonio_ID28.tar.bz2__content/gonio_ID28/det130_g45_0001p.cbf angle: 45.0
[7]:
#Define the transformation of the geometry as function of the goniometrer position.
# by default scale1 = pi/180 (convert deg to rad) and scale2 = 0.
from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer
goniotrans2d = GeometryTransformation(param_names = ["dist", "poni1", "poni2",
"rot1", "rot2",
"scale1", "scale2"],
dist_expr="dist",
poni1_expr="poni1",
poni2_expr="poni2",
rot1_expr="scale1 * pos + rot1",
rot2_expr="scale2 * pos + rot2",
rot3_expr="0.0")
[8]:
#Definition of the parameters start values and the bounds
param = {"dist":0.30,
"poni1":0.08,
"poni2":0.08,
"rot1":0,
"rot2":0,
"scale1": numpy.pi/180., # rot2 is in radians, while the motor position is in degrees
"scale2": 0
}
#Defines the bounds for some variables. We start with very strict bounds
bounds = {"dist": (0.25, 0.31),
"poni1": (0.07, 0.1),
"poni2": (0.07, 0.1),
"rot1": (-0.01, 0.01),
"rot2": (-0.01, 0.01),
"scale1": (numpy.pi/180., numpy.pi/180.), #strict bounds on the scale: we expect the gonio to be precise
"scale2": (0, 0) #strictly bound to 0
}
[9]:
gonioref2d = GoniometerRefinement(param, #initial guess
bounds=bounds,
pos_function=get_angle,
trans_function=goniotrans2d,
detector=detector,
wavelength=wavelength)
print("Empty goniometer refinement object:")
print(gonioref2d)
Empty goniometer refinement object:
GoniometerRefinement with 0 geometries labeled: .
[10]:
# Populate with the images and the control points
for fn in images:
base = os.path.splitext(fn)[0]
bname = os.path.basename(base)
fimg = fabio.open(fn)
sg =gonioref2d.new_geometry(bname, image=fimg.data, metadata=bname,
control_points=base+".npt",
calibrant=calibrant)
print(sg.label, "Angle:", sg.get_position())
print("Filled refinement object:")
print(gonioref2d)
det130_g0_0001p Angle: 0.0
det130_g17_0001p Angle: 17.0
det130_g45_0001p Angle: 45.0
Filled refinement object:
GoniometerRefinement with 3 geometries labeled: det130_g0_0001p, det130_g17_0001p, det130_g45_0001p.
[11]:
# Initial refinement of the goniometer model with 5 dof
gonioref2d.refine2()
Cost function before refinement: 0.0007718180621418161
[0.3 0.08 0.08 0. 0. 0.01745329
0. ]
fun: 4.4188913962221164e-08
jac: array([ 8.45390424e-09, 8.92671967e-08, 3.37939728e-07, -9.31373960e-08,
3.55469623e-08, -3.52535668e-03, 9.75932426e-04])
message: 'Optimization terminated successfully.'
nfev: 191
nit: 21
njev: 21
status: 0
success: True
x: array([0.28456254, 0.08940175, 0.08892749, 0.00442128, 0.0029766 ,
0.01745329, 0. ])
Cost function after refinement: 4.4188913962221164e-08
GonioParam(dist=0.28456254071696413, poni1=0.08940174752967277, poni2=0.08892749449109176, rot1=0.004421281197936293, rot2=0.0029766044781369553, scale1=0.017453292519943295, scale2=0.0)
maxdelta on: dist (0) 0.3 --> 0.28456254071696413
[11]:
array([0.28456254, 0.08940175, 0.08892749, 0.00442128, 0.0029766 ,
0.01745329, 0. ])
[12]:
# Remove constrains on the refinement:
gonioref2d.bounds=None
gonioref2d.refine2()
Cost function before refinement: 4.4188913962221164e-08
[0.28456254 0.08940175 0.08892749 0.00442128 0.0029766 0.01745329
0. ]
fun: 1.872456502089939e-08
jac: array([-1.85415310e-08, -5.84912461e-08, -3.53206487e-07, -4.17128709e-07,
2.09677793e-07, 2.08230218e-06, 1.22462904e-07])
message: 'Optimization terminated successfully.'
nfev: 124
nit: 13
njev: 13
status: 0
success: True
x: array([ 2.84547605e-01, 8.86540301e-02, 8.93070277e-02, 5.48717450e-03,
5.54901577e-03, 1.74619285e-02, -2.09889821e-05])
Cost function after refinement: 1.872456502089939e-08
GonioParam(dist=0.2845476046521533, poni1=0.0886540301310109, poni2=0.08930702773696293, rot1=0.005487174500206645, rot2=0.005549015769013916, scale1=0.017461928498376328, scale2=-2.0988982115318873e-05)
maxdelta on: rot2 (4) 0.0029766044781369553 --> 0.005549015769013916
[12]:
array([ 2.84547605e-01, 8.86540301e-02, 8.93070277e-02, 5.48717450e-03,
5.54901577e-03, 1.74619285e-02, -2.09889821e-05])
[13]:
# Check the calibration on all 3 images
fig, ax = subplots(1, 3, figsize=(9, 3) )
for idx,lbl in enumerate(gonioref2d.single_geometries):
sg = gonioref2d.single_geometries[lbl]
if sg.control_points.get_labels():
sg.geometry_refinement.set_param(gonioref2d.get_ai(sg.get_position()).param)
jupyter.display(sg=sg, ax=ax[idx])
WARNING:matplotlib.legend:No handles with labels found to put in legend.
WARNING:matplotlib.legend:No handles with labels found to put in legend.
[14]:
#Create a MultiGeometry integrator from the refined geometry:
angles = []
images = []
for sg in gonioref2d.single_geometries.values():
angles.append(sg.get_position())
images.append(sg.image)
multigeo = gonioref2d.get_mg(angles)
multigeo.radial_range=(0, 63)
print(multigeo)
# Integrate the whole set of images in a single run:
res_mg = multigeo.integrate1d(images, 10000)
ax = jupyter.plot1d(res_mg, label="multigeo")
for lbl, sg in gonioref2d.single_geometries.items():
ai = gonioref2d.get_ai(sg.get_position())
img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
ax.plot(*res, "--", label=lbl)
ax.legend()
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.19.
File "/users/kieffer/.venv/py37/lib/python3.7/site-packages/pyFAI/multi_geometry.py", line 149, in integrate1d
mask=mask, flat=flat)
MultiGeometry integrator with 3 geometries on (0, 63) radial range (2th_deg) and (-180, 180) azimuthal range (deg)
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.19.
File "<ipython-input-14-d90507300a41>", line 19, in <module>
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
[14]:
<matplotlib.legend.Legend at 0x7f539ecd2780>
[15]:
#Let's focus on the inner most ring on the image taken at 45°:
#ax.set_xlim(21.5, 21.7)
ax.set_xlim(29.0, 29.2)
ax.set_ylim(0, 5e11)
[15]:
(0, 500000000000.0)
On all three imges, the rings on the outer side of the detector are shifted in compatison with the average signal comming from the other two images. This phenomenon could be related to volumetric absorption of the photon in the thickness of the detector.
To be able to investigate this phenomenon further, the goniometer geometry is saved in a JSON file:
[16]:
gonioref2d.save("id28.json")
with open("id28.json") as f:
print(f.read())
{
"content": "Goniometer calibration v2",
"detector": "Pilatus 1M",
"detector_config": {},
"wavelength": 6.968e-11,
"param": [
0.2845476046521533,
0.0886540301310109,
0.08930702773696293,
0.005487174500206645,
0.005549015769013916,
0.017461928498376328,
-2.0988982115318873e-05
],
"param_names": [
"dist",
"poni1",
"poni2",
"rot1",
"rot2",
"scale1",
"scale2"
],
"pos_names": [
"pos"
],
"trans_function": {
"content": "GeometryTransformation",
"param_names": [
"dist",
"poni1",
"poni2",
"rot1",
"rot2",
"scale1",
"scale2"
],
"pos_names": [
"pos"
],
"dist_expr": "dist",
"poni1_expr": "poni1",
"poni2_expr": "poni2",
"rot1_expr": "scale1 * pos + rot1",
"rot2_expr": "scale2 * pos + rot2",
"rot3_expr": "0.0",
"constants": {
"pi": 3.141592653589793
}
}
}
Peak profile¶
Let’s plot the full-width at half maximum for every peak in the different intergated profiles:
[17]:
#Peak profile
from scipy.interpolate import interp1d
from scipy.optimize import bisect
def calc_fwhm(integrate_result, calibrant):
"calculate the tth position and FWHM for each peak"
delta = integrate_result.intensity[1:] - integrate_result.intensity[:-1]
maxima = numpy.where(numpy.logical_and(delta[:-1]>0, delta[1:]<0))[0]
minima = numpy.where(numpy.logical_and(delta[:-1]<0, delta[1:]>0))[0]
maxima += 1
minima += 1
tth = []
FWHM = []
for tth_rad in calibrant.get_2th():
tth_deg = tth_rad*integrate_result.unit.scale
if (tth_deg<=integrate_result.radial[0]) or (tth_deg>=integrate_result.radial[-1]):
continue
idx_theo = abs(integrate_result.radial-tth_deg).argmin()
id0_max = abs(maxima-idx_theo).argmin()
id0_min = abs(minima-idx_theo).argmin()
I_max = integrate_result.intensity[maxima[id0_max]]
I_min = integrate_result.intensity[minima[id0_min]]
tth_maxi = integrate_result.radial[maxima[id0_max]]
I_thres = (I_max + I_min)/2.0
if minima[id0_min]>maxima[id0_max]:
if id0_min == 0:
min_lo = integrate_result.radial[0]
else:
min_lo = integrate_result.radial[minima[id0_min-1]]
min_hi = integrate_result.radial[minima[id0_min]]
else:
if id0_min == len(minima) -1:
min_hi = integrate_result.radial[-1]
else:
min_hi = integrate_result.radial[minima[id0_min+1]]
min_lo = integrate_result.radial[minima[id0_min]]
f = interp1d(integrate_result.radial, integrate_result.intensity-I_thres)
tth_lo = bisect(f, min_lo, tth_maxi)
tth_hi = bisect(f, tth_maxi, min_hi)
FWHM.append(tth_hi-tth_lo)
tth.append(tth_deg)
return tth, FWHM
fig, ax = subplots()
ax.plot(*calc_fwhm(res_mg, calibrant), "o", label="multi")
for lbl, sg in gonioref2d.single_geometries.items():
ai = gonioref2d.get_ai(sg.get_position())
img = sg.image * ai.dist * ai.dist / ai.pixel1 / ai.pixel2
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
t,w = calc_fwhm(res, calibrant=calibrant)
ax.plot(t, w,"-o", label=lbl)
ax.set_title("Peak shape as function of the angle")
ax.set_xlabel(res_mg.unit.label)
ax.legend()
WARNING:pyFAI.DEPRECATION:Function _integrate1d_legacy is deprecated since pyFAI version 0.19.
File "<ipython-input-17-ae89462c1e23>", line 51, in <module>
res = ai.integrate1d(img, 5000, unit="2th_deg", method="splitpixel")
[17]:
<matplotlib.legend.Legend at 0x7f539ec4f3c8>
Conclusion:¶
Can the FWHM and peak position be corrected using raytracing and deconvolution ?