LABO4 - tutorial analisi dati di spettroscopia¶
Obiettivo
Acquisire le conoscenze di base per analizzare i dati con python 3 (usanto le più note librerie in ambito scientifico quali numpy, matplotlib e scipy). In particolare serve a prendere confidenza con:
- leggere diversi file di dati
- creare/manipolare istrogrammi
- eseguire fit e usarne i risultati
- scrivere/usare funzioni
- gestire la grafica
Attività di laboratorio
L'attività di laboratorio descritta in questo tutorial consiste nel calibrare la risposta dell'apparato di acquisizione dati utilizzando sorgenti con spettri di emissione noti. In questo esempio verranno utilizzati i dati acquisiti con due rivelatori a scintillazione (NaI, LrBr) e una sorgente di Rd-226 (gamma). L'informazione salvata nei dati è l'ampiezza del segnale della luce di scintillazione emessa dal rivelatore (letta da fotomoltiplicatori e amplificata di fattori noti) quando questo viene colpito dalla radiazione emessa dalla sorgente.
Lo spettro dei dati acquisiti va ricondotto a quello noto della sorgente per caratterizzare il funzionamento di tutto l'apparato (linearità della risposta del rivelatore, risoluzione energetica, dipendenza dall'alimentazione del fotomoltiplicatore).
Formato dei dati
I dati sono file di testo con una riga per ogni evento acquisito e più colonne. La prima contiene il tempo in cui il segnale è stato registrato (unix time), le successive le ampiezze registrate sui canali dello strumento di digitalizzazione del segnale (un Waveform Digitizer):
...
001480064955 249 84 66 106
001480064955 3717 85 68 106
001480064955 2498 75 69 107
001480064955 276 86 68 106
001480064955 450 94 67 106
001480064955 388 75 67 105
001480064955 277 87 67 106
...In questo caso analizzeremo solo un segnale (la seconda colonna).
Costruire il programma di analisi
Per prima cosa vanno richiamati (import) i moduli necessari. Teoricamente lo si può fare in qualsiasi punto del programma (purché precedente al primo utilizzo) ma è consigliabile farlo all'inizio per questioni di leggibilità .
Il modulo (o "libreria") viene caricato e subito gli viene assegnato un alias per rendere più facile il suo utilizzo all'interno del programma. E per aumentarne la leggibilità . Ogni funzione del modulo numpy verrà ad esempio richiamata precedendola da "np" (anziché scrivere per interno "numpy").
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt
# Questo è un "magic command" per avere i grafici all'interno del notebook.
# Non è necessario se usate un qualsiasi altro editor/IDE
%matplotlib inline
Lettura del file di dati (per comodità scriviamo il nome del file in una stringa) usando la funzione loadtxt del modulo numpy. L'opzione unpack serve a "spacchettare" l'output della lettura nei diversi array che predisponiamo prima del segno "=". Altrimenti il comando loadtxt di default scriverebbe tutto in un'unica grande matrice. In questo caso le colonne sono poche e vale la pena importarle in variabili diverse.
datafile = './datadir/run020752.dat'
time, a, b, c, d = np.loadtxt(datafile, unpack=True)
I dati che ci interessano sono ora contenuti nell'array (o "vettore") di numpy a. Tecnicamente si tratta di un "array n-dimensionale". Per comodità possiamo anche rinominarlo e stamparne il contenuto per controllo.
dati = a
print('Dati è un oggetto del tipo ', type(dati))
print('Il numero totale di elementi che contiene è ',len(dati))
print(dati)
Usando le librerie grafiche matplotlib possiamo guardare qual è il loro andamento evento per evento (cioè riga per riga) usando la funzione della matplotlib plot. Il comando è preceduto dalla stringa "plt" perché è l'alias con cui all'inizio abbiamo definito il modulo matplotlib.pyplot.
plt.plot(dati)
plt.show() # <-- mettetelo ogni volta che volete mostrare il plot appena costruito
Di default il comando plot utilizza un marker molto grande (tutti i marker), per rimpicciolirlo ad esempio si può usare un "pixel" (puntino) e rendere la colorazione parzialmente trasparente:
plt.plot(dati, marker=',',alpha=0.3)
plt.show()
Quello che interessa è comunque la distribuzione cumulativa di tutti gli eventi, cioè un istogramma:
plt.hist(dati)
plt.show()
Di default il numero di bin è basso (10), per aumentarlo:
plt.hist(dati,500)
plt.show()
Rendiamo la figura un po' più grande e la distribuzione più visibile. Aggiungiamo gli assi
plt.figure(figsize=(15,5))
plt.hist(dati,1000, fc='orange')
plt.xlabel('ADC', fontsize=14)
plt.ylabel('eventi', fontsize=14)
plt.show()
Probabilmente il numero intorno a 8000 è una saturazione che non ci interessa. Ci concentriamo sui dati che indicano un valore fino a 4000 ADC. Per farlo usiamo una funzione utile degli array di numpy che tecnicamente prende il nome di boolean/mask index of arrays. In sostanza è un comodo modo per filtrare il contenuto di un array di dati sulla base di alcune condizini. In questo caso ci interessa limitare il set di dati a quelli con un valore di ADC inferiore a 4000. Possiamo quindi creare un secondo array di dati filtrati e poi rifare il grafico:
datif = dati[(dati<4000)]
plt.figure(figsize=(15,5))
plt.hist(datif,500, fc='orange')
plt.xlabel('ADC', fontsize=14)
plt.ylabel('eventi', fontsize=14)
plt.show()
A questo punto vogliamo misurare la posizione dei diversi picchi. Per farlo possiamo cercare di fare un fit dei dati usando una generica funzione gaussiana e limitando il range del fit all'intervallo contenente il picco che ci interessa. La funzione gaussiana (normale) è presente in moltissime librerie, ma è facile da definire:
def gaus(x, A, mu, sigma):
return A*np.exp(-(x-mu)**2/(2*sigma**2))
Per quanto riguarda il fit, usiamo il comando curve_fit del modulo scipy.optimize importato inizialmente con l'alias opt. Come si legge nella documentazione la funzione prende come argomenti la curva da fittare, la funzione con cui fittarla e un array di parametri di inizializzazione. La curva da fittare è disponibile come output del plot, occorre però memorizzarne il contenuto in opportuni oggetti.
h, bins, _ = plt.hist(datif,500, fc='orange') # "_" è un nome di variabile valido e tipicamente sta ad indicare una variabile "dummy" che non si utilizzerÃ
print('h è un array di dimensione ',len(h),"e contiene tutti i valori dell'istogramma generato")
print('Stampo primi 3 elementi:')
for i in range(3):
print(h[i])
print('Stampo ultimi 3 elementi:')
for i in range(3):
print(h[len(h)-i-1])
Vediamo il contenuto di bins. Attenzione che se h ha dimensioni N bins avrà dimensione N+1. Sono gli edge dei bin dell'istrogramma, quindi sempre uno in più:
print('bins è un array di dimensione ',len(bins),"e contiene tutti i valori dell'istogramma generato")
print('Stampo primi 3 elementi:')
for i in range(3):
print(bins[i])
print('Stampo ultimi 3 elementi:')
for i in range(3):
print(bins[len(bins)-i-1])
Tornando al nostro spettro, inizializziamo i parametri e facciamo un fit se una certa regione, ad esempio quella intorno a 1400 ADC.
def gaus(x, A, mu, sigma):
return A*np.exp(-(x-mu)**2/(2*sigma**2))
# definisco sotto-insieme dei dati
picco = datif[(1300<datif) & (datif<1500)]
# histo come prima, ma centrato sul picco. Memorizzo bin e valori
h, bins, _ = plt.hist(picco,100, fc='orange')
# definizio parametri sensati per inizializzare il fit
initPars = [200, 1420, 50]
# eseguo il fit e memorizzo i risultati (parametri e matrice di covarianza degli errori)
# attenzione che "x" e "y" del fit devono avere uguali dimensioni. Quindi dell'array di bin prendiamo
# l'edge di sinistra, in totale quindi tutti i componenti tranne l'ultimo.
optimizedPars, pCov = opt.curve_fit(gaus,bins[:-1],h, p0=initPars)
# aggiungo al plot la funzione con i parametri di inizializzazione (per controllo) e dopo il fit
plt.plot(bins,gaus(bins,*initPars),':r', linewidth=3, label='inizializzazione fit')
plt.plot(bins,gaus(bins,*optimizedPars),'--k', linewidth=3, label='risultato fit')
plt.xlabel('ADC', fontsize=14)
plt.ylabel('eventi', fontsize=14)
plt.legend()
plt.show()
Il simbolo "*" che precede una lista/array tecnicamente si chiama unpacking-argument-lists.
Gli errori relativi ai parametri del fit sono ottenibili dalla radice quadrata dei valori sulla diagonale della matrice di covarianza (in questo caso 3x3):
print('La forma della matrice di covarianza è ',pCov.shape)
print('Il contenuto\n',pCov)
print('La diagonale ',np.diag(pCov))
print(100*'-')
print('Il risultato del fit con errori è quindi: ')
for i in range(3):
print('Parametro %d: %.f +- %.f'%(i, optimizedPars[i], np.sqrt(np.diag(pCov))[i] ) )
A questo punto mettiamo insieme le cose e proviamo a fare fit dei picchi più evidenti inizializzando i valori in modo opportuno. Prima conviene creare un unico istrogramma e ricavarne x e y per poi lavorare su quelle.
# Usiamo solo i dati sotto i 4000 ADC
plt.figure(figsize=(15,5))
h, bins, _ = plt.hist(dati[(dati<4000)],500, fc='orange')
# memorizziamo bin e valori più comodi (per le x prendiamo il centro dei bin)
hy = np.array(h)
hx = np.array([(bins[i+1]+bins[i])/2 for i in range(len(bins)-1)])
plt.plot(hx,hy,'-k')
plt.show()
Ora creiamo a occhio i differenti range di valori per inizializzare i fit dei vari picchi.
# definisco una lista vuota e "appendo" liste di array per i vari picchi.
# idem per i parametri di inizializzazione
range_picchi = []
initPars_picchi = []
# picco 1
range_picchi.append((260, 310))
initPars_picchi.append([5000, 300, 50])
# picco 2
range_picchi.append((480, 560))
initPars_picchi.append([1000, 550, 50])
# picco 3
range_picchi.append((600, 680))
initPars_picchi.append([1000, 650, 50])
# picco 4
range_picchi.append((720, 800))
initPars_picchi.append([1000, 750, 50])
# picco 5
range_picchi.append((820, 920))
initPars_picchi.append([2000, 900, 50])
# picco 6
range_picchi.append((1350, 1500))
initPars_picchi.append([1000,1400, 50])
initPars_picchi
plt.figure(figsize=(15,5))
h, bins, _ = plt.hist(dati[(dati<2000)],400, histtype='step', ec='k')
# memorizziamo bin e valori più comodi (per le x prendiamo il centro dei bin)
hy = np.array(h)
hx = np.array([(bins[i+1]+bins[i])/2 for i in range(len(bins)-1)])
# eseguo il fit più volte memorizzando gli output
optPars_picchi=[]
pCov_picchi=[]
for i in range(6):
# definisco la zona da fittare
hxfit = hx[ (hx> range_picchi[i][0]) & (hx<range_picchi[i][1]) ]
hyfit = hy[ (hx>range_picchi[i][0]) & (hx<range_picchi[i][1]) ]
optimizedPars, pCov = opt.curve_fit(gaus,hxfit,hyfit, p0=initPars_picchi[i])
optPars_picchi.append(optimizedPars)
pCov_picchi.append(pCov)
plt.plot(hxfit,gaus(hxfit,*optimizedPars),'-',color='C%d'%i, linewidth=3, label='%d: $\mu$ %.1f, $\sigma$ %.1f'%(i,optimizedPars[1],optimizedPars[2]))
plt.grid()
plt.legend(fontsize=16)
plt.show()
Ora possiamo raffinare i fit, magari aggiungendo il background e così che i valori dei fit gaussiani siano più affidabili. Ad esempio una retta e una gaussiana, picco per picco:
plt.figure(figsize=(20,10))
def fitpicco(x,m,q,A,mu,sigma):
return m*x+q+gaus(x,A,mu,sigma)
optPars_picchi2=[]
pCov_picchi2=[]
for i in range(6):
plt.subplot(2,3,1+i)
hxfit = hx[ (hx> range_picchi[i][0]-30) & (hx<range_picchi[i][1]+30) ]
hyfit = hy[ (hx>range_picchi[i][0]-30) & (hx<range_picchi[i][1]+30) ]
optimizedPars, pCov = opt.curve_fit(fitpicco,hxfit,hyfit, p0=[0,0]+list(initPars_picchi[i]))
optPars_picchi2.append(optimizedPars[2:]) # memorizzo solo dal 3 in poi (i primi due sono la retta)
pCov_picchi2.append(pCov[2:])
plt.plot(hx,hy,'-k', label='dati')
plt.plot(hxfit,fitpicco(hxfit,*optimizedPars),'-r', linewidth=3, label='fit ($\mu=%.1f$, $\sigma=%.1f$)'%(optimizedPars[3],np.abs(optimizedPars[4])))
plt.legend(fontsize=16)
plt.show()
print(optPars_picchi2)
Ora vediamo di fare una calibrazione (fit lineare ADC vs Energia) usando mean e sigma
picchi_KeV = [77, 185, 241, 295, 351, 609]
picchi_ADC = np.array(optPars_picchi2)[:,1]
errori_ADC = np.absolute(np.array(optPars_picchi2)[:,2])
plt.errorbar(picchi_KeV, picchi_ADC, errori_ADC, marker='o', color='k', label='Dati', linewidth=0.5)
m, q = np.polyfit(picchi_KeV, picchi_ADC, 1, w=errori_ADC)
x = np.linspace(picchi_KeV[0], picchi_KeV[-1],100)
plt.plot(x,m*x+q,'--r',label='Fit y={0:.1f} x + {1:.1f}'.format(m,q), linewidth=2)
plt.legend()
Ora cerchiamo di condensare tutto in un unico programma. Proviamo a creare un'unica funzione che dopo aver aperto un file ci restituisca i valori ottenuti dalla calibrazione.