Tutorial analisi spettri e correlazioni

Introduzione

Questo quaderno si inserisce nel materiale didattico che verrà fornito agli studenti del secondo modulo del corso di laboratorio di fisica 3A (per gli amici Labo 4).

Abbiamo deciso di iniziare analizzando lo spettro del $^{226}$Ra,che è uno degli esempi più delicati e completi con cui avrete a che fare.

Nel caso durante il corso nascessero esigenze particolari riguardo a qualche dubbio nell'effettuare una certa procedura, oppure vi trovaste a dover lavorare con dati in formati strani, esattamente come durante lo scorso semestre, potete contattarmi e d'accordo con gli altri tutor e con i docenti, vedrò di venirvi incontro nel miglior modo possibile. La mia mail è scarsi@studenti.uninsubria.it

Scopo

In questo quaderno verranno presentati alcuni spunti necessari per poter costruire istogrammi relativi a spettri, fittarne i picchi e costruire istogrammi bidimensionali per studiare le correlazioni.

Il dataset utilizzato uno di quelli relativi all'analisi della spettroscopia gamma, ma le medesime tecniche si applicano ad un generico contesto in cui si hanno dei rivelatori da cui, sotto certe condizioni di trigger, viene estratto il massimo della waveform e si deve valutare la distribuzione dei massimi (come per esempio nel contesto della spettroscopia alpha).

Strumenti

L'analisi viene effettuata utilizzando Python 3 ed in particolare moduli numpy, matplotlib e scipy per la gestione e la visualizzazione di dati numerici.

Setup sperimentale e contesto

Il setup sperimentale è costituito da due tubi fotomoltiplicatori (PMT) connessi ad una catena elettronica di amplificazione e shaping. Ogni qualvolta almeno uno dei due PMT è sopra una soglia prestabilita, vengono acquisiti i massimi delle forme d'onda di ciascun canale.

Un'altra possibilità su cui siamo al lavoro proprio di questi tempo è di acquisire direttamente il segnale proveniente dai PMT ed effettuare la parte di amplificazione, shaping e ricerca del massimo via software. Nella parte superiore si vede il segnale originario relativo ad un LaBr$_3$, mentre nella parte inferiore il segnale shapeato, su cui andare a ricercare il massimo.

Formato dei dati

I dati sono dei files di testo (in formato ASCII): ogni riga corrisponde ad un evento (ovvero ogni volta che almeno uno dei due PMT ha avuto valore sopra soglia), mentre ogni colonna contiene una certa quantità. Aprendo un file, esso ha la forma

...
001480064955 591 87 66 106 
001480064955 130 6 4095 69 
001480064955 235 92 66 106 
001480064955 304 90 66 106 
001480064955 858 87 66 107 
001480064955 257 90 66 106 
001480064955 249 84 66 106 
001480064955 3717 85 68 106 
001480064955 2498 75 69 107 
...

dove la prima colonna rappresenta il tempo unix) (ovvero il numero di secondi trascorsi dalla mezzanotte del 1° Gennaio 1970) mentre le successive quattro colonne rappresentano il massimo del segnale per ciascun canale. I due PMT sono connessi ai primi due canali. Poiché l'analisi è perfettamente equivalente, noi ci soffermeremo solamente sul primo canale.

Analisi dei dati

Iniziamo caricando (import) i moduli che verranno utulizzati durante l'analisi e definendo un alias per comodità nel richiamarli.

Load dei dati

Carichiamo i dati utilizzando il metodo np.loadtxt. Il carattere "_", pur essendo un nome lecito per una variabile, viene convenzionalmente utilizzato per indicare una variabile di nessuna importanza. Infatti per caricare una matrice con N colonne ed utilizzare l'opzione unpack=True (che permette di memorizzare ciascuna colonna in un diverso vettore) è necessario passare tante variabili quante sono le colonne, indipendentemente dal fatto che alcune non sono interessanti.

Possiamo subito vedere che forma abbiano i dati caricati

Visualizzazione dei dati

La prima cosa da fare quando si hanno dei dati è sicuramente guardarli. Uno potrebbe essere tentato di plottarli direttamente

Ma purtroppo non si vede un granchè. Possiamo provare a plottare singoli punti non uniti da linee, eventualmente con marker più picolo

Il problema risiede nel fatto che questi dati non sono da guardare in questo modo

Costruzione degli spettri

A noi infatti interessa studiare la distribuzione dei massimi, sapere in un certo intervallo di valori quanti eventi sono stati registrati: dobbiamo costruire un istogramma

Matplotlib mette a disposizione il metodo hist che permette di effettuare l'istogramma di un vettore

Da questo istogramma si vede poco e nulla. Proviamo ad aumentare il numero di bin

C'è ancora un po' da lavorare per ottenere un bello spettro, anche se qualcosa già si inizia ad intravedere. Vi segnalo che anche numpy consente di effettuare istogrammi. Personalmente preferisco mantenere distinta la parte in cui vado a costruire l'istogramma e la parte in cui vado a plottare l'istogramma costruito, pertanto di seguito tutti gli istogrammi saranno creati con il metodo np.histogram.

Notiamo come i due metodi (matplotlib e numpy) accettino gli stessi argomenti, ad eccezione dei parametri relativi alla grafica. Infatti il metodo di matplotlib va a richiamare proprio quello di numpy.

Istogrammi 1D e spettri

Il metodo np.histogram accetta come argomento obbligatorio il dataset da istogrammare. Possiamo poi specificargli alcuni parametri opzionali, come per esempio il numero di bin o il range. Gli argomenti che ritorna invece sono il numero di eventi per ciascun bin e il bordo dei bin.

Il secondo vettore ha un elemento in più del primo e, per ottenere le coordinate $xy$ da plottare, è necessario sfasare tutti i punti, ad eccezione dell'ultimo, di mezzo bin.

Quando devo disegnare più curve in uno stesso plot con impostazioni analoghe, per comodità definisco un dizionario che le contiene e le passo alla funzione plt.plot effettuando l'unpacking dei dizionari mediante il doppio asterisco.

Per poter disegnare un istogramma con il metodo plt.plot è necessario specificare l'opzione drawstyle (o la sua abbreviazione ds) a steps-mid, che indica di tracciare il grafico a gradini e che il gradino debba essere centrato rispetto alla corrispondente ascissa.

Da questo istogramma si vede poco e nulla. Prima di tutto, poiché l'ADC, ovvero lo strumento che campiona la forma d'onda ha una risoluzione a 12 bit, il maggior numero che ci può restituire è $2^{12}-1=4095$, ed al di sopra potrebbe esserci una saturazione. Abbiamo anche provato a mettere l'asse $y$ in scala logaritmica (ax.set_yscale("log"))

Una riflessione finale: Abbiamo fatto il taglio sull'istogramma già costruito, in poche parole visualizziamo solamente una porzione dell'istogramma costruito. Poiché ci sono fondati motivi teorici per fissare un limite massimo, sarebbe buona norma mettere il range sull'istogramma da costruire. In questo esempio ho aggirato il problema facendo mooolti bin

Versione corretta

Et voilà! Ecco il nostro spettro!

Che cosa mostra questo plot? Questo plot ci mostra la distribuzione del massimo misurato da un PMT, ovvero ci dice in ogni intervallo di energie (misurate in unità arbitrarie di ADC) quanti eventi sono stati registrati. L'area sottesa a questo grafico è invece il numero totale di eventi utilizzati.

Da questo plot si vedono emergere alcuni picchi, che dovrebbero corrispondere ai fotoni che vengono emessi.

Fit dei picchi

Una volta costruito l'istogramma, vogliamo effettuare il fit dei picchi: in particolare vogliamo trovare la $x$ del massimo e il corrispondente errore: effettuiamo quindi il fit con una funzione gaussiana.

Notare in particolare come i picchi non siano centrati sullo 0 ma abbiano un offset verticale: per tenerne conto possiamo sommare un polinomio di primo grado oppure un esponenziale decrescente. Se il passaggio precedente è stato effettuato correttamente può sembrare che siano montati su una retta decrescente: attenzione a non farsi ingannare dalle impostazioni degli assi.

I picchi da fittare sono i primi 7, che hanno energie 77, 185, 241, 295, 351, 609 e 1120 keV. (Il secondo picco nell'immagine sottostante compare a 186/187, nel logbook fornito era dato a 185. In ogni caso non fa tutta sta gran differenza...)

Si tratta di fit un po' delicati da far convergere, serve un po' di pazienza: per aiutarsi a stimare gli starting points può essere utile ricorrere al seguente magic command %matplotlib qt che consente di aprire il plot in una finestra grafica esterna, e dove muovendosi col puntatore si possono visualizzare le coordinate. Per visualizzare i plot nel modo usuale, esiste il magic command %matplotlib inline.

Per comodità ho definito una funzione che mi aiuta a mantenermi compatto il codice, in modo tale che se facessi anche sono una picoola modifica, non debba propagarla per 7 volte, una per ciascun picco.

Notiamo in particolare come il parametro absolute_sigma = True permetta di utilizzare nel fit gli errori veri e non degli errori riscalati di un certo fattore tale per cui il $\bar\chi^2\approx 1$

Qualcuno potrebbe pensare che io sia un pazzo ad aver stimato a mano con precisione fino alla n cifra decimale gli starting points. In realtà mi sono servito di uno strumento molto comodo che ho leggermente customizzato e che permette in modo visuale di aggiustare a mano i parametri di fit. In particolare, rispetto alla versione disponibile online, ho aggiunto la possibilità di passare errori al fit (cosa di vitale importanza per un fisico). Appena avrò modo, metterò un link sul mio sito.

Per poter essere dinamico è necessario aprirlo in una finestra grafica esterna (%matplotlib qt), mentre in linea (%matplotlib inline) viene stampata solo un'immagine statica.

Ottenimento dei valori estratti dal fit

Dopo aver effettuato il fit è importante ottenere i parametri risultanti dal fit, in particolare ci interessa ottenere $\mu$ e $\sigma$ delle gaussiane con i relativi errori. I parametri si trovano direttamente in popt, mentre pcov costituisce la matrice delle covarianze. Sulla diagonale si trovano le varianze, ovvero il quadrato della deviazione standard.

Inoltre, essendo la gaussiana $$\dfrac{1}{\sigma\,\sqrt{2\pi}}\,e^{-\dfrac{(x-\mu)^2}{2\sigma^2}}$$ normalizzata, ovvero di area unitaria, il prefattore $a$ di $$\dfrac{a}{\sigma\,\sqrt{2\pi}}\,e^{-\dfrac{(x-\mu)^2}{2\sigma^2}}$$ rappresenta l'area della gaussiana.

Notiamo che nella funzione gaussiana, la "forma" viene data dal secondo termine (l'esponenziale) mentre il prefattore prima è semplicemente una costante. Talvolta (per comodità o per esempio nella parametrizzazione di PAW) la gaussiana viene indicata come $$A\,e^{-\dfrac{(x-\mu)^2}{2\sigma^2}}$$ tale formulazione è assolutamente lecita, in poche parole si è definito $A=\dfrac{a}{\sigma\,\sqrt{2\pi}}$. In questo caso bisogna però stare attenti che il parametro $A$ non rappresenta direttamente l'area.

Siccome io sono pigro, questo è un modo di farsi aiutare dal pitone a formattare le tabelle $\LaTeX $, cosicchè sia sufficiente copiare questo output ed incollarlo nella relazione, senza dover inserire ogni numero a mano.

Costruzione della retta di calibrazione

Una volta che sono stati individuati i picchi con il fit gaussiano, è possibile convertire gli spettri da unità arbitrarie di ADC ad energie in keV. Sapendo l'energia vera dei picchi (vedi sopra) ed avendo stabilito i picchi che si vedono quali energie dovrebbero avere, è possibile plottare l'energia in unità di ADC (la $x$ del picco, la $\mu$ estratta dal fit) in funzione dell'energia vera. Questa retta costituisce la retta di calibrazione (più precisamente l'inverso della retta di calibrazione, ovvero quella funzione che accetta come argomento un'energia in ADC e restituisce l'energia in keV. Per stabilire chi mettere sulle $x$ e chi sulle $y$ chiediamoci chi ha errore e chi assumiamo come esatto.

Inizio quindi definendo la funzione di fit. È vero che per una retta esiste la funzione polifit, ma ci sono state delle volte in cui mi ha dato problemi, e comunque non si possono specificare gli errori.

Successivamente definisco dei vettori contenenti i valori veri dei picchi, i valori estratti dal fit e gli errori.

Procedo quindi effettuando il fit ed andando ad effettuare il plot

Plot dell'istogramma in energia

A questo punto per plottare l'istogramma in energia non ci resta che copiare il codice usato in precedenza, applicando una certa trasformazione alle $x$ per ottenere lo spettro della nostra sorgente di $^{226}$Ra.

Volendo si può colorare tutta l'area al di sotto dello spettro, in modo tale da ottenere un grafico a barre. Si tratta di una semplice operazione estetica. Il metodo plt.fill_between permette di riempire l'area tra due curve, nel nostro caso tra lo 0 e l'istogramma tracciato.

Pro tip

Durante l'attività di spettroscopia gamma acquisirete un file per ogni sorgente: piuttosto che aprire ciascun file salvandosi i dati in differenti variabili, può essere comodo automatizzare la procedura. Per ogni evento vengono salvati due differenti files: uno ".dat" che contiene i dati e uno ".log" che contiene alcune informazioni sull'acquisizione.

In una certa cartella ho due sottocartelle:

# Questa lista contiene le varie matrici di dati, ogni matrice è un file
listaMatrici = []
# Questa lista invece contiene 
listaNum = []

# Ciclo su tutti i files nella cartella dati, indietro di un livello ("..")
for e in os.scandir(r"..\DATI"):
    if e.name.split(".")[1] != "dat": continue

    tmp = np.loadtxt(e.path)
    listaMatrici.append(tmp.copy())
    listaNum.append(int(e.name.split(".")[0][-2:]))

zipIt = zip(listaNum, range(len(listaNum)))
# Questo dizionario accetta come input il numero di run (in realtà gli ultimi due numeri) e dà come output l'indice progressivo con cui è stato caricato. Non è automaticamente garantito che i dati vengano caricati in ordine alfabetico
myDict = dict(zipIt)

Attenzione! È fodnamentale appendere tmp.copy() e non semplicemente tmp: nel primo caso si appende un nuovo oggetto per ogni file, mentre nel secondo caso si appende sempre lo stesso oggetto (e quindi tutte le matrici sarebbero uguali all'ultimo file aperto).

Per effettuare la procedura di analisi e fit, anziché costruire manualmente plot ed istogrammi per ciascuna sorgente, ho implementato una funzione la cui firma è

def makeFit(numRun, sogliaMin, sogliaMax, startPars, rangeFit, valVero, sorgente,\
            fig, ax, primoPicco = True):

che può essere un utile spunto. Ognuno di voi sceglierà comunque la strategia che più ritiene opportuna. In particolare l'opzione primoPicco serve a tener conto di sorgenti che abbiano più picchi (come il caso del $^{226}$Ra che abbiamo visto in questa esercitazione.

Istogrammi bidimensionali

I biplot o istogrammi bidimensionali possono essere costruiti facilmente sia usando matplotlib che numpy utilizzando funzioni analoghe. In questo caso utilizziamo il metodo plt.hist2d, che accetta come argomenti due vettori, oltre a parametri analoghi al caso monodimensionale (numero di bin, range, ...). Inoltre specifichiamo la mappa di colore che intendiamo utilizzare.

Il biplot effettua un binnaggio lungo i due assi ed in ogni quadratino viene riportato il numero delle volte che il primo canale rientrava in quel bin lungo x ed il secondo canale il quel bin lungo y, ovvero il numero di volte che il primo PMT ha visto una certa energia ed il secondo PMT un'altra certa energia. Il numero di volte viene codificato con un colore.

Questo tipo di plot sono molto potenti in quanto permettono di studiare le correlazioni tra i due canali. Nel nostro caso abbiamo due fotoni che possono venire emessi in coincidenza (a distanza di pochissimo tempo) e venir rivelati uno per ciascun PMT: questo giustifica gli addensamenti che si notano.

Iniziamo quindi a caricare i dati relativi al $^{22}$Na. In questo caso è necessario utilizzare i dati di entrambi i PMT.

Costruisco quindi l'istogramma bidimensionale

Spesso può essere utile consultare gli istogrammi in scala logaritmica per accentuare piccole differenze. Per farlo è sufficientemente specificare tra gli argomenti norm = mpl.colors.LogNorm()

Notiamo che tale grafico ha qualcosa di visualmente fastidioso, in particolare introduce una discontinuità, alternando bianco e blu. Ciò è legato al fatto che dove si ha il bianco si hanno 0 conteggi, il cui logaritmo non è definito. Poiché i nostri dati sono sempre riferiti a conteggi, quindi sono 0, 1, 2.... può essere utile ridefinire a 0 il caso in cui vi siano stati 0 conteggi, mantenendo così l'uniformità. Il modulo copy ci permette di copiare un oggetto in modo completo, in modo tale che noi possiamo modificare una versione clonata di un oggetto di libreria.

Il significato di questo plot vi sarà sicuramente più chiaro dalla teoria vista a lezione e dall'esperienza relativa in laboratorio. Senza la pretesa di essere esaustivo, dalla mia relazione:

Si possono distinguere due bande verticali/orizzontali in corrispondenza dei picchi fotoelettrici a 511 keV e 1274 keV: il primo corrisponde al picco di annichilazione, mentre il secondo al fotone emesso dal decadimento della sorgente. Non si nota un significativo addensamento in prossimità dell’intersezione a 511 keV in quanto la sorgente era posta davanti ai due rivelatori, mentre i due fotoni di annichilazione vengono emessi in direzioni opposte: non risulta quindi possibile rivelare nei due rivelatori i due fotoni emessi dall’annichilazione, e quando entrambi i rivelatori hanno rivelato un fotone da 511 keV non può che essere un evento casuale, dato dalla rivelazione di un fotone ciascuno proveniente da due annichilazioni differenti. Le due bande diagonali rappresentano la situazione in cui un fotone all’interno di un rivelatore ceda parte della sua energia per effetto Compton e il fotone Compton, di bassa energia, ceda la rimanente energia all’altro rivelatore: in questo caso tutta l’energia è stata ceduta, ma ai due rivelatori vengono di volta in volta ceduti differenti valori. Sotto la diagonale a 511 keV si può ben distinguere il Compton continuum. Si osservi in particolare come le bande diagonali non raggiungano mai gli assi: questo si verifica poiché il massimo trasferimento di energia mediante effetto Compton non coincide col picco fotoelettrico, ed infatti negli spettri la spalla Compton si trova sempre prima del corrispondente picco fotoelettrico.

Talvolta può essere comodo (soprattutto per sovrapporre due istogrammi bidimensionali, nel cui caso servirebbero "4 coordinate") costruire un "boxplot", ovvero un istogramma 2D, dove l'informazione sul numero di conteggi non è codificata nel colore ma nella larghezza del marker. In questo senso, utilizzando differenti colori si possono sovrapporre più istogrammi bidimensionali.

Un'altra cosa utile può essere di tracciare delle curve di livello, ovvero delle curve che uniscono punti alla stessa quota. Anche questa operazione può essere fatta comodamente sfruttando plt.contour, oppure plt.contourf per colorare tra le curve di livello.

Appendice A: l'importanza di guardare i dati

Questi dati hanno la particolarità di avere alcuni indici statistici uguali fino alla seconda cifra decimale (media x/y, std x/y e indice di correlazione di Pearson). Questo semplice esempio vuole mettere in luce l'importanza di guardare i dati, non limitarsi a farseli riassumere da qualche indice statistico

Appendice B: il problema della funzione spigolosa

Talvolta, per esempio quando si effettua il fit gaussiano di un picco, può accadere di avere pochi punti a disposizione. Ricalcando quanto descritto nel tutorial del prof. Mascagna, generiamo un esempio di dati e vediamo come approcciare il problema

Proviamo ora ad eseguire il fit con la gaussiana

Per risolvere il problema creiamo un vettore "più denso", che abbia i medesimi estremi, ma con più punti in mezzo, in modo tale che la funzione venga disegnata più smussata. Tale problema in linea teorica l'abbiamo già incontrato con la retta di calibrazione ma, trattandosi di una retta, l'avere pochi punti in mezzo non ha alcuna rilevanza: bastano infatti due punti per definire una retta nello spazio.