diff options
| -rwxr-xr-x | utils/plot-energy | 64 |
1 files changed, 64 insertions, 0 deletions
diff --git a/utils/plot-energy b/utils/plot-energy index 35e1e89..292ae7b 100755 --- a/utils/plot-energy +++ b/utils/plot-energy @@ -39,6 +39,18 @@ SNOMAN_MASS = { AV_RADIUS = 600.0 +# Data cleaning bitmasks. +DC_MUON = 0x1 +DC_JUNK = 0x2 +DC_CRATE_ISOTROPY = 0x4 +DC_QVNHIT = 0x8 +DC_NECK = 0x10 +DC_FLASHER = 0x20 +DC_ESUM = 0x40 +DC_OWL = 0x80 +DC_OWL_TRIGGER = 0x100 +DC_FTS = 0x200 + 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 @@ -57,6 +69,7 @@ if __name__ == '__main__': import matplotlib.pyplot as plt import numpy as np import pandas as pd + import sys parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") @@ -72,6 +85,20 @@ if __name__ == '__main__': for i, event in enumerate(data['data']): for ev in event['ev']: if 'fit' not in ev: + fit_results.append(( + ev['run'], + ev['gtr'], + ev['nhit'], + ev['gtid'], + ev['dc'], + 0, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan)) continue 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 @@ -114,7 +141,10 @@ if __name__ == '__main__': fit_results.append(( ev['run'], + ev['gtr'], + ev['nhit'], ev['gtid'], + ev['dc'], id, fit_result['posx'], fit_result['posy'], @@ -130,7 +160,10 @@ if __name__ == '__main__': # https://github.com/pandas-dev/pandas/issues/4464 array = np.array(fit_results, 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 ('id',np.int), # particle id ('x', np.double), # x ('y',np.double), # y @@ -142,6 +175,37 @@ if __name__ == '__main__': ) df = pd.DataFrame.from_records(array) + # remove events 200 microseconds after a muon + muons = df[(df.dc & DC_MUON) != 0] + + print(len(df)) + print("nmuons = %i" % len(muons)) + + df = df[(df.dc & DC_MUON) == 0] + + if muons.size: + # FIXME: need to deal with 50 MHz clock rollover + df = df[~np.any((df.gtr.values > muons.gtr.values[:,np.newaxis]) & (df.gtr.values <= (muons.gtr.values[:,np.newaxis] + 200e3)),axis=0)] + + print(len(df)) + + # perform prompt event data cleaning + df = df[df.dc & (DC_JUNK | DC_CRATE_ISOTROPY | DC_QVNHIT | DC_FLASHER | DC_NECK) == 0] + + print(len(df)) + + # apply prompt event selection + df = df[df.nhit >= 100] + + print(len(df)) + + # get rid of events which don't have a fit + nan = np.isnan(df.fmin.values) + df = df[~nan] + + 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() |