Pitone per analizzare tracce del rivelatore bete basato su acquisizione tramite schede audio¶

Inizio a partire da quanto mi ha spedito Erik

In [1]:
import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl

from scipy.optimize import curve_fit
from scipy.signal import find_peaks

import soundfile as sf

In questa cella definisco se voglio fare il plot dei picchi trovati (per controllo) e quanti campionamenti utilizzare per il plot. Nota: i picchi vengono comunque cercati su tutto il set di dati

In [2]:
plot = True
N = 1000000

Carico il file di dati, disponibili qui. Il questo quaderno si può invece trovare qui

In [3]:
%%time
path = 'co60.flac'                                                  
data, samplerate = sf.read(path)   
print(data.shape)

(255777024, 2)
CPU times: total: 9.83 s
Wall time: 16.8 s

Ottengo il vettore relativo ai dati. Non ho letto la documentazione in prima persona, ma mi fido che i dati siano nella seconda colonna

In [4]:
%%time
DATI = data[:,1]
DATI.shape
dati = np.abs(DATI)

#DATI = DATI[:1000000]

CPU times: total: 11.8 s
Wall time: 49.7 s

Inizio a plottare la wf, così per guardarla in faccia

In [5]:
# %matplotlib qt
%matplotlib inline


fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(DATI[:N], c = "tab:green")

plt.show()

Visto che a noi interessa l'ampiezza, immagino che ne vada considerato il valore assoluto

In [6]:
# qt per avere la finestra grafica
%matplotlib inline 

dati = np.abs(DATI) # Valore assoluto

fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(dati[:N], c = "tab:green")
ax.axhline(y = .25, ls = "--", c = "tab:orange")

plt.show()

Inizio del mio lavoro¶

Noto che mettendo una soglia a 0.25 i picchi risultano ben al di sopra della baseline...

Ricerca dei picchi¶

Senza volermi troppo male, uso find peaks. Per un riferimento con dei disegnini chiari, si veda la documentazione di matlab per la medesima funzione. Si tratta di una funzione molto comoda che permette di individuare picchi in un vettore, impostando se necessario vari parametri.

Vado dunque ad effettuare la ricerca dei picchi, e li disegno sopra alla wf originale, per assocurarmi della bontà del lavoro

In [7]:
%%time
%matplotlib inline


# Settings - eventualmente provare a giocarci un po'
threshold = .1 # <-- Set threshold here
minWidth = 3 # <-- Set minimum width of a single peak
minDist = 50 # <-- Set minimum distance between two consecutive peaks

peaks, _ = find_peaks(dati, height = threshold, width = minWidth, distance = minDist) # <-- x dei picchi
print(f"I have found {peaks.shape} peaks")

PHvect = dati[peaks] # <-- y dei picchi = Pulse Height


# Plot
if plot:
    fig, ax = plt.subplots()
    fig.set_size_inches(12,5)

    ax.plot(dati[:N], c = "tab:green", label = "Original abs wf")
    ax.axhline(y = threshold, ls = "--", c = "tab:orange", label = "Threshold")
    ax.plot(peaks[peaks<N], PHvect[peaks<N], ls = "", marker = "o", color = "tab:red", label = "Found peaks")

    ax.set_xlabel("Number of sample", fontsize = 14)
    ax.set_ylabel("Ampiezza", fontsize = 14)
    ax.set_title("", fontsize = 16)

    #ax.grid()
    ax.legend(loc = "upper right")

    plt.show()

I have found (1548,) peaks

CPU times: total: 11.7 s
Wall time: 21.8 s

Istogramma delle PH¶

Vado ad effettuare e plottare l'istogramma delle Pulse Height. Nota: sembra che ci sia un po' di saturazione

In [8]:
h, bins = np.histogram(PHvect, bins = 30)
binc = bins[:-1] + (bins[1]-bins[0])/2


fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Spectrum")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)

#ax.plot(binc, myFun(binc, *popt))

ax.grid()
ax.set_xlabel("Energy", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("232-Th Spectrum", fontsize = 16)
ax.legend()
#ax.set_yscale("log")

plt.show()
In [9]:
dt = np.diff(peaks)

h, bins = np.histogram(dt, bins = 30)
binc = bins[:-1] + (bins[1]-bins[0])/2


fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Tempi")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)

#ax.plot(binc, myFun(binc, *popt))

ax.grid()
ax.set_xlabel("$\Delta T$", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("Tempi di interarrivo", fontsize = 16)
ax.legend()
#ax.set_yscale("log")

plt.show()
In [ ]: