diff options
| -rwxr-xr-x | utils/plot-energy | 241 |
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() |