update plot.py to plot the angular resolution

1 files changed, 18 insertions, 0 deletions

diff --git a/utils/plot.py b/utils/plot.py
index 44728ea..e786d0e 100755
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -51,6 +51,7 @@ if __name__ == '__main__':
         dy = []
         dz = []
         dT = []
+        thetas = []
         for event in data['data']:
             # get the particle ID
             id = event['mctk'][-1]['id']
@@ -72,6 +73,16 @@ if __name__ == '__main__':
             true_posz = event['mcvx'][0]['posz']
             posz = event['ev'][0]['fit'][0]['posz']
             dz.append(posz-true_posz)
+            dirx = event['mctk'][-1]['dirx']
+            diry = event['mctk'][-1]['diry']
+            dirz = event['mctk'][-1]['dirz']
+            true_dir = [dirx,diry,dirz]
+            true_dir = np.array(true_dir)/np.linalg.norm(true_dir)
+            theta = event['ev'][0]['fit'][0]['theta']
+            phi = event['ev'][0]['fit'][0]['phi']
+            dir = [np.sin(theta)*np.cos(phi),np.sin(theta)*np.sin(phi),np.cos(theta)]
+            dir = np.array(dir)/np.linalg.norm(dir)
+            thetas.append(np.degrees(np.arccos(np.dot(true_dir,dir))))
 
         mean, mean_error, std, std_error = get_stats(dT)
         print("dT      = %.2g +/- %.2g" % (mean, mean_error))
@@ -85,6 +96,8 @@ if __name__ == '__main__':
         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))
+        mean, mean_error, std, std_error = get_stats(thetas)
+        print("std(theta) = %4.2g +/- %.2g" % (std, std_error))
 
         plt.figure(1)
         plot_hist(dT, label=filename)
@@ -98,6 +111,9 @@ if __name__ == '__main__':
         plt.figure(4)
         plot_hist(dz, label=filename)
         plt.xlabel("Z Position difference (cm)")
+        plt.figure(5)
+        plot_hist(thetas, label=filename)
+        plt.xlabel(r"$\theta$ (deg)")
 
     plt.figure(1)
     plt.legend()
@@ -107,4 +123,6 @@ if __name__ == '__main__':
     plt.legend()
     plt.figure(4)
     plt.legend()
+    plt.figure(5)
+    plt.legend()
     plt.show()