update plot.py to plot the log likelihood ratio

1 files changed, 13 insertions, 1 deletions

diff --git a/utils/plot.py b/utils/plot.py
index a7dfa6d..1f2019b 100755
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -10,6 +10,9 @@ from scipy.stats import iqr
 import matplotlib
 matplotlib.use("Qt5Agg")
 
+IDP_E_MINUS  =    20
+IDP_MU_MINUS =    22
+
 SNOMAN_MASS = {
     20: 0.511,
     21: 0.511,
@@ -21,7 +24,7 @@ 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)
+    n = max(int((np.max(x)-np.min(x))/h),10)
     bins = np.linspace(np.min(x),np.max(x),n)
     plt.hist(x, bins=bins, histtype='step', label=label)
 
@@ -58,6 +61,7 @@ if __name__ == '__main__':
         dz = []
         dT = []
         thetas = []
+        likelihood_ratio = []
         for event in data['data']:
             # get the particle ID
             id = event['mctk'][-1]['id']
@@ -89,6 +93,9 @@ if __name__ == '__main__':
             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))))
+            fmin_electron = event['ev'][0]['fit'][IDP_E_MINUS]['fmin']
+            fmin_muon = event['ev'][0]['fit'][IDP_MU_MINUS]['fmin']
+            likelihood_ratio.append(fmin_muon-fmin_electron)
 
         mean, mean_error, std, std_error = get_stats(dT)
         print("dT      = %.2g +/- %.2g" % (mean, mean_error))
@@ -120,6 +127,9 @@ if __name__ == '__main__':
         plt.figure(5)
         plot_hist(thetas, label=filename)
         plt.xlabel(r"$\theta$ (deg)")
+        plt.figure(6)
+        plot_hist(likelihood_ratio, label=filename)
+        plt.xlabel(r"Log Likelihood Ratio ($e/\mu$)")
 
     plt.figure(1)
     plt.legend()
@@ -131,4 +141,6 @@ if __name__ == '__main__':
     plt.legend()
     plt.figure(5)
     plt.legend()
+    plt.figure(6)
+    plt.legend()
     plt.show()