diff options
| -rwxr-xr-x | utils/sddm/plot_energy.py | 31 | ||||
| -rwxr-xr-x | utils/submit-grid-jobs | 60 |
2 files changed, 33 insertions, 58 deletions
diff --git a/utils/sddm/plot_energy.py b/utils/sddm/plot_energy.py index 1c1b7e3..b4882d9 100755 --- a/utils/sddm/plot_energy.py +++ b/utils/sddm/plot_energy.py @@ -141,6 +141,20 @@ def michel_cut(ev): michel['muon_gtid'] = -1 return michel +def tag_michels(ev): + """ + Tags michel events with the michel column and also adds a column linking to + the muon GTID. The muon also gets labelled with the stopping_muon column. + """ + michels = michel_cut(ev) + ev['michel'] = np.zeros(len(ev),dtype=np.bool) + ev.loc[michels.index,'michel'] = True + ev['muon_gtid'] = -np.ones(len(ev),dtype=np.int) + ev.loc[michels.index,'muon_gtid'] = michels['muon_gtid'] + ev['stopping_muon'] = np.zeros(len(ev),dtype=np.bool) + ev.loc[ev.gtid.isin(ev.muon_gtid[ev.muon_gtid > 0].values),'stopping_muon'] = 1 + return ev + def atmospheric_events(ev): """ Tags atmospheric events which have a neutron follower. @@ -452,6 +466,12 @@ def get_events(filenames, merge_fits=False): fits = pd.concat([pd.read_hdf(filename, "fits") for filename in filenames],ignore_index=True) rhdr = pd.concat([pd.read_hdf(filename, "rhdr") for filename in filenames],ignore_index=True) + if len(ev) == 0: + if merge_fits: + return ev + else: + return ev, fits + first_gtid = rhdr.set_index('run').to_dict()['first_gtid'] # First, remove junk events since orphans won't have a 50 MHz clock and so @@ -566,6 +586,17 @@ def get_events(filenames, merge_fits=False): # retrigger cut ev = ev.groupby('run',group_keys=False).apply(retrigger_cut) + # Label instrumentals + ev['noise'] = ev.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_ITC | DC_ESUM) != 0 + ev['neck'] = ((ev.dc & DC_NECK) != 0) & ~ev.noise + ev['flasher'] = ((ev.dc & DC_FLASHER) != 0) & ~(ev.noise | ev.neck) & (ev.nhit < 1000) + ev['breakdown'] = ((ev.dc & (DC_FLASHER | DC_BREAKDOWN)) != 0) & ~(ev.noise | ev.neck) & (ev.nhit >= 1000) + ev['muon'] = ((ev.dc & DC_MUON) != 0) & ~(ev.noise | ev.neck | ev.flasher | ev.breakdown) + ev['signal'] = ~(ev.noise | ev.neck | ev.flasher | ev.breakdown | ev.muon) + ev['instrumental'] = ~ev.signal + + ev = ev.groupby('run',group_keys=False).apply(tag_michels) + if merge_fits: ev = pd.merge(fits,ev,how='inner',on=['run','gtid']) # Set the index to (run, gtid) so we can set columns from the single particle results diff --git a/utils/submit-grid-jobs b/utils/submit-grid-jobs index f73a57b..c5f72b0 100755 --- a/utils/submit-grid-jobs +++ b/utils/submit-grid-jobs @@ -444,6 +444,7 @@ if __name__ == '__main__': import pandas as pd from sddm.dc import DC_MUON, DC_JUNK, DC_CRATE_ISOTROPY, DC_QVNHIT, DC_NECK, DC_FLASHER, DC_ESUM, DC_OWL, DC_OWL_TRIGGER, DC_FTS, DC_ITC, DC_BREAKDOWN from sddm import read_hdf + from sddm.plot_energy import get_events parser = argparse.ArgumentParser("submit grid jobs", formatter_class=argparse.ArgumentDefaultsHelpFormatter) parser.add_argument("filenames", nargs='*', help="input files") @@ -556,8 +557,7 @@ if __name__ == '__main__': else: check_call([zdab_cat,filename,"-o",output.name],stderr=f) - ev = read_hdf(output.name, "ev") - rhdr = read_hdf(output.name, "rhdr") + ev, fits = get_events([output.name], merge_fits=False) if len(ev) == 0: continue @@ -567,62 +567,6 @@ if __name__ == '__main__': log.notice("Skipping %s because it's already in the database" % tail) continue - first_gtid = rhdr.set_index('run').to_dict()['first_gtid'] - - # First, remove junk events since orphans won't have a 50 MHz clock and so - # could screw up the 50 MHz clock unwrapping - ev = ev[ev.dc & DC_JUNK == 0] - - # We need the events to be in time order here in order to calculate the - # delta t between events. It's not obvious exactly how to do this. You - # could sort by GTID, but that wraps around. Similarly we can't sort by the - # 50 MHz clock because it also wraps around. Finally, I'm hesitant to sort - # by the 10 MHz clock since it can be unreliable. - # - # Update: Phil proposed a clever way to get the events in order using the - # GTID: - # - # > The GTID rollover should be easy to handle because there should never - # > be two identical GTID's in a run. So if you order the events by GTID, - # > you can assume that events with GTID's that come logically before the - # > first GTID in the run must have occurred after the other events. - # - # Therefore, we can just add 0x1000000 to all GTIDs before the first GTID - # in the event and sort on that. We get the first GTID from the RHDR bank. - ev['gtid_sort'] = ev['gtid'].copy() - - ev = ev.groupby('run',as_index=False).apply(gtid_sort,first_gtid=first_gtid).reset_index(level=0,drop=True) - - ev = ev.sort_values(by=['run','gtid_sort'],kind='mergesort') - - for run, ev_run in ev.groupby('run'): - # Warn about 50 MHz clock jumps since they could indicate that the - # events aren't in order. - dt = np.diff(ev_run.gtr) - if np.count_nonzero((np.abs(dt) > 1e9) & (dt > -0x7ffffffffff*20.0/2)): - log.warn("Warning: %i 50 MHz clock jumps in run %i. Are the events in order?" % \ - (np.count_nonzero((np.abs(dt) > 1e9) & (dt > -0x7ffffffffff*20.0/2)),run)) - - # unwrap the 50 MHz clock within each run - ev.gtr = ev.groupby(['run'],group_keys=False)['gtr'].transform(unwrap_50_mhz_clock) - - for run, ev_run in ev.groupby('run'): - # Warn about GTID jumps since we could be missing a potential flasher - # and/or breakdown, and we need all the events in order to do a - # retrigger cut - if np.count_nonzero(np.diff(ev_run.gtid) != 1): - log.warn("Warning: %i GTID jumps in run %i" % (np.count_nonzero(np.diff(ev_run.gtid) != 1),run)) - - # calculate the time difference between each event and the previous event - # so we can tag retrigger events - ev['dt'] = ev.groupby(['run'],group_keys=False)['gtr'].transform(lambda x: np.concatenate(([1e9],np.diff(x.values)))) - - # Now, select prompt events. - ev = ev.groupby('run',group_keys=False).apply(prompt_event) - - ev['instrumental'] = False - ev.loc[ev.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK | DC_ITC | DC_BREAKDOWN | DC_MUON) != 0,'instrumental'] = True - nevents = 0 njobs = 0 for index, event in ev.iterrows(): |