update plot.py

- determine bin width using the Freedman Diaconis rule
- print out error on the standard deviation

1 files changed, 42 insertions, 10 deletions

diff --git a/utils/plot.py b/utils/plot.py
index 6a2fe73..44728ea 100755
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -1,6 +1,8 @@
 #!/usr/bin/env python
 from __future__ import print_function, division
 import yaml
+import numpy as np
+from scipy.stats import iqr
 
 SNOMAN_MASS = {
     20: 0.511,
@@ -9,6 +11,28 @@ SNOMAN_MASS = {
     23: 105.658
 }
 
+def plot_hist(x, label=None):
+    # determine the bin width using the Freedman Diaconis rule
+    # see https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule
+    h = 2*iqr(x)/len(x)**(1/3)
+    n = int((np.max(x)-np.min(x))/h)
+    bins = np.linspace(np.min(x),np.max(x),n)
+    plt.hist(x, bins=bins, histtype='step', label=label)
+
+def get_stats(x):
+    """
+    Returns a tuple (mean, error mean, std, error std) for the values in x.
+
+    The formula for the standard error on the standard deviation comes from
+    https://stats.stackexchange.com/questions/156518.
+    """
+    mean = np.mean(x)
+    std = np.std(x)
+    n = len(x)
+    u4 = np.mean((x-mean)**4)
+    error = np.sqrt((u4-(n-3)*std**4/(n-1))/n)/(2*std)
+    return mean, std/np.sqrt(n), std, error
+
 if __name__ == '__main__':
     import argparse
     import matplotlib.pyplot as plt
@@ -26,7 +50,7 @@ if __name__ == '__main__':
         dx = []
         dy = []
         dz = []
-        e = []
+        dT = []
         for event in data['data']:
             # get the particle ID
             id = event['mctk'][-1]['id']
@@ -38,7 +62,7 @@ if __name__ == '__main__':
             # so we need to convert it into kinetic energy
             ke = true_energy - mass
             energy = event['ev'][0]['fit'][0]['energy']
-            e.append(energy-ke)
+            dT.append(energy-ke)
             true_posx = event['mcvx'][0]['posx']
             posx = event['ev'][0]['fit'][0]['posx']
             dx.append(posx-true_posx)
@@ -49,22 +73,30 @@ if __name__ == '__main__':
             posz = event['ev'][0]['fit'][0]['posz']
             dz.append(posz-true_posz)
 
-        print("dT = %.2g +/- %.2g" % (np.mean(e), np.std(e)))
-        print("dx = %4.2g +/- %.2g" % (np.mean(dx), np.std(dx)))
-        print("dy = %4.2g +/- %.2g" % (np.mean(dy), np.std(dy)))
-        print("dz = %4.2g +/- %.2g" % (np.mean(dz), np.std(dz)))
+        mean, mean_error, std, std_error = get_stats(dT)
+        print("dT      = %.2g +/- %.2g" % (mean, mean_error))
+        print("std(dT) = %.2g +/- %.2g" % (std, std_error))
+        mean, mean_error, std, std_error = get_stats(dx)
+        print("dx      = %4.2g +/- %.2g" % (mean, mean_error))
+        print("std(dx) = %4.2g +/- %.2g" % (std, std_error))
+        mean, mean_error, std, std_error = get_stats(dy)
+        print("dy      = %4.2g +/- %.2g" % (mean, mean_error))
+        print("std(dy) = %4.2g +/- %.2g" % (std, std_error))
+        mean, mean_error, std, std_error = get_stats(dz)
+        print("dz      = %4.2g +/- %.2g" % (mean, mean_error))
+        print("std(dz) = %4.2g +/- %.2g" % (std, std_error))
 
         plt.figure(1)
-        plt.hist(e, histtype='step', label=filename)
+        plot_hist(dT, label=filename)
         plt.xlabel("Kinetic Energy difference (MeV)")
         plt.figure(2)
-        plt.hist(dx, histtype='step', label=filename)
+        plot_hist(dx, label=filename)
         plt.xlabel("X Position difference (cm)")
         plt.figure(3)
-        plt.hist(dy, histtype='step', label=filename)
+        plot_hist(dy, label=filename)
         plt.xlabel("Y Position difference (cm)")
         plt.figure(4)
-        plt.hist(dz, histtype='step', label=filename)
+        plot_hist(dz, label=filename)
         plt.xlabel("Z Position difference (cm)")
 
     plt.figure(1)