aboutsummaryrefslogtreecommitdiff
path: root/utils/convert-ratdb
blob: 5dd6ee42ec3656ccaaf1567dd3d6074194957f68 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#!/usr/bin/env python
"""
Script to convert the PMT database in RAT to a CSV file (with spaces instead of commas).
"""
from __future__ import print_function, division
import yaml

if __name__ == '__main__':
    import argparse
    import numpy as np

    parser = argparse.ArgumentParser("convert a RATDB file -> CSV file")
    parser.add_argument("filename", help="RATDB filename")
    parser.add_argument("-o", "--output", help="output filename", required=True)
    args = parser.parse_args()

    s = ''
    with open(args.filename) as f:
        for line in f:
            if line.startswith('//'):
                continue
            s += line

    data = yaml.load(s)

    csv = np.zeros((len(data['x']),7),dtype=np.float64)
    csv[:,0] = data['x']
    csv[:,1] = data['y']
    csv[:,2] = data['z']
    csv[:,3] = data['u']
    csv[:,4] = data['v']
    csv[:,5] = data['w']
    csv[:,6] = data['pmt_type']

    # reverse PMT normals
    csv[:,3] = -csv[:,3]
    csv[:,4] = -csv[:,4]
    csv[:,5] = -csv[:,5]

    header = "Format: x, y, z, u, v, w, pmt_type"

    np.savetxt(args.output, csv, fmt=['%10.2f']*6 + ['%i'],
               header=header)
>.groupby('id')): if id == 20: plt.subplot(2,3,1) elif id == 22: plt.subplot(2,3,2) elif id == 2020: plt.subplot(2,3,4) elif id == 2022: plt.subplot(2,3,5) elif id == 2222: plt.subplot(2,3,6) if muons: plt.hist(np.log10(df_id.ke.values/1000), bins=np.linspace(0,4.5,100), histtype='step') plt.xlabel("log10(Energy (GeV))") else: bins = np.logspace(np.log10(20),np.log10(10e3),21) plt.hist(df_id.ke.values, bins=bins, histtype='step') plt.gca().set_xscale("log") plt.xlabel("Energy (MeV)") plt.title('$' + ''.join([particle_id[int(''.join(x))] for x in grouper(str(id),2)]) + '$') if len(df): plt.tight_layout() def plot_hist(df, muons=False): 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) if muons: plt.hist(np.log10(df_id.ke.values/1000), bins=np.linspace(0,4.5,100), histtype='step') plt.xlabel("log10(Energy (GeV))") else: plt.hist(df_id.ke.values, bins=np.linspace(20,10e3,100), histtype='step') plt.xlabel("Energy (MeV)") plt.title(str(id)) if len(df): plt.tight_layout() if __name__ == '__main__': import argparse import numpy as np import pandas as pd import sys import h5py from sddm.plot_energy import * from sddm.plot import * from sddm import setup_matplotlib parser = argparse.ArgumentParser("plot fit results") parser.add_argument("filenames", nargs='+', help="input files") parser.add_argument("--save", action='store_true', default=False, help="save corner plots for backgrounds") parser.add_argument("--nhit-thresh", type=int, default=None, help="nhit threshold to apply to events before processing (should only be used for testing to speed things up)") args = parser.parse_args() setup_matplotlib(args.save) import matplotlib.pyplot as plt # Loop over runs to prevent using too much memory evs = [] data_filenames = [filename for filename in args.filenames if 'SNO' in filename] mc_filenames = [filename for filename in args.filenames if 'SNO' not in filename] if len(data_filenames): rhdr = pd.concat([read_hdf(filename, "rhdr").assign(filename=filename) for filename in data_filenames],ignore_index=True) for run, df in rhdr.groupby('run'): evs.append(get_events(df.filename.values, merge_fits=True, nhit_thresh=args.nhit_thresh)) if len(mc_filenames): evs.append(get_events(mc_filenames, merge_fits=True, nhit_thresh=args.nhit_thresh)) ev = pd.concat(evs) ev = ev[ev.prompt & ~np.isnan(ev.fmin)] ev = ev[ev.ke > 20] with pd.option_context('display.max_rows', None, 'display.max_columns', None): print("Noise events") print(ev[ev.noise][['psi','x','y','z','id1','id2']]) print("Muons") print(ev[ev.muon][['psi','r','id1','id2','id3','energy1','energy2','energy3']]) print("Neck") print(ev[ev.neck & ev.psi < 6][['psi','r','id1','cos_theta']]) print("Flashers") print(ev[ev.flasher & ev.udotr > 0]) print("Signal") print(ev[ev.signal]) # save as PDF b/c EPS doesn't support alpha values if args.save: plot_corner_plot(ev[ev.breakdown],"Breakdowns",save="breakdown_corner_plot") plot_corner_plot(ev[ev.muon],"Muons",save="muon_corner_plot") plot_corner_plot(ev[ev.flasher],"Flashers",save="flashers_corner_plot") plot_corner_plot(ev[ev.neck],"Neck",save="neck_corner_plot") plot_corner_plot(ev[ev.noise],"Noise",save="noise_corner_plot") plot_corner_plot(ev[ev.signal],"Signal",save="signal_corner_plot") else: plot_corner_plot(ev[ev.breakdown],"Breakdowns") plot_corner_plot(ev[ev.muon],"Muons") plot_corner_plot(ev[ev.flasher],"Flashers") plot_corner_plot(ev[ev.neck],"Neck") plot_corner_plot(ev[ev.noise],"Noise") plot_corner_plot(ev[ev.signal],"Signal") fig = plt.figure() plot_hist2(ev[ev.flasher]) despine(fig,trim=True) plt.suptitle("Flashers") fig = plt.figure() plot_hist2(ev[ev.muon],muons=True) despine(fig,trim=True) plt.suptitle("Muons") plt.show()