From e64ddcb421ffe44b276e7f68cc35f25ffdb8e25e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 25 Sep 2019 16:59:27 +0200 Subject: Add peakogram/detector-shift --- crystfel-demo.c | 41 +++++++++++ data/crystfel-demo.glade | 10 +-- files/detector-shift | 163 +++++++++++++++++++++++++++++++++++++++++++ files/peakogram-stream | 175 +++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 385 insertions(+), 4 deletions(-) create mode 100755 files/detector-shift create mode 100755 files/peakogram-stream 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 @@ True False - indexamajig -i files.lst -o output.stream --indexing=none --peaks=zaef + indexamajig -i files.lst -o 5HT2B-peaks.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8 --indexing=none True True @@ -215,7 +215,7 @@ True False - indexamajig -i files.lst -o 5HT2B-all.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8 + indexamajig -i files.lst -o 5HT2B-nocell.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8 True True @@ -307,7 +307,7 @@ True False - indexamajig -i files.lst -o output.stream --peaks=zaef -p lysozyme.cell + indexamajig -i files.lst -o 5HT2B-cell.stream -j `nproc` --peaks=peakfinder8 --threshold=400 --min-pix-count=2 --min-snr=8 -p 5HT2B.cell True True @@ -422,6 +422,7 @@ True True True + True @@ -435,6 +436,7 @@ True True True + True @@ -491,7 +493,7 @@ True False - partialator -i indexing.stream -o merged.hkl --model=unity + partialator -i 5HT2B-cell.stream -o merged.hkl --model=unity True True 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 +# 2016 Mamoru Suzuki +# 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 +# 2014-2016 Thomas Grant +# +# 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 . + +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() + -- cgit v1.2.3