add a analysis script

scripts/analyze_waveform.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('module_file', help='a json database of modules and channels')
+parser.add_argument('output_file', help='Output csv file')
+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()
+
+
+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)])
+
+for iev in np.arange(0, tree.GetEntries()):
+    tree.GetEntry(iev)
+    if iev % 100 == 0:
+        print('processed {}'.format(iev), end='\r')
+    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)
+

scripts/draw_events.py

@@ -25,7 +25,7 @@ def fit_pedestal(bk, outlier_thres=2.0, min_lengths=10):
     return ped_mean, ped_err
 
 
-def get_peak(vals, thres=8.0, twindow=(0, 192), zero_thres=5.0):
+def get_peak(vals, thres=8.0, twindow=(0, 63), zero_thres=5.0):
     if len(vals) < 1:
         return 0, 0, 0, 0
     # find peaks
@@ -70,13 +70,21 @@ 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_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',
+    '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',
     'B', 'B', 'B', 'B',
@@ -117,6 +125,7 @@ for iev in evs:
         )
 
     fig, axs = plt.subplots(nrows, ncols, figsize=figsize)
+    lines, labels = [], []
     for (i, ax) in enumerate(axs.flat):
         box = dict(boxstyle='round', facecolor='white', alpha=0.5)
         ax.text(0.68, 0.83, ch[i], transform=ax.transAxes, fontsize=16, verticalalignment='top', bbox=box)
@@ -124,7 +133,6 @@ for iev in evs:
         ax.patch.set_alpha(pmt_colors[pmt[i]][1])
 
         sp = samples[ch[i]]
-        lines, labels = [], []
         if len(sp) > 0:
             pp = pprops[ch[i]]
             ax.plot(np.arange(len(sp)), sp, label='Waveform Data')

src/esb_analyze.cpp

@@ -315,7 +315,7 @@ void waveform_analysis(const std::vector<uint32_t> &raw, BranchData &res)
     // find peaks
     TSpectrum s;
     s.SetResolution(0.5);
-    int npeaks = s.SearchHighRes(wfbuf, bkbuf, res.nraw, 2.0, 5, false, 3, true, 3);
+    int npeaks = s.SearchHighRes(wfbuf, bkbuf, res.nraw, 3.0, 10, false, 5, true, 3);
 
     // fill branch data
     double *pos = s.GetPositionX();