aboutsummaryrefslogtreecommitdiff
path: root/scripts/detector-shift
blob: d9a3fab465b12657dbfdd0b14f8a3cb0b1a5a2f2 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#!/usr/bin/env python

#
# Determine mean detector shift based on prediction refinement results
#
# Copyright (c) 2015 Deutsches Elektronen-Synchrotron DESY,
#                    a research centre of the Helmholtz Association.
#
# Author:
#    2015 Thomas White <taw@physics.org>
#

import sys
import os
import re
import matplotlib.pyplot as plt

f = open(sys.argv[1], 'r')
if len(sys.argv) > 2:
    geom = sys.argv[2]
    have_geom = 1
else:
    have_geom = 0

# Determine the mean shifts
x_shifts = []
y_shifts = []
z_shifts = []

prog1 = re.compile("^predict_refine/det_shift\sx\s=\s([0-9\.\-]+)\sy\s=\s([0-9\.\-]+)\smm$")
prog2 = re.compile("^predict_refine/clen_shift\s=\s([0-9\.\-]+)\smm$")

while True:

    fline = f.readline()
    if not fline:
        break

    match = prog1.match(fline)
    if match:
        xshift = float(match.group(1))
        yshift = float(match.group(2))
        x_shifts.append(xshift)
        y_shifts.append(yshift)

    match = prog2.match(fline)
    if match:
        zshift = float(match.group(1))
        z_shifts.append(zshift)

f.close()

if len(x_shifts) != len(z_shifts):
    print 'Warning: number of xy shifts not equal to number of z shifts'

n = len(z_shifts)
z_shifts = [t for t in z_shifts if abs(t) < 2.0]
if len(z_shifts) < n:
    print 'Warning: %i camera length shifts of more than 2 mm ignored' \
           % (n-len(z_shifts))

mean_x = sum(x_shifts) / len(x_shifts)
mean_y = sum(y_shifts) / len(y_shifts)
mean_z = sum(z_shifts) / len(z_shifts)
print 'From %i measurements:' % len(x_shifts)
print 'Mean shifts: dx = %.2f mm,  dy = %.2f mm;  dz = %.2f mm' \
      % (mean_x,mean_y,mean_z)

# Apply shifts to geometry
if have_geom:

    out = os.path.splitext(geom)[0]+'-predrefine.geom'
    print 'Applying corrections to %s, output filename %s' % (geom,out)
    g = open(geom, 'r')
    h = open(out, 'w')

    prog1 = re.compile("^\s*res\s+=\s+([0-9\.]+)\s")
    prog2 = re.compile("^\s*(.*)\/res\s+=\s+([0-9\.]+)\s")
    prog3 = re.compile("^\s*(.*)\/corner_x\s+=\s+([0-9\.\-]+)\s")
    prog4 = re.compile("^\s*(.*)\/corner_y\s+=\s+([0-9\.\-]+)\s")
    default_res = 0
    while True:

        fline = g.readline()
        if not fline:
            break

        match = prog1.match(fline)
        if match:
            default_res = float(match.group(1))
            print 'default res %f' % (default_res)
            h.write(fline)
            continue

        match = prog2.match(fline)
        if match:
            panel = match.group(1)
            panel_res = float(match.group(2))
            print 'panel res %s / %f' % (panel, panel_res)
            h.write(fline)
            continue

        match = prog3.match(fline)
        if match:
            panel = match.group(1)
            panel_cnx = float(match.group(2))
            res = default_res     # FIXME!
            h.write('%s/corner_x = %f\n' % (panel,panel_cnx+(mean_x*res*1e-3)))
            continue

        match = prog4.match(fline)
        if match:
            panel = match.group(1)
            panel_cny = float(match.group(2))
            res = default_res     # FIXME!
            h.write('%s/corner_y = %f\n' % (panel,panel_cny+(mean_y*res*1e-3)))
            continue

        h.write(fline)

    g.close()
    h.close()

plt.suptitle('Detector shifts according to prediction refinement')
plt.subplot(121, autoscale_on=False, aspect='equal')
plt.plot(x_shifts, y_shifts, 'rx')
plt.plot(0, 0, 'bo')
plt.axis([-2,2,-2,2])
plt.title('In-plane')
plt.xlabel('x shift / mm')
plt.ylabel('y shift / mm')
plt.grid(True)
plt.subplot(122)
plt.hist(z_shifts,bins=30,range=[-2,2])
plt.title('Camera length')
plt.xlabel('z shift / mm')
plt.ylabel('Frequency')
plt.grid(True)
plt.show()