MX Calibrate

Calibrate a translation table from a set of powder diffraction images taken at various sample-detector distances.

This is a notebook replacement of the MX-calibrate tool from pyFAI with advanced features.

Start with some constant definition:

[1]:
calibrant_name = "CeO2"
detector_name = "Pilatus 2M"
file_pattern = "massif1/test-powder*.cbf"
result_file = "MX-calibrate.json"
wavelength = None # set a value to override the one in the headers
[2]:
!wget http://www.silx.org/pub/pyFAI/massif1.tar.bz2
!tar -xvjf massif1.tar.bz2
--2024-09-12 12:30:57--  http://www.silx.org/pub/pyFAI/massif1.tar.bz2
Resolving www.silx.org (www.silx.org)... 195.154.237.27
Connecting to www.silx.org (www.silx.org)|195.154.237.27|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 6784503 (6.5M) [application/x-bzip2]
Saving to: ‘massif1.tar.bz2.1’

massif1.tar.bz2.1   100%[===================>]   6.47M  22.4MB/s    in 0.3s

2024-09-12 12:30:57 (22.4 MB/s) - ‘massif1.tar.bz2.1’ saved [6784503/6784503]

massif1/
massif1/test-powder_5_0001.poni
massif1/test-powder_4_0001.cbf
massif1/test-powder_6_0001.poni
massif1/test-powder_8_0001.cbf
massif1/test-powder_7_0001.cbf
massif1/test-powder_7_0001.poni
massif1/test-powder_3_0001.cbf
massif1/test-powder_1_0001.poni
massif1/test-powder_4_0001.poni
massif1/test-powder_8_0001.poni
massif1/test-powder_2_0001.poni
massif1/test-powder_6_0001.cbf
massif1/test-powder_2_0001.cbf
massif1/test-powder_5_0001.cbf
massif1/test-powder_3_0001.poni
massif1/test-powder_1_0001.cbf
[3]:
%matplotlib widget
#inline2
[4]:
import os
import glob
import logging
import numpy
from matplotlib.pyplot import subplots
from scipy.stats import linregress
import fabio
import pyFAI
from pyFAI.gui import jupyter
import pyFAI.calibrant
from pyFAI.gui.jupyter.calib import Calibration
from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement
from pyFAI.gui.cli_calibration import AbstractCalibration
import pyFAI.gui.mpl_calib
pyFAI.gui.mpl_calib.logger.setLevel(logging.ERROR)
print(f"Running pyFAI version {pyFAI.version}")
WARNING:pyFAI.gui.matplotlib:matplotlib already loaded, setting its backend may not work
Running pyFAI version 2024.9.0-dev0
[5]:
detector = pyFAI.detector_factory(detector_name)
calibrant = pyFAI.calibrant.get_calibrant(calibrant_name)
files = sorted(glob.glob(file_pattern))
print("input files: "+" ".join(files))
input files: massif1/test-powder_1_0001.cbf massif1/test-powder_2_0001.cbf massif1/test-powder_3_0001.cbf massif1/test-powder_4_0001.cbf massif1/test-powder_5_0001.cbf massif1/test-powder_6_0001.cbf massif1/test-powder_7_0001.cbf massif1/test-powder_8_0001.cbf
[6]:
first = fabio.open(files[0])

def get_dectris_headers(fimg):
    """return the dectris headers from a Pilatus detector"""
    res = {}
    for line in fimg.header.get("_array_data.header_contents", "").split("\n"):
        words = line.split()
        if len(words)>=3:
            key = words[1]
            for v in words[2:]:
                try:
                    vf = float(v)
                except:
                    continue
                if not("." in v or "e" in v):
                    vf = int(v)
                res[key] = vf
    return res

get_dectris_headers(first)
[6]:
{'Silicon': 0.00045,
 'Pixel_size': 0.000172,
 'N_oscillations': 1,
 'Chi': 0.0,
 'Phi': 0.0,
 'Kappa': 0.0,
 'Alpha': 0.0,
 'Polarization': 0.99,
 'Detector_2theta': 0.0,
 'Angle_increment': 1.0,
 'Transmission': 100.0,
 'Flux': 436215830143.2828,
 'Detector_Voffset': 0.0,
 'Detector_distance': 0.126474,
 'Wavelength': 0.965459,
 'N_excluded_pixels:': 321,
 'Threshold_setting': 6421,
 'Count_cutoff': 1048500,
 'Tau': 0,
 'Exposure_period': 0.02115,
 'Exposure_time': 0.02,
 'Start_angle': 0.0}
[7]:
if wavelength is None:
    wavelength = get_dectris_headers(first)["Wavelength"] * 1e-10
calibrant.wavelength = wavelength
[8]:
#apply mask to the detector
mask = numpy.logical_or(detector.mask, first.data<0)
detector.mask = mask

Manual calibration of the first image

[9]:
# Important: select the ring number before right-click on the ring. Finally click on the refine button
calib = Calibration(img=first.data,
                    mask=mask,
                    detector=detector,
                    wavelength=wavelength,
                    calibrant=calibrant)
[10]:
calib.extract_cpt()
# calib.geoRef.rot1 = calib.geoRef.rot2 = calib.geoRef.rot3 = 0
calib.refine(fixed=["wavelength", "rot3"])
Before refinement, the geometry is:
Detector Pilatus 2M      PixelSize= 172µm, 172µm         BottomRight (3)
Wavelength= 9.654590e-11 m
SampleDetDist= 1.003811e-01 m   PONI= -3.757365e-04, -1.204715e-05 m    rot1=-0.000164  rot2=-0.000283  rot3=0.000000 rad
DirectBeamDist= 100.381 mm      Center: x=0.026, y=-2.349 pix   Tilt= 0.019° tiltPlanRotation= -59.866° 𝛌= 0.965Å
Detector Pilatus 2M      PixelSize= 172µm, 172µm         BottomRight (3)
Wavelength= 9.654590e-11 m
SampleDetDist= 1.003811e-01 m   PONI= -3.757304e-04, -1.204024e-05 m    rot1=-0.000164  rot2=-0.000283  rot3=0.000000 rad
DirectBeamDist= 100.381 mm      Center: x=0.026, y=-2.349 pix   Tilt= 0.019° tiltPlanRotation= -59.866° 𝛌= 0.965Å
Detector Pilatus 2M      PixelSize= 172µm, 172µm         BottomRight (3)
Wavelength= 9.654590e-11 m
SampleDetDist= 1.003811e-01 m   PONI= -3.757330e-04, -1.204277e-05 m    rot1=-0.000164  rot2=-0.000283  rot3=0.000000 rad
DirectBeamDist= 100.381 mm      Center: x=0.026, y=-2.349 pix   Tilt= 0.019° tiltPlanRotation= -59.866° 𝛌= 0.965Å
Detector Pilatus 2M      PixelSize= 172µm, 172µm         BottomRight (3)
Wavelength= 9.654590e-11 m
SampleDetDist= 1.003811e-01 m   PONI= -3.757330e-04, -1.204277e-05 m    rot1=-0.000164  rot2=-0.000283  rot3=0.000000 rad
DirectBeamDist= 100.381 mm      Center: x=0.026, y=-2.349 pix   Tilt= 0.019° tiltPlanRotation= -59.866° 𝛌= 0.965Å
[11]:
#Return to `inline` mode
%matplotlib inline
calib.peakPicker.widget.fig.show()

Check that the beam-center and the distance is correct and how much they are off.

Calibration of all frames in automatic mode:

[12]:
# Definition of the geometry translation function:
get_distance = lambda fimg: get_dectris_headers(fimg)["Detector_distance"]

geotrans = GeometryTransformation(param_names = ["dist_offset",
                                                 "poni1", "poni2", "rot1","rot2",
                                                "dist_scale", "poni1_scale", "poni2_scale"],
                                  dist_expr="pos * dist_scale + dist_offset",
                                  poni1_expr="pos * poni1_scale + poni1",
                                  poni2_expr="pos * poni2_scale + poni2",
                                  rot1_expr="rot1",
                                  rot2_expr="rot2",
                                  rot3_expr="0.0")

param = {
         "dist_offset": calib.geoRef.dist-get_distance(first),
         "poni1": calib.geoRef.poni1,
         "poni2": calib.geoRef.poni2,
         "rot1": calib.geoRef.rot1,
         "rot2": calib.geoRef.rot2,
         "dist_scale": 1.0,
         "poni1_scale": 0.0,
         "poni2_scale": 0.0,
}

print(param)
{'dist_offset': np.float64(-0.026092892938312504), 'poni1': np.float64(-0.00037573300834645016), 'poni2': np.float64(-1.2042769862867359e-05), 'rot1': np.float64(-0.00016408204750500895), 'rot2': np.float64(-0.00028267253920143603), 'dist_scale': 1.0, 'poni1_scale': 0.0, 'poni2_scale': 0.0}
[13]:
# Definition of the geometry refinement: the parameter order is the same as the param_names



gonioref = GoniometerRefinement(param, #initial guess
                                pos_function=get_distance,
                                trans_function=geotrans,
                                detector=detector,
                                wavelength=wavelength)
print("Empty refinement object:")
print(gonioref)
Empty refinement object:
GoniometerRefinement with 0 geometries labeled: .
[14]:
# Let's populate the goniometer refinement object with all geometries:
for fn in files:
    base = os.path.splitext(fn)[0]
    fimg = fabio.open(fn)
    local_calib = AbstractCalibration(img=fimg.data, mask=mask,
                                      detector=detector,
                                      wavelength=wavelength,
                                      calibrant=calibrant)
    local_calib.preprocess()
    local_calib.fixed = ["wavelength", "rot3"]
    local_calib.ai = gonioref.get_ai(get_distance(fimg))
    local_calib.extract_cpt()
    sg = gonioref.new_geometry(os.path.basename(base), image=fimg.data, metadata=fimg,
                              control_points=local_calib.peakPicker.points,
                              geometry=local_calib.ai,
                              calibrant=calibrant)

print("Filled refinement object:")
print(gonioref)
print(os.linesep+"\tLabel \t Distance")
for k, v in gonioref.single_geometries.items():
    print(k,v.get_position())
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
ERROR:root:No diffraction image available => not showing the contour
Filled refinement object:
GoniometerRefinement with 8 geometries labeled: test-powder_1_0001, test-powder_2_0001, test-powder_3_0001, test-powder_4_0001, test-powder_5_0001, test-powder_6_0001, test-powder_7_0001, test-powder_8_0001.

        Label    Distance
test-powder_1_0001 0.126474
test-powder_2_0001 0.141749
test-powder_3_0001 0.199249
test-powder_4_0001 0.171074
test-powder_5_0001 0.226674
test-powder_6_0001 0.293162
test-powder_7_0001 0.357899
test-powder_8_0001 0.484611
[15]:
fig, ax = subplots(len(files), figsize=(10, 10*len(files)))
for sp, sg in zip(ax, gonioref.single_geometries.values()):
    jupyter.display(sg=sg, ax=sp, label=sg.label)
../../../../_images/usage_tutorial_Goniometer_MX-calibrate_MX-calibrate_17_0.png
../../../../_images/usage_tutorial_Goniometer_MX-calibrate_MX-calibrate_17_1.png
[16]:
gonioref.refine3(fix=["dist_scale", "poni1_scale", "poni2_scale"])
Free parameters: ['dist_offset', 'poni1', 'poni2', 'rot1', 'rot2']
Fixed: {'dist_scale': 1.0, 'poni1_scale': 0.0, 'poni2_scale': 0.0}
 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 0.00013393119308234097
       x: [-2.690e-02 -5.431e-04  3.443e-04 -9.480e-04  2.735e-03]
     nit: 12
     jac: [-2.467e-08 -6.754e-08  2.978e-08  3.201e-08 -2.547e-10]
    nfev: 76
    njev: 12
Constrained Least square 0.00013515746371649427 --> 0.00013393119308234097
maxdelta on rot2: -0.00028267253920143603 --> 0.0027349029831441
[16]:
np.float64(0.00013393119308234097)
[17]:
gonioref.refine3(fix=[])
Free parameters: ['dist_offset', 'poni1', 'poni2', 'rot1', 'rot2', 'dist_scale', 'poni1_scale', 'poni2_scale']
Fixed: {}
 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 0.00013309356727363595
       x: [-2.465e-02 -1.991e-03 -6.307e-04  4.048e-04  2.098e-04
            9.896e-01  1.030e-02  5.912e-03]
     nit: 19
     jac: [ 1.849e-07  1.323e-07  2.279e-07 -3.294e-07 -6.368e-08
           -2.767e-08  3.281e-07 -5.561e-07]
    nfev: 174
    njev: 19
Constrained Least square 0.00013393119308234097 --> 0.00013309356727363595
maxdelta on dist_scale: 1.0 --> 0.9895731490421706
[17]:
np.float64(0.00013309356727363595)

Interpretation of this fit:

[18]:
gonioref.get_ai(0.2)
[18]:
Detector Pilatus 2M      PixelSize= 172µm, 172µm         BottomRight (3)
Wavelength= 9.654590e-11 m
SampleDetDist= 1.732631e-01 m   PONI= 6.831004e-05, 5.517154e-04 m      rot1=0.000405  rot2=0.000210  rot3=0.000000 rad
DirectBeamDist= 173.263 mm      Center: x=2.800, y=0.609 pix    Tilt= 0.026° tiltPlanRotation= 152.598° 𝛌= 0.965Å
[19]:
gonioref.get_ai(0.3)
[19]:
Detector Pilatus 2M      PixelSize= 172µm, 172µm         BottomRight (3)
Wavelength= 9.654590e-11 m
SampleDetDist= 2.722204e-01 m   PONI= 1.098207e-03, 1.142930e-03 m      rot1=0.000405  rot2=0.000210  rot3=0.000000 rad
DirectBeamDist= 272.220 mm      Center: x=6.004, y=6.717 pix    Tilt= 0.026° tiltPlanRotation= 152.598° 𝛌= 0.965Å

Persistence of this fit

[20]:
gonioref.save(result_file)
[21]:
with open(result_file) as r:
    print(r.read())
{
  "content": "Goniometer calibration v2",
  "detector": "Pilatus 2M",
  "detector_config": {
    "orientation": 3
  },
  "wavelength": 9.65459e-11,
  "param": [
    -0.024651575968886893,
    -0.0019914832674963775,
    -0.0006307139722537774,
    0.00040475111186017863,
    0.00020981746873407635,
    0.9895731490421706,
    0.010298966520882262,
    0.005912146676571252
  ],
  "param_names": [
    "dist_offset",
    "poni1",
    "poni2",
    "rot1",
    "rot2",
    "dist_scale",
    "poni1_scale",
    "poni2_scale"
  ],
  "pos_names": [
    "pos"
  ],
  "trans_function": {
    "content": "GeometryTransformation",
    "param_names": [
      "dist_offset",
      "poni1",
      "poni2",
      "rot1",
      "rot2",
      "dist_scale",
      "poni1_scale",
      "poni2_scale"
    ],
    "pos_names": [
      "pos"
    ],
    "dist_expr": "pos * dist_scale + dist_offset",
    "poni1_expr": "pos * poni1_scale + poni1",
    "poni2_expr": "pos * poni2_scale + poni2",
    "rot1_expr": "rot1",
    "rot2_expr": "rot2",
    "rot3_expr": "0.0",
    "constants": {
      "pi": 3.141592653589793
    }
  }
}

Interpretation of the fit:

[22]:
distances = []
f_distances = []
f_poni1 = []
f_poni2 = []
g_distances = []
g_poni1 = []
g_poni2 = []
for sg in gonioref.single_geometries.values():
    distance = sg.get_position()
    distances.append(distance)
    sg.geometry_refinement.refine3(fix=["wavelength", "rot3"])
    f_distances.append(sg.geometry_refinement.dist)
    f_poni1.append(sg.geometry_refinement.poni1)
    f_poni2.append(sg.geometry_refinement.poni2)
    ai = gonioref.get_ai(distance)
    g_distances.append(ai.dist)
    g_poni1.append(ai.poni1)
    g_poni2.append(ai.poni2)
[23]:
fig,ax = subplots(2)
ax[0].plot(distances, f_distances, label="fitted individually")
ax[0].plot(distances, g_distances, label="fitted table")
ax[0].set_title("Observed deviations:")
ax[1].set_xlabel("Encoder distance (m)")
ax[1].plot(distances, f_poni1, label="poni1 fitted individually")
ax[1].plot(distances, g_poni1, label="poni1 fitted table")
ax[1].plot(distances, f_poni2, label="poni2 fitted individually")
ax[1].plot(distances, g_poni2, label="poni2 fitted table")
ax[0].set_ylabel("Fitted distance (m)")
ax[1].set_ylabel("Fitted PONIs (m)")
ax[0].legend()
ax[1].legend()
pass
../../../../_images/usage_tutorial_Goniometer_MX-calibrate_MX-calibrate_28_0.png
[24]:
obtained_x = []
obtained_y = []

for dst in distances:
    geo = gonioref.get_ai(dst)
    fit2d = geo.getFit2D()
    obtained_x.append(fit2d["centerX"])
    obtained_y.append(fit2d["centerY"])
fig,ax = subplots()
ax.plot(distances, obtained_x, label="X")
ax.plot(distances, obtained_y, label="Y")
ax.set_title("Observed variation:")
ax.set_ylabel("Beam center (pixels)")
ax.set_xlabel("Encoder distance (m)")
ax.legend()
pass
../../../../_images/usage_tutorial_Goniometer_MX-calibrate_MX-calibrate_29_0.png
[25]:
#Simply print out the result:
lrx = linregress(distances, obtained_x)
lry = linregress(distances, obtained_y)
print(f"Beam-center X: {lrx}")
print(f"Beam-center y: {lry}")
print()
print(f"beamcenter_x = {lrx.intercept:.3f} {lrx.slope:+.3f} * encoder_value")
print(f"beamcenter_y = {lry.intercept:.3f} {lry.slope:+.3f} * encoder_value")
Beam-center X: LinregressResult(slope=np.float64(32.044278036961266), intercept=np.float64(-3.608931505386055), rvalue=np.float64(1.0), pvalue=np.float64(2.5000000000000343e-60), stderr=np.float64(0.0), intercept_stderr=np.float64(0.0))
Beam-center y: LinregressResult(slope=np.float64(61.08486205927713), intercept=np.float64(-11.608462786430533), rvalue=np.float64(1.0), pvalue=np.float64(2.5000000000000343e-60), stderr=np.float64(0.0), intercept_stderr=np.float64(0.0))

beamcenter_x = -3.609 +32.044 * encoder_value
beamcenter_y = -11.608 +61.085 * encoder_value

Nota:

The degradation between 0.3 and 0.5m correspond to the image 6->7 and is related to the disparition of the third ring!

Conclusion:

This notebook demonstrates: * The usage of the geometry calibration in Jupyter-lab to calibrate the first image * The creation of a goniometer-refinement * The population of this goniometer-refinement with automatic control-point extraction * The fit of the table, first with the constrains of a perfecty aligned table, then with a mis-aligned table

In our case the table is miss-aligned in the horizontal direction by 2.3mm/meter (i.e. 2.3 mradian). This should be taken into account when calculating the beam-center at different distances.