Commit 3c220ee9 authored by Remy Moll's avatar Remy Moll
Browse files
parents 31510004 aad4ba33
......@@ -14,7 +14,9 @@
"from scipy import signal\n",
"import decode\n",
"\n",
"DECODER = decode.ArrayLoader(\"data/\")"
"DECODER = decode.ArrayLoader(\"data/\")\n",
"\n",
"# Just a push"
]
},
{
......
%% Cell type:code id: tags:
```
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [20, 10]
import scipy as sc
from scipy import signal
import decode
DECODER = decode.ArrayLoader("data/")
# Just a push
```
%% Cell type:markdown id: tags:
# Useful functions
%% Cell type:code id: tags:
```
def basic_noise_filter(events):
list_events = []
for e in events:
keep = True
for ch in e:
abs_ch = np.abs(ch[...,1])
if np.average(abs_ch) > 0.5 * np.max(abs_ch):
# if the average of the value grows toward the peak value, either the signal is nearly constant (empty noise) or the signal fluctuates a lot and has a high average (amplification noise)
keep = False
if keep:
list_events.append(e)
print("basic_noise_filter kept {} events out of {} (-> {}% kept.)".format(len(list_events), len(events), int(100 * len(list_events) / len(events))))
return list_events
def peak_noise_filter(events, threshold=10, ret_peaks=False, peak_t=False):
list_events = []
list_peaks = []
for e in events:
keep = True
peaks_ind = [] # peak of each channel for one particular event
for ch in e:
abs_ch = np.abs(ch[...,1])
# avg = np.average(ch[...,1])
# extreme = np.max(ch[...,1])
# if np.abs(extreme - avg) <= 0.04 or avg >= 0.02:
# # a high difference means it's likely more than just noise
# # a low difference across at least one channel means that there is noise
# keep = False
peak_ind, _ = signal.find_peaks(abs_ch, height=0.07, width=5)
# TODO set height relative to average
peaks_y = abs_ch[peak_ind]
if len(peak_ind) == 0 or len(peak_ind) > threshold:
# no peak found or too many peaks. in both cases the signal is likely bad
# keep = False
# break
pass
elif len(peak_ind) < int(0.5 * threshold): # few actual peaks, keep the first
peaks_ind.append(peak_ind[0])
else:
# threshold/2 < len(peak_ind) < threshold
# -> relatively few events, we filter out the noise by taking the maximum
if np.max(peaks_y) <= 2 * np.average(peaks_y):
# no clear global maximum -> noise? (clear as in outstanding, implying a fluctuating signal again)
# TODO check how much impact this has
# keep = False
# break
pass
else:
peaks_ind.append(peak_ind[np.argmax(peaks_y)])
# keep the index of the global maximum, the other ones must have been background noise since the highest peak was relatively outstanding
if keep:
list_events.append(e)
if ret_peaks:
if peak_t:
list_peaks.append(ch[...,0][peaks_ind]) # return actual timestamps
else:
list_peaks.append(peaks_ind) # return raw indices
print("peak_noise_filter kept {} events out of {} (-> {}% kept.)".format(len(list_events), len(events), int(100 * len(list_events) / len(events))))
if ret_peaks:
return list_events, np.array(list_peaks)
else:
return list_events
def rebin_signal(to_bin, nbins=3):
rebin = to_bin
for i in range(nbins):
rebin = (rebin[:-1:2] + rebin[1::2])/2
return rebin
# Data originally in units: length in mm, time in ns
def calculate_drift_velocity(min_time, max_time):
scale_length = 1e-1 # corresponds to cm
scale_time = 1e-3 # corresponds to microseconds
chamber_width = 94.3 * scale_length
drift_time = (max_time - min_time)*scale_time
drift_velocity = chamber_width/drift_time
return drift_velocity
# Data originally in units: voltage in V, length in mm
def calculate_electric_field(cathode_voltage, anode_voltage):
scale_length = 1e-1 # corresponds to cm
scale_potential = 1e-3 # corresponds to kV
chamber_width = 94.3 * scale_length
potential = (cathode_voltage + anode_voltage)*scale_potential
electric_field = potential/chamber_width
return electric_field
```
%% Cell type:code id: tags:
```
events = DECODER.get_data(2000, events=5000, delay=603)
base_filtered = basic_noise_filter(events)
complex_filtered, channel_peak_indexes = peak_noise_filter(base_filtered, threshold=15, ret_peaks=True, peak_t=False)
plt.figure()
counter = 0
for i, event in enumerate(complex_filtered[:5]):
peak_ind = channel_peak_indexes[i]
for j, ch in enumerate(event[:3]):
x = rebin_signal(ch[...,0], 5)
y = rebin_signal(ch[...,1], 5)
# ind = peak_ind[j]
plt.plot(x,y)
plt.plot(ch[...,0], ch[...,1])
# plt.plot(x[ind], y[ind], "+")
counter += 1
plt.show()
print(counter)
```
%%%% Output: stream
0
Warning: No data found.
%%%% Output: error
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_17428/1669732101.py in <module>
1 events = DECODER.get_data(2000, events=5000, delay=603)
----> 2 base_filtered = basic_noise_filter(events)
3 complex_filtered, channel_peak_indexes = peak_noise_filter(base_filtered, threshold=15, ret_peaks=True, peak_t=False)
4 plt.figure()
5 counter = 0
~\AppData\Local\Temp/ipykernel_17428/1075906179.py in basic_noise_filter(events)
1 def basic_noise_filter(events):
2 list_events = []
----> 3 for e in events:
4 keep = True
5 for ch in e:
TypeError: 'NoneType' object is not iterable
%% Cell type:code id: tags:
```
def gaussian(x, *p):
A, mu, sigma = p
return A * np.exp(-(x -mu)**2/(2*sigma**2))
```
%% Cell type:code id: tags:
```
anode_voltage = 2400
cathode_voltage = 1000
events=5000
delay=758 #0 for default
events = DECODER.get_data(cathode_voltage, events=events, delay=delay)
base_filtered = basic_noise_filter(events)
rebins = np.array([[rebin_signal(ch, nbins=3) for ch in base_filtered[...,i]] for i in range(len(base_filtered[1]))])
complex_filtered, channel_peak_indexes = peak_noise_filter(rebins, threshold=15, ret_peaks=True, peak_t=True)
drift_velocity_array = []
nbins = 20
fig, axs = plt.subplots(4, 2, sharex=True)
subplots = axs.flat
fig.suptitle("Voltage = ${}$V".format(cathode_voltage))
print(channel_peak_indexes.shape)
for i in range(channel_peak_indexes.shape[1]): # usually 8
subplt = subplots[i]
subplt.set(xlim=[0,3000])
subplt.set_title("Channel {}".format(i+1))
subplt.set(ylabel='Number of events', xlabel='Time (ns)')
hist, bin_edges = np.histogram(channel_peak_indexes[...,i], bins=nbins)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
guess = [80, 1500, 40]
# guess = [300, 1800, 400]
# coeff, var_matrix = sc.optimize.curve_fit(gaussian, bin_centers, hist, p0=guess)
coeff = guess # no actual fitting
subplt.hist(channel_peak_indexes[...,i], bins=nbins)
subplt.plot(bin_centers, gaussian(bin_centers, *coeff))
subplt.axvline(coeff[1], linestyle='dashed', linewidth=1, label="$\mu$")
nsigma = 1.9
max_time = coeff[1] + nsigma * coeff[2]
min_time = coeff[1] - nsigma * coeff[2]
subplt.axvline(max_time, linewidth=1, label="$\mu + {} \sigma$".format(nsigma))
subplt.axvline(min_time, linewidth=1, label="$\mu - {} \sigma$".format(nsigma))
subplt.legend()
drift_velocity = calculate_drift_velocity(min_time, max_time)
print(f"Drift Velocity Channel {i+1}: {drift_velocity} cm * us^-1")
drift_velocity_array.append(drift_velocity)
plt.show()
drift_velocity_array = np.array(drift_velocity_array[3:-1])
mean_drift_velocity = np.mean(drift_velocity_array)
std_drift_velocity = np.std(drift_velocity_array)
electric_field = calculate_electric_field(cathode_voltage, anode_voltage)
print(f"Electric Field: {electric_field}" +'\n' +
f"Mean Drift Velocity: {mean_drift_velocity}" + '\n' +
f"Standard Deviation DV: {std_drift_velocity}" + '\n')
```
%%%% Output: stream
basic_noise_filter kept 2845 events out of 4999 (-> 56% kept.)
%%%% Output: error
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_13400/4254197204.py in <module>
7
8 base_filtered = basic_noise_filter(events)
----> 9 rebins = np.array([[rebin_signal(ch, nbins=3) for ch in base_filtered[...,i]] for i in range(len(base_filtered[1]))])
10 complex_filtered, channel_peak_indexes = peak_noise_filter(rebins, threshold=15, ret_peaks=True, peak_t=True)
11 drift_velocity_array = []
~\AppData\Local\Temp/ipykernel_13400/4254197204.py in <listcomp>(.0)
7
8 base_filtered = basic_noise_filter(events)
----> 9 rebins = np.array([[rebin_signal(ch, nbins=3) for ch in base_filtered[...,i]] for i in range(len(base_filtered[1]))])
10 complex_filtered, channel_peak_indexes = peak_noise_filter(rebins, threshold=15, ret_peaks=True, peak_t=True)
11 drift_velocity_array = []
TypeError: list indices must be integers or slices, not tuple
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment