update plot-energy to include ockham factor

1 files changed, 67 insertions, 6 deletions

diff --git a/utils/plot-energy b/utils/plot-energy
index 3c2b587..35e1e89 100755
--- a/utils/plot-energy
+++ b/utils/plot-energy
@@ -37,6 +37,8 @@ SNOMAN_MASS = {
     23: 105.658
 }
 
+AV_RADIUS = 600.0
+
 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
@@ -71,14 +73,45 @@ if __name__ == '__main__':
             for ev in event['ev']:
                 if 'fit' not in ev:
                     continue
-                for id, fit_result in ev['fit'].iteritems():
+                for id, fit_result in [x for x in ev['fit'].iteritems() if isinstance(x[0],int)]:
                     # FIXME: Should I just store the particle ids in the YAML
                     # output as a list of particle ids instead of a single
                     # integer?
                     ids = map(int,chunks(str(id),2))
                     energy = 0.0
+                    skip = False
                     for i, ke in zip(ids,np.atleast_1d(fit_result['energy'])):
                         energy += ke + SNOMAN_MASS[i]
+
+                        # This is a bit of a hack. It appears that many times
+                        # the fit will actually do much better by including a
+                        # very low energy electron or muon. I believe the
+                        # reason for this is that of course my likelihood
+                        # function is not perfect (for example, I don't include
+                        # the correct angular distribution for Rayleigh
+                        # scattered light), and so the fitter often wants to
+                        # add a very low energy electron or muon to fix things.
+                        #
+                        # Ideally I would fix the likelihood function, but for
+                        # now we just discard any fit results which have a very
+                        # low energy electron or muon.
+                        if len(ids) > 1 and i == 20 and ke < 20.0:
+                            skip = True
+
+                        if len(ids) > 1 and i == 22 and ke < 200.0:
+                            skip = True
+
+                    if skip:
+                        continue
+
+                    # Calculate the approximate Ockham factor.
+                    # See Chapter 20 in "Probability Theory: The Logic of Science" by Jaynes
+                    #
+                    # Note: This is a really approximate form by assuming that
+                    # the shape of the likelihood space is equal to the average
+                    # uncertainty in the different parameters.
+                    w = len(ids)*np.log(0.1*0.001) + np.sum(np.log(fit_result['energy'])) + len(ids)*np.log(1e-4/(4*np.pi))
+
                     fit_results.append((
                         ev['run'],
                         ev['gtid'],
@@ -88,7 +121,8 @@ if __name__ == '__main__':
                         fit_result['posz'],
                         fit_result['t0'],
                         energy,
-                        fit_result['fmin']))
+                        fit_result['fmin'] - w,
+                        fit_result['psi']/ev['nhit']))
 
     # create a dataframe
     # note: we have to first create a numpy structured array since there is no
@@ -103,16 +137,43 @@ if __name__ == '__main__':
                             ('z',np.double),     # z
                             ('t0',np.double),    # t0
                             ('ke',np.double),    # kinetic energy
-                            ('fmin',np.double)]  # negative log likelihood
+                            ('fmin',np.double),  # negative log likelihood
+                            ('psi',np.double)]   # goodness of fit parameter
                    )
     df = pd.DataFrame.from_records(array)
 
     # get the best fit
     df = df.sort_values('fmin').groupby(['run','gtid']).first()
 
-    for id, df_id in df.groupby('id'):
-        plt.figure()
-        plot_hist(df_id.ke.values)
+    # require r < 6 meters
+    df = df[np.sqrt(df.x.values**2 + df.y.values**2 + df.z.values**2) < AV_RADIUS]
+
+    # Note: Need to design and apply a psi based cut here, and apply the muon
+    # and neutron follower cuts.
+
+    for id, df_id in sorted(df.groupby('id')):
+        if id == 20:
+            plt.subplot(3,4,1)
+        elif id == 22:
+            plt.subplot(3,4,2)
+        elif id == 2020:
+            plt.subplot(3,4,5)
+        elif id == 2022:
+            plt.subplot(3,4,6)
+        elif id == 2222:
+            plt.subplot(3,4,7)
+        elif id == 202020:
+            plt.subplot(3,4,9)
+        elif id == 202022:
+            plt.subplot(3,4,10)
+        elif id == 202222:
+            plt.subplot(3,4,11)
+        elif id == 222222:
+            plt.subplot(3,4,12)
+
+        plt.hist(df_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step')
         plt.xlabel("Energy (MeV)")
         plt.title(str(id))
+
+    plt.tight_layout()
     plt.show()