Commit mes

analysis.ipynb

@@ -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