-
Notifications
You must be signed in to change notification settings - Fork 3
/
analyze-waveforms
executable file
·163 lines (137 loc) · 6.32 KB
/
analyze-waveforms
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
#!/usr/bin/env python3
from __future__ import print_function, division
import h5py
import numpy as np
def iqr(x):
return np.percentile(x,75) - np.percentile(x,25)
def get_times(x, data, baseline=10):
"""
Returns the times at which the waveforms in `x` cross 40% of their minimum
value.
"""
data = np.asarray(data)
# Get the first 10 ns of every waveform to calculate the noise level
noise = iqr(data[:,np.where(x < x[0] + baseline)[0]])
# Get events with a pulse
pulses = np.min(data,axis=-1) < -noise*5
# Select events with pulses. If there are no matches (which might be the
# case for the triggering channel), then don't apply the selection.
if np.count_nonzero(pulses):
data = data[pulses]
argmin = np.argmin(data,axis=-1)
threshold = 0.4*data[np.arange(data.shape[0]),argmin]
return x[data.shape[1]-np.argmax(((np.arange(data.shape[1]) < argmin[:,np.newaxis]) & (data > threshold[:,np.newaxis]))[:,::-1],axis=-1)-1]
def get_window(x, data, left=1, right=10):
"""
Returns the indices start and stop over which you should integrate the
waveforms in `x`. The window is found by calculating the median hit time
for all pulses in `x` and then going back 10 ns and forward 125 ns for sodium and 10 ns otherwise.
"""
data = np.asarray(data)
t = get_times(x,data)
mean_hit_time = np.median(t)
a, b = np.searchsorted(x,[mean_hit_time-left,mean_hit_time+right])
if a < 0:
a = 0
if b > len(x) - 1:
b = len(x) - 1
return a, b
def integrate(x, data, a, b):
"""
Integrate all waveforms in `data` with times `x`.
"""
# i = v/r
# divide by 50 ohms to convert to a charge
return -np.trapz(data[:,a:b],x=x[a:b])*1000/50.0
def get_bins(x):
"""
Returns bins for the data `x` using the Freedman Diaconis rule. See
https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule.
The equation below returns bins using said rule. A smaller bin value than this was chosen after finetuning for better fitting results.
bin_width = 2*iqr(x)/(len(x)**(1/3.0))
"""
bin_width = 0.5*iqr(x)/(len(x)**(1/3.0))
return np.arange(np.min(x),np.max(x),bin_width)
def chunks(lst, n):
"""Yield successive n-sized chunks from lst."""
for i in range(0, len(lst), n):
yield (i,i + n)
def convert_data(f, channel, start, stop):
if 'yinc' in dict(f[channel].attrs):
x = f[channel].attrs['xorg'] + np.linspace(0,f[channel].attrs['xinc']*f[channel].attrs['points'],int(f[channel].attrs['points']))
if ':WAVeform:FORMat' in dict(f['settings'].attrs) and f['settings'].attrs[':WAVeform:FORMat'] != 'ASC':
# convert word values -> voltages if the data was saved in a non-ascii format
y = f[channel][start:stop]*f[channel].attrs['yinc'] + f[channel].attrs['yorg']
else:
y = f[channel][start:stop]
else:
# In older versions of the code, I stored xorg, xinc, etc.
# in the main HDF5 group and not on a per channel basis
x = f.attrs['xorg'] + np.linspace(0,f.attrs['xinc']*f.attrs['points'],int(f.attrs['points']))
if ':WAVeform:FORMat' in dict(f['settings'].attrs) and f['settings'].attrs[':WAVeform:FORMat'] != 'ASC':
# convert word values -> voltages if the data was saved in a non-ascii format
y = f[channel][start:stop]*f.attrs['yinc'] + f.attrs['yorg']
else:
y = f[channel][start:stop]
return x*1e9, y
if __name__ == '__main__':
from argparse import ArgumentParser
import ROOT
import matplotlib.pyplot as plt
parser = ArgumentParser(description='Analyze data from the Agilent scope')
parser.add_argument('filenames',nargs='+',help='input filenames (hdf5 format)')
parser.add_argument('-o','--output', default=None, help='output file name', required=True)
parser.add_argument('--sodium', default=False, action='store_true', help='flag to indicate data is from a sodium source')
parser.add_argument('--laser', default=False, action='store_true', help='flag to indicate data is from a laser source')
parser.add_argument('--plot', default=False, action='store_true', help='plot the waveforms and charge integral')
parser.add_argument('--chunks', default=10000, type=int, help='number of waveforms to process at a time')
args = parser.parse_args()
charge = {}
for filename in args.filenames:
with h5py.File(filename) as f:
for channel in f:
if channel == 'settings':
continue
N = len(f[channel])
for i in range(0, N, args.chunks):
x, y = convert_data(f,channel,i,i+args.chunks)
if i == 0:
if args.sodium:
a, b = get_window(x,y,left=50,right=300)
else:
a, b = get_window(x,y,left=50,right=300)
if channel in charge:
charge[channel] = np.concatenate((charge[channel],integrate(x,y, a, b)))
else:
charge[channel] = integrate(x,y, a, b)
if args.plot:
plt.figure()
plt.subplot(2,1,1)
plt.plot(x,y[:100].T)
plt.xlabel("Time (ns)")
plt.ylabel("Voltage (V)")
plt.axvline(x[a])
plt.axvline(x[b])
plt.subplot(2,1,2)
plt.plot(x,np.median(y,axis=0))
plt.xlabel("Time (ns)")
plt.ylabel("Voltage (V)")
plt.axvline(x[a])
plt.axvline(x[b])
plt.suptitle(channel)
f = ROOT.TFile(args.output,"recreate")
for channel in charge:
bins = get_bins(charge[channel])
h = ROOT.TH1D(channel,"Charge Integral for %s" % channel,len(bins),bins[0],bins[-1])
for x in charge[channel]:
h.Fill(x)
h.GetXaxis().SetTitle("Charge (pC)")
if args.plot:
plt.figure()
plt.hist(charge[channel],bins=bins,histtype='step',label=channel)
plt.xlabel("Charge (pC)")
plt.legend()
h.Write()
f.Close()
if args.plot:
plt.show()