Detector alignment simulation (WIP)

4 files changed, 96 insertions, 27 deletions

diff --git a/julia/CrystFEL/src/CrystFEL.jl b/julia/CrystFEL/src/CrystFEL.jl
index aa1b4518..47f09d41 100644
--- a/julia/CrystFEL/src/CrystFEL.jl
+++ b/julia/CrystFEL/src/CrystFEL.jl
@@ -21,6 +21,7 @@ export TetragonalLattice, HexagonalLattice, RhombohedralLattice, CubicLattice
 export PrimitiveCell, ACenteredCell, BCenteredCell, CCenteredCell
 export BodyCenteredCell, FaceCenteredCell, RhombohedralCell, RhombohedralCellOnHexagonalAxes
 export NoUniqueAxis, UnknownUniqueAxis, UniqueAxisA, UniqueAxisB, UniqueAxisC
+export rotatecell
 
 include("detgeom.jl")
 using .DetGeoms
diff --git a/julia/CrystFEL/src/cell.jl b/julia/CrystFEL/src/cell.jl
index b094770a..80250221 100644
--- a/julia/CrystFEL/src/cell.jl
+++ b/julia/CrystFEL/src/cell.jl
@@ -7,6 +7,7 @@ export TetragonalLattice, HexagonalLattice, RhombohedralLattice, CubicLattice
 export PrimitiveCell, ACenteredCell, BCenteredCell, CCenteredCell
 export BodyCenteredCell, FaceCenteredCell, RhombohedralCell, RhombohedralCellOnHexagonalAxes
 export NoUniqueAxis, UnknownUniqueAxis, UniqueAxisA, UniqueAxisB, UniqueAxisC
+export rotatecell
 
 
 # Represents the real C-side (opaque) structure.
@@ -305,4 +306,21 @@ function Base.show(io::IO, uc::UnitCell)
 end
 
 
+mutable struct CrystFELQuaternion
+    w::Cdouble
+    x::Cdouble
+    y::Cdouble
+    z::Cdouble
+end
+
+function rotatecell(uc, quat)
+    q = CrystFELQuaternion(quat...)
+    out = @ccall libcrystfel.cell_rotate(uc.internalptr::Ptr{InternalUnitCell},
+                                         q::CrystFELQuaternion)::Ptr{InternalUnitCell}
+    if out == C_NULL
+        throw(ErrorException("Failed to rotate unit cell"))
+    end
+    UnitCell(out)
+end
+
 end   # of module
diff --git a/julia/alignment-test-moved.geom b/julia/alignment-test-moved.geom
new file mode 100644
index 00000000..c1943056
--- /dev/null
+++ b/julia/alignment-test-moved.geom
@@ -0,0 +1,49 @@
+adu_per_photon = 1
+res = 10000
+clen = 100.0 mm
+photon_energy = 9000 eV
+
+dim0 = %
+dim1 = ss
+dim2 = fs
+data = /data/data
+
+q0/dim0 = 0
+q0/min_fs = 0
+q0/min_ss = 0
+q0/max_fs = 1024
+q0/max_ss = 1024
+q0/fs = x
+q0/ss = y
+q0/corner_x = -1056
+q0/corner_y = 30
+
+q1/dim0 = 1
+q1/min_fs = 0
+q1/min_ss = 0
+q1/max_fs = 1024
+q1/max_ss = 1024
+q1/fs = x
+q1/ss = y
+q1/corner_x = 30
+q1/corner_y = 30
+
+q2/dim0 = 2
+q2/min_fs = 0
+q2/min_ss = 0
+q2/max_fs = 1024
+q2/max_ss = 1024
+q2/fs = x
+q2/ss = y
+q2/corner_x = 30
+q2/corner_y = -1054
+
+q3/dim0 = 3
+q3/min_fs = 0
+q3/min_ss = 0
+q3/max_fs = 1024
+q3/max_ss = 1024
+q3/fs = x
+q3/ss = y
+q3/corner_x = -1054
+q3/corner_y = -1054
diff --git a/julia/alignment-test.jl b/julia/alignment-test.jl
index 345c7d77..e6e1ab63 100644
--- a/julia/alignment-test.jl
+++ b/julia/alignment-test.jl
@@ -1,44 +1,45 @@
 using CrystFEL
 using Random
 
-# Create empty image for simulation purposes
-dtempl = loaddatatemplate("julia/alignment-test.geom")
-image = Image(dtempl)
-
-# Create a crystal and calculate its reflections
-cell = UnitCell(MonoclinicLattice, PrimitiveCell, 123, 45, 80, 90, 97, 90)
-cr = Crystal(cell)  # FIXME: Random rotation
-truth = predictreflections(cr, image)
-
 # "Simulate" a diffraction pattern from the reflections
-function sketch_pattern(reflist)
-    peaks = PeakList()
+function sketch_pattern(image, cr)
+    reflist = predictreflections(cr, image)
+    peaklist = PeakList()
     for refl in reflist
         if randn() > 0
             let dpos = refl.detectorposition
-                push!(peaks, dpos.fs, dpos.ss, dpos.panelnumber, 100.0)
+                push!(peaklist, dpos.fs, dpos.ss, dpos.panelnumber, 100.0)
             end
         end
     end
-    return peaks
+    return peaklist
 end
 
-image.peaklist = sketch_pattern(truth)
 
-# Index the pattern
-indexer = Indexer("asdf", dtempl, cell, retry=false, multilattice=false, refine=true)
-index(image, indexer)
+function randomrotation(cell)
+    rotatecell(cell, (1,0,0,0))
+end
 
 
-# Utility routine for visualising peaks
-function plotpanel(image, pn)
-    x = []
-    y = []
-    for pk in image.peaklist
-        if pk.panelnumber == pn
-            push!(x, pk.fs)
-            push!(y, pk.ss)
-        end
+function simulate_and_index(cell, image_true, dtempl_moved, n)
+
+    indexer = Indexer("asdf", dtempl_moved, cell, retry=false, multilattice=false, refine=true)
+
+    for _ in 1:n
+
+        cr = Crystal(randomrotation(cell))
+        peaklist = sketch_pattern(image_true, cr)
+        image_moved = Image(dtempl_moved)
+
+        image_moved.peaklist = peaklist
+        index(image_moved, indexer)
+
     end
-    scatter(x, y)
 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-moved.geom")
+simulate_and_index(cell, image_true, dtempl_moved, 100)