add analysis scripts

scripts/analyze_waveform.py

@@ -26,8 +26,8 @@ def get_channels(json_path):
 parser = argparse.ArgumentParser('Raw waveform analysis')
 
 parser.add_argument('root_file', help='a root file of waveform data')
-parser.add_argument('module_file', help='a json database of modules and channels')
-parser.add_argument('output_file', help='Output csv file')
+parser.add_argument('output_dir', help='output directory')
+parser.add_argument('--prefix', dest='prefix', help='prefix for the output files', default='')
 parser.add_argument('--peak-thres', dest='pthres', help='sigma threshold for finding the peak, default 8', type=float, default=8)
 parser.add_argument('--peak-min', dest='pmin', help='lower limit of the peak adc value, default 10', type=int, default=10)
 parser.add_argument('--peak-max', dest='pmax', help='upper limit of the peak adc value, default 8000', type=int, default=8000)
@@ -36,29 +36,79 @@ parser.add_argument('--time-max', dest='tmax', help='upper limit of the time win
 
 args = parser.parse_args()
 
+# figsize = (24, 16)
+# nrows, ncols = 3, 4
+# ch = [
+#     'S2_1', 'S2_2', 'S2_3', 'S2_4',
+#     'S2_5', 'S2_6', 'S2_7', 'S2_8',
+#     'S2_9', 'S2_10', 'S2_11',
+# ]
+
+# figsize = (16, 16)
+# nrows, ncols = 3, 3
+# ch = [
+#     'C4', 'C6_4', 'C7_4',
+#     'C5_4', 'C9_4', 'C8_4',
+#     'C1', 'C2', 'C3',
+# ]
+
+figsize = (16, 16)
+nrows, ncols = 4, 4
+ch = [
+    'Cer11_5', 'Cer12_5', 'Cer13_5', 'Cer14_5',
+    'Cer21_5', 'Cer22_5', 'Cer23_5', 'Cer24_5',
+    'Cer31_5', 'Cer32_5', 'Cer33_5', 'Cer34_5',
+    'Cer41_5', 'Cer42_5', 'Cer43_5', 'Cer44_5',
+]
 
 f = ROOT.TFile.Open(args.root_file, 'read')
 tree = f.EvTree
-ch = get_channels(args.module_file)
-res = pd.DataFrame(columns=[c + '_peak' for c in list(ch)] + [c + '_pos' for c in list(ch)])
+
+props = {c: [np.array([]), np.array([])] for c in ch}
 
 for iev in np.arange(0, tree.GetEntries()):
     tree.GetEntry(iev)
-    if iev % 100 == 0:
+    if iev % 1000 == 0:
         print('processed {}'.format(iev), end='\r')
+
+    # channels
     evpeaks, evposes = [], []
-    for c in ch:
-        br_peak = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
-        br_pos = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
-
-        ch_peak, ch_pos = 0, 0
-        for peak, pos in zip(br_peak, br_pos):
-            if pos < args.tmax and pos > args.tmin and peak > ch_peak:
-                ch_peak = peak
-                ch_pos = pos
-        evpeaks.append(ch_peak)
-        evposes.append(ch_pos)
-    res.loc[iev, :] = evpeaks + evposes
-
-res.to_csv(args.output_file)
+    for c, prop in props.items():
+        prop[0] = np.concatenate((prop[0], branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)))
+        prop[1] = np.concatenate((prop[1], branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)))
+print('processed {}'.format(iev))
+
+
+bins = np.arange(0, 5000, step=1)
+indices = (bins[1:] + bins[:-1])/2.
+peaks = pd.DataFrame(index=indices, data={c: np.histogram(prop[0], bins)[0] for c, prop in props.items()})
+
+bins = np.arange(0, 64, step=1)
+indices = (bins[1:] + bins[:-1])/2.
+poses = pd.DataFrame(index=indices, data={c: np.histogram(prop[1], bins)[0] for c, prop in props.items()})
+
+
+# plot
+def plot_hist(df, ny, nx, x_label, y_label, fs=(16, 16), fontsize=18):
+    box = dict(boxstyle='round', facecolor='wheat', alpha=0.3)
+    fig, axs = plt.subplots(ny, nx, figsize=fs)
+    for i, ax in enumerate(axs.flat):
+        if i >= len(df.columns):
+            continue
+        ax.text(0.75, 0.90, df.columns[i], transform=ax.transAxes, fontsize=fontsize - 4, verticalalignment='top', bbox=box)
+        ax.bar(df.index.values, df.iloc[:, i].values, width=pd.Series(data=df.index).diff(1).mean())
+        # ax.patch.set_facecolor('wheat')
+        # ax.patch.set_alpha(0.05)
+        # ax.set_yscale('log')
+    fig.tight_layout(rect=(0.03, 0.05, 0.98, 0.95))
+    fig.text(0.5, 0.02, x_label, ha='center', fontsize=fontsize)
+    fig.text(0.02, 0.5, y_label, ha='center', rotation=90, fontsize=fontsize)
+    return fig
+
+
+plot_hist(peaks, nrows, ncols, 'ADC Value', 'Counts', figsize).savefig(os.path.join(args.output_dir, args.prefix + 'peaks.png'))
+plot_hist(poses, nrows, ncols, 'Sample Number', 'Counts', figsize).savefig(os.path.join(args.output_dir, args.prefix + 'timings.png'))
+
+peaks.to_csv(os.path.join(args.output_dir, args.prefix + 'peaks.csv'))
+poses.to_csv(os.path.join(args.output_dir, args.prefix + 'timings.csv'))
 

scripts/analyze_waveform_time.py

+import ROOT
+import os
+import json
+import numpy as np
+import pandas as pd
+import argparse
+from matplotlib import pyplot as plt
+from scipy import signal
+from collections.abc import Iterable
+
+
+def branch_to_array1d(br, t):
+    return np.ndarray((len(br), ), dtype=t, buffer=br)
+
+
+def get_channels(json_path):
+    dbf = open(json_path)
+    db = json.load(dbf)
+    channels = dict()
+    for mod, mod_prop in db.items():
+        for ch in mod_prop['channels']:
+            channels[ch['name']] = ch['type']
+    return channels
+
+
+parser = argparse.ArgumentParser('Raw waveform analysis')
+
+parser.add_argument('root_file', help='a root file of waveform data')
+parser.add_argument('output_dir', help='output directory')
+parser.add_argument('--prefix', dest='prefix', help='prefix for the output files', default='')
+parser.add_argument('--peak-thres', dest='pthres', help='sigma threshold for finding the peak, default 8', type=float, default=8)
+parser.add_argument('--peak-min', dest='pmin', help='lower limit of the peak adc value, default 10', type=int, default=10)
+parser.add_argument('--peak-max', dest='pmax', help='upper limit of the peak adc value, default 8000', type=int, default=8000)
+parser.add_argument('--time-min', dest='tmin', help='lower limit of the time window, default 0', type=int, default=0)
+parser.add_argument('--time-max', dest='tmax', help='upper limit of the time window, default 63', type=int, default=63)
+
+args = parser.parse_args()
+
+# figsize = (24, 16)
+# nrows, ncols = 3, 4
+# ch = [
+#     'S2_1', 'S2_2', 'S2_3', 'S2_4',
+#     'S2_5', 'S2_6', 'S2_7', 'S2_8',
+#     'S2_9', 'S2_10', 'S2_11',
+# ]
+
+# figsize = (16, 16)
+# nrows, ncols = 3, 3
+# ch = [
+#     'C4', 'C6_4', 'C7_4',
+#     'C5_4', 'C9_4', 'C8_4',
+#     'C1', 'C2', 'C3',
+# ]
+
+figsize = (16, 16)
+nrows, ncols = 4, 4
+ch = [
+    'Cer11_5', 'Cer12_5', 'Cer13_5', 'Cer14_5',
+    'Cer21_5', 'Cer22_5', 'Cer23_5', 'Cer24_5',
+    'Cer31_5', 'Cer32_5', 'Cer33_5', 'Cer34_5',
+    'Cer41_5', 'Cer42_5', 'Cer43_5', 'Cer44_5',
+]
+
+ch_pos = [
+    (20, 30), (20, 30), (20, 30), (20, 30),
+    (25, 35), (25, 35), (25, 35), (23, 33),
+    (23, 33), (23, 33), (23, 33), (23, 33),
+    (20, 30), (20, 30), (15, 25), (25, 35),
+]
+
+f = ROOT.TFile.Open(args.root_file, 'read')
+tree = f.EvTree
+
+props = {c: [np.array([]), np.array([])] for c in ch}
+cher_sums = pd.DataFrame(index=np.arange(0, tree.GetEntries()), columns=['sum', 'nhits'], dtype=('float64', 'int64'))
+
+for iev in np.arange(0, tree.GetEntries()):
+    tree.GetEntry(iev)
+    if iev > 1000:
+        break
+    if iev % 1000 == 0:
+        print('processed {}'.format(iev), end='\r')
+
+    csum = 0
+    cnhits = 0
+    # channels
+    evpeaks, evposes = [], []
+    for p, (c, prop) in zip(ch_pos, props.items()):
+        epeak, epos = 0, 0
+        peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
+        poses = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
+        for pos, peak in zip(poses, peaks):
+            if pos <= p[1] and pos >= p[0] and peak > epeak:
+                epeak = peak
+                epos = pos
+
+        if epeak >= 1:
+            cnhits += 1
+            csum += epeak
+            prop[0] = np.append(prop[0], epeak)
+            prop[1] = np.append(prop[1], epos)
+
+    cher_sums.loc[iev] = (csum, cnhits)
+print('processed {}'.format(iev))
+
+
+bins = np.arange(0, 5000, step=1)
+indices = (bins[1:] + bins[:-1])/2.
+peaks = pd.DataFrame(index=indices, data={c: np.histogram(prop[0], bins)[0] for c, prop in props.items()})
+
+bins = np.arange(0, 64, step=1)
+indices = (bins[1:] + bins[:-1])/2.
+poses = pd.DataFrame(index=indices, data={c: np.histogram(prop[1], bins)[0] for c, prop in props.items()})
+
+
+# plot
+def plot_hist(df, ny, nx, x_label, y_label, fs=(16, 16), fontsize=18):
+    box = dict(boxstyle='round', facecolor='wheat', alpha=0.3)
+    fig, axs = plt.subplots(ny, nx, figsize=fs)
+    for i, ax in enumerate(axs.flat):
+        if i >= len(df.columns):
+            continue
+        ax.text(0.75, 0.90, df.columns[i], transform=ax.transAxes, fontsize=fontsize - 4, verticalalignment='top', bbox=box)
+        ax.bar(df.index.values, df.iloc[:, i].values, width=pd.Series(data=df.index).diff(1).mean())
+        # ax.patch.set_facecolor('wheat')
+        # ax.patch.set_alpha(0.05)
+        # ax.set_yscale('log')
+    fig.tight_layout(rect=(0.03, 0.05, 0.98, 0.95))
+    fig.text(0.5, 0.02, x_label, ha='center', fontsize=fontsize)
+    fig.text(0.02, 0.5, y_label, ha='center', rotation=90, fontsize=fontsize)
+    return fig
+
+
+plot_hist(peaks, nrows, ncols, 'ADC Value', 'Counts', figsize).savefig(os.path.join(args.output_dir, args.prefix + 'peaks.png'))
+plot_hist(poses, nrows, ncols, 'Sample Number', 'Counts', figsize).savefig(os.path.join(args.output_dir, args.prefix + 'timings.png'))
+
+peaks.to_csv(os.path.join(args.output_dir, args.prefix + 'peaks.csv'))
+poses.to_csv(os.path.join(args.output_dir, args.prefix + 'timings.csv'))
+
+cher_sums[cher_sums.notnull().all(axis=1)].to_csv(os.path.join(args.output_dir, args.prefix + 'sum.csv'))
+

scripts/draw_events.py

@@ -70,20 +70,27 @@ tree = f.EvTree
 # alignment of channels
 figsize = (16, 16)
 nrows, ncols = 4, 4
-# ch = [
-#     'Cer11_5', 'Cer12_5', 'Cer13_5', 'Cer14_5',
-#     'Cer21_5', 'Cer22_5', 'Cer23_5', 'Cer24_5',
-#     'Cer31_5', 'Cer32_5', 'Cer33_5', 'Cer34_5',
-#     'Cer41_5', 'Cer42_5', 'Cer43_5', 'Cer44_5',
-# ]
-
 ch = [
-    'Cer11_1', 'Cer11_2', 'Cer11_3', 'Cer11_4',
-    'Cer12_1', 'Cer12_2', 'Cer12_3', 'Cer12_4',
-    'Cer13_1', 'Cer13_2', 'Cer13_3', 'Cer13_4',
-    'Cer14_1', 'Cer14_2', 'Cer14_3', 'Cer14_4',
+    'Cer11_5', 'Cer12_5', 'Cer13_5', 'Cer14_5',
+    'Cer21_5', 'Cer22_5', 'Cer23_5', 'Cer24_5',
+    'Cer31_5', 'Cer32_5', 'Cer33_5', 'Cer34_5',
+    'Cer41_5', 'Cer42_5', 'Cer43_5', 'Cer44_5',
+]
+
+ch_pos = [
+    (20, 30), (20, 30), (20, 30), (20, 30),
+    (25, 35), (25, 35), (25, 35), (23, 33),
+    (23, 33), (23, 33), (23, 33), (23, 33),
+    (20, 30), (20, 30), (15, 25), (25, 35),
 ]
 
+# ch = [
+#     'Cer11_1', 'Cer11_2', 'Cer11_3', 'Cer11_4',
+#     'Cer12_1', 'Cer12_2', 'Cer12_3', 'Cer12_4',
+#     'Cer13_1', 'Cer13_2', 'Cer13_3', 'Cer13_4',
+#     'Cer14_1', 'Cer14_2', 'Cer14_3', 'Cer14_4',
+# ]
+
 
 pmt = [
     'A', 'A', 'A', 'A',
@@ -99,6 +106,16 @@ pmt_colors = {
     'D': ('lightskyblue', 0.05),
 }
 
+# figsize = (16, 16)
+# nrows, ncols = 6, 4
+# ch = [
+#     'C1', 'C2', 'C3', 'C4',
+#     'C5_1', 'C5_2', 'C5_3', 'C5_4',
+#     'C6_1', 'C6_2', 'C6_3', 'C6_4',
+#     'C7_1', 'C7_2', 'C7_3', 'C7_4',
+#     'C8_1', 'C8_2', 'C8_3', 'C8_4',
+#     'C9_1', 'C9_2', 'C9_3', 'C9_4',
+# ]
 
 def branch_to_array1d(br, t):
     return np.ndarray((len(br), ), dtype=t, buffer=br)
@@ -140,10 +157,12 @@ for iev in evs:
             xr = ax.get_xlim()
             ax.fill_between(xr, pp[2]-pp[3], pp[2]+pp[3], color='red', alpha=0.25, label='Pedestal Error')
             ax.set_xlim(xr[0], xr[1])
-            ax.axvline(x=args.tmin, color='k', linestyle=':', label='Time Window')
-            ax.axvline(x=args.tmax, color='k', linestyle=':')
+            time = ch_pos[i]
+            ax.axvline(x=time[0], color='k', linestyle=':', label='Time Window')
+            ax.axvline(x=time[1], color='k', linestyle=':')
             for peak, pos in zip(pp[0], pp[1]):
-                if pos <= args.tmax and pos >= args.tmin:
+                # if pos <= args.tmax and pos >= args.tmin:
+                if pos <= time[1] and pos >= time[0]:
                     ax.plot(pos, peak + pp[2], 'kv', label='Peak')
             lines, labels = ax.get_legend_handles_labels()
 

scripts/extract_calosum.py

+import ROOT
+import os
+import json
+import numpy as np
+import pandas as pd
+import argparse
+from matplotlib import pyplot as plt
+from scipy import signal
+from collections.abc import Iterable
+
+
+def branch_to_array1d(br, t):
+    return np.ndarray((len(br), ), dtype=t, buffer=br)
+
+
+def get_channels(json_path):
+    dbf = open(json_path)
+    db = json.load(dbf)
+    channels = dict()
+    for mod, mod_prop in db.items():
+        for ch in mod_prop['channels']:
+            channels[ch['name']] = ch['type']
+    return channels
+
+
+parser = argparse.ArgumentParser('Raw waveform analysis')
+
+parser.add_argument('root_file', help='a root file of waveform data')
+parser.add_argument('output_dir', help='output directory')
+parser.add_argument('--prefix', dest='prefix', help='prefix for the output files', default='')
+parser.add_argument('--peak-thres', dest='pthres', help='sigma threshold for finding the peak, default 8', type=float, default=8)
+parser.add_argument('--peak-min', dest='pmin', help='lower limit of the peak adc value, default 10', type=int, default=10)
+parser.add_argument('--peak-max', dest='pmax', help='upper limit of the peak adc value, default 8000', type=int, default=8000)
+parser.add_argument('--time-min', dest='tmin', help='lower limit of the time window, default 0', type=int, default=0)
+parser.add_argument('--time-max', dest='tmax', help='upper limit of the time window, default 63', type=int, default=63)
+
+args = parser.parse_args()
+
+calo_chs = [
+    'C1', 'C2', 'C3', 'C4',
+    'C5_1', 'C5_2', 'C5_3', 'C5_4',
+    'C6_1', 'C6_2', 'C6_3', 'C6_4',
+    'C7_1', 'C7_2', 'C7_3', 'C7_4',
+    'C8_1', 'C8_2', 'C8_3', 'C8_4',
+    'C9_1', 'C9_2', 'C9_3', 'C9_4',
+]
+calo_time = (25, 35)
+calo_thres = 1
+
+f = ROOT.TFile.Open(args.root_file, 'read')
+tree = f.EvTree
+
+props = {c: [np.array([]), np.array([])] for c in ch}
+calo_sums = pd.DataFrame(columns=['calo sum', 'calo nhits'], dtype='float64')
+
+for iev in np.arange(0, tree.GetEntries()):
+    tree.GetEntry(iev)
+    if iev % 1000 == 0:
+        print('processed {}'.format(iev), end='\r')
+
+    # calorimeter trigger
+    calo_sum = 0
+    calo_fires = 0
+    for c in calo_chs:
+        peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
+        times = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
+        cp, ct = 0, 0
+        for peak, time in zip(peaks, times):
+            if time <= calo_time[1] and time >= calo_time[0] and peak > cp:
+                cp = peak
+        calo_sum += cp
+        if cp > 1:
+            calo_fires += 1
+    if calo_sum < calo_thres:
+        continue
+    calo_sums.loc[iev] = (calo_sum, calo_fires)
+print('processed {}'.format(iev))
+
+calo_sums.to_csv(os.path.join(args.output_dir, args.prefix + 'calo_sum.csv'))
+

scripts/old/fit_npe.C

@@ -21,3 +21,4 @@ void fit_npe(const char *path, int nhit = 5)
     gaus1->Draw("same");
     gaus2->Draw("same");
 }
+