update analysis scripts

scripts/analyze_waveform.py

@@ -26,32 +26,10 @@ 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('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)
+parser.add_argument('output', help='path to the output csv file')
 
 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 = [
@@ -61,9 +39,27 @@ ch = [
     'Cer41_5', 'Cer42_5', 'Cer43_5', 'Cer44_5',
 ]
 
+ch_pos = [
+    -15.5, -14.5, -15.5, -15.5,
+    -10.5, -9.5, -10.5, -14.5,
+    -14.5, -14.5, -14.5, -14.5,
+    -14.5, -14.5, -21.5, -9.5
+]
+pos_width = 3
+
+# trigger channels
+# tr = ['S2_1', 'S2_2', 'S2_3', 'S2_4', 'S2_5', 'S2_6', 'S2_7', 'S2_8', 'S2_9', 'S2_10', 'S2_11']
+# tr_time = (17, 23)
+tr = ['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']
+tr_time = (28, 32)
+
 f = ROOT.TFile.Open(args.root_file, 'read')
 tree = f.EvTree
 
+trg_ch = np.ndarray(shape=(tree.GetEntries(), ), dtype=object)
+trg_val = np.ndarray(shape=(tree.GetEntries(), 2), dtype='float64')
+ch_val = np.ndarray(shape=(tree.GetEntries(), len(ch)*2), dtype='float64')
 props = {c: [np.array([]), np.array([])] for c in ch}
 
 for iev in np.arange(0, tree.GetEntries()):
@@ -71,44 +67,58 @@ for iev in np.arange(0, tree.GetEntries()):
     if iev % 1000 == 0:
         print('processed {}'.format(iev), end='\r')
 
+    # triggers
+    ev_trg_peaks, ev_trg_poses = [], []
+    for c in tr:
+        peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
+        poses = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
+        tpeak, tpos = 0, 0
+        for peak, pos in zip(peaks, poses):
+            if pos <= tr_time[1] and pos >= tr_time[0] and peak > tpeak:
+                tpeak = peak
+                tpos = pos
+        ev_trg_peaks.append(tpeak)
+        ev_trg_poses.append(tpos)
+
+    # no triggers
+    if not len(ev_trg_peaks):
+        continue
+
+    # get the maximum peak from all trigger channels
+    ich = np.argmax(ev_trg_peaks)
+    trg_ch[iev] = tr[ich]
+    trg_val[iev] = (ev_trg_peaks[ich], ev_trg_poses[ich])
+
     # channels
-    evpeaks, evposes = [], []
-    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)))
+    for i, c in enumerate(ch):
+        peaks = branch_to_array1d(tree.__getattr__(c + '_Ppeak'), np.float32)
+        poses = branch_to_array1d(tree.__getattr__(c + '_Ptime'), np.float32)
+        props[c][0] = np.concatenate((props[c][0], peaks))
+        props[c][1] = np.concatenate((props[c][1], poses - ev_trg_poses[ich]))
+
+        # check timing
+        cpeak, cpos = 0, 0
+        for peak, pos in zip(peaks, poses):
+            pos -= ev_trg_poses[ich]
+            if pos <= ch_pos[i] + pos_width and pos >= ch_pos[i] - pos_width and peak > cpeak:
+                cpeak = peak
+                cpos = pos
+        ch_val[iev][[i, len(ch) + i]] = (cpeak, cpos)
+
 print('processed {}'.format(iev))
 
 
+result = pd.DataFrame(index=np.arange(0, tree.GetEntries()),
+                      columns=['trg_peak', 'trg_pos'] + [c + '_peak' for c in ch] + [c + '_pos' for c in ch],
+                      data=np.concatenate((trg_val, ch_val), axis=1))
+result.loc[:, 'trg_ch'] = trg_ch
+result.to_csv(args.output)
+
 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()})
+pd.DataFrame(index=indices, data={c: np.histogram(prop[0], bins)[0] for c, prop in props.items()}).to_csv('peaks.csv')
 
-bins = np.arange(0, 64, step=1)
+bins = np.arange(-64, 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'))
+pd.DataFrame(index=indices, data={c: np.histogram(prop[1], bins)[0] for c, prop in props.items()}).to_csv('timings.csv')
 

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, 3.0, 10, false, 5, true, 3);
+    int npeaks = s.SearchHighRes(wfbuf, bkbuf, res.nraw, 3.0, 20, false, 5, true, 3);
 
     // fill branch data
     double *pos = s.GetPositionX();