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
|
using CrystFEL
using Random
using Plots
using MillepedeII
# "Simulate" a diffraction pattern
function sketch_pattern(image, cr)
reflist = predictreflections(cr, image)
peaklist = PeakList()
for refl in reflist
if randn() > 0
let dpos = refl.detectorposition
push!(peaklist, round(dpos.fs), round(dpos.ss), dpos.panelnumber, 100.0)
end
end
end
return peaklist
end
function simulate_and_index(cell, image_true, dtempl_moved, mille, n)
indexer = Indexer("asdf", dtempl_moved, cell, retry=false, multilattice=false, refine=true)
nidx = 0
for i in 1:n
# Create a diffraction pattern for a random orientation
cr = Crystal(rotatecell(cell))
peaklist = sketch_pattern(image_true, cr)
# Make an image with the correct spot positions,
# but with an incorrect geometry
image_moved = Image(dtempl_moved)
image_moved.peaklist = peaklist
# Index the pattern (and store Mille data),
# based on the incorrect geometry
index(image_moved, indexer, mille=mille)
if image_moved.n_crystals > 0
nidx += 1
end
if i % 100 == 0
print("*")
else
print(".")
end
end
println("")
println("Indexed ", nidx, " out of ", n, " frames")
end
dtempl_true = loaddatatemplate("julia/alignment-test.geom")
image_true = Image(dtempl_true)
cell = UnitCell(MonoclinicLattice, PrimitiveCell, 123, 45, 80, 90, 97, 90)
dtempl_moved = loaddatatemplate("julia/alignment-test.geom")
translategroup!(dtempl_moved, "q1", 200e-6, 0, 0)
let mille = Mille("mille.dat")
simulate_and_index(cell, image_true, dtempl_moved, mille, 100)
close(mille)
end
function ploth(a, n, offs, label)
histogram(map(x->x.residual/x.esd,
filter(x->in(n, keys(x.globalgradients)),
a[offs:3:end])),
label=label)
end
function plotfs(filename)
t = loadmille(filename)
l = @layout([q0 q1; q2 q3])
a = collect(Iterators.flatten(t))
q0 = ploth(a, 101, 1, "q0")
q1 = ploth(a, 201, 1, "q1")
q2 = ploth(a, 301, 1, "q2")
q3 = ploth(a, 401, 1, "q3")
plot(q0, q1, q2, q3, layout=l, plot_title="Fast scan residual")
end
function plotss(filename)
t = loadmille(filename)
l = @layout([q0 q1; q2 q3])
a = collect(Iterators.flatten(t))
q0 = ploth(a, 102, 2, "q0")
q1 = ploth(a, 202, 2, "q1")
q2 = ploth(a, 302, 2, "q2")
q3 = ploth(a, 402, 2, "q3")
plot(q0, q1, q2, q3, layout=l, plot_title="Slow scan residual")
end
function plotr(filename)
t = loadmille(filename)
l = @layout([q0 q1; q2 q3])
a = collect(Iterators.flatten(t))
meas = histogram(map(x->x.residual, a[3:3:end]))
plot(meas, plot_title="Excitation error residual")
end
|
