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

Superior plotting

parent 461ab165
This diff is collapsed.
This diff is collapsed.
......@@ -21,8 +21,13 @@ class Channel:
@property
def peak_index(self):
return self.peaks[0]
# TODO Check if the first value is the most accurate or if the index of maximum voltage is better
if self.peaks.size != 0:
absolute_peak = np.argmin(self.voltage[self.peaks])
return self.peaks[absolute_peak]
return self.peaks[0] # don't actually use the first value blindly
else:
return 0
@property
def peak_time(self):
......@@ -56,6 +61,10 @@ class Channel:
if len(self.peaks) > 5:
return False
if np.std(self.voltage[self.peaks]) >= 0.15:
# if all peaks have similar height, how would we know, which one is the right one?
return False
return True
......@@ -74,8 +83,8 @@ class Event:
@property
def is_usable(self):
good = sum([ch.is_usable for ch in self.channels])
ends_good = self.channels[0].is_usable and self.channels[-1].is_usable
return (good >= self.min_good_channels) and ends_good
# ends_good = self.channels[0].is_usable and self.channels[-1].is_usable
return (good >= self.min_good_channels)# and ends_good
def recreate_muon_path(self, drift_velocity):
......
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"plt.rcParams['figure.figsize'] = [20, 10]\n",
"\n",
"import scipy as sc\n",
"from scipy import signal\n",
"import decode\n",
"\n",
"DECODER = decode.ArrayLoader(\"data/\")\n",
"\n",
"# Just a push"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Useful functions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def basic_noise_filter(events):\n",
" list_events = []\n",
" for e in events:\n",
" keep = True\n",
" for ch in e:\n",
" abs_ch = np.abs(ch[...,1])\n",
" if np.average(abs_ch) > 0.5 * np.max(abs_ch):\n",
" # 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)\n",
" keep = False\n",
" if keep:\n",
" list_events.append(e)\n",
"\n",
" print(\"basic_noise_filter kept {} events out of {} (-> {}% kept.)\".format(len(list_events), len(events), int(100 * len(list_events) / len(events))))\n",
" return list_events\n",
"\n",
"\n",
"def peak_noise_filter(events, threshold=10, ret_peaks=False, peak_t=False):\n",
" list_events = []\n",
" list_peaks = []\n",
" for e in events:\n",
" keep = True\n",
" peaks_ind = [] # peak of each channel for one particular event\n",
" for ch in e:\n",
" abs_ch = np.abs(ch[...,1])\n",
" # avg = np.average(ch[...,1])\n",
" # extreme = np.max(ch[...,1])\n",
" # if np.abs(extreme - avg) <= 0.04 or avg >= 0.02:\n",
" # # a high difference means it's likely more than just noise\n",
" # # a low difference across at least one channel means that there is noise\n",
" # keep = False\n",
" peak_ind, _ = signal.find_peaks(abs_ch, height=0.07, width=5)\n",
" # TODO set height relative to average\n",
" peaks_y = abs_ch[peak_ind]\n",
"\n",
" if len(peak_ind) == 0 or len(peak_ind) > threshold:\n",
" # no peak found or too many peaks. in both cases the signal is likely bad\n",
" # keep = False\n",
" # break\n",
" pass\n",
"\n",
" elif len(peak_ind) < int(0.5 * threshold): # few actual peaks, keep the first\n",
" peaks_ind.append(peak_ind[0])\n",
"\n",
" else: \n",
" # threshold/2 < len(peak_ind) < threshold \n",
" # -> relatively few events, we filter out the noise by taking the maximum\n",
" if np.max(peaks_y) <= 2 * np.average(peaks_y):\n",
" # no clear global maximum -> noise? (clear as in outstanding, implying a fluctuating signal again)\n",
" # TODO check how much impact this has\n",
" # keep = False\n",
" # break\n",
" pass\n",
" else:\n",
" peaks_ind.append(peak_ind[np.argmax(peaks_y)])\n",
" # keep the index of the global maximum, the other ones must have been background noise since the highest peak was relatively outstanding\n",
"\n",
" if keep:\n",
" list_events.append(e)\n",
" if ret_peaks:\n",
" if peak_t:\n",
" list_peaks.append(ch[...,0][peaks_ind]) # return actual timestamps\n",
" else:\n",
" list_peaks.append(peaks_ind) # return raw indices\n",
"\n",
" print(\"peak_noise_filter kept {} events out of {} (-> {}% kept.)\".format(len(list_events), len(events), int(100 * len(list_events) / len(events))))\n",
"\n",
" if ret_peaks:\n",
" return list_events, np.array(list_peaks)\n",
" else:\n",
" return list_events\n",
"\n",
"\n",
"def rebin_signal(to_bin, nbins=3):\n",
" rebin = to_bin\n",
" for i in range(nbins):\n",
" rebin = (rebin[:-1:2] + rebin[1::2])/2\n",
" return rebin\n",
"\n",
"# Data originally in units: length in mm, time in ns\n",
"def calculate_drift_velocity(min_time, max_time):\n",
"\n",
" scale_length = 1e-1 # corresponds to cm\n",
" scale_time = 1e-3 # corresponds to microseconds\n",
"\n",
" chamber_width = 94.3 * scale_length\n",
" drift_time = (max_time - min_time)*scale_time\n",
" \n",
" drift_velocity = chamber_width/drift_time\n",
" return drift_velocity\n",
"\n",
"\n",
"# Data originally in units: voltage in V, length in mm\n",
"def calculate_electric_field(cathode_voltage, anode_voltage):\n",
" \n",
" scale_length = 1e-1 # corresponds to cm\n",
" scale_potential = 1e-3 # corresponds to kV\n",
" \n",
" chamber_width = 94.3 * scale_length\n",
" potential = (cathode_voltage + anode_voltage)*scale_potential\n",
"\n",
" electric_field = potential/chamber_width\n",
" return electric_field\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0\n",
"Warning: No data found.\n"
]
},
{
"ename": "TypeError",
"evalue": "'NoneType' object is not iterable",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_17428/1669732101.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mevents\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mDECODER\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_data\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m2000\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mevents\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m5000\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdelay\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m603\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mbase_filtered\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbasic_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mevents\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 3\u001b[0m \u001b[0mcomplex_filtered\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mchannel_peak_indexes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpeak_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbase_filtered\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthreshold\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mret_peaks\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpeak_t\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 4\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[0mcounter\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_17428/1075906179.py\u001b[0m in \u001b[0;36mbasic_noise_filter\u001b[1;34m(events)\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mbasic_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mevents\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0mlist_events\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[1;32mfor\u001b[0m \u001b[0me\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mevents\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[0mkeep\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mch\u001b[0m \u001b[1;32min\u001b[0m \u001b[0me\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mTypeError\u001b[0m: 'NoneType' object is not iterable"
]
}
],
"source": [
"events = DECODER.get_data(2000, events=5000, delay=603)\n",
"base_filtered = basic_noise_filter(events)\n",
"complex_filtered, channel_peak_indexes = peak_noise_filter(base_filtered, threshold=15, ret_peaks=True, peak_t=False)\n",
"plt.figure()\n",
"counter = 0\n",
"for i, event in enumerate(complex_filtered[:5]):\n",
" peak_ind = channel_peak_indexes[i]\n",
" for j, ch in enumerate(event[:3]):\n",
" x = rebin_signal(ch[...,0], 5)\n",
" y = rebin_signal(ch[...,1], 5)\n",
" # ind = peak_ind[j]\n",
" \n",
" plt.plot(x,y)\n",
" plt.plot(ch[...,0], ch[...,1])\n",
" # plt.plot(x[ind], y[ind], \"+\")\n",
" counter += 1\n",
" plt.show()\n",
"\n",
"\n",
"print(counter)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def gaussian(x, *p):\n",
" A, mu, sigma = p\n",
" return A * np.exp(-(x -mu)**2/(2*sigma**2))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"basic_noise_filter kept 2845 events out of 4999 (-> 56% kept.)\n"
]
},
{
"ename": "TypeError",
"evalue": "list indices must be integers or slices, not tuple",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_13400/4254197204.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[0mbase_filtered\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbasic_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mevents\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m \u001b[0mrebins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mrebin_signal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mch\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnbins\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mch\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mbase_filtered\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m...\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbase_filtered\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 10\u001b[0m \u001b[0mcomplex_filtered\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mchannel_peak_indexes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpeak_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrebins\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthreshold\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mret_peaks\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpeak_t\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 11\u001b[0m \u001b[0mdrift_velocity_array\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_13400/4254197204.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m 7\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[0mbase_filtered\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbasic_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mevents\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m \u001b[0mrebins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mrebin_signal\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mch\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnbins\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mch\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mbase_filtered\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m...\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbase_filtered\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 10\u001b[0m \u001b[0mcomplex_filtered\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mchannel_peak_indexes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpeak_noise_filter\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrebins\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthreshold\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m15\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mret_peaks\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpeak_t\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 11\u001b[0m \u001b[0mdrift_velocity_array\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mTypeError\u001b[0m: list indices must be integers or slices, not tuple"
]
}
],
"source": [
"anode_voltage = 2400\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",
"base_filtered = basic_noise_filter(events)\n",
"rebins = np.array([[rebin_signal(ch, nbins=3) for ch in base_filtered[...,i]] for i in range(len(base_filtered[1]))])\n",
"complex_filtered, channel_peak_indexes = peak_noise_filter(rebins, threshold=15, ret_peaks=True, peak_t=True)\n",
"drift_velocity_array = []\n",
"\n",
"nbins = 20\n",
"fig, axs = plt.subplots(4, 2, sharex=True)\n",
"subplots = axs.flat\n",
"fig.suptitle(\"Voltage = ${}$V\".format(cathode_voltage))\n",
"print(channel_peak_indexes.shape)\n",
"for i in range(channel_peak_indexes.shape[1]): # usually 8\n",
" subplt = subplots[i]\n",
" subplt.set(xlim=[0,3000])\n",
" subplt.set_title(\"Channel {}\".format(i+1))\n",
" subplt.set(ylabel='Number of events', xlabel='Time (ns)')\n",
"\n",
" hist, bin_edges = np.histogram(channel_peak_indexes[...,i], bins=nbins)\n",
" bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2\n",
" guess = [80, 1500, 40]\n",
" # guess = [300, 1800, 400]\n",
" # coeff, var_matrix = sc.optimize.curve_fit(gaussian, bin_centers, hist, p0=guess)\n",
" coeff = guess # no actual fitting\n",
"\n",
" subplt.hist(channel_peak_indexes[...,i], bins=nbins)\n",
" subplt.plot(bin_centers, gaussian(bin_centers, *coeff))\n",
" subplt.axvline(coeff[1], linestyle='dashed', linewidth=1, label=\"$\\mu$\")\n",
"\n",
" nsigma = 1.9\n",
" max_time = coeff[1] + nsigma * coeff[2]\n",
" min_time = coeff[1] - nsigma * coeff[2]\n",
"\n",
" subplt.axvline(max_time, linewidth=1, label=\"$\\mu + {} \\sigma$\".format(nsigma))\n",
" subplt.axvline(min_time, linewidth=1, label=\"$\\mu - {} \\sigma$\".format(nsigma))\n",
"\n",
" subplt.legend()\n",
"\n",
" drift_velocity = calculate_drift_velocity(min_time, max_time)\n",
" print(f\"Drift Velocity Channel {i+1}: {drift_velocity} cm * us^-1\")\n",
"\n",
" drift_velocity_array.append(drift_velocity)\n",
"plt.show()\n",
"\n",
"drift_velocity_array = np.array(drift_velocity_array[3:-1])\n",
"mean_drift_velocity = np.mean(drift_velocity_array)\n",
"std_drift_velocity = np.std(drift_velocity_array)\n",
"\n",
"electric_field = calculate_electric_field(cathode_voltage, anode_voltage)\n",
"\n",
"print(f\"Electric Field: {electric_field}\" +'\\n' +\n",
"f\"Mean Drift Velocity: {mean_drift_velocity}\" + '\\n' +\n",
"f\"Standard Deviation DV: {std_drift_velocity}\" + '\\n')"
]
}
],
"metadata": {
"interpreter": {
"hash": "63fd5069d213b44bf678585dea6b12cceca9941eaf7f819626cde1f2670de90d"
},
"kernelspec": {
"display_name": "Python 3.9.7 64-bit",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
},
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
}
%% 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
This diff is collapsed.
This diff is collapsed.
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