Fitting wavelength with multiple sample-detector distances
Wavelength and detector distance are correlated parameters when fitting the geometry using a single detector position. One should fix either one, or the other. The simultaneous fitting of several images taken at various detector distances has proven to alleviate this limitation (https://doi.org/10.1107/S1600577519013328).
This tutorial explains how to perform a multi-position geometry refinement, this within pyFAI using a notebook.
The dataset was recorded at the DanMAX beamline at MaxIV (Lund, Sweden) and made available by Mads Ry Jørgensen. It represents LaB6 reference material collected at 20keV at with a Pilatus detector at various distances from the sample.
Loading data
All the data are stored into a single HDF5 file following the Nexus convention. This file has been reprocessed and differs from what is acquired at the DanMAX beamline.
[1]:
%pylab inline
import numpy
from matplotlib.pyplot import subplots
from pyFAI.goniometer import GeometryTransformation, ExtendedTransformation, SingleGeometry,\
GoniometerRefinement, Goniometer
from pyFAI.calibrant import get_calibrant
import h5py
import pyFAI
from pyFAI.gui import jupyter
import time
start_time = time.perf_counter()
import logging
logging.basicConfig(level=logging.WARNING)
print(f"Running pyFAI version {pyFAI.version}")
Populating the interactive namespace from numpy and matplotlib
Running pyFAI version 0.21.0-dev1
[2]:
#Download data from internet
from silx.resources import ExternalResources
#Comment out and configure the proxy if you are behind a firewall
#os.environ["http_proxy"] = "http://proxy.company.com:3128"
downloader = ExternalResources("pyFAI", "http://www.silx.org/pub/pyFAI/gonio/", "PYFAI_DATA")
data_file = downloader.getfile("LaB6_20keV.h5")
print("file downloaded:", data_file)
file downloaded: /tmp/pyFAI_testdata_kieffer/LaB6_20keV.h5
[3]:
h5 = h5py.File(data_file)
images = h5["entry_0000/DanMAX/Pilatus/data"][()]
distances = h5["entry_0000/DanMAX/sdd/value"][()]
energy = h5["entry_0000/DanMAX/monochromator/energy"][()]
print("Distances: ", distances)
print("Energy:", energy)
LaB6 = get_calibrant("LaB6")
wavelength_guess = pyFAI.units.hc/energy*1e-10
print("Wavelength:", wavelength_guess)
LaB6.wavelength = wavelength_guess
Distances: [174. 274. 374. 474. 574. 674.]
Energy: 20
Wavelength: 6.199209921660013e-11
[4]:
fig, ax = subplots(2,3, figsize=(18,9))
for a,i in zip(ax.ravel(), images):
jupyter.display(i, ax=a)
[5]:
detector = pyFAI.detector_factory("Pilatus2MCdTe")
detector.mask = numpy.min(images, axis=0)<0
This dataset is composed of 6 images collected between 150 and 650 mm with a Pilatus 2M CdTe detector. Debye-Scherrer rings are very nicely visible and a fully automated calibration will be performed.
Automatic calibration
Since those images are pretty nice, one can read the beam-center position at (x=749, y=1573, ). The tilt and other distortion will be neglected in this first stage. We will now perform the automatic ring extraction
[6]:
geometries = {}
for dist, img in zip(distances, images):
ai = pyFAI.azimuthalIntegrator.AzimuthalIntegrator(detector=detector, wavelength=wavelength_guess)
ai.setFit2D(dist, 749, 1573)
geo = SingleGeometry(dist, img, metadata=dist, calibrant=LaB6, detector=detector, geometry=ai)
geometries[dist] = geo
[7]:
# Process the last image, the one with fewer rings:
# First extract some control points:
geo.control_points = geo.extract_cp()
# Visualization
ax = jupyter.display(sg=geo)
[8]:
print("Optimization of the geometry ... residual error is:", geo.geometry_refinement.refine2())
print(geo.geometry_refinement)
Optimization of the geometry ... residual error is: 4.946354913323504e-09
Detector Pilatus CdTe 2M PixelSize= 1.720e-04, 1.720e-04 m
Wavelength= 6.199210e-11m
SampleDetDist= 6.744518e-01m PONI= 2.710143e-01, 1.309248e-01m rot1=0.002969 rot2= 0.000522 rot3= 0.000000 rad
DirectBeamDist= 674.455mm Center: x=749.550, y=1577.712 pix Tilt=0.173 deg tiltPlanRotation= 170.022 deg
[9]:
# Extraction of the control points for all geometries:
fig, ax = subplots(2,3, figsize=(15, 10))
for a, lbl in zip(ax.ravel(), geometries):
geo = geometries[lbl]
geo.control_points = geo.extract_cp()
# print(geo.control_points)
print(f"Optimization of the geometry {lbl}, residual error is: {geo.geometry_refinement.refine2()}")
jupyter.display(sg=geo, ax=a)
a.get_legend().remove()
Optimization of the geometry 174.0, residual error is: 5.0379114709395904e-08
WARNING:matplotlib.legend:No handles with labels found to put in legend.
Optimization of the geometry 274.0, residual error is: 2.8761882965883173e-08
WARNING:matplotlib.legend:No handles with labels found to put in legend.
Optimization of the geometry 374.0, residual error is: 1.5386732553563293e-08
WARNING:matplotlib.legend:No handles with labels found to put in legend.
Optimization of the geometry 474.0, residual error is: 1.0845675708554825e-08
WARNING:matplotlib.legend:No handles with labels found to put in legend.
Optimization of the geometry 574.0, residual error is: 7.808604900480222e-09
WARNING:matplotlib.legend:No handles with labels found to put in legend.
Optimization of the geometry 674.0, residual error is: 4.313329049235334e-09
At this stage, we have 6 images and between 38 and 6 rings per image which is enough to start calibrating them all-together.
Sample stage setup
We will optimize the energy in addition to all other parameters except the rotation along the beam (rot3).
dist0represents the offset of the sample-detector stage. The associatedscale0is expected to be 1e-3 to convert milimeters in meters.poni1represents the vertical position of the center and the associatedscale1is expected to be null.poni2represents the horizontal position of the center and the associatedscale2is expected to be null.All other rotation are expected to be null as well.
[10]:
goniotrans = ExtendedTransformation(param_names = ["dist0", "scale0",
"poni1", "scale1",
"poni2", "scale2",
"rot1", "rot2",
"energy"],
dist_expr="dist0 + pos*scale0",
poni1_expr="poni1 + pos*scale1",
poni2_expr="poni2 + pos*scale2",
rot1_expr="rot1",
rot2_expr="rot2",
rot3_expr="0",
wavelength_expr="hc/energy*1e-10")
[11]:
# Starting parameters ...
param = {"dist0": 0.0,
"poni1": geo.geometry_refinement.poni1,
"poni2": geo.geometry_refinement.poni2,
"rot1": 0.0,
"rot2": 0.0,
"scale0": 0.001,
"scale1": 0.0,
"scale2": 0.0,
"energy": energy,
}
[12]:
#Defines the bounds for some variables
bounds = {"dist0": ( -0.1, 0.1),
"poni1": ( 0.0, 0.4),
"poni2": ( 0.0, 0.3),
"rot1": (-1.0, 1.0),
"rot2": (-1.0, 1.0),
"scale0": (-1.1, 1.1),
"scale1": (-1.1, 1.1),
"scale2": (-1.1, 1.1),
"energy": (energy-1,energy+1)
}
[13]:
def distance(param):
"""Since the label is directly the distance ..."""
return float(param)
assert 152.0 == distance(152)
[14]:
gonioref = GoniometerRefinement(param, # Initial guess
bounds=bounds, # Enforces constrains
pos_function=distance,
trans_function=goniotrans,
detector=detector,
wavelength=wavelength_guess)
print("Empty refinement object:", gonioref)
for lbl, geo in geometries.items():
sg = gonioref.new_geometry(str(lbl), image=geo.image, metadata=lbl, control_points=geo.control_points, calibrant=LaB6)
print(lbl, sg.get_position())
print("Populated refinement object:", gonioref)
Empty refinement object: GoniometerRefinement with 0 geometries labeled: .
174.0 174.0
274.0 274.0
374.0 374.0
474.0 474.0
574.0 574.0
674.0 674.0
Populated refinement object: GoniometerRefinement with 6 geometries labeled: 174.0, 274.0, 374.0, 474.0, 574.0, 674.0.
At this stage, the GoniometerRefinement is fully populated and can directly be optimzied:
## Optimization of all parameters (including the energy)
All optimizer available in scipy are exposed in pyFAI, the default one is slsqp which takes into account bounds and other constrains. It is very robust but not the most precise. This is why we finish with a simplex step (without bounds).
[15]:
print("Global cost after refinement:", gonioref.refine3())
Free parameters: ['dist0', 'scale0', 'poni1', 'scale1', 'poni2', 'scale2', 'rot1', 'rot2', 'energy']
Fixed: {}
fun: 5.067921390992141e-08
jac: array([-5.05846979e-06, -5.51428730e-04, 5.04780779e-06, 3.78535004e-04,
3.96370630e-08, 7.49089263e-04, -5.40181355e-07, 4.97459019e-07,
1.92552112e-07])
message: 'Optimization terminated successfully'
nfev: 251
nit: 22
njev: 22
status: 0
success: True
x: array([-3.25546066e-04, 1.00169179e-03, 2.70557111e-01, -1.67112565e-06,
1.28772296e-01, 2.15655486e-06, 1.96998599e-03, 2.96385779e-03,
2.00000888e+01])
Constrained Least square 3.514045281574116e-05 --> 5.067921390992141e-08
maxdelta on rot2: 0.0 --> 0.002963857793809446
Global cost after refinement: 5.067921390992141e-08
[16]:
# The `simplex` provides a refinement without bonds but more percise
gonioref.refine3(method="simplex")
WARNING:pyFAI.goniometer:No bounds for optimization method Nelder-Mead
Free parameters: ['dist0', 'scale0', 'poni1', 'scale1', 'poni2', 'scale2', 'rot1', 'rot2', 'energy']
Fixed: {}
/usr/lib/python3/dist-packages/scipy/optimize/_minimize.py:535: RuntimeWarning: Method Nelder-Mead cannot handle constraints nor bounds.
warn('Method %s cannot handle constraints nor bounds.' % method,
final_simplex: (array([[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913351e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913351e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01],
[-3.85949158e-04, 1.00278822e-03, 2.70525242e-01,
-7.83913352e-07, 1.28776088e-01, 2.56852858e-06,
2.33723757e-03, 2.32373650e-03, 2.00150826e+01]]), array([4.79705028e-08, 4.79705028e-08, 4.79705028e-08, 4.79705028e-08,
4.79705028e-08, 4.79705028e-08, 4.79705028e-08, 4.79705028e-08,
4.79705028e-08, 4.79705028e-08]))
fun: 4.79705028183527e-08
message: 'Maximum number of function evaluations has been exceeded.'
nfev: 1800
nit: 1144
status: 1
success: False
x: array([-3.85949158e-04, 1.00278822e-03, 2.70525242e-01, -7.83913352e-07,
1.28776088e-01, 2.56852858e-06, 2.33723757e-03, 2.32373650e-03,
2.00150826e+01])
Constrained Least square 5.067921390992141e-08 --> 4.79705028183527e-08
maxdelta on energy: 20.000088772130283 --> 20.01508263171841
[16]:
4.79705028183527e-08
[17]:
gonioref.save("table.json")
[18]:
print(open("table.json").read())
{
"content": "Goniometer calibration v2",
"detector": "Pilatus CdTe 2M",
"detector_config": {},
"wavelength": 6.194538424573843e-11,
"param": [
-0.0003859491583822855,
0.0010027882171541542,
0.27052524195268257,
-7.839133519131975e-07,
0.12877608790110912,
2.5685285779882493e-06,
0.002337237566171065,
0.002323736501912757,
20.01508263171841
],
"param_names": [
"dist0",
"scale0",
"poni1",
"scale1",
"poni2",
"scale2",
"rot1",
"rot2",
"energy"
],
"pos_names": [
"pos"
],
"trans_function": {
"content": "ExtendedTransformation",
"param_names": [
"dist0",
"scale0",
"poni1",
"scale1",
"poni2",
"scale2",
"rot1",
"rot2",
"energy"
],
"pos_names": [
"pos"
],
"dist_expr": "dist0 + pos*scale0",
"poni1_expr": "poni1 + pos*scale1",
"poni2_expr": "poni2 + pos*scale2",
"rot1_expr": "rot1",
"rot2_expr": "rot2",
"rot3_expr": "0",
"wavelength_expr": "hc/energy*1e-10",
"constants": {
"pi": 3.141592653589793,
"hc": 12.398419843320026,
"q": 1.602176634e-19
}
}
}
[19]:
gonioref.wavelength, gonioref._wavelength, 1e-10*pyFAI.units.hc/gonioref.nt_param(*gonioref.param).energy
[19]:
(6.194538424573843e-11, 6.199209921660013e-11, 6.194538424573844e-11)
[20]:
print("Ensure calibrant's wavelength has been updated:")
LaB6
Ensure calibrant's wavelength has been updated:
[20]:
LaB6 Calibrant with 152 reflections at wavelength 6.194538424573843e-11
[21]:
print(f"Total execution time: {time.perf_counter()-start_time:.3f}s")
Total execution time: 81.837s
This method ensures the wavelength has been updated in all objects important for the calibration.
At this stage, calibration has been performed on a set of points extracted automatically. Each individual image should be controled to ensure control points are homogeneously distributed along the ring. If not, those images should see their control-points re-extracted (starting from the latest/best model) and refined again. Basically, this comes down to running cells 8 and after again.
Conclusions
When the dataset is “clean”, auto-extraction of control points works out of the box
Multi-position calibration to be performed in a minute once the model is known
Energy can be refined with this methodology.