diff --git a/scripts/analyze_waveform.py b/scripts/analyze_waveform.py
index 8c63f4c53617f2f6d2b79f569db90fb7d72c5c80..1597b4ea18410e49f01ca834aaca5be70d1289cb 100644
--- a/scripts/analyze_waveform.py
+++ b/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')
diff --git a/src/esb_analyze.cpp b/src/esb_analyze.cpp
index ed1334ff42cd7845218c75646431d220b9c373a8..f69262601e5381b2d15d32f419d528f800ac1f2a 100644
--- a/src/esb_analyze.cpp
+++ b/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();