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
|
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, dpos.fs, 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)
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 i % 100 == 0
print("*")
else
print(".")
end
end
println("")
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 plotresiduals(filename)
t = loadmille(filename)
l = @layout([q0 q1; q2 q3])
a = collect(Iterators.flatten(t))
function ploth(n, offs, label)
histogram(map(x->x.residual,filter(x->in(n, keys(x.globalgradients)), a[offs:3:end])), label=label)
end
q0 = ploth(101, 1, "q0")
q1 = ploth(201, 1, "q1")
q2 = ploth(301, 1, "q2")
q3 = ploth(401, 1, "q3")
plot(q0, q1, q2, q3, layout=l, plot_title="Fast scan residual")
#q0 = ploth(102, 2, "q0")
#q1 = ploth(202, 2, "q1")
#q2 = ploth(302, 2, "q2")
#q3 = ploth(402, 2, "q3")
#plot(q0, q1, q2, q3, layout=l, plot_title="Slow scan residual", reuse=false)
end
|