diff --git a/scripts/analyze_waveform.py b/scripts/analyze_waveform.py new file mode 100644 index 0000000000000000000000000000000000000000..0298e5f037c002835446a209d9484e62c03be096 --- /dev/null +++ b/scripts/analyze_waveform.py @@ -0,0 +1,64 @@ +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) + diff --git a/scripts/draw_events.py b/scripts/draw_events.py index c7b63b46ccdb70a7a97d8f8c98a5d31319fad963..4ef063886dd2c8eb0517c2b3f7d23878d14a30d5 100644 --- a/scripts/draw_events.py +++ b/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') diff --git a/src/esb_analyze.cpp b/src/esb_analyze.cpp index b6fb1a1ae7ae5bab70e70ac5da6c36cb6c663066..ed1334ff42cd7845218c75646431d220b9c373a8 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, 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();