Add peakogram/detector-shift

4 files changed, 385 insertions, 4 deletions

diff --git a/crystfel-demo.c b/crystfel-demo.c
index 3bde9d8..b4e1025 100644
--- a/crystfel-demo.c
+++ b/crystfel-demo.c
@@ -228,6 +228,42 @@ gint check_near_bragg(GtkWidget *widget, struct crystfeldemo *demo)
 }
 
 
+gint peakogram(GtkWidget *widget, struct crystfeldemo *demo)
+{
+	GError *error = NULL;
+	GSubprocess *sub;
+
+	sub = g_subprocess_new(G_SUBPROCESS_FLAGS_NONE,
+	                       &error, "sh", "-c",
+	                       "${CRYSTFEL_DEMO_FILES}/peakogram-stream -i "
+	                       "${CRYSTFEL_DEMO_FILES}/cell.stream",
+	                       NULL);
+
+	if ( sub == NULL ) {
+		printf("Failed to start demo process\n");
+	}
+	return 0;
+}
+
+
+gint detector_shift(GtkWidget *widget, struct crystfeldemo *demo)
+{
+	GError *error = NULL;
+	GSubprocess *sub;
+
+	sub = g_subprocess_new(G_SUBPROCESS_FLAGS_NONE,
+	                       &error, "sh", "-c",
+	                       "${CRYSTFEL_DEMO_FILES}/detector-shift "
+	                       "${CRYSTFEL_DEMO_FILES}/cell.stream",
+	                       NULL);
+
+	if ( sub == NULL ) {
+		printf("Failed to start demo process\n");
+	}
+	return 0;
+}
+
+
 static int change_to_tempdir()
 {
 	char tmpdir[64];
@@ -302,6 +338,11 @@ int main(int argc, char *argv[])
 	gtk_builder_add_callback_symbol(builder, "stop_near_bragg",
 	                                G_CALLBACK(stop_near_bragg));
 
+	gtk_builder_add_callback_symbol(builder, "peakogram",
+	                                G_CALLBACK(peakogram));
+	gtk_builder_add_callback_symbol(builder, "detector_shift",
+	                                G_CALLBACK(detector_shift));
+
 	gtk_builder_connect_signals(builder, &demo);
 	gtk_widget_show_all(window);
 	gtk_main();
diff --git a/data/crystfel-demo.glade b/data/crystfel-demo.glade
index f214fd6..942cd52 100644
--- a/data/crystfel-demo.glade
+++ b/data/crystfel-demo.glade
@@ -121,7 +121,7 @@
                       <object class="GtkLabel">
                         <property name="visible">True</property>
                         <property name="can_focus">False</property>
-                        <property name="label" translatable="yes">indexamajig -i files.lst -o output.stream --indexing=none --peaks=zaef</property>
+                        <property name="label" translatable="yes">indexamajig -i files.lst -o 5HT2B-peaks.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8 --indexing=none</property>
                         <property name="wrap">True</property>
                         <property name="selectable">True</property>
                       </object>
@@ -215,7 +215,7 @@
                       <object class="GtkLabel">
                         <property name="visible">True</property>
                         <property name="can_focus">False</property>
-                        <property name="label" translatable="yes">indexamajig -i files.lst -o 5HT2B-all.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8</property>
+                        <property name="label" translatable="yes">indexamajig -i files.lst -o 5HT2B-nocell.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8</property>
                         <property name="wrap">True</property>
                         <property name="selectable">True</property>
                       </object>
@@ -307,7 +307,7 @@
                       <object class="GtkLabel">
                         <property name="visible">True</property>
                         <property name="can_focus">False</property>
-                        <property name="label" translatable="yes">indexamajig -i files.lst -o output.stream --peaks=zaef -p lysozyme.cell</property>
+                        <property name="label" translatable="yes">indexamajig -i files.lst -o 5HT2B-cell.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8 -p 5HT2B.cell</property>
                         <property name="wrap">True</property>
                         <property name="selectable">True</property>
                       </object>
@@ -422,6 +422,7 @@
                             <property name="visible">True</property>
                             <property name="can_focus">True</property>
                             <property name="receives_default">True</property>
+                            <signal name="clicked" handler="peakogram" swapped="no"/>
                           </object>
                           <packing>
                             <property name="expand">True</property>
@@ -435,6 +436,7 @@
                             <property name="visible">True</property>
                             <property name="can_focus">True</property>
                             <property name="receives_default">True</property>
+                            <signal name="clicked" handler="detector_shift" swapped="no"/>
                           </object>
                           <packing>
                             <property name="expand">True</property>
@@ -491,7 +493,7 @@
                       <object class="GtkLabel">
                         <property name="visible">True</property>
                         <property name="can_focus">False</property>
-                        <property name="label" translatable="yes">partialator -i indexing.stream -o merged.hkl --model=unity</property>
+                        <property name="label" translatable="yes">partialator -i 5HT2B-cell.stream -o merged.hkl --model=unity</property>
                         <property name="wrap">True</property>
                         <property name="selectable">True</property>
                       </object>
diff --git a/files/detector-shift b/files/detector-shift
new file mode 100755
index 0000000..ce3bf59
--- /dev/null
+++ b/files/detector-shift
@@ -0,0 +1,163 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Determine mean detector shift based on prediction refinement results
+#
+# Copyright © 2015-2018 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
+#
+# Author:
+#    2015-2018 Thomas White <taw@physics.org>
+#    2016      Mamoru Suzuki <mamoru.suzuki@protein.osaka-u.ac.jp>
+#    2018      Chun Hong Yoon
+#
+
+import sys
+import os
+import re
+import numpy as np
+import matplotlib.pyplot as plt
+
+if sys.argv[1] == "-":
+    f = sys.stdin
+else:
+    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()
+
+mean_x = sum(x_shifts) / len(x_shifts)
+mean_y = sum(y_shifts) / len(y_shifts)
+print('Mean shifts: dx = {:.2} mm,  dy = {:.2} mm'.format(mean_x,mean_y))
+print('Shifts will be applied to geometry file when you close the graph window')
+print('Click anywhere on the graph to override the detector shift')
+
+def plotNewCentre(x, y):
+    circle1 = plt.Circle((x,y),.1,color='r',fill=False)
+    fig.gca().add_artist(circle1)
+    plt.plot(x, y, 'b8', color='m')
+    plt.grid(True)
+
+def onclick(event):
+    print('New shifts: dx = {:.2} mm,  dy = {:.2} mm'.format(event.xdata, event.ydata))
+    print('Shifts will be applied to geometry file when you close the graph window')
+    mean_x = event.xdata
+    mean_y = event.ydata
+    plotNewCentre(mean_x, mean_y)
+
+nbins = 200
+H, xedges, yedges = np.histogram2d(x_shifts,y_shifts,bins=nbins)
+H = np.rot90(H)
+H = np.flipud(H)
+Hmasked = np.ma.masked_where(H==0,H)
+
+# Plot 2D histogram using pcolor
+plt.ion()
+fig2 = plt.figure()
+cid = fig2.canvas.mpl_connect('button_press_event', onclick)
+plt.pcolormesh(xedges,yedges,Hmasked)
+plt.title('Detector shifts according to prediction refinement')
+plt.xlabel('x shift / mm')
+plt.ylabel('y shift / mm')
+plt.plot(0, 0, 'bH', color='c')
+fig = plt.gcf()
+cbar = plt.colorbar()
+cbar.ax.set_ylabel('Counts')
+plotNewCentre(mean_x, mean_y)
+plt.show(block=True)
+
+# Apply shifts to geometry
+if have_geom:
+
+    out = os.path.splitext(geom)[0]+'-predrefine.geom'
+    print('Applying corrections to {}, output filename {}'.format(geom,out))
+    g = open(geom, 'r')
+    h = open(out, 'w')
+    panel_resolutions = {}
+
+    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))
+            h.write(fline)
+            continue
+
+        match = prog2.match(fline)
+        if match:
+            panel = match.group(1)
+            panel_res = float(match.group(2))
+            default_res =  panel_res
+            panel_resolutions[panel] = panel_res
+            h.write(fline)
+            continue
+
+        match = prog3.match(fline)
+        if match:
+            panel = match.group(1)
+            panel_cnx = float(match.group(2))
+            if panel in panel_resolutions:
+                res = panel_resolutions[panel]
+            else:
+                res = default_res
+                print('Using default resolution ({} px/m) for panel {}'.format(res, panel))
+            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))
+            if panel in panel_resolutions:
+                res = panel_resolutions[panel]
+            else:
+                res = default_res
+                print('Using default resolution ({} px/m) for panel {}'.format(res, panel))
+            h.write('%s/corner_y = %f\n' % (panel,panel_cny+(mean_y*res*1e-3)))
+            continue
+
+        h.write(fline)
+
+    g.close()
+    h.close()
+
diff --git a/files/peakogram-stream b/files/peakogram-stream
new file mode 100755
index 0000000..874f6a9
--- /dev/null
+++ b/files/peakogram-stream
@@ -0,0 +1,175 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Check a stream for saturation
+#
+# Copyright © 2016-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
+# Copyright © 2016      The Research Foundation for SUNY
+#
+# Authors:
+#   2016-2017 Thomas White <taw@physics.org>
+#   2014-2016 Thomas Grant <tgrant@hwi.buffalo.edu>
+#
+# This file is part of CrystFEL.
+#
+# CrystFEL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# CrystFEL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+
+import sys
+import argparse
+import math as m
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.colors import LogNorm
+
+def c2(a):
+    return m.cos(a) * m.cos(a)
+
+def s2(a):
+    return m.sin(a) * m.sin(a)
+
+# Return 1/d for hkl in cell, in 1/Angstroms
+def resolution(scell, shkl):
+
+    a = float(scell[0])*10.0
+    b = float(scell[1])*10.0
+    c = float(scell[2])*10.0  # nm -> Angstroms
+
+    al = m.radians(float(scell[3]))
+    be = m.radians(float(scell[4]))
+    ga = m.radians(float(scell[5])) # in degrees
+
+    h = int(shkl[0])
+    k = int(shkl[1])
+    l = int(shkl[2])
+
+    pf = 1.0 - c2(al) - c2(be) - c2(ga) + 2.0*m.cos(al)*m.cos(be)*m.cos(ga)
+    n1 = h*h*s2(al)/(a*a) + k*k*s2(be)/(b*b) + l*l*s2(ga)/(c*c)
+    n2a = 2.0*k*l*(m.cos(be)*m.cos(ga) - m.cos(al))/(b*c)
+    n2b = 2.0*l*h*(m.cos(ga)*m.cos(al) - m.cos(be))/(c*a)
+    n2c = 2.0*h*k*(m.cos(al)*m.cos(be) - m.cos(ga))/(a*b)
+
+    return m.sqrt((n1 + n2a + n2b + n2c) / pf)
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument("-i", default="my.stream", help="stream filename")
+parser.add_argument("-l", action="store_true", help="log scale y-axis")
+parser.add_argument("--rmin", type=float, help="minimum resolution cutoff (1/d in Angstroms^-1)")
+parser.add_argument("--rmax", type=float, help="maximum resolution cutoff (1/d in Angstroms^-1)")
+parser.add_argument("--imin", type=float, help="minimum peak intensity cutoff")
+parser.add_argument("--imax", type=float, help="maximum peak intensity cutoff")
+parser.add_argument("--nmax", default=np.inf, type=int, help="maximum number of peaks to read")
+parser.add_argument("-o", default="peakogram", help="output file prefix")
+args = parser.parse_args()
+
+data = []
+n=0
+in_list = 0
+cell = []
+
+if args.i == "-":
+    f = sys.stdin
+else:
+    f = open(args.i)
+
+if f:
+    for line in f:
+
+        if line.find("Cell parameters") != -1:
+            cell[0:3] = line.split()[2:5]
+            cell[3:6] = line.split()[6:9]
+            continue
+        if line.find("Reflections measured after indexing") != -1:
+            in_list = 1
+            continue
+        if line.find("End of reflections") != -1:
+            in_list = 0
+        if in_list == 1:
+            in_list = 2
+            continue
+        elif in_list != 2:
+            continue
+
+        # From here, we are definitely handling a reflection line
+
+        # Add reflection to list
+        columns = line.split()
+        n += 1
+        try:
+            data.append([resolution(cell, columns[0:3]),columns[5]])
+        except:
+            print("Error with line: "+line.rstrip("\r\n"))
+            print("Cell: "+str(cell))
+
+        if n%1000==0:
+            sys.stdout.write("\r%i predicted reflections found" % n)
+            sys.stdout.flush()
+
+        if n >= args.nmax:
+            break
+
+
+
+data = np.asarray(data,dtype=float)
+
+sys.stdout.write("\r%i predicted reflections found" % n)
+sys.stdout.flush()
+
+print("")
+
+x = data[:,0]
+y = data[:,1]
+
+xmin = np.min(x[x>0])
+xmax = np.max(x)
+ymin = np.min(y[y>0])
+ymax = np.max(y)
+
+if args.rmin is not None:
+    xmin = args.rmin
+if args.rmax is not None:
+    xmax = args.rmax
+if args.imin is not None:
+    ymin = args.imin
+if args.imax is not None:
+    ymax = args.imax
+
+keepers = np.where((x>=xmin) & (x<=xmax) & (y>=ymin) & (y<=ymax))
+
+x = x[keepers]
+y = y[keepers]
+
+if args.l:
+    y = np.log10(y)
+    ymin = np.log10(ymin)
+    ymax = np.log10(ymax)
+
+bins=300
+H,xedges,yedges = np.histogram2d(y,x,bins=bins)
+
+fig = plt.figure()
+ax1 = plt.subplot(111)
+plot = ax1.pcolormesh(yedges,xedges,H, norm=LogNorm())
+cbar = plt.colorbar(plot)
+plt.xlim([xmin,xmax])
+plt.ylim([ymin,ymax])
+plt.xlabel("1/d (A^-1)")
+if args.l:
+    plt.ylabel("Log(Reflection max intensity)")
+else:
+    plt.ylabel("Reflection max intensity")
+plt.title(args.i)
+plt.show()
+