diff --git a/scripts/analyze_waveform.py b/scripts/analyze_waveform.py
index 0298e5f037c002835446a209d9484e62c03be096..8c63f4c53617f2f6d2b79f569db90fb7d72c5c80 100644
--- a/scripts/analyze_waveform.py
+++ b/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'))
diff --git a/scripts/analyze_waveform_time.py b/scripts/analyze_waveform_time.py
new file mode 100644
index 0000000000000000000000000000000000000000..f179ea9d2df7cdd24dc88410e426157d068c023f
--- /dev/null
+++ b/scripts/analyze_waveform_time.py
@@ -0,0 +1,141 @@
+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'))
+
diff --git a/scripts/draw_events.py b/scripts/draw_events.py
index 4ef063886dd2c8eb0517c2b3f7d23878d14a30d5..b751c9688e4fe0156a72f2c9f7acdc92087bedc7 100644
--- a/scripts/draw_events.py
+++ b/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()
diff --git a/scripts/extract_calosum.py b/scripts/extract_calosum.py
new file mode 100644
index 0000000000000000000000000000000000000000..389d34f6caad5ddc47bb71522736271ab5974c59
--- /dev/null
+++ b/scripts/extract_calosum.py
@@ -0,0 +1,80 @@
+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'))
+
diff --git a/scripts/old/fit_npe.C b/scripts/old/fit_npe.C
index da494bb5725346d7f47cc11a8af43116820b6c07..c261b5941ae0815515545badf198c76c9e6ed365 100644
--- a/scripts/old/fit_npe.C
+++ b/scripts/old/fit_npe.C
@@ -21,3 +21,4 @@ void fit_npe(const char *path, int nhit = 5)
gaus1->Draw("same");
gaus2->Draw("same");
}
+