ImXPAD S540 detector at D2AMΒΆ

This tutorial corresponds to the calibration the goniometer an ImXPAD detector composed of 8 stripes of 7 modules, many of which are defective, on a goniometer.

This detector is mounted on the goniometer 2theta arm at the D2AM beam-line, French CRG at the ESRF synchrotron.

The raw data files are available at: http://www.silx.org/pub/pyFAI/gonio/D2AM-15/

%pylab nbagg
Populating the interactive namespace from numpy and matplotlib
import glob, os, time
start_time = time.time()
import fabio, pyFAI
from pyFAI.goniometer import GeometryTransformation, GoniometerRefinement, Goniometer
from pyFAI.gui import jupyter
#Definition of the detector and deplay of an image and its mask:

d5 = pyFAI.detector_factory("D5Geom.h5")
print(d5.shape)

images = glob.glob("*.edf.gz")
fimg = fabio.open(images[-1])

for k,v in fimg.header.items():
    print(k, ": ", v)

f,ax=subplots(1,2)
ax[0].imshow(d5.mask, origin="lower")
ax[0].set_title("D5 Mask")
ax[1].imshow(numpy.arcsinh(fimg.data), cmap="inferno", origin="lower")
ax[1].set_title(fimg.filename)
(960, 578)
EDF_DataBlockID :  0.Image.Psd
EDF_BinarySize :  4439040
EDF_HeaderSize :  1536
ByteOrder :  LowByteFirst
DataType :  DoubleValue
Dim_1 :  578
Dim_2 :  960
Image :  0
HeaderID :  EH:000000:000000:000000
Size :  4439040
VersionNumber :  1
Epoch :  1481332406.5250968933
det_sample_dist :  0
y_beam :  0
x_beam :  0
Lambda :  0.495938
offset :  0
count_time :  120
point_no :  66
scan_no :  906
preset :  0
col_end :  559
col_beg :  0
row_end :  959
row_beg :  0
counter_pos :  120 2647 175 14.9682 0 88.8742 14.9682 0 25 25 1791 9.11272e+09 0 0 0 174.849 0 0 14.9682
counter_mne :  sec vct1 vct2 vct3 vct4 Imach pseudoC pfoil Emono Ecod img roi1 roi2 roi3 roi4 pico1 pico2 pico3 pico4
motor_pos :  65.9999 0.081186 89.9909 -89.9917 -0.0027 0.0021 57.1227 134.747 -32.9503 0.16656 -5 0.47558 -1.5 0 0 0 4.53604 0.1417 1.04 1.04022 1.04022 -4.4 -1.10211 -0.543725 -9.962 -14.038 -16.865 -7.195 24 -2.038 24.06 4.835 1
motor_mne :  del eta chi phi nu mu keta kap kphi tsx tsy tsz rox roy tox toy mono inc1 courb courbb courbf omega khimono gamma su6 sd6 sf6 sb6 vg6 vo6 hg6 ho6 rien
suffix :  .edf
prefix :  16Dec08D5_
dir :  /users/opd02/raw
run :  1791
title :  CCD Image
<IPython.core.display.Javascript object>
<matplotlib.text.Text at 0x7ff16c4be780>
# Define wavelength and create our "large" LaB6 calibrant

wavelength = 0.495938 * 1e-10
from pyFAI.calibrant import Cell, Calibrant
c = Cell.cubic(4.1568260)
c.save("LaB6", dmin=0.2)
LaB6 = Calibrant("LaB6.D")
LaB6.wavelength = wavelength
print("2theta max: ", numpy.degrees(LaB6.get_2th()[-1]))
print("Number of reflections: ", len(LaB6.get_2th()))
2theta max:  179.173497672
Number of reflections:  236
#Use a few manually calibrated images:

npt_files = [ i for i in glob.glob("*.npt") if "new" not in i]
npt_files.sort()
npt_files[0]
print("Number of hand-calibrated images :",len(npt_files))
Number of hand-calibrated images : 15
# Definition of the goniometer translation function:
# The detector rotates vertically, around the horizontal axis, i.e. rot2.
# Rotation both around axis 1 and axis 2 are allowed

goniotrans = GeometryTransformation(param_names = ["dist", "poni1", "poni2",
                                                   "rot1", "rot2", "rot3", "scale1", "scale2" ],
                                    dist_expr="dist",
                                    poni1_expr="poni1",
                                    poni2_expr="poni2",
                                    rot1_expr="scale1 * pos +rot1",
                                    rot2_expr="scale2 * pos + rot2",
                                    rot3_expr="rot3")


#Definition of the function reading the goniometer angle from the filename of the image.

def get_angle(metadata):
    """Takes the angle from the first motor position and returns the angle of the goniometer arm"""
    return float(metadata["motor_pos"].split()[0])

print('filename', fimg.filename, "angle:",get_angle(fimg.header))
filename 16Dec08D5_1791-rsz.edf.gz angle: 65.9999
# Definition of the geometry refinement: the parameter order is the same as the param_names

rot3 = numpy.pi/2
scale1 = -numpy.pi/180
scale2 = 0
param = {"dist":0.5,
         "poni1":0.05,
         "poni2":0.05,
         "rot1":0,
         "rot2":0,
         "rot3": rot3,
         "scale1": scale1,
         "scale2": scale2,
        }
#Defines the bounds for some variables
bounds = {"dist": (0.2, 0.8),
          "poni1": (0, 0.1),
          "poni2": (0, 0.1),
          "rot1": (-0.1, 0.1),
          "rot2": (-0.1, 0.1),
          "rot3": (rot3, rot3), #strict bounds on rot3
          #"scale1": (scale1, scale1),
          #"scale2": (scale2, scale2),
         }
gonioref = GoniometerRefinement(param, #initial guess
                                bounds=bounds,
                                pos_function=get_angle,
                                trans_function=goniotrans,
                                detector=d5, wavelength=wavelength)
print("Empty refinement object:", gonioref)

#Let's populate the goniometer refinement object with all control point files:

for fn in npt_files[:]:
    base = os.path.splitext(fn)[0]
    fimg = fabio.open(base + ".edf.gz")
    sg =gonioref.new_geometry(base, image=fimg.data, metadata=fimg.header, control_points=fn, calibrant=LaB6)
    print(base, "Angle:", sg.get_position())


print("Filled refinement object:")
print(gonioref)
Empty refinement object: GoniometerRefinement with 0 geometries labeled: .
16Dec08D5_1725-rsz Angle: -0.003
16Dec08D5_1726-rsz Angle: 0.9998
16Dec08D5_1727-rsz Angle: 2.0
16Dec08D5_1728-rsz Angle: 2.9998
16Dec08D5_1729-rsz Angle: 4.0002
16Dec08D5_1730-rsz Angle: 4.9998
16Dec08D5_1735-rsz Angle: 10.0001
16Dec08D5_1742-rsz Angle: 16.9996
16Dec08D5_1749-rsz Angle: 24.0001
16Dec08D5_1756-rsz Angle: 30.9997
16Dec08D5_1763-rsz Angle: 37.9999
16Dec08D5_1770-rsz Angle: 44.9997
16Dec08D5_1777-rsz Angle: 52.0
16Dec08D5_1784-rsz Angle: 58.9995
16Dec08D5_1791-rsz Angle: 65.9999
Filled refinement object:
GoniometerRefinement with 15 geometries labeled: 16Dec08D5_1725-rsz, 16Dec08D5_1726-rsz, 16Dec08D5_1727-rsz, 16Dec08D5_1728-rsz, 16Dec08D5_1729-rsz, 16Dec08D5_1730-rsz, 16Dec08D5_1735-rsz, 16Dec08D5_1742-rsz, 16Dec08D5_1749-rsz, 16Dec08D5_1756-rsz, 16Dec08D5_1763-rsz, 16Dec08D5_1770-rsz, 16Dec08D5_1777-rsz, 16Dec08D5_1784-rsz, 16Dec08D5_1791-rsz.
# Initial refinement of the goniometer model with 5 dof

gonioref.refine2()
Cost function before refinement: 0.000307661216084
     fun: 1.985522552167619e-08
     jac: array([  2.40897930e-08,   3.04889987e-07,   4.13706377e-07,
        -1.75572721e-07,  -9.57396855e-07,   2.39808173e-14,
        -5.24588331e-06,  -7.15055308e-07])
 message: 'Optimization terminated successfully.'
    nfev: 178
     nit: 17
    njev: 17
  status: 0
 success: True
       x: array([  5.20464974e-01,   6.34779274e-02,   4.46077057e-02,
         3.03313916e-03,   7.46909752e-03,   1.57079633e+00,
        -1.74719084e-02,  -3.84262988e-04])
Cost function after refinement: 1.985522552167619e-08
GonioParam(dist=0.52046497447469431, poni1=0.063477927365486408, poni2=0.044607705706853802, rot1=0.0030331391589748357, rot2=0.0074690975201246379, rot3=1.5707963267948966, scale1=-0.017471908444320044, scale2=-0.00038426298826876057)
maxdelta on: dist (0) 0.5 --> 0.520464974475
array([  5.20464974e-01,   6.34779274e-02,   4.46077057e-02,
         3.03313916e-03,   7.46909752e-03,   1.57079633e+00,
        -1.74719084e-02,  -3.84262988e-04])
width=3
height=int(ceil(len(gonioref.single_geometries)/width))
fig,ax = subplots(height, width,figsize=(10,15))
for idx, sg in enumerate(gonioref.single_geometries.values()):
    sg.geometry_refinement.set_param(gonioref.get_ai(sg.get_position()).param)
    jupyter.display(sg=sg, ax=ax[idx//width, idx%width])
<IPython.core.display.Javascript object>
/scisoft/users/jupyter/jupy34/lib/python3.4/site-packages/matplotlib/axes/_axes.py:545: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "
# Final pass of refinement with all constrains removed, very fine refinement

gonioref.bounds = None
gonioref.refine2("slsqp", eps=1e-13, maxiter=10000, ftol=1e-12)
Cost function before refinement: 1.98552255217e-08
     fun: 1.9855093621939147e-08
     jac: array([  8.40415502e-09,   0.00000000e+00,   0.00000000e+00,
        -1.50215999e-08,  -8.69763870e-07,  -3.11019910e-09,
        -3.47206415e-05,  -5.44781151e-07])
 message: 'Optimization terminated successfully.'
    nfev: 15
     nit: 1
    njev: 1
  status: 0
 success: True
       x: array([  5.20464974e-01,   6.34779274e-02,   4.46077057e-02,
         3.03313916e-03,   7.46909771e-03,   1.57079633e+00,
        -1.74719008e-02,  -3.84262869e-04])
Cost function after refinement: 1.9855093621939147e-08
GonioParam(dist=0.52046497447285478, poni1=0.063477927365486408, poni2=0.044607705706853802, rot1=0.0030331391622627612, rot2=0.0074690977104984162, rot3=1.5707963267955773, scale1=-0.017471900844671508, scale2=-0.00038426286902714811)
maxdelta on: scale1 (6) -0.0174719084443 --> -0.0174719008447
array([  5.20464974e-01,   6.34779274e-02,   4.46077057e-02,
         3.03313916e-03,   7.46909771e-03,   1.57079633e+00,
        -1.74719008e-02,  -3.84262869e-04])
#Create a MultiGeometry integrator from the refined geometry:

angles = []
images = []
for sg in gonioref.single_geometries.values():
    angles.append(sg.get_position())
    images.append(sg.image)

multigeo = gonioref.get_mg(angles)
multigeo.radial_range=(0, 80)
print(multigeo)
MultiGeometry integrator with 15 geometries on (0, 80) radial range (2th_deg) and (-180, 180) azimuthal range (deg)
# Integrate the whole set of images in a single run:

res = multigeo.integrate1d(images, 10000)
fig, ax = subplots()
ax.plot(*res)
ax.set_xlabel(res.unit.label)
ax.set_ylabel("Intensity")

#Note the large number of peaks due to hot pixels ....
<IPython.core.display.Javascript object>
<matplotlib.text.Text at 0x7ff11472d2e8>
#Add hot pixels to the mask: pixel which are 15x more intense than the median in their ring.

thres = 15

old_mask = d5.mask.astype("bool", copy=True)
new_mask = d5.mask.astype("bool", copy=True)

for ai,img in zip(multigeo.ais,images):
    b,a = ai.separate(img, 1000, restore_mask=0)
    b[old_mask] = 0
    b[b<0] = 0
    print(sum(b>thres*a))
    new_mask = numpy.logical_or(new_mask, (b>thres*a))

print(sum(old_mask), sum(new_mask), sum(new_mask)-sum(old_mask))
4
4
4
4
4
4
4
4
4
4
4
3
3
3
3
198678 198683 5
# Update the mask
for ai in multigeo.ais:
    ai.detector.mask = new_mask

# Integrate the whole set of images in a single run:
res2 = multigeo.integrate1d(images, 10000)
fig, ax = subplots()
ax.plot(*res, label="Before hot-pixel removal")
ax.plot(*res2, label="After hot-pixel removal")
ax.legend()
ax.set_xlabel(res.unit.label)
ax.set_ylabel("Intensity")
<IPython.core.display.Javascript object>
<matplotlib.text.Text at 0x7ff11470c6a0>
# Integrate the whole set of images in 2D:

res2d = multigeo.integrate2d(images, 1000, 360)
fig, ax = subplots()
ax.imshow(numpy.arcsinh(res2d[0]), cmap="inferno", origin="lower", aspect="auto",
          extent=[res2d[1][0], res2d[1][-1], res2d[2][0], res2d[2][-1]])
ax.set_xlabel(res.unit.label)
ax.set_ylabel(r"$\chi$ angle")
<IPython.core.display.Javascript object>
<matplotlib.text.Text at 0x7ff114664b00>
print("Total execution time", time.time()-start_time)
Total execution time 18.806222677230835

Previous topic

Calibration of a 2 theta arm with a Pilatus 100k detector

Next topic

pyFAI scripts manual

This Page