plot Michel electrons and muons in plot-energy

This commit is a major update to the plot-energy script. The most significant
changes are:

- new prompt event selection

    We now define prompt events as any event >= 100 nhit which is at least 250
    ms away from the last 100 nhit event.

- add Michel electron event selection.

    Michel electrons are now selected as any event between 800 ns and 20
    microseconds after a muon or prompt event. Additionally they are required
    to be >= 100 nhit and pass an extra set of data cleaning cuts (compared to
    the prompt events).

- add an atmospheric "sideband"

    I also updated the script to plot the energy distribution and particle ID
    for all events that *do* have a neutron follower as a kind of sideband
    which should contain only atmospheric events.

- plot muon energy spectrum and angular distribution

1 files changed, 211 insertions, 30 deletions

diff --git a/utils/plot-energy b/utils/plot-energy
index 6867983..f62ac0e 100755
--- a/utils/plot-energy
+++ b/utils/plot-energy
@@ -148,6 +148,8 @@ if __name__ == '__main__':
                         fit_result['posz'],
                         fit_result['t0'],
                         energy,
+                        np.atleast_1d(fit_result['theta'])[0],
+                        np.atleast_1d(fit_result['phi'])[0],
                         fit_result['fmin'] - w,
                         fit_result['psi']/ev['nhit']))
 
@@ -164,87 +166,154 @@ if __name__ == '__main__':
                             ('z',np.double),     # z
                             ('t0',np.double),    # t0
                             ('ke',np.double),    # kinetic energy
+                            ('theta',np.double), # direction of 1st particle
+                            ('phi',np.double),   # direction of 2nd particle
                             ('fmin',np.double),  # negative log likelihood
                             ('psi',np.double)]   # goodness of fit parameter
                    )
     df = pd.DataFrame.from_records(array)
 
     array = np.array(events,
-                     dtype=[('run',np.int),      # run number
-                            ('gtr',np.double),   # 50 MHz clock in ns
-                            ('nhit',np.int),     # number of PMTs hit
-                            ('gtid',np.int),     # gtid
-                            ('dc',np.int),       # data cleaning word
-                            ('ftpx',np.double),       # data cleaning word
-                            ('ftpy',np.double),       # data cleaning word
-                            ('ftpz',np.double),       # data cleaning word
-                            ('rsp_energy',np.double)]       # data cleaning word
+                     dtype=[('run',np.int),           # run number
+                            ('gtr',np.double),        # 50 MHz clock in ns
+                            ('nhit',np.int),          # number of PMTs hit
+                            ('gtid',np.int),          # gtid
+                            ('dc',np.int),            # data cleaning word
+                            ('ftpx',np.double),       # FTP fitter X position
+                            ('ftpy',np.double),       # FTP fitter Y position
+                            ('ftpz',np.double),       # FTP fitter Z position
+                            ('rsp_energy',np.double)] # RSP energy
                    )
     df_ev = pd.DataFrame.from_records(array)
 
-    # remove events 200 microseconds after a muon
-    muons = df_ev[(df_ev.dc & DC_MUON) != 0]
+    # Make sure events are in order. We use run number and GTID here which
+    # should be in time order as well except for bit flips in the GTID
+    # This is important for later when we look at time differences in the 50
+    # MHz clock. We should perhaps do a check here based on the 10 MHz clock to
+    # make sure things are in order
+    df_ev = df_ev.sort_values(by=['run','gtid'])
 
     print("number of events = %i" % len(df_ev))
 
-    print("number of muons = %i" % len(muons))
-
-    df_ev = df_ev[(df_ev.dc & DC_MUON) == 0]
-
-    print("number of events after muon cut = %i" % len(df_ev))
-
-    if muons.size:
-        # FIXME: need to deal with 50 MHz clock rollover
-        df_ev = df_ev[~np.any((df_ev.gtr.values > muons.gtr.values[:,np.newaxis]) & (df_ev.gtr.values <= (muons.gtr.values[:,np.newaxis] + 200e3)),axis=0)]
-
-    print("number of events after muon follower cut = %i" % len(df_ev))
-
-    # perform prompt event data cleaning
+    # First, do basic data cleaning which is done for all events
     df_ev = df_ev[df_ev.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK) == 0]
 
     print("number of events after data cleaning = %i" % len(df_ev))
 
+    # Now, we select events tagged by the muon tag which should tag only external muons
+    # We keep the sample of muons since it's needed later to identify Michel
+    # electrons and to apply the muon follower cut
+    muons = df_ev[(df_ev.dc & DC_MUON) != 0]
+
+    print("number of muons = %i" % len(muons))
+
     # select prompt events
     # FIXME: how to deal with two prompt events one after another
     prompt = df_ev[df_ev.nhit >= 100]
 
+    # We define a prompt event here as any event with an NHIT > 100 and whose
+    # previous > 100 nhit event was more than 250 ms ago
+    #
+    # Should I use 10 MHz clock here since the time isn't critical and then I
+    # don't have to worry about 50 MHz clock rollover?
+    prompt_mask = np.concatenate(([True],np.diff(prompt.gtr.values) > 250e6))
+
+    # Michel electrons and neutrons can be any event which is not a prompt
+    # event
+    follower = pd.concat([df_ev[df_ev.nhit < 100],prompt[~prompt_mask]])
+
+    # Apply the prompt mask
+    prompt = prompt[prompt_mask]
+
+    # Try to identify Michel electrons. Currenly, the event selection is based
+    # on Richie's thesis. Here, we do the following:
+    #
+    # 1. Apply more data cleaning cuts to potential Michel electrons
+    # 2. Nhit >= 100
+    # 3. It must be > 800 ns and less than 20 microseconds from a prompt event
+    #    or a muon
+    prompt_plus_muons = pd.concat([prompt,muons])
+
     print("number of events after prompt nhit cut = %i" % len(prompt))
 
-    if prompt.size:
+    if prompt_plus_muons.size and follower.size:
+        # require Michel events to pass more of the SNO data cleaning cuts
+        michel = follower[follower.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK | DC_ESUM | DC_OWL | DC_OWL_TRIGGER | DC_FTS) == DC_FTS]
+        michel = michel[michel.nhit >= 100]
+        # accept events which had a muon more than 800 nanoseconds but less than 20 microseconds before them
+        # The 800 nanoseconds cut comes from Richie's thesis. He also mentions
+        # that the In Time Channel Spread Cut is very effective at cutting
+        # electrons events caused by muons, so I should implement that
+        michel = michel[np.any((michel.gtr.values > prompt_plus_muons.gtr.values[:,np.newaxis] + 800) & (michel.gtr.values < prompt_plus_muons.gtr.values[:,np.newaxis] + 20e3),axis=0)]
+    else:
+        michel = prompt_plus_muons.copy()
+
+    if prompt.size and follower.size:
         # FIXME: need to deal with 50 MHz clock rollover
         # neutron followers have to obey stricter set of data cleaning cuts
-        neutron = df_ev[df_ev.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK | DC_ESUM | DC_OWL | DC_OWL_TRIGGER | DC_FTS) == DC_FTS]
+        neutron = follower[follower.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK | DC_ESUM | DC_OWL | DC_OWL_TRIGGER | DC_FTS) == DC_FTS]
         neutron = neutron[~np.isnan(neutron.ftpx) & ~np.isnan(neutron.rsp_energy)]
         r = np.sqrt(neutron.ftpx**2 + neutron.ftpy**2 + neutron.ftpz**2)
         neutron = neutron[r < AV_RADIUS]
         neutron = neutron[neutron.rsp_energy > 4.0]
         # neutron events accepted after 20 microseconds and before 250 ms (50 ms during salt)
+        df_atm = prompt[np.any((neutron.gtr.values > prompt.gtr.values[:,np.newaxis] + 20e3) & (neutron.gtr.values < prompt.gtr.values[:,np.newaxis] + 250e6),axis=1)]
         df_ev = prompt[~np.any((neutron.gtr.values > prompt.gtr.values[:,np.newaxis] + 20e3) & (neutron.gtr.values < prompt.gtr.values[:,np.newaxis] + 250e6),axis=1)]
     else:
         df_ev = prompt
+        df_atm = prompt
 
     print("number of events after neutron follower cut = %i" % len(df_ev))
 
-    df = pd.merge(df,df_ev,how='inner',on=['run','gtid'])
+    # Get rid of muon events in our main event sample
+    prompt = prompt[(prompt.dc & DC_MUON) == 0]
+
+    print("number of events after muon cut = %i" % len(df_ev))
+
+    if muons.size:
+        # FIXME: need to deal with 50 MHz clock rollover
+        # remove events 200 microseconds after a muon
+        prompt = prompt[~np.any((prompt.gtr.values > muons.gtr.values[:,np.newaxis]) & (prompt.gtr.values <= (muons.gtr.values[:,np.newaxis] + 200e3)),axis=0)]
+
+    print("number of events after muon follower cut = %i" % len(df_ev))
+
+    df_atm = pd.merge(df,df_atm,how='inner',on=['run','gtid'])
+    muons = pd.merge(df,muons,how='inner',on=['run','gtid'])
+    michel = pd.merge(df,michel,how='inner',on=['run','gtid'])
+    df = pd.merge(df,prompt,how='inner',on=['run','gtid'])
 
     # get rid of events which don't have a fit
     nan = np.isnan(df.fmin.values)
     df = df[~nan]
 
+    nan_atm = np.isnan(df_atm.fmin.values)
+    df_atm = df_atm[~nan_atm]
+
+    nan_muon = np.isnan(muons.fmin.values)
+    muons = muons[~nan_muon]
+
+    nan_michel = np.isnan(michel.fmin.values)
+    michel = michel[~nan_michel]
+
     if np.count_nonzero(nan):
         print("skipping %i events because they are missing fit information!" % np.count_nonzero(nan),file=sys.stderr)
 
     # get the best fit
     df = df.sort_values('fmin').groupby(['run','gtid']).first()
+    df_atm = df_atm.sort_values('fmin').groupby(['run','gtid']).first()
+    michel_best_fit = michel.sort_values('fmin').groupby(['run','gtid']).first()
+    muon_best_fit = muons.sort_values('fmin').groupby(['run','gtid']).first()
+    muons = muons[muons.id == 22]
 
     # require r < 6 meters
     df = df[np.sqrt(df.x.values**2 + df.y.values**2 + df.z.values**2) < AV_RADIUS]
+    df_atm = df_atm[np.sqrt(df_atm.x.values**2 + df_atm.y.values**2 + df_atm.z.values**2) < AV_RADIUS]
 
     print("number of events after radius cut = %i" % len(df))
 
-    # Note: Need to design and apply a psi based cut here, and apply the muon
-    # and neutron follower cuts.
+    # Note: Need to design and apply a psi based cut here
 
+    plt.figure()
     for id, df_id in sorted(df.groupby('id')):
         if id == 20:
             plt.subplot(3,4,1)
@@ -269,5 +338,117 @@ if __name__ == '__main__':
         plt.xlabel("Energy (MeV)")
         plt.title(str(id))
 
-    plt.tight_layout()
+    plt.suptitle("Without Neutron Follower")
+
+    if len(df):
+        plt.tight_layout()
+
+    plt.figure()
+    for id, df_atm_id in sorted(df_atm.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_atm_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step')
+        plt.xlabel("Energy (MeV)")
+        plt.title(str(id))
+
+    plt.suptitle("With Neutron Follower")
+
+    if len(df_atm):
+        plt.tight_layout()
+
+    plt.figure()
+    for id, michel_id in sorted(michel_best_fit.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(michel_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step')
+        plt.xlabel("Energy (MeV)")
+        plt.title(str(id))
+
+    plt.suptitle("Michel Electrons")
+
+    if len(michel):
+        plt.tight_layout()
+
+    plt.figure()
+    for id, muon_id in sorted(muon_best_fit.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(muon_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step')
+        plt.xlabel("Energy (MeV)")
+        plt.title(str(id))
+
+    plt.suptitle("External Muons")
+
+    if len(muon_best_fit):
+        plt.tight_layout()
+
+    plt.figure()
+    plt.subplot(2,1,1)
+    plt.hist(muons.ke.values, bins=np.logspace(3,7,100), histtype='step')
+    plt.xlabel("Energy (MeV)")
+    plt.gca().set_xscale("log")
+    plt.subplot(2,1,2)
+    plt.hist(np.cos(muons.theta.values), bins=np.linspace(-1,1,100), histtype='step')
+    plt.xlabel(r"$\cos(\theta)$")
+    plt.suptitle("Muons")
+
+    michel = michel[michel.id == 20]
+
+    plt.figure()
+    plt.hist(michel.ke.values, bins=np.linspace(20,100,100), histtype='step', label="My fitter")
+    plt.hist(michel[~np.isnan(michel.rsp_energy.values)].rsp_energy.values, bins=np.linspace(20,100,100), histtype='step',label="RSP")
+    plt.xlabel("Energy (MeV)")
+    plt.title("Michel Electrons")
+    plt.legend()
     plt.show()