indexamajig: Estimate resolution based on peaks only

1 files changed, 46 insertions, 0 deletions

diff --git a/scripts/ave-peak-resolution b/scripts/ave-peak-resolution
new file mode 100755
index 00000000..468cffd9
--- /dev/null
+++ b/scripts/ave-peak-resolution
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Find mean diffracting resolution
+#
+# Copyright © 2014-2017 Deutsches Elektronen-Synchrotron DESY,
+#                       a research centre of the Helmholtz Association.
+#
+# Author:
+#    2014-2017 Thomas White <taw@physics.org>
+#
+
+import sys
+import numpy
+import matplotlib.pyplot as plt
+
+f = open(sys.argv[1])
+a = []
+pulseId = 0
+
+while True:
+    fline = f.readline()
+    if not fline:
+        break
+    if fline.find("hdf5/XFEL/pulseId") != -1:
+        pulseId = int(fline.split('= ')[1])
+    if fline.find("peak_resolution") != -1:
+        res = float(fline.split('= ')[1].split(' ')[0].rstrip("\r\n"))
+        if pulseId != 0 and pulseId != 488 and pulseId != 496:
+            a.append(res)
+        continue
+
+f.close()
+
+b = numpy.array(a)
+print(" Mean: {:.2} nm^-1 = {:.2} A".format(numpy.mean(b),10.0/numpy.mean(b)))
+print(" Best: {:.2} nm^-1 = {:.2} A".format(numpy.max(b),10.0/numpy.max(b)))
+print("Worst: {:.2} nm^-1 = {:.2} A".format(numpy.min(b),10.0/numpy.min(b)))
+print("Std deviation: {:.2} nm^-1".format(numpy.std(b)))
+
+plt.hist(a,bins=30)
+plt.title('Resolution based on peak search')
+plt.xlabel('Resolution / nm^-1')
+plt.ylabel('Frequency')
+plt.grid(True)
+plt.show()