Commit c219b5c7 authored by Remy Moll's avatar Remy Moll
Browse files

More data

parent 78202c50
......@@ -234,9 +234,9 @@
],
"source": [
"anode_voltage = 2400\n",
"cathode_voltage = 4000\n",
"events=15000\n",
"delay=0 #0 for default\n",
"cathode_voltage = 1000\n",
"events=5000\n",
"delay=758 #0 for default\n",
"\n",
"events = DECODER.get_data(cathode_voltage, events=events, delay=delay)\n",
"\n",
%% Cell type:code id: tags:
 
```
import numpy as np
import matplotlib
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/")
```
 
%% 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
 
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
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
 
 
# 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(4000, events=15000)#, delay=504)
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[:50]):
peak_ind = channel_peak_indexes[i]
for j, ch in enumerate(event):
x = ch[...,0]
y = ch[...,1]
ind = peak_ind[j]
 
plt.plot(x,y)
plt.plot(x[ind], y[ind], "+")
counter += 1
 
print(counter)
 
 
plt.show()
```
 
%%%% Output: stream
 
basic_noise_filter kept 11000 events out of 14999 (-> 73% kept.)
peak_noise_filter kept 4763 events out of 11000 (-> 43% kept.)
400
 
%%%% Output: display_data
 
[Hidden Image Output]
 
%% 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 = 4000
events=15000
delay=0 #0 for default
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)
complex_filtered, channel_peak_indexes = peak_noise_filter(base_filtered, threshold=15, ret_peaks=True, peak_t=True)
 
drift_velocity_array = []
 
nbins = 25
fig, axs = plt.subplots(4, 2, sharex=True)
subplots = axs.flat
fig.suptitle("Voltage = ${}$V".format(cathode_voltage))
 
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)
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 11000 events out of 14999 (-> 73% kept.)
peak_noise_filter kept 4763 events out of 11000 (-> 43% kept.)
Drift Velocity Channel 1: 4.624873275491098 cm * us^-1
Drift Velocity Channel 2: 4.781531436711113 cm * us^-1
Drift Velocity Channel 3: 4.8733604591528845 cm * us^-1
Drift Velocity Channel 4: 4.853026146421106 cm * us^-1
Drift Velocity Channel 5: 4.819756533991151 cm * us^-1
Drift Velocity Channel 6: 4.691955468690552 cm * us^-1
Drift Velocity Channel 7: 4.544244333425682 cm * us^-1
Drift Velocity Channel 8: -4.140187085683373 cm * us^-1
 
%%%% Output: display_data
 
[Hidden Image Output]
 
%%%% Output: stream
 
Electric Field: 0.6786850477200425
Mean Drift Velocity: 3.631070071025027
Standard Deviation DV: 2.9392270792080044
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
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