Calibration of a very large Pilatus detector with overlapping grid position

This tutorial presents the calibration of the Pilatus 900kw CdTe which is a very large 2D detector (4500x200) running at ESRF ID06-LVP . The detector is so large that the grid needs to be displaced in front of the detector and the operation needs to be performed several times.

The overall strategy is very similar to the ID15 detector calibration, excepts that all needs to be done 3 times, one for each of the grid position: left, center and right:

1. Image preprocessing

2. Peak picking

3. Grid assignment

4. Displacement fitting

5. Reconstruction of the pixel position

6. Saving into a detector definition file

7. Validation of the geometry with a 2D integration

Each module being made by lithographic processes, the error within a module will be assumeed to be constant. We will use the name “displacement of the module” to describe the rigide movement of the module.

This tutorial uses data acquired by Marie Ruat from the ESRF detector group during the commissionning of the detector. The ID06-LVP is acknowledged for commissionning beam-time and fruitful discussion.

This detector contains 18 half-modules, each bound to a single CdTe monocrystal sensor and is designed for high energy X-ray radiation detection. Due to the construction procedure, these half-modules could show a misalignment within the detector plane. While the manufacturer (Dectris) garanties a precision within a pixel (172µm), the miss-alignment of certain modules can be seen while calibrating Debye-Scherrer ring using refereance sample. So the aim of this work is to provide a detector description with a better precision better than the original detector.

This work will be performed on the image of a grid available: http://www.silx.org/pub/pyFAI/detector_calibration

It is a good exercise to calibrate all rings of the later image using the pyFAI-calib2 tool. A calibration close to perfection is needed to visualize the module miss-alignement we aim at correcting.

%matplotlib inline
# %matplotlib nbagg

#many imports which will be used all along the notebook
import time
start_time = time.perf_counter()
import os
import pyFAI
import fabio
import glob
import numpy
from numpy.lib.stride_tricks import as_strided
from collections import namedtuple
from math import sin, cos, sqrt
from scipy.ndimage import convolve, binary_dilation
from scipy.spatial import distance_matrix
from scipy.optimize import minimize
from matplotlib.pyplot import subplots
from pyFAI.ext.bilinear import Bilinear
from pyFAI.ext.watershed import InverseWatershed
from silx.resources import ExternalResources

print("Using pyFAI verison: ", pyFAI.version)

Triplet = namedtuple("Triplet", "left center right")

from matplotlib import colors
logcolor = colors.LogNorm(1e5, 3e5)
normcolor = colors.LogNorm(1, 2)

# Some compound types ...
dt = numpy.dtype([('y', numpy.float64),
                  ('x', numpy.float64),
                  ('i', numpy.int64)])

dl = numpy.dtype([('y', numpy.float64),
                  ('x', numpy.float64),
                  ('i', numpy.int64),
                  ('Y', numpy.int64),
                  ('X', numpy.int64)])

Using pyFAI verison:  0.20.0-beta1

#Download all data:
#Nota: Configure here your proxy if you are behind a firewall
#os.environ["http_proxy"] = "http://proxy.comany.com:3128"
downloader = ExternalResources("detector_calibration", "http://www.silx.org/pub/pyFAI/detector_calibration/")
median21_left = downloader.getfile("Pilatus900kwID06_median21_left.npy")
median21_center = downloader.getfile("Pilatus900kwID06_median21_center.npy")
median21_right = downloader.getfile("Pilatus900kwID06_median21_right.npy")
mask_left = downloader.getfile("Pilatus900kwID06_mask_left.npy")
mask_center = downloader.getfile("Pilatus900kwID06_mask_center.npy")
mask_right = downloader.getfile("Pilatus900kwID06_mask_right.npy")
minimum = downloader.getfile("Pilatus900kwID06_minimum.npy")

img_left = fabio.open(median21_left).data
img_center = fabio.open(median21_center).data
img_right = fabio.open(median21_right).data
fig,ax = subplots(3, figsize=(20,6))
ax[0].set_title("Grid points acquired (after median filter)")
ax[0].imshow(img_left, interpolation="bilinear", norm=logcolor, cmap="inferno")
ax[1].imshow(img_center, interpolation="bilinear", norm=logcolor, cmap="inferno")
ax[2].imshow(img_right, interpolation="bilinear", norm=logcolor, cmap="inferno")
pass

def display(triplet, **kwargs):
    fig,ax = subplots(3, figsize=(20, 9))
    ax[0].set_title("left")
    ax[0].imshow(triplet.left, **kwargs)
    ax[0].set_xlim(0, 1580)
    ax[1].set_title("center")
    ax[1].imshow(triplet.center, **kwargs)
    ax[1].set_xlim(1400, 3010)
    ax[2].set_title("right")
    ax[2].imshow(triplet.right, **kwargs)
    ax[2].set_xlim(2830, triplet.right.shape[-1]+5)
    return ax
data = Triplet(img_left, img_center, img_right)
display(data, interpolation="bilinear", norm=logcolor, cmap="inferno")
pass

Image processing

There are 4 pre-processing steps which are needed.

• Define for each module a unique identifier which will be used later on during the fitting procedure

• Define the proper mask: each module is the assembly of 4x2 sub-modules and there are (3) interpolated pixels between each sub-module, such “unreliable pixels should be masked out as well

• Correct the grid image by the smoothed image to have a constant background.

• Convolve the raw image with a typical hole shape to allow a precise spotting of the hole center.

pilatus = pyFAI.detector_factory("Pilatus_900kw_CdTe")
print(pilatus)
print(pilatus.shape)
mask1 = pilatus.mask
module_size = pilatus.MODULE_SIZE
module_gap = pilatus.MODULE_GAP
submodule_size = (96,60)

Detector Pilatus CdTe 900kw      PixelSize= 1.720e-04, 1.720e-04 m
(195, 4439)

#1 + 2 Calculation of the module_id and the interpolated-mask:
mid = numpy.zeros(pilatus.shape, dtype=int)
mask2 = numpy.zeros(pilatus.shape, dtype=int)
idx = 1
for i in range(1):
    y_start = i*(module_gap[0] + module_size[0])
    y_stop = y_start + module_size[0]
    for j in range(9):
        x_start = j*(module_gap[1] + module_size[1])
        x_stop = x_start + module_size[1]
        mid[y_start:y_stop,x_start: x_start+module_size[1]//2] = idx
        idx+=1
        mid[y_start:y_stop,x_start+module_size[1]//2: x_stop] = idx
        idx+=1
        mask2[y_start+submodule_size[0]-1:y_start+submodule_size[0]+2,
              x_start:x_stop] = 1
        for k in range(1,8):
            mask2[y_start:y_stop,
              x_start+k*(submodule_size[1]+1)-1:x_start+k*(submodule_size[1]+1)+2] = 1

fix, ax = subplots(2, figsize=(20,4))
ax[0].set_title("Module Id and inter-module mask")
ax[0].imshow(mid)
ax[1].imshow(mask2)
pass

#Extra masking: bad pixels marked by the detector
mask0 = fabio.open(minimum).data<0

bad = numpy.where(mask0 | mask1 | mask2 | fabio.open(mask_left).data)
data.left[bad] = numpy.nan
bad = numpy.where(mask0 | mask1 | mask2 | fabio.open(mask_center).data)
data.center[bad] = numpy.nan
bad = numpy.where(mask0 | mask1 | mask2 | fabio.open(mask_right).data)
data.right[bad] = numpy.nan

display(data, interpolation="bilinear", norm=logcolor)
pass

# The Nan-masked image contains now only valid values (and Nan elsewhere). We will make a large median filter to
# build up a smooth image without gaps.
#
# This function is backported from future version of numpy ... it allows to expose a winbowed view
# to perform the nanmedian-filter

def sliding_window_view(x, shape, subok=False, readonly=True):
    """
    Creates sliding window views of the N dimensional array with the given window
    shape. Window slides across each dimension of `x` and extract subsets of `x`
    at any window position.
    Parameters
    ----------
    x : array_like
        Array to create sliding window views of.
    shape : sequence of int
        The shape of the window. Must have same length as the number of input array dimensions.
    subok : bool, optional
        If True, then sub-classes will be passed-through, otherwise the returned
        array will be forced to be a base-class array (default).
    readonly : bool, optional
        If set to True, the returned array will always be readonly view.
        Otherwise it will return writable copies(see Notes).
    Returns
    -------
    view : ndarray
        Sliding window views (or copies) of `x`. view.shape = x.shape - shape + 1
    See also
    --------
    as_strided: Create a view into the array with the given shape and strides.
    broadcast_to: broadcast an array to a given shape.
    Notes
    -----
    ``sliding_window_view`` create sliding window views of the N dimensions array
    with the given window shape and its implementation based on ``as_strided``.
    Please note that if readonly set to True, views are returned, not copies
    of array. In this case, write operations could be unpredictable, so the returned
    views are readonly. Bear in mind that returned copies (readonly=False) will
    take more memory than the original array, due to overlapping windows.
    For some cases there may be more efficient approaches to calculate transformations
    across multi-dimensional arrays, for instance `scipy.signal.fftconvolve`, where combining
    the iterating step with the calculation itself while storing partial results can result
    in significant speedups.
    Examples
    --------
    >>> i, j = np.ogrid[:3,:4]
    >>> x = 10*i + j
    >>> shape = (2,2)
    >>> np.lib.stride_tricks.sliding_window_view(x, shape)
    array([[[[ 0,  1],
             [10, 11]],
            [[ 1,  2],
             [11, 12]],
            [[ 2,  3],
             [12, 13]]],
           [[[10, 11],
             [20, 21]],
            [[11, 12],
             [21, 22]],
            [[12, 13],
             [22, 23]]]])
    """
    np = numpy
    # first convert input to array, possibly keeping subclass
    x = np.array(x, copy=False, subok=subok)

    try:
        shape = np.array(shape, np.int)
    except:
        raise TypeError('`shape` must be a sequence of integer')
    else:
        if shape.ndim > 1:
            raise ValueError('`shape` must be one-dimensional sequence of integer')
        if len(x.shape) != len(shape):
            raise ValueError("`shape` length doesn't match with input array dimensions")
        if np.any(shape <= 0):
            raise ValueError('`shape` cannot contain non-positive value')

    o = np.array(x.shape) - shape  + 1 # output shape
    if np.any(o <= 0):
        raise ValueError('window shape cannot larger than input array shape')

    if type(readonly) != bool:
        raise TypeError('readonly must be a boolean')

    strides = x.strides
    view_strides = strides

    view_shape = np.concatenate((o, shape), axis=0)
    view_strides = np.concatenate((view_strides, strides), axis=0)
    view = as_strided(x, view_shape, view_strides, subok=subok, writeable=not readonly)

    if not readonly:
        return view.copy()
    else:
        return view

%%time
#Calculate a background image using a large median filter ... takes a while
shape = (13,13)
padded = Triplet(*(numpy.pad(i, tuple((i//2,) for i in shape), mode="edge") for i in data))
print(padded.left.shape)

(207, 4451)
CPU times: user 2.95 ms, sys: 0 ns, total: 2.95 ms
Wall time: 2.67 ms

%%time
background = Triplet(*[numpy.nanmedian(sliding_window_view(i, shape), axis = (-2,-1)) for i in padded])
print(background.left.shape)

/usr/lib/python3/dist-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: All-NaN slice encountered
  r, k = function_base._ureduce(a, func=_nanmedian, axis=axis, out=out,

(195, 4439)
CPU times: user 13.7 s, sys: 1.08 s, total: 14.8 s
Wall time: 14.8 s

display(background, norm=logcolor, interpolation="bilinear")
pass

normalized = Triplet(*(i/j for i,j in zip(data, background)))
display(normalized, interpolation="nearest", norm=normcolor)
pass

fig,ax = subplots(1,3, figsize=(9,5))

ax[0].hist(normalized.left.ravel(), 100, range=(0,2))
ax[1].hist(normalized.center.ravel(), 100, range=(0,2))
ax[2].hist(normalized.right.ravel(), 100, range=(0,2))
pass

For a precise measurement of the peak position, one trick is to convolve the image with a pattern which looks like a hole of the grid.

#Definition of the convolution kernel
ksize = 5
y,x = numpy.ogrid[-(ksize-1)//2:ksize//2+1,-(ksize-1)//2:ksize//2+1]
d = numpy.sqrt(y*y+x*x)

#Fade out curve definition
fadeout = lambda x: 1/(1+numpy.exp(3*(x-2.2)))

kernel = fadeout(d)
mini=kernel.sum()
print(mini)

fig,ax = subplots(1,3)
ax[0].imshow(d)
ax[0].set_title("Distance array")

ax[1].plot(numpy.linspace(0,5,100),fadeout(numpy.linspace(0,5,100)))
ax[1].set_title("fade-out curve")

ax[2].imshow(kernel)
ax[2].set_title("Convolution kernel")
pass

15.439885086158014

my_smooth = Triplet(*(convolve(i, kernel, mode="constant", cval=0)/mini for i in normalized))
print(my_smooth.center.shape)
display(my_smooth)
pass

(195, 4439)

all_masks = mask0 | mask1 | mask2
big_mask = binary_dilation(all_masks, iterations=ksize//2+1+1)
print(all_masks.sum(), big_mask.sum())

62453 208997

Peak picking

We use the watershed module from pyFAI to retrieve all peak positions. Those regions are sieved out respectively for:

• their size, it should be larger than the kernel itself

• the peaks too close to masked regions are removed

• the intensity of the peak

%%time
tmp = []
for i in my_smooth:
    iw = InverseWatershed(i)
    iw.init()
    iw.merge_singleton()
    all_regions = set(iw.regions.values())
    regions = [i for i in all_regions if i.size>mini]
    tmp.append(regions)
regions = Triplet(*tmp)

CPU times: user 5.3 s, sys: 189 ms, total: 5.49 s
Wall time: 5.48 s

#Remove peaks on masked region
sieved_region = Triplet(*([i for i in j if not big_mask[(i.index//pilatus.shape[-1], i.index%pilatus.shape[-1])]]
                          for j in regions))
print("Number of peaks not on masked areea : %s %s %s"%
      (len(sieved_region[0]),len(sieved_region[1]),len(sieved_region[2])))

Number of peaks not on masked areea : 4773 4772 4885

# Histogram of peak height:
s = Triplet(*(numpy.array([i.maxi for i in j]) for j in sieved_region))

fig, ax = subplots(3, figsize=(15,6))
[ax[i].hist(s[i], 100) for i in range(3)]
[ax[i].set_yscale("log") for i in range(3)]
pass

#sieve-out for peak intensity
int_mini = 1.2
peaks = Triplet(*([(i.index//pilatus.shape[-1], i.index%pilatus.shape[-1]) for i in j if (i.maxi)>int_mini]
                  for j in sieved_region))
print("Number of remaining peaks with I>%s: %s"%(int_mini, [len(i) for i in peaks]))

peaks_raw = Triplet(*(numpy.array(i) for i in peaks))

Number of remaining peaks with I>1.2: [248, 242, 252]

# Finally the peak positions are interpolated using a second order taylor expansion
# in thevinicy of the maximum value of the signal:

#Create bilinear interpolator
bl = [Bilinear(i) for i in my_smooth]

#Overlay raw peak coordinate and refined peak positions

ref_peaks = [[b.local_maxi(p) for p in peaki] for b, peaki in zip(bl, peaks)]
ax = display(data)
peaks_ref = [numpy.array(i) for i in ref_peaks]
for i in range(3):
    ax[i].plot(peaks_raw[i][:,1], peaks_raw[i][:, 0], ".r")
    ax[i].plot(peaks_ref[i][:,1],peaks_ref[i][:, 0], ".b")
ax[0].set_title("Extracted peak position (red: raw, blue: refined)")
pass

display(Triplet(mid, mid, mid))
pass

At this stage we have about 3x250 peaks (with sub-pixel precision) which are visually distributed on all modules. Some modules have their peaks located along sub-module boundaries which are masked out, hence they have fewer ontrol point for the calculation. Let’s assign each peak to a module identifier. This allows to print out the number of peaks per module:

yxi = Triplet(*[numpy.array([i+(mid[round(i[0]),round(i[1])],)
                   for i in j], dtype=dt)
                for j in ref_peaks])
print("Number of keypoint per module:")
cp = Triplet([numpy.NaN], [numpy.NaN], [numpy.NaN])
for i in range(1,mid.max()+1):
    cp.left.append((yxi.left[:]["i"] == i).sum())
    cp.center.append((yxi.center[:]["i"] == i).sum())
    cp.right.append((yxi.right[:]["i"] == i).sum())
    print("Module id:",i,
          "left cp:", (yxi.left[:]["i"] == i).sum(),
          "center cp:", (yxi.center[:]["i"] == i).sum(),
          "right cp:", (yxi.right[:]["i"] == i).sum())

fig, ax = subplots(figsize=(8,6))
ax.plot(cp.left, "-or", label="left")
ax.plot(cp.center,"-og",  label="center")
ax.plot(cp.right,"-ob",  label="right")
ax.set_ylabel("Number of control point")
ax.set_xlabel("Module id")
ax.set_xticks(numpy.arange(1, 19))
ax.legend()
ax.set_title("Overlapping of the 3 grid positions")
pass

Number of keypoint per module:
Module id: 1 left cp: 33 center cp: 0 right cp: 0
Module id: 2 left cp: 48 center cp: 0 right cp: 0
Module id: 3 left cp: 29 center cp: 0 right cp: 0
Module id: 4 left cp: 42 center cp: 0 right cp: 0
Module id: 5 left cp: 42 center cp: 0 right cp: 0
Module id: 6 left cp: 36 center cp: 12 right cp: 0
Module id: 7 left cp: 18 center cp: 27 right cp: 0
Module id: 8 left cp: 0 center cp: 48 right cp: 0
Module id: 9 left cp: 0 center cp: 36 right cp: 0
Module id: 10 left cp: 0 center cp: 41 right cp: 0
Module id: 11 left cp: 0 center cp: 48 right cp: 0
Module id: 12 left cp: 0 center cp: 24 right cp: 24
Module id: 13 left cp: 0 center cp: 6 right cp: 30
Module id: 14 left cp: 0 center cp: 0 right cp: 48
Module id: 15 left cp: 0 center cp: 0 right cp: 42
Module id: 16 left cp: 0 center cp: 0 right cp: 30
Module id: 17 left cp: 0 center cp: 0 right cp: 48
Module id: 18 left cp: 0 center cp: 0 right cp: 30

Grid assignment

The calibration is performed using a regular grid, the idea is to assign to each peak of coordinates (x,y) the integer value (X, Y) which correspond to the grid corrdinate system.

The first step is to measure the grid pitch which correspond to the distance (in pixels) from one peak to the next. This is easily obtained from a pair-wise distribution function.

# pairwise distance calculation using scipy.spatial.distance_matrix

dist = [distance_matrix(i, i) for i in peaks_ref]

fig, ax = subplots(3,figsize=(15,6))
for i in range(3):
    ax[i].hist(dist[i].ravel(), 250, range=(0,100))
ax[0].set_title("Pair-wise distribution function")
ax[2].set_xlabel("distance in pixel")
pass

The histogram of the pair-distribution function has a first peak at 0 and the second peak between 29 and 30. Let’s start the fit with this value

Two other parameters correspond to the offset, in pixel for the grid index (X,Y) = (0,0). The easiest is to measure the smallest x and y for the first module.

#from pair-wise distribution histogram
step = 29.3

def index_module(module, position, step):
    "Return the peak position for the given module at the grid position. The guess parameter is also provided"
    indexed_module = yxi[position][yxi[position][:]["i"] == module]
    y_min = indexed_module["y"].min()
    x_min = indexed_module["x"].min()
    print("offset for the first peak: ", x_min, y_min)

    indexed = numpy.zeros(indexed_module.size, dtype=dl)
    delta_max = 0

    for i,v in enumerate(indexed_module):
        Y = int(round((v["y"] - y_min)/step))
        X = int(round((v["x"] - x_min)/step))
        indexed[i]["y"] = v["y"]
        indexed[i]["x"] = v["x"]
        indexed[i]["i"] = v["i"]
        indexed[i]["Y"] = Y
        indexed[i]["X"] = X
        #print(delta_max)
        delta_max = max(delta_max, sqrt((v["y"]-Y*step-y_min)**2 + (v["x"]-X*step-x_min)**2)/step)
        print("peak id: %2s %35s Y:%d (Δ=%.3f) X:%s (Δ=%.3f)"%
              (i,v, Y, (v["y"]-Y*step-y_min)/step, X, (v["x"]-X*step-x_min)/step))
    if delta_max>0.1:
        print("Assignment is error prone with delta_max (in steps) =", delta_max)
    guess = [step, y_min, x_min, 0]
    return indexed, guess

#work with the first module and fit the peak positions
indexed1, guess1 = index_module(1, 0, step)
guess1

offset for the first peak:  2.052476167678833 6.0
peak id:  0        (64.98052979, 2.05247617, 1) Y:2 (Δ=0.013) X:0 (Δ=0.000)
peak id:  1       (64.64517257, 31.64133275, 1) Y:2 (Δ=0.002) X:1 (Δ=0.010)
peak id:  2       (64.75793037, 90.28502041, 1) Y:2 (Δ=0.005) X:3 (Δ=0.011)
peak id:  3      (64.85140824, 148.84530333, 1) Y:2 (Δ=0.009) X:5 (Δ=0.010)
peak id:  4      (64.86348662, 207.42443478, 1) Y:2 (Δ=0.009) X:7 (Δ=0.009)
peak id:  5      (64.95814403, 236.73580122, 1) Y:2 (Δ=0.012) X:8 (Δ=0.010)
peak id:  6         (35.02235405, 2.8385089, 1) Y:1 (Δ=-0.009) X:0 (Δ=0.027)
peak id:  7       (35.39691702, 31.71216562, 1) Y:1 (Δ=0.003) X:1 (Δ=0.012)
peak id:  8       (6.34152079, 236.75626817, 1) Y:0 (Δ=0.012) X:8 (Δ=0.010)
peak id:  9        (6.31276822, 207.3469238, 1) Y:0 (Δ=0.011) X:7 (Δ=0.007)
peak id: 10       (6.25261039, 148.84848826, 1) Y:0 (Δ=0.009) X:5 (Δ=0.010)
peak id: 11                        (6., 90., 1) Y:0 (Δ=0.000) X:3 (Δ=0.002)
peak id: 12        (6.10884568, 31.65180257, 1) Y:0 (Δ=0.004) X:1 (Δ=0.010)
peak id: 13     (123.48410541, 207.45542562, 1) Y:4 (Δ=0.010) X:7 (Δ=0.010)
peak id: 14      (123.59747732, 236.6557506, 1) Y:4 (Δ=0.014) X:8 (Δ=0.007)
peak id: 15       (152.97377014, 2.05581188, 1) Y:5 (Δ=0.016) X:0 (Δ=0.000)
peak id: 16      (152.51041356, 31.59106046, 1) Y:5 (Δ=0.000) X:1 (Δ=0.008)
peak id: 17      (152.63233629, 90.28906977, 1) Y:5 (Δ=0.005) X:3 (Δ=0.011)
peak id: 18     (152.80133259, 148.82351807, 1) Y:5 (Δ=0.010) X:5 (Δ=0.009)
peak id: 19     (152.91486572, 236.64688599, 1) Y:5 (Δ=0.014) X:8 (Δ=0.007)
peak id: 20       (123.17985493, 2.57148761, 1) Y:4 (Δ=-0.001) X:0 (Δ=0.018)
peak id: 21      (123.25800872, 31.71460831, 1) Y:4 (Δ=0.002) X:1 (Δ=0.012)
peak id: 22      (123.35842642, 90.25260955, 1) Y:4 (Δ=0.005) X:3 (Δ=0.010)
peak id: 23     (123.43735978, 148.85604818, 1) Y:4 (Δ=0.008) X:5 (Δ=0.010)
peak id: 24       (181.66608471, 2.49644381, 1) Y:6 (Δ=-0.005) X:0 (Δ=0.015)
peak id: 25      (181.81992741, 31.69740948, 1) Y:6 (Δ=0.001) X:1 (Δ=0.012)
peak id: 26      (181.94562443, 90.26021093, 1) Y:6 (Δ=0.005) X:3 (Δ=0.011)
peak id: 27      (182.11482991, 148.8560147, 1) Y:6 (Δ=0.011) X:5 (Δ=0.010)
peak id: 28     (182.17538108, 207.34319276, 1) Y:6 (Δ=0.013) X:7 (Δ=0.007)
peak id: 29      (182.24330163, 236.7036297, 1) Y:6 (Δ=0.015) X:8 (Δ=0.009)
peak id: 30       (35.49848831, 90.35593623, 1) Y:1 (Δ=0.007) X:3 (Δ=0.014)
peak id: 31      (35.52367088, 207.50462353, 1) Y:1 (Δ=0.008) X:7 (Δ=0.012)
peak id: 32      (35.64302954, 236.68204087, 1) Y:1 (Δ=0.012) X:8 (Δ=0.008)

[29.3, 6.0, 2.052476167678833, 0]

The grid looks very well aligned with the axes which makes this step easier but nothing garanties it is perfect, so the rotation of the grid has to be measured as well.

The default rotation will be zero and will be fitted later on.

Once the indexes X,Y determined for eack peak, one can fit the parameter to properly align the grid with the first module. Those 4 parameters are step-size, x_min, y_min and angle

# Align center grid on module #9
indexed9, guess9 = index_module(9, 1, step)
guess9

offset for the first peak:  1983.78669282794 10.031553290784359
peak id:  0    (126.94783784, 1984.12777565, 9) Y:4 (Δ=-0.010) X:0 (Δ=0.012)
peak id:  1     (126.8751576, 2013.45500252, 9) Y:4 (Δ=-0.012) X:1 (Δ=0.013)
peak id:  2    (126.95567599, 2042.78409229, 9) Y:4 (Δ=-0.009) X:2 (Δ=0.014)
peak id:  3    (127.07517209, 2072.16480435, 9) Y:4 (Δ=-0.005) X:3 (Δ=0.016)
peak id:  4    (127.15928899, 2130.88592764, 9) Y:4 (Δ=-0.002) X:5 (Δ=0.020)
peak id:  5    (127.15745774, 2189.65744391, 9) Y:4 (Δ=-0.003) X:7 (Δ=0.026)
peak id:  6    (156.13855915, 1984.19626783, 9) Y:5 (Δ=-0.013) X:0 (Δ=0.014)
peak id:  7     (156.11295008, 2013.5757513, 9) Y:5 (Δ=-0.014) X:1 (Δ=0.017)
peak id:  8    (156.25047338, 2042.90473287, 9) Y:5 (Δ=-0.010) X:2 (Δ=0.018)
peak id:  9     (156.26182449, 2072.2062602, 9) Y:5 (Δ=-0.009) X:3 (Δ=0.018)
peak id: 10     (156.3414048, 2130.94891568, 9) Y:5 (Δ=-0.006) X:5 (Δ=0.023)
peak id: 11    (156.41063583, 2189.68413225, 9) Y:5 (Δ=-0.004) X:7 (Δ=0.027)
peak id: 12     (10.04249137, 1983.78669283, 9) Y:0 (Δ=0.000) X:0 (Δ=0.000)
peak id: 13      (10.0950248, 2013.15991376, 9) Y:0 (Δ=0.002) X:1 (Δ=0.002)
peak id: 14     (10.03155329, 2042.56048054, 9) Y:0 (Δ=0.000) X:2 (Δ=0.006)
peak id: 15     (10.15170847, 2071.91228989, 9) Y:0 (Δ=0.004) X:3 (Δ=0.008)
peak id: 16     (10.18692389, 2130.68992156, 9) Y:0 (Δ=0.005) X:5 (Δ=0.014)
peak id: 17      (10.17274876, 2189.4159292, 9) Y:0 (Δ=0.005) X:7 (Δ=0.018)
peak id: 18    (185.58678433, 2131.00290171, 9) Y:6 (Δ=-0.008) X:5 (Δ=0.024)
peak id: 19     (185.62742698, 2189.6441454, 9) Y:6 (Δ=-0.007) X:7 (Δ=0.026)
peak id: 20    (185.32707188, 1984.23386431, 9) Y:6 (Δ=-0.017) X:0 (Δ=0.015)
peak id: 21    (185.31160745, 2013.65180957, 9) Y:6 (Δ=-0.018) X:1 (Δ=0.019)
peak id: 22    (185.40640566, 2042.94130556, 9) Y:6 (Δ=-0.015) X:2 (Δ=0.019)
peak id: 23    (185.43761739, 2072.24295565, 9) Y:6 (Δ=-0.013) X:3 (Δ=0.019)
peak id: 24      (39.3746087, 2189.38625699, 9) Y:1 (Δ=0.001) X:7 (Δ=0.017)
peak id: 25     (39.22990538, 1983.90723486, 9) Y:1 (Δ=-0.003) X:0 (Δ=0.004)
peak id: 26     (39.34481928, 2013.20239243, 9) Y:1 (Δ=0.000) X:1 (Δ=0.004)
peak id: 27     (39.29185569, 2042.61863551, 9) Y:1 (Δ=-0.001) X:2 (Δ=0.008)
peak id: 28     (39.30902126, 2071.98339974, 9) Y:1 (Δ=-0.001) X:3 (Δ=0.010)
peak id: 29     (39.39395407, 2130.70946553, 9) Y:1 (Δ=0.002) X:5 (Δ=0.014)
peak id: 30      (68.53318956, 1983.9473334, 9) Y:2 (Δ=-0.003) X:0 (Δ=0.005)
peak id: 31     (68.52492085, 2013.39607638, 9) Y:2 (Δ=-0.004) X:1 (Δ=0.011)
peak id: 32     (68.56618777, 2042.62876877, 9) Y:2 (Δ=-0.002) X:2 (Δ=0.008)
peak id: 33      (68.6289348, 2072.04879676, 9) Y:2 (Δ=-0.000) X:3 (Δ=0.012)
peak id: 34      (68.67246765, 2130.7804866, 9) Y:2 (Δ=0.001) X:5 (Δ=0.017)
peak id: 35     (68.60224375, 2189.44976753, 9) Y:2 (Δ=-0.001) X:7 (Δ=0.019)

[29.3, 10.031553290784359, 1983.78669282794, 0]

The error in positionning each of the pixel is less than 0.1 pixel which is already excellent and will allow a straight forward fit.

Let’s do th same for the left and right grid positions:

# Align left grid on module #6
indexed6, guess6 = index_module(6, 0, step)
guess6

offset for the first peak:  1262.76739102602 6.5021456480026245
peak id:  0     (65.23745903, 1262.83622795, 6) Y:2 (Δ=0.005) X:0 (Δ=0.002)
peak id:  1     (65.13842112, 1321.38845757, 6) Y:2 (Δ=0.001) X:2 (Δ=0.001)
peak id:  2     (65.14430301, 1380.07568739, 6) Y:2 (Δ=0.001) X:4 (Δ=0.004)
peak id:  3     (65.07097046, 1409.40630868, 6) Y:2 (Δ=-0.001) X:5 (Δ=0.005)
peak id:  4     (65.06450427, 1438.74234235, 6) Y:2 (Δ=-0.001) X:6 (Δ=0.006)
peak id:  5     (65.12864372, 1468.09552926, 6) Y:2 (Δ=0.001) X:7 (Δ=0.008)
peak id:  6      (6.67313325, 1262.76739103, 6) Y:0 (Δ=0.006) X:0 (Δ=0.000)
peak id:  7      (6.55283955, 1321.53733033, 6) Y:0 (Δ=0.002) X:2 (Δ=0.006)
peak id:  8      (6.62240306, 1380.14368123, 6) Y:0 (Δ=0.004) X:4 (Δ=0.006)
peak id:  9      (6.50214565, 1409.53155875, 6) Y:0 (Δ=0.000) X:5 (Δ=0.009)
peak id: 10      (6.61714479, 1438.72696936, 6) Y:0 (Δ=0.004) X:6 (Δ=0.005)
peak id: 11      (6.62297457, 1468.09986952, 6) Y:0 (Δ=0.004) X:7 (Δ=0.008)
peak id: 12    (123.86860578, 1262.77460632, 6) Y:4 (Δ=0.006) X:0 (Δ=0.000)
peak id: 13    (123.73915467, 1321.47262728, 6) Y:4 (Δ=0.001) X:2 (Δ=0.004)
peak id: 14    (123.74013492, 1380.08734463, 6) Y:4 (Δ=0.001) X:4 (Δ=0.004)
peak id: 15    (123.61760235, 1409.49391153, 6) Y:4 (Δ=-0.003) X:5 (Δ=0.008)
peak id: 16    (123.64603654, 1438.67660466, 6) Y:4 (Δ=-0.002) X:6 (Δ=0.004)
peak id: 17    (123.68760154, 1468.00486562, 6) Y:4 (Δ=-0.000) X:7 (Δ=0.005)
peak id: 18    (153.12439226, 1262.79544121, 6) Y:5 (Δ=0.004) X:0 (Δ=0.001)
peak id: 19    (153.01401874, 1321.43544179, 6) Y:5 (Δ=0.000) X:2 (Δ=0.002)
peak id: 20    (153.00659179, 1380.08258364, 6) Y:5 (Δ=0.000) X:4 (Δ=0.004)
peak id: 21    (152.92952084, 1409.37969038, 6) Y:5 (Δ=-0.002) X:5 (Δ=0.004)
peak id: 22    (152.92040372, 1438.72945455, 6) Y:5 (Δ=-0.003) X:6 (Δ=0.006)
peak id: 23    (152.92236792, 1468.03059884, 6) Y:5 (Δ=-0.003) X:7 (Δ=0.006)
peak id: 24    (182.43571654, 1262.84406386, 6) Y:6 (Δ=0.005) X:0 (Δ=0.003)
peak id: 25    (182.29259774, 1321.38074964, 6) Y:6 (Δ=-0.000) X:2 (Δ=0.000)
peak id: 26    (182.29060692, 1380.04957829, 6) Y:6 (Δ=-0.000) X:4 (Δ=0.003)
peak id: 27     (182.18869503, 1409.3573423, 6) Y:6 (Δ=-0.004) X:5 (Δ=0.003)
peak id: 28    (182.18435012, 1438.72177231, 6) Y:6 (Δ=-0.004) X:6 (Δ=0.005)
peak id: 29    (182.15922907, 1467.99772772, 6) Y:6 (Δ=-0.005) X:7 (Δ=0.004)
peak id: 30     (35.95338714, 1262.80419599, 6) Y:1 (Δ=0.005) X:0 (Δ=0.001)
peak id: 31     (35.82896177, 1321.46761093, 6) Y:1 (Δ=0.001) X:2 (Δ=0.003)
peak id: 32      (35.8987331, 1380.10332788, 6) Y:1 (Δ=0.003) X:4 (Δ=0.005)
peak id: 33     (35.81407647, 1409.42496902, 6) Y:1 (Δ=0.000) X:5 (Δ=0.005)
peak id: 34      (35.8450772, 1438.75288776, 6) Y:1 (Δ=0.001) X:6 (Δ=0.006)
peak id: 35      (35.82283981, 1468.1203842, 6) Y:1 (Δ=0.001) X:7 (Δ=0.009)

[29.3, 6.5021456480026245, 1262.76739102602, 0]

# Align right grid on module #12
indexed12, guess12 = index_module(12, 2, step)
guess12

offset for the first peak:  2852.8019351810217 8.436505198478699
peak id:  0    (66.94455156, 2852.87055804, 12) Y:2 (Δ=-0.003) X:0 (Δ=0.002)
peak id:  1    (66.97839728, 2882.18239576, 12) Y:2 (Δ=-0.002) X:1 (Δ=0.003)
peak id:  2    (66.97670252, 2911.47377473, 12) Y:2 (Δ=-0.002) X:2 (Δ=0.002)
peak id:  3    (67.10628749, 2940.83689539, 12) Y:2 (Δ=0.002) X:3 (Δ=0.005)
peak id:  4   (154.74125573, 2911.58230206, 12) Y:5 (Δ=-0.007) X:2 (Δ=0.006)
peak id:  5   (154.89283001, 2940.90798742, 12) Y:5 (Δ=-0.001) X:3 (Δ=0.007)
peak id:  6   (154.68613148, 2853.03283477, 12) Y:5 (Δ=-0.009) X:0 (Δ=0.008)
peak id:  7   (154.71002036, 2882.29963639, 12) Y:5 (Δ=-0.008) X:1 (Δ=0.007)
peak id:  8     (8.64040324, 2940.72405878, 12) Y:0 (Δ=0.007) X:3 (Δ=0.001)
peak id:  9     (8.53160283, 2911.47199836, 12) Y:0 (Δ=0.003) X:2 (Δ=0.002)
peak id: 10      (8.56915843, 2882.1179322, 12) Y:0 (Δ=0.005) X:1 (Δ=0.001)
peak id: 11   (183.88110682, 2853.06453107, 12) Y:6 (Δ=-0.012) X:0 (Δ=0.009)
peak id: 12   (183.92152359, 2882.34221697, 12) Y:6 (Δ=-0.011) X:1 (Δ=0.008)
peak id: 13   (183.98796783, 2911.64271811, 12) Y:6 (Δ=-0.008) X:2 (Δ=0.008)
peak id: 14   (184.08316714, 2940.95591466, 12) Y:6 (Δ=-0.005) X:3 (Δ=0.009)
peak id: 15      (8.4365052, 2852.80193518, 12) Y:0 (Δ=0.000) X:0 (Δ=0.000)
peak id: 16   (125.41962573, 2852.95906276, 12) Y:4 (Δ=-0.007) X:0 (Δ=0.005)
peak id: 17   (125.43723008, 2882.23844126, 12) Y:4 (Δ=-0.007) X:1 (Δ=0.005)
peak id: 18   (125.44167745, 2911.44268954, 12) Y:4 (Δ=-0.007) X:2 (Δ=0.001)
peak id: 19   (125.64018717, 2940.84172335, 12) Y:4 (Δ=0.000) X:3 (Δ=0.005)
peak id: 20    (37.70449063, 2852.81877436, 12) Y:1 (Δ=-0.001) X:0 (Δ=0.001)
peak id: 21    (37.76720728, 2882.15749501, 12) Y:1 (Δ=0.001) X:1 (Δ=0.002)
peak id: 22     (37.73162082, 2911.4721604, 12) Y:1 (Δ=-0.000) X:2 (Δ=0.002)
peak id: 23    (37.85726659, 2940.76938577, 12) Y:1 (Δ=0.004) X:3 (Δ=0.002)

[29.3, 8.436505198478699, 2852.8019351810217, 0]

The error in positionning each of the pixel is less than 0.1 pixel which is already excellent and will allow a straight forward fit.

All rotations will be performed around the center of each module:

#Calculate the center of every single module for rotation around this center.
centers = {i: numpy.array([[numpy.where(mid == i)[1].mean()], [numpy.where(mid == i)[0].mean()]]) for i in range(1, 19)}
for k,v in centers.items():
    print(k,v.ravel())

1 [121.  97.]
2 [364.5  97. ]
3 [615.  97.]
4 [858.5  97. ]
5 [1109.   97.]
6 [1352.5   97. ]
7 [1603.   97.]
8 [1846.5   97. ]
9 [2097.   97.]
10 [2340.5   97. ]
11 [2591.   97.]
12 [2834.5   97. ]
13 [3085.   97.]
14 [3328.5   97. ]
15 [3579.   97.]
16 [3822.5   97. ]
17 [4073.   97.]
18 [4316.5   97. ]

# Define a rotation of a module around the center of the module ...

def rotate(angle, xy, module):
    "Perform the rotation of the xy points around the center of the given module"
    rot = [[cos(angle),-sin(angle)],
           [sin(angle), cos(angle)]]
    center = centers[module]
    return numpy.dot(rot, xy - center) + center

Fit the grid on a reference module for each position

The cost function for the reference module for each grid position is calculated as the sum of distances squared in pixel space. It uses 4 parameters which are step-size, y_min, x_min, and angle

def cost_grid(param, module, indexed):
    """Cost function for moving the grid on a reference module
    contains: step, y_min, x_min, angle
    returns the sum of distance squared in pixel space
    """
    step = param[0]
    y_min = param[1]
    x_min = param[2]
    angle = param[3]
    XY = numpy.vstack((indexed["X"], indexed["Y"]))
    xy_min = [[x_min], [y_min]]
    xy_guess = rotate(angle, step * XY + xy_min, module)
    delta = xy_guess - numpy.vstack((indexed["x"], indexed["y"]))
    return (delta*delta).sum()

def fit_grid(module, position, guess, indexed):
    """
    :param module:

    :return: fit result of the grid aligned on the module
    """
    where = "left center right".split()[position]
    print(f"Align {where} position on module #{module}")
    print(f"Before optimization {guess} cost=", cost_grid(guess, module, indexed))
    res = minimize(cost_grid, guess, (module, indexed), method = "slsqp")
    print(res)
    print("Average displacement (pixels): ",sqrt(res.fun/len(indexed)))
    return res

# Alignment on the left side
res6 = fit_grid(6, 0, guess6, indexed6)
print("#"*50)
# Alignment on the center
res9 = fit_grid(9, 1, guess9, indexed9)
print("#"*50)
# Alignment on the right side
res12 = fit_grid(12, 2, guess12, indexed12)

Align left position on module #6
Before optimization [29.3, 6.5021456480026245, 1262.76739102602, 0] cost= 1.0084329241894014
     fun: 0.40058707550013467
     jac: array([-0.0014098 , -0.00021778, -0.00010217, -0.04366066])
 message: 'Optimization terminated successfully'
    nfev: 39
     nit: 5
    njev: 5
  status: 0
 success: True
       x: array([ 2.93026713e+01,  6.52410821e+00,  1.26288272e+03, -2.52541590e-04])
Average displacement (pixels):  0.10548658096598178
##################################################
Align center position on module #9
Before optimization [29.3, 10.031553290784359, 1983.78669282794, 0] cost= 9.755912514600535
     fun: 2.6138054098617345
     jac: array([-0.00210074, -0.00039294, -0.00011373, -0.00774959])
 message: 'Optimization terminated successfully'
    nfev: 43
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([ 2.92981305e+01,  9.88812013e+00,  1.98421340e+03, -4.07651924e-04])
Average displacement (pixels):  0.26945445965782905
##################################################
Align right position on module #12
Before optimization [29.3, 8.436505198478699, 2852.8019351810217, 0] cost= 1.228755947570468
     fun: 0.3619580972987666
     jac: array([-3.33767384e-04,  6.27711415e-06, -7.71209598e-05, -1.86840594e-02])
 message: 'Optimization terminated successfully'
    nfev: 43
     nit: 6
    njev: 6
  status: 0
 success: True
       x: array([ 2.92539932e+01,  8.52560255e+00,  2.85299471e+03, -5.44452215e-04])
Average displacement (pixels):  0.12280711456364655

At this step, the grid is perfectly aligned with the first module. This module is used as the reference one and all other are aligned along it, using this first fit:

def index_positon(ref_module, position, guess, thres=1):
    """Return the indexed peak position for all the module for the grid position

    """
    where = "left center right".split()[position]
    print(f"Aligning all peak positions of {where} position on module #{ref_module}")
    step, y_min, x_min, angle = guess
    xy = numpy.vstack((yxi[position]["x"], yxi[position]["y"]))
    indexed = numpy.zeros(xy.shape[1], dtype=dl)
    indexed["y"] = yxi[position]["y"]
    indexed["x"] = yxi[position]["x"]
    indexed["i"] = yxi[position]["i"]
    xy_min = [[x_min], [y_min]]
    XY_app = (rotate(-angle, xy, module=ref_module) - xy_min) / step
    XY_int = numpy.round((XY_app)).astype("int")
    indexed["X"] = XY_int[0]
    indexed["Y"] = XY_int[1]
    xy_guess = rotate(angle, step * XY_int + xy_min, module=ref_module)
    err = (XY_app - XY_int)*step
    delta = numpy.sqrt((err**2).sum(axis=0))
    print(f"suspicious: {(delta>thres).sum()} / {delta.size} with threshold {thres}:")
    suspicious = indexed[numpy.where(abs(delta>thres))]
    print(suspicious)
    return indexed

print(len(index_positon(9, 1, guess=res9.x, thres=3)))
print(len(index_positon(6, 0, guess=res6.x, thres=1)))
print(len(index_positon(12, 2, guess=res12.x, thres=6)))

Aligning all peak positions of center position on module #9
suspicious: 9 / 242 with threshold 3:
[(11.02989381, 2866.03500895, 12, 0, 30)
 (11.06291354, 2924.90754639, 12, 0, 32)
 (11.30030179, 2983.61285341, 13, 0, 34)
 (40.42757976, 2983.54515862, 13, 1, 34)
 (40.2171582 , 2865.98276701, 12, 1, 30)
 (40.24599957, 2924.85079978, 12, 1, 32)
 (69.6720044 , 2983.31674147, 13, 2, 34)
 (69.36312002, 2865.95335668, 12, 2, 30)
 (69.39535499, 2924.79366741, 12, 2, 32)]
242
Aligning all peak positions of left position on module #6
suspicious: 5 / 248 with threshold 1:
[( 35.02235405,   2.8385089 , 1, 1, -43)
 (  7.        , 412.        , 2, 0, -29)
 (  6.        ,  90.        , 1, 0, -40)
 (152.51041356,  31.59106046, 1, 5, -42)
 (181.66608471,   2.49644381, 1, 6, -43)]
248
Aligning all peak positions of right position on module #12
suspicious: 7 / 252 with threshold 6:
[( 40.99528235, 4349.52190217, 18, 1, 51)
 ( 41.14180143, 4408.21720226, 18, 1, 53)
 ( 11.92475192, 4408.41050243, 18, 0, 53)
 ( 11.79140341, 4349.55060229, 18, 0, 51)
 ( 11.77533133, 4290.84790458, 18, 0, 49)
 ( 70.3750506 , 4408.0958024 , 18, 2, 53)
 (128.8388529 , 4407.83846572, 18, 4, 53)]
252

Only 7 peaks have an initial displacement of more than 6 pixels, all located in module 18, which is the furthest away from module 12. All other assignment are in 3 pixel on the left side and 1 pixel in the central position. The visual inspection confirms all localizations are valid.

There are 18 (half-)modules which have each of them 2 translations and one rotation. Only 7 of the are fitted in the first step. In addition to the step size, this represents 22 degrees of freedom for the fit. The first module is used to align the grid, all other modules are then aligned along this grid.

#his contains all peaks with their index
indexed = Triplet(index_positon(6, 0, guess=res6.x, thres=7),
                  index_positon(9, 1, guess=res9.x, thres=7),
                  index_positon(12, 2, guess=res12.x, thres=7))

Aligning all peak positions of left position on module #6
suspicious: 0 / 248 with threshold 7:
[]
Aligning all peak positions of center position on module #9
suspicious: 0 / 242 with threshold 7:
[]
Aligning all peak positions of right position on module #12
suspicious: 0 / 252 with threshold 7:
[]

The submodule cost funtion is the sum of the squares of the difference in the pixel space:

#here are defined the piot point for switching reference:
pivots = [6, 9, 12]

def submodule_cost(param, module, position):
    """contains: step,
                 y_min_6L, x_min_6L, angle_6L,
                 y_min_9C, x_min_9C, angle_9C,
                 y_min_12R, x_min_12R, angle_12R,

                 y_min_1, x_min_1, angle_1,
                 y_min_2, x_min_2, angle_2,
                 y_min_3, x_min_3, angle_3, ...

    :param: array with 64 parameters
    :param module: module number from 1 to 18
    :param position: 0, 1 or 2
    :returns: the sum of distance squared in pixel space for the given module with the given grid.

    """

    step = param[0]
    mask = indexed[position]["i"] == module
    if mask.sum() == 0 :
        return 0
    substack = indexed[position][mask]
    pivot = pivots[position]

    y_min_grid, x_min_grid, angle_grid = param[1+3*(position): 4+3*(position)]

    XY = numpy.vstack((substack["X"], substack["Y"]))

    xy_min_grid = numpy.array([[x_min_grid], [y_min_grid]])
    xy_guess1 = rotate(angle_grid, step * XY + xy_min_grid, pivot)

    #print(y_min_grid, x_min_grid, angle_grid)

    if module == pivot:
        #print("Not much to do for reference module as it is naturally alligned")
        delta = xy_guess1 - numpy.vstack((substack["x"], substack["y"]))
    else:
        "perform the correction for given module"
        y_min, x_min, angle = param[7+3*module: 10+3*module]
        xy_min = numpy.array([[x_min], [y_min]])
        xy_guess = rotate(angle, xy_guess1 + xy_min, module)
        delta = xy_guess - numpy.vstack((substack["x"], substack["y"]))
        #print(y_min, x_min, angle)
    return (delta*delta).sum()

def print_res64(param):
    res = ["step: %.3f"%param[0]]
    f=lambda p:"Δx: %8.3f, Δy: %6.3f rot: %6.3f°"%(p[1], p[0], numpy.rad2deg(p[2]))
    res.append(f" {pivots[0]}L: {f(param[1:4])}")
    res.append(f" {pivots[1]}C: {f(param[4:7])}")
    res.append(f"{pivots[2]}R: {f(param[7:10])}")
    for i in range(1,19):
        w = "L" if i<min(pivots) else "R" if i>max(pivots) else "C"
        res.append("%2s%s: "%(i,w)+f(param[7+3*i:10+3*i]))
    print("\n".join(res))

submodule_cost(numpy.zeros(64), 1, 0)

1155395.4186643849

# Evaluated the guess and print the cost:
guess64 = numpy.zeros(64)
guess64[:4] = res6.x
guess64[4:7] = res9.x[1:]
guess64[7:10] = res12.x[1:]
print_res64(guess64)

step: 29.303
 6L: Δx: 1262.883, Δy:  6.524 rot: -0.014°
 9C: Δx: 1984.213, Δy:  9.888 rot: -0.023°
12R: Δx: 2852.995, Δy:  8.526 rot: -0.031°
 1L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 2L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 3L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 4L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 5L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 6C: Δx:    0.000, Δy:  0.000 rot:  0.000°
 7C: Δx:    0.000, Δy:  0.000 rot:  0.000°
 8C: Δx:    0.000, Δy:  0.000 rot:  0.000°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.000, Δy:  0.000 rot:  0.000°
11C: Δx:    0.000, Δy:  0.000 rot:  0.000°
12C: Δx:    0.000, Δy:  0.000 rot:  0.000°
13R: Δx:    0.000, Δy:  0.000 rot:  0.000°
14R: Δx:    0.000, Δy:  0.000 rot:  0.000°
15R: Δx:    0.000, Δy:  0.000 rot:  0.000°
16R: Δx:    0.000, Δy:  0.000 rot:  0.000°
17R: Δx:    0.000, Δy:  0.000 rot:  0.000°
18R: Δx:    0.000, Δy:  0.000 rot:  0.000°

print("\nContribution to the total cost of each module/position")
for m in range(1,19):
    print(m, "%10f\t %10f\t %10f"%tuple(submodule_cost(guess64, m, i) for i in range(3)))

Contribution to the total cost of each module/position
1  24.826181       0.000000        0.000000
2  30.830874       0.000000        0.000000
3   7.344046       0.000000        0.000000
4   7.663464       0.000000        0.000000
5   2.550972       0.000000        0.000000
6   0.400587      46.645569        0.000000
7   6.197929      54.642847        0.000000
8   0.000000      26.282742        0.000000
9   0.000000       2.634828        0.000000
10   0.000000     23.952726        0.000000
11   0.000000    145.809223        0.000000
12   0.000000    172.924027        1.338200
13   0.000000     56.798298       20.032284
14   0.000000      0.000000      138.498766
15   0.000000      0.000000      223.065279
16   0.000000      0.000000      274.254628
17   0.000000      0.000000      639.247608
18   0.000000      0.000000      527.538219

Fit all modules of the central grid position:

# Fit the center position:

def cost_all_center(param):
    """contains: step, y_min_1, x_min_1, angle_1, ...
    returns the sum of distance squared in pixel space
    """
    return sum(submodule_cost(param, module=i, position=1) for i in range(6,13))

print_res64(guess64)
print(cost_all_center(guess64))
res_center = minimize(cost_all_center, guess64, method = "slsqp")
print_res64(res_center.x)
print(res_center)

step: 29.303
 6L: Δx: 1262.883, Δy:  6.524 rot: -0.014°
 9C: Δx: 1984.213, Δy:  9.888 rot: -0.023°
12R: Δx: 2852.995, Δy:  8.526 rot: -0.031°
 1L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 2L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 3L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 4L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 5L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 6C: Δx:    0.000, Δy:  0.000 rot:  0.000°
 7C: Δx:    0.000, Δy:  0.000 rot:  0.000°
 8C: Δx:    0.000, Δy:  0.000 rot:  0.000°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.000, Δy:  0.000 rot:  0.000°
11C: Δx:    0.000, Δy:  0.000 rot:  0.000°
12C: Δx:    0.000, Δy:  0.000 rot:  0.000°
13R: Δx:    0.000, Δy:  0.000 rot:  0.000°
14R: Δx:    0.000, Δy:  0.000 rot:  0.000°
15R: Δx:    0.000, Δy:  0.000 rot:  0.000°
16R: Δx:    0.000, Δy:  0.000 rot:  0.000°
17R: Δx:    0.000, Δy:  0.000 rot:  0.000°
18R: Δx:    0.000, Δy:  0.000 rot:  0.000°
472.8919615302689
step: 29.287
 6L: Δx: 1262.883, Δy:  6.524 rot: -0.014°
 9C: Δx: 1984.248, Δy:  9.923 rot: -0.023°
12R: Δx: 2852.995, Δy:  8.526 rot: -0.031°
 1L: Δx:   -0.000, Δy:  0.000 rot: -0.000°
 2L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 3L: Δx:   -0.000, Δy:  0.000 rot:  0.000°
 4L: Δx:   -0.000, Δy:  0.000 rot:  0.000°
 5L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 6C: Δx:   -1.038, Δy: -1.262 rot: -0.319°
 7C: Δx:   -0.455, Δy: -1.318 rot:  0.016°
 8C: Δx:   -0.136, Δy: -0.641 rot:  0.004°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.798, Δy:  0.357 rot:  0.010°
11C: Δx:    1.866, Δy:  0.705 rot:  0.062°
12C: Δx:    2.872, Δy:  1.060 rot:  0.095°
13R: Δx:    0.000, Δy:  0.000 rot:  0.000°
14R: Δx:    0.000, Δy:  0.000 rot:  0.000°
15R: Δx:    0.000, Δy:  0.000 rot:  0.000°
16R: Δx:    0.000, Δy:  0.000 rot:  0.000°
17R: Δx:    0.000, Δy:  0.000 rot:  0.000°
18R: Δx:    0.000, Δy:  0.000 rot:  0.000°
     fun: 24.99843343683122
     jac: array([-1.03936195e-01,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        1.02663040e-03, -1.89185143e-03,  7.82108307e-01,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00, -6.29425049e-05,  6.07967377e-04, -6.62517548e-03,
       -4.48226929e-05,  7.52449036e-04, -6.42776489e-04,  8.05854797e-05,
        5.07831573e-04,  7.46250153e-04,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  2.06470490e-04, -8.49723816e-04,  1.12357140e-02,
        5.77926636e-04, -1.83582306e-03, -4.41551208e-04,  3.47137451e-04,
       -1.29771233e-03,  7.09772110e-04,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00])
 message: 'Optimization terminated successfully'
    nfev: 1701
     nit: 25
    njev: 25
  status: 0
 success: True
       x: array([ 2.92866114e+01,  6.52410821e+00,  1.26288272e+03, -2.52541590e-04,
        9.92266288e+00,  1.98424796e+03, -4.08389445e-04,  8.52560255e+00,
        2.85299471e+03, -5.44452215e-04,  4.43141471e-16, -3.56524178e-15,
       -4.42484061e-16,  7.65225719e-17,  4.90227298e-16,  1.05097212e-15,
        1.00843762e-16, -1.63522698e-16,  3.98694271e-16,  8.90225607e-17,
       -2.29036169e-16,  1.17762056e-16,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00, -1.26206052e+00, -1.03837867e+00, -5.56579208e-03,
       -1.31779372e+00, -4.55166603e-01,  2.70894309e-04, -6.41391927e-01,
       -1.36493114e-01,  6.63013003e-05,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  3.57324229e-01,  7.97695112e-01,  1.69995912e-04,
        7.05362252e-01,  1.86623546e+00,  1.08590970e-03,  1.06025620e+00,
        2.87224614e+00,  1.65628899e-03,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00])

Fit all modules of the left grid position:

def cost_all_left(param):
    """contains: step, y_min_1, x_min_1, angle_1, ...
    returns the sum of distance squared in pixel space
    """
    return sum(submodule_cost(param, module=i, position=0) for i in range(1,6))

print_res64(res_center.x)
print(cost_all_left(res_center.x))
res_left = minimize(cost_all_left, res_center.x, method = "slsqp")
print_res64(res_left.x)
print(res_left)

step: 29.287
 6L: Δx: 1262.883, Δy:  6.524 rot: -0.014°
 9C: Δx: 1984.248, Δy:  9.923 rot: -0.023°
12R: Δx: 2852.995, Δy:  8.526 rot: -0.031°
 1L: Δx:   -0.000, Δy:  0.000 rot: -0.000°
 2L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 3L: Δx:   -0.000, Δy:  0.000 rot:  0.000°
 4L: Δx:   -0.000, Δy:  0.000 rot:  0.000°
 5L: Δx:    0.000, Δy:  0.000 rot:  0.000°
 6C: Δx:   -1.038, Δy: -1.262 rot: -0.319°
 7C: Δx:   -0.455, Δy: -1.318 rot:  0.016°
 8C: Δx:   -0.136, Δy: -0.641 rot:  0.004°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.798, Δy:  0.357 rot:  0.010°
11C: Δx:    1.866, Δy:  0.705 rot:  0.062°
12C: Δx:    2.872, Δy:  1.060 rot:  0.095°
13R: Δx:    0.000, Δy:  0.000 rot:  0.000°
14R: Δx:    0.000, Δy:  0.000 rot:  0.000°
15R: Δx:    0.000, Δy:  0.000 rot:  0.000°
16R: Δx:    0.000, Δy:  0.000 rot:  0.000°
17R: Δx:    0.000, Δy:  0.000 rot:  0.000°
18R: Δx:    0.000, Δy:  0.000 rot:  0.000°
177.2287870711433
step: 29.318
 6L: Δx: 1262.815, Δy:  6.800 rot:  0.018°
 9C: Δx: 1984.248, Δy:  9.923 rot: -0.023°
12R: Δx: 2852.995, Δy:  8.526 rot: -0.031°
 1L: Δx:    0.128, Δy: -0.273 rot:  0.035°
 2L: Δx:   -0.215, Δy: -0.023 rot:  0.005°
 3L: Δx:   -0.014, Δy:  0.313 rot:  0.048°
 4L: Δx:   -0.023, Δy:  0.242 rot: -0.044°
 5L: Δx:    0.057, Δy:  0.017 rot: -0.023°
 6C: Δx:   -1.038, Δy: -1.262 rot: -0.319°
 7C: Δx:   -0.455, Δy: -1.318 rot:  0.016°
 8C: Δx:   -0.136, Δy: -0.641 rot:  0.004°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.798, Δy:  0.357 rot:  0.010°
11C: Δx:    1.866, Δy:  0.705 rot:  0.062°
12C: Δx:    2.872, Δy:  1.060 rot:  0.095°
13R: Δx:    0.000, Δy:  0.000 rot:  0.000°
14R: Δx:    0.000, Δy:  0.000 rot:  0.000°
15R: Δx:    0.000, Δy:  0.000 rot:  0.000°
16R: Δx:    0.000, Δy:  0.000 rot:  0.000°
17R: Δx:    0.000, Δy:  0.000 rot:  0.000°
18R: Δx:    0.000, Δy:  0.000 rot:  0.000°
     fun: 3.1468682956116614
     jac: array([-2.55170465e-03,  1.39370561e-03, -4.24772501e-04, -3.12974602e-01,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  1.19021535e-03,  2.53379345e-04,
       -9.42820311e-03,  8.86440277e-04, -6.56247139e-04, -1.65006518e-03,
        4.60982323e-04,  9.81301069e-04, -5.31286001e-03, -8.94248486e-04,
       -1.64186954e-03,  7.08445907e-03, -3.13997269e-04,  6.37471676e-04,
       -1.17188096e-02,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00])
 message: 'Optimization terminated successfully'
    nfev: 1309
     nit: 19
    njev: 19
  status: 0
 success: True
       x: array([ 2.93179468e+01,  6.80042099e+00,  1.26281542e+03,  3.07983783e-04,
        9.92266288e+00,  1.98424796e+03, -4.08389445e-04,  8.52560255e+00,
        2.85299471e+03, -5.44452215e-04, -2.73083912e-01,  1.28385037e-01,
        6.18372111e-04, -2.29987461e-02, -2.15395245e-01,  9.08149490e-05,
        3.12562343e-01, -1.42701163e-02,  8.33332628e-04,  2.42427908e-01,
       -2.29928859e-02, -7.62558422e-04,  1.74204378e-02,  5.73986978e-02,
       -4.04518418e-04, -1.26206052e+00, -1.03837867e+00, -5.56579208e-03,
       -1.31779372e+00, -4.55166603e-01,  2.70894309e-04, -6.41391927e-01,
       -1.36493114e-01,  6.63013003e-05,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  3.57324229e-01,  7.97695112e-01,  1.69995912e-04,
        7.05362252e-01,  1.86623546e+00,  1.08590970e-03,  1.06025620e+00,
        2.87224614e+00,  1.65628899e-03,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00])

Fit all modules of the right grid position:

def cost_all_right(param):
    """contains: step, y_min_1, x_min_1, angle_1, ...
    returns the sum of distance squared in pixel space
    """
    return sum(submodule_cost(param, module=i, position=2) for i in range(13,19))

bounds = [(None, None) for i in range(64)]
bounds[0] = (29, 30)
for i in range(3, 64, 3):
    bounds[i] = (-.1, .1) #limit rotations


print_res64(res_left.x)
print(cost_all_right(res_left.x))
res_right = minimize(cost_all_right, res_left.x, method = "slsqp", bounds=bounds)
print_res64(res_right.x)
print(res_right)

step: 29.318
 6L: Δx: 1262.815, Δy:  6.800 rot:  0.018°
 9C: Δx: 1984.248, Δy:  9.923 rot: -0.023°
12R: Δx: 2852.995, Δy:  8.526 rot: -0.031°
 1L: Δx:    0.128, Δy: -0.273 rot:  0.035°
 2L: Δx:   -0.215, Δy: -0.023 rot:  0.005°
 3L: Δx:   -0.014, Δy:  0.313 rot:  0.048°
 4L: Δx:   -0.023, Δy:  0.242 rot: -0.044°
 5L: Δx:    0.057, Δy:  0.017 rot: -0.023°
 6C: Δx:   -1.038, Δy: -1.262 rot: -0.319°
 7C: Δx:   -0.455, Δy: -1.318 rot:  0.016°
 8C: Δx:   -0.136, Δy: -0.641 rot:  0.004°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.798, Δy:  0.357 rot:  0.010°
11C: Δx:    1.866, Δy:  0.705 rot:  0.062°
12C: Δx:    2.872, Δy:  1.060 rot:  0.095°
13R: Δx:    0.000, Δy:  0.000 rot:  0.000°
14R: Δx:    0.000, Δy:  0.000 rot:  0.000°
15R: Δx:    0.000, Δy:  0.000 rot:  0.000°
16R: Δx:    0.000, Δy:  0.000 rot:  0.000°
17R: Δx:    0.000, Δy:  0.000 rot:  0.000°
18R: Δx:    0.000, Δy:  0.000 rot:  0.000°
1624.4976047642317
step: 29.317
 6L: Δx: 1262.815, Δy:  6.800 rot:  0.018°
 9C: Δx: 1984.248, Δy:  9.923 rot: -0.023°
12R: Δx: 2853.232, Δy:  9.558 rot:  0.048°
 1L: Δx:    0.128, Δy: -0.273 rot:  0.035°
 2L: Δx:   -0.215, Δy: -0.023 rot:  0.005°
 3L: Δx:   -0.014, Δy:  0.313 rot:  0.048°
 4L: Δx:   -0.023, Δy:  0.242 rot: -0.044°
 5L: Δx:    0.057, Δy:  0.017 rot: -0.023°
 6C: Δx:   -1.038, Δy: -1.262 rot: -0.319°
 7C: Δx:   -0.455, Δy: -1.318 rot:  0.016°
 8C: Δx:   -0.136, Δy: -0.641 rot:  0.004°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.798, Δy:  0.357 rot:  0.010°
11C: Δx:    1.866, Δy:  0.705 rot:  0.062°
12C: Δx:    2.872, Δy:  1.060 rot:  0.095°
13R: Δx:   -0.687, Δy: -0.704 rot:  0.089°
14R: Δx:   -0.555, Δy: -0.101 rot:  0.044°
15R: Δx:   -0.073, Δy:  0.165 rot:  0.048°
16R: Δx:    0.307, Δy:  0.348 rot:  0.012°
17R: Δx:    0.481, Δy:  0.600 rot:  0.074°
18R: Δx:    0.777, Δy:  0.727 rot:  0.078°
     fun: 7.2981677337707005
     jac: array([-1.36465335e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00, -8.65423679e-02,
        2.14300752e-02, -2.13366965e+02,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  8.15318227e-02,  2.25723982e-02,
        4.11015803e+00,  5.12102246e-02,  2.25861073e-02, -6.15992308e-01,
        5.80726266e-02, -1.39374137e-02,  1.13379014e+01, -1.39135778e-01,
        5.60290217e-02, -4.29959089e+00, -1.00755453e-01, -3.62920761e-03,
       -3.02239978e+00, -3.75187397e-02, -6.21220469e-02,  4.88940239e+00])
 message: 'Optimization terminated successfully'
    nfev: 2314
     nit: 35
    njev: 34
  status: 0
 success: True
       x: array([ 2.93170572e+01,  6.80042099e+00,  1.26281542e+03,  3.07983783e-04,
        9.92266288e+00,  1.98424796e+03, -4.08389445e-04,  9.55770601e+00,
        2.85323223e+03,  8.45691522e-04, -2.73083912e-01,  1.28385037e-01,
        6.18372111e-04, -2.29987461e-02, -2.15395245e-01,  9.08149490e-05,
        3.12562343e-01, -1.42701163e-02,  8.33332628e-04,  2.42427908e-01,
       -2.29928859e-02, -7.62558422e-04,  1.74204378e-02,  5.73986978e-02,
       -4.04518418e-04, -1.26206052e+00, -1.03837867e+00, -5.56579208e-03,
       -1.31779372e+00, -4.55166603e-01,  2.70894309e-04, -6.41391927e-01,
       -1.36493114e-01,  6.63013003e-05,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  3.57324229e-01,  7.97695112e-01,  1.69995912e-04,
        7.05362252e-01,  1.86623546e+00,  1.08590970e-03,  1.06025620e+00,
        2.87224614e+00,  1.65628899e-03, -7.03833830e-01, -6.86912213e-01,
        1.54985730e-03, -1.00918737e-01, -5.55084013e-01,  7.60968791e-04,
        1.64985794e-01, -7.31128161e-02,  8.40076237e-04,  3.48200888e-01,
        3.07113044e-01,  2.01699756e-04,  6.00277182e-01,  4.80659625e-01,
        1.28857896e-03,  7.27422679e-01,  7.77118206e-01,  1.36064347e-03])

Fit all modules togeather

def cost_all(param):
    """contains: step, y_min_1, x_min_1, angle_1, ...
    returns the sum of distance squared in pixel space
    """
    return (sum(submodule_cost(param, module=i, position=0) for i in range(1, 19))+
            sum(submodule_cost(param, module=i, position=1) for i in range(1, 19))+
            sum(submodule_cost(param, module=i, position=2) for i in range(1, 19)))

print_res64(res_right.x)
print(cost_all(res_right.x))
res_all = minimize(cost_all, res_right.x, method = "slsqp")
print_res64(res_all.x)
print(res_all)

step: 29.317
 6L: Δx: 1262.815, Δy:  6.800 rot:  0.018°
 9C: Δx: 1984.248, Δy:  9.923 rot: -0.023°
12R: Δx: 2853.232, Δy:  9.558 rot:  0.048°
 1L: Δx:    0.128, Δy: -0.273 rot:  0.035°
 2L: Δx:   -0.215, Δy: -0.023 rot:  0.005°
 3L: Δx:   -0.014, Δy:  0.313 rot:  0.048°
 4L: Δx:   -0.023, Δy:  0.242 rot: -0.044°
 5L: Δx:    0.057, Δy:  0.017 rot: -0.023°
 6C: Δx:   -1.038, Δy: -1.262 rot: -0.319°
 7C: Δx:   -0.455, Δy: -1.318 rot:  0.016°
 8C: Δx:   -0.136, Δy: -0.641 rot:  0.004°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.798, Δy:  0.357 rot:  0.010°
11C: Δx:    1.866, Δy:  0.705 rot:  0.062°
12C: Δx:    2.872, Δy:  1.060 rot:  0.095°
13R: Δx:   -0.687, Δy: -0.704 rot:  0.089°
14R: Δx:   -0.555, Δy: -0.101 rot:  0.044°
15R: Δx:   -0.073, Δy:  0.165 rot:  0.048°
16R: Δx:    0.307, Δy:  0.348 rot:  0.012°
17R: Δx:    0.481, Δy:  0.600 rot:  0.074°
18R: Δx:    0.777, Δy:  0.727 rot:  0.078°
253.63552675118063
step: 29.337
 6L: Δx: 1262.787, Δy:  6.450 rot:  0.065°
 9C: Δx: 1984.276, Δy:  9.712 rot:  0.054°
12R: Δx: 2852.533, Δy:  8.280 rot:  0.001°
 1L: Δx:    0.904, Δy:  1.034 rot: -0.012°
 2L: Δx:    0.397, Δy:  1.085 rot: -0.042°
 3L: Δx:    0.448, Δy:  1.212 rot:  0.001°
 4L: Δx:    0.256, Δy:  0.943 rot: -0.091°
 5L: Δx:    0.197, Δy:  0.511 rot: -0.070°
 6C: Δx:   -0.131, Δy: -0.197 rot: -0.396°
 7C: Δx:    0.260, Δy: -0.507 rot: -0.075°
 8C: Δx:    0.064, Δy: -0.243 rot: -0.074°
 9C: Δx:    0.000, Δy:  0.000 rot:  0.000°
10C: Δx:    0.120, Δy:  0.087 rot: -0.068°
11C: Δx:    0.802, Δy:  0.098 rot: -0.015°
12C: Δx:    1.379, Δy:  0.125 rot:  0.017°
13R: Δx:    0.136, Δy:  0.698 rot:  0.153°
14R: Δx:   -0.187, Δy:  1.530 rot:  0.092°
15R: Δx:    0.142, Δy:  2.004 rot:  0.095°
16R: Δx:    0.328, Δy:  2.394 rot:  0.060°
17R: Δx:    0.348, Δy:  2.854 rot:  0.122°
18R: Δx:    0.491, Δy:  3.184 rot:  0.125°
     fun: 66.45961832136827
     jac: array([ 1.35029793e-01, -9.63211060e-05,  1.64031982e-04, -4.82873917e-02,
        1.04713440e-03,  3.05175781e-05,  3.52234840e-01,  9.56535339e-04,
        3.71265411e-03,  2.41114616e-01,  1.55448914e-04, -2.86102295e-06,
        5.71308136e-02,  1.30939484e-03,  1.62792206e-03,  1.27172470e-02,
       -2.12001801e-03, -1.89590454e-03, -5.01747131e-02,  7.29560852e-04,
        1.00135803e-04,  1.78995132e-02,  4.86373901e-05,  2.43186951e-04,
        9.49859619e-03, -4.97817993e-04, -1.21116638e-04, -7.12776184e-03,
        7.77244568e-04,  1.78337097e-04, -2.01988220e-03, -3.50952148e-04,
       -8.16345215e-04,  1.49402618e-02,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  5.28335571e-04,  1.20162964e-03, -3.67984772e-02,
       -1.09481812e-03, -1.13105774e-03, -1.35173798e-02, -8.10623169e-05,
       -3.06129456e-04,  2.24399567e-03,  4.03404236e-04, -7.90596008e-04,
        9.23824310e-03, -1.07765198e-04,  1.30939484e-03, -5.34019470e-02,
       -9.72747803e-05, -1.52969360e-03, -8.30364227e-03, -8.97407532e-04,
        1.76715851e-03, -4.48608398e-02,  1.56307220e-03,  2.12574005e-03,
       -8.47148895e-03, -4.23431396e-04,  7.79151917e-04,  2.05459595e-02])
 message: 'Optimization terminated successfully'
    nfev: 4443
     nit: 65
    njev: 65
  status: 0
 success: True
       x: array([ 2.93371266e+01,  6.45023148e+00,  1.26278747e+03,  1.13239278e-03,
        9.71162045e+00,  1.98427565e+03,  9.43710563e-04,  8.27953350e+00,
        2.85253315e+03,  9.74544593e-06,  1.03405136e+00,  9.03817728e-01,
       -2.05950984e-04,  1.08468692e+00,  3.96913106e-01, -7.33585395e-04,
        1.21235726e+00,  4.47974357e-01,  8.76479738e-06,  9.42655736e-01,
        2.56352273e-01, -1.58700475e-03,  5.11024696e-01,  1.97256463e-01,
       -1.22888839e-03, -1.96523436e-01, -1.30596734e-01, -6.91752024e-03,
       -5.06999541e-01,  2.60446600e-01, -1.31291505e-03, -2.43189788e-01,
        6.37054395e-02, -1.28578196e-03,  9.60351311e-16,  4.52603901e-15,
        1.28453865e-15,  8.67808165e-02,  1.20408008e-01, -1.18213812e-03,
        9.84346285e-02,  8.02432359e-01, -2.66208434e-04,  1.25040560e-01,
        1.37881899e+00,  3.04199645e-04,  6.98363043e-01,  1.35938283e-01,
        2.66226477e-03,  1.52992179e+00, -1.87296236e-01,  1.59809116e-03,
        2.00437136e+00,  1.41943838e-01,  1.65825005e-03,  2.39408611e+00,
        3.27570959e-01,  1.03957700e-03,  2.85422874e+00,  3.47848594e-01,
        2.12854841e-03,  3.18408863e+00,  4.90926824e-01,  2.18972897e-03])

print("\nContribution to the total cost of each module/position")
for m in range(1,19):
    print(m, "%10f\t %10f\t %10f"%tuple(submodule_cost(res_all.x, m, i) for i in range(3)))

Contribution to the total cost of each module/position
1   2.039591       0.000000        0.000000
2   1.210239       0.000000        0.000000
3   0.190345       0.000000        0.000000
4   0.183150       0.000000        0.000000
5   0.275581       0.000000        0.000000
6   1.697196       2.102496        0.000000
7   1.282112       7.021678        0.000000
8   0.000000       7.179887        0.000000
9   0.000000       5.218395        0.000000
10   0.000000      3.490917        0.000000
11   0.000000      4.629866        0.000000
12   0.000000      3.613325        4.132154
13   0.000000     11.716505        3.812564
14   0.000000      0.000000        1.124001
15   0.000000      0.000000        0.521402
16   0.000000      0.000000        0.611840
17   0.000000      0.000000        2.221689
18   0.000000      0.000000        2.184686

Reconstruction of the pixel position

The pixel position can be obtained from the standard Pilatus detector. Each module is then displaced according to the fitted values, except the first one which is left where it is.

pixel_coord = pyFAI.detector_factory("Pilatus900kwCdTe").get_pixel_corners()
pixel_coord_raw = pixel_coord.copy()

def correct(x, y, dx, dy, angle, module):
    "apply the correction dx, dy and angle to those pixels ..."
    trans = numpy.array([[dx],
                         [dy]])
    xy_guess = numpy.vstack((x.ravel(),
                             y.ravel()))
    xy_cor = rotate(-angle, xy_guess, module) - trans
    xy_cor.shape = ((2,)+x.shape)
    return xy_cor[0], xy_cor[1]


for module in range(1, 19):
    # Extract the pixel corners for one module
    module_idx = numpy.where(mid == module)
    one_module = pixel_coord_raw[module_idx]


    #retrieve the fitted values
    dy, dx, angle = res_all.x[7+3*module: 10+3*module]

    #z = one_module[...,0]
    y = one_module[...,1]/pilatus.pixel1
    x = one_module[...,2]/pilatus.pixel2

    x_cor, y_cor = correct(x, y, dx, dy, angle, module)
    if i>12:
        dy, dx, angle = res_all.x[7+3*12: 10+3*12]
        x_cor, y_cor = correct(x_cor, y_cor, dx, dy, angle, 12)
    if i<6:
        dy, dx, angle = res_all.x[7+3*6: 10+3*6]
        x_cor, y_cor = correct(x_cor, y_cor, dx, dy, angle, 6)

    one_module[...,1] = y_cor * pilatus.pixel1
    one_module[...,2] = x_cor * pilatus.pixel2
    #Update the array
    pixel_coord_raw[module_idx] = one_module

pilatus.set_pixel_corners(pixel_coord_raw)
pilatus.mask = mask0 | mask1 | mask2
pilatus.save("Pilatus_ID06_raw.h5")

displ = numpy.sqrt(((pixel_coord - pixel_coord_raw)**2).sum(axis=-1))
displ /= pilatus.pixel1 #convert in pixel units
fig, ax = subplots(figsize=(8,6))
ax.hist(displ.ravel(), 100)
ax.set_title("Displacement of pixels versus the reference representation")
ax.set_xlabel("Error in pixel size (172µm)")
pass

#Kabsch alignment of the whole detector ...

unmasked = numpy.logical_not(all_masks)
misaligned = numpy.vstack((pixel_coord_raw[..., 2].ravel(), #x
                           pixel_coord_raw[..., 1].ravel())) #y

reference = numpy.vstack((pixel_coord[..., 2].ravel(), #x
                          pixel_coord[..., 1].ravel())) #y

def kabsch(P, R):
    "Align P on R"
    centroid_P = P.mean(axis=0)
    centroid_R = R.mean(axis=0)
    centered_P = P - centroid_P
    centered_R = R - centroid_R
    C = numpy.dot(centered_P.T, centered_R)
    V, S, W = numpy.linalg.svd(C)
    d = (numpy.linalg.det(V) * numpy.linalg.det(W)) < 0.0
    if d:
        S[-1] = -S[-1]
        V[:, -1] = -V[:, -1]
    # Create Rotation matrix U
    U = numpy.dot(V, W)
    P = numpy.dot(centered_P, U)
    return P + centroid_R

%time aligned = kabsch(misaligned.T, reference.T).T

CPU times: user 57 ms, sys: 15.9 ms, total: 72.9 ms
Wall time: 72.8 ms

displ = numpy.sqrt(((aligned-reference)**2).sum(axis=0))
displ /= pilatus.pixel1 #convert in pixel units
fig, ax = subplots(figsize=(8,6))
ax.hist(displ.ravel(), 100)
ax.set_title("Displacement of pixels versus the reference representation")
ax.set_xlabel("Pixel size (172µm)")
pass

pixel_coord_aligned = pixel_coord.copy()
pixel_coord_aligned[...,1] = aligned[1,:].reshape(pixel_coord.shape[:-1])
pixel_coord_aligned[...,2] = aligned[0,:].reshape(pixel_coord.shape[:-1])

pilatus.set_pixel_corners(pixel_coord_aligned)
pilatus.mask = all_masks
pilatus.save("Pilatus_ID06_final.h5")

fig, ax = subplots(2, figsize=(20, 4))
i0 = ax[0].imshow((pixel_coord_aligned[...,2].mean(axis=-1) - pixel_coord[...,2].mean(axis=-1))/pilatus.pixel2)
ax[0].set_title("Displacement x (in pixel)")
i1 = ax[1].imshow((pixel_coord_aligned[...,1].mean(axis=-1) - pixel_coord[...,1].mean(axis=-1))/pilatus.pixel1)
ax[1].set_title("Displacement y (in pixel)")
fig.colorbar(i0, ax=ax[0])
fig.colorbar(i1, ax=ax[1])
pass

Conclusion

This tutorial presents the way to calibrate a large module based detector using a small grid. The HDF5 file generated is directly useable by any parts of pyFAI, the reader is invited in calibrating the rings images with the default definition and with this optimized definition and check the residual error is almost divided by a factor two.

To come back on the precision of the localization of the pixel: not all the pixel are within the specifications provided by Dectris which claims the misaliment of the modules is within one pixel.

Nota: There is not validation yet of this modelization of the detector. There has been no parallax effect corrections so far.

print(f"Total execution time: {time.perf_counter() - start_time:.3f} s")

Total execution time: 43.272 s