Tutorial analisi DEVA

Scopo di questo quaderno è illustrare come aprire i dati relativi a DEVA raccolti al PS del CERN (giugno 2022) e alla BTF di Frascati (aprile 2023).

DISCLAMER: Quando carico i dati, i miei percorsi dei file iniziano con /eos/project...., chiaramente quando si fa l'analisi sul vostro PC bisogna mettere il corretto percorso del file.

I dati sono disponibili in due formati

• HDF5, un formato matriciale, per il quale può essere necessario installare l'apposita libreria conda install -c anaconda h5py
• ASCII, in questo caso ad una data run corrispondono una serie di files, uno per ciascuna spill (CERN) o uno ogni 500 eventi (BTF). Si possono caricare con una funzione tipo questa

import numpy as np
import glob

def loadAsciiSpill(numRun):
    lst = glob.glob(rf".\ascii_seldom\*{numRun}*.dat")

    print(f"loading {lst[0]}")
    dataMatrix = np.loadtxt(lst[0])

    for i, file in enumerate(lst):
        if i==0: continue

        print(f"loading {file}")

        # Leggo il file
        tmpMat = np.loadtxt(file)

        # Lo appendo alla matriciona
        dataMatrix = np.vstack((dataMatrix.copy(), tmpMat.copy()))


    return dataMatrix

Scopo dell'analisi

In tutti i casi, almeno come primo step, suggerisco di iniziare con lo studio del fascio (profilo, divergenza, verificare che vi sia correlazione tra le x e le y dei canali dei silici...). Successivamente si può passare allo studio delle distribuzioni delle pulse height e dei tempi, quindi si possono ricostruire le mappe dei punti di impatto sulla faccia di DEVA. Consultando il logbook, in base a come è posizionato il calorimetro si dovrebbe studiare cosa succede, verificare che si riescono a vedere le varie tiles di scintillatore. Infine si potrebbe studiare il cosiddetto DEVA effect, quel fenomeno per cui una particella da 2 GeV o due particelle da 1 GeV (o qualsiasi altra combinazione) non rilasciano la stessa quantità di energia, per via del non contenimento... Lo puoi vedere quando hai per esempio un fascio da X energia chiedendo di avere un solo cluster, o quando hai un fascio da X/2 energia, ma chiedendo di avere due cluster

CERN 2022

Il logbook di riferimento è disponibile al seguente link, se non avete accesso, scrivetemi. Le run relative a DEVA iniziano in particolare a pagina 6 (Cercare "Added Deva to the setup" nel logbook), dove subito dopo sono illustrati anche i canali come sono collegati nei digitizer. Ricordo che sia le ampiezze che i tempi sono espressi in unità arbitrarie. Entrambi i digitizer hanno una risoluzione a 14 bit (range $0..2^{14}$) mentre i tempi sono espressi in unità di Tick: il 5730 ha una frequenza di campionamento di 500 MHz, mentre il 1724 di 100 MHz.

Subito dopo, nel logbook si trova uno schema di come è fatta DEVA (si tratta di un calorimetro a campionamento, cosituito da un'alternanza di piombo e scintillatore plastico) ed una possibile tesi di riferimento.

import sys, os, re, glob
import time
import numpy as np
import h5py
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit

Come descritto nel logbook, i dati sono disponibili a questo indirizzo. Consultando il logbook mi sono reso conto che non esistono gli HDF5, ma solo i numpy o i files ASCII, si vede che ho iniziato a farli dal test beam successivo. I file compact sono in particolare un unico ASCII con tutti i vari ASCII appesi uno sotto l'altro

dat = np.load("/eos/project/i/insulab-como/testBeam/giugno2022/npz/run501078.npz")["data"]
dat.shape

(9322, 72)

dat

array([[5.082 , 5.9246, 4.8279, ..., 0.    , 0.    , 0.    ],
       [4.6222, 3.8962, 4.3318, ..., 0.    , 0.    , 0.    ],
       [4.8525, 4.66  , 4.5852, ..., 0.    , 0.    , 0.    ],
       ...,
       [4.961 , 0.242 , 4.598 , ..., 0.    , 0.    , 0.    ],
       [8.9899, 1.1494, 8.6575, ..., 0.    , 0.    , 0.    ],
       [6.2539, 8.8939, 5.8447, ..., 0.    , 0.    , 0.    ]])

I dati in questo caso sono delle tabelle dove ogni riga corrisponde ad un evento, mentre ogni colonna ad una diversa grandezza fisica, descritta nel logbook (edit: nel logbook non sono riuscito a trovarlo, ma le posizioni e le quantità del digitizer, ph e tempi, sono corrette e nella posizione corretta). In particolare

• 4 xpos
• 4+4 (Credo num strip e num cluster, ma va verificato, comunque non importantissimo per ora)
• 16 baseline
• 16 ph
• 16 tempi di picco
• 5 xgonio (da ignorare, non c'era nessun goniometro)
• anche gli ultimi 7 non riesco a trovare cosa siano, probabilmente c'è il tempo, il tempo assoluto, il numero di evento, il numero di step di scansione del goniometro. Almeno per ora lasciamo perdere

# Una posizione negativa significa che c'è stato un errore, e pertanto deve essere saltata
(dat[:,:4]<0).sum()

9

xpos = dat[:,:4]

print(xpos.min())
condToExclude = np.all(((xpos > 0) & (xpos < 10)), axis = 1)
print(condToExclude.shape)
xpos = xpos[condToExclude]

# Creo 4 vettori con le coordinate
y1 = xpos[:,0]
x1 = xpos[:,1]
y2 = xpos[:,2]
x2 = xpos[:,3]

# 12 posizioni + numcluster + numstrip
# 16 base 16 ph 16 tempi
baseline = dat[:, 12:12+16][condToExclude]
ph = dat[:, 12+16:12+16+16][condToExclude]
time = dat[:, 12+16+16:12+16+16+16][condToExclude]

print(xpos.shape)

-4000.0
(9322,)
(9313, 4)

Profilo del fascio

fig, ax = plt.subplots(2,2)
fig.set_size_inches(12,12)
fig.subplots_adjust(hspace = .2)

fig.suptitle("Profilo del fascio", fontsize = 16)

ax = ax.flatten()

for i in range(len(ax)):
    
    # Faccio l'istogramma
    h, bins = np.histogram(xpos[:,i], bins = 384, range = (0, 0.0242*384)) # Passo strip 242 um - 384 bin perchè ho 384 strip
    
    # Ottengo il centro del bin
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    # Lo plotto
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green", lw = .7,)
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Pos [cm]",fontsize = 14)
    
    

plt.show()

A questo punto come già visto per la tua tesi, puoi vedere se c'è correlazione tra le coordinate e calcolare la divergenza del fascio. Link degli esempi che avevo già realizzato 1 e 2

Distribuzione delle PH

ph.shape

(9313, 16)

fig, ax = plt.subplots(4,4)
fig.set_size_inches(40, 40)
fig.subplots_adjust(hspace = .2)


ax = ax.flatten()

for i in range(len(ax)):
    
    # Faccio l'istogramma
    h, bins = np.histogram(ph[:,i], bins = 100,)
    
    # Ottengo il centro del bin
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    # Lo plotto
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green", )
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Energy [ADC]",fontsize = 14)
    
    ax[i].set_yscale("log")
    
    

plt.show()

Distribuzione dei tempi

fig, ax = plt.subplots(4,4)
fig.set_size_inches(40, 40)
fig.subplots_adjust(hspace = .2)


ax = ax.flatten()

for i in range(len(ax)):
    
    # Faccio l'istogramma
    h, bins = np.histogram(time[:,i], bins = 100,)
    
    # Ottengo il centro del bin
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    # Lo plotto
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green", )
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Time [Tick]",fontsize = 14)
    
    

plt.show()

Baseline

fig, ax = plt.subplots(4,4)
fig.set_size_inches(40, 40)
fig.subplots_adjust(hspace = .2)


ax = ax.flatten()

for i in range(len(ax)):
    
    # Faccio l'istogramma
    h, bins = np.histogram(baseline[:,i], bins = 100,)
    
    # Ottengo il centro del bin
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    # Lo plotto
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green", )
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Energy [ADC]",fontsize = 14)
    
    ax[i].set_yscale("log")

    

plt.show()

Ho fatto i tre plot precedenti per verificare che effettivamente le cose fossero come mi aspetto. La baseline serve solamente perchè la pulse height è calcolata come differenza tra il massimo del picco rispetto alla baseline. Sui grafici di PH e tempo si possono mettere i tagli, esattamente come descritto negli altri esempi. Si possono anche costruire gli istogrammi di correlazione PH/Tempo, dove è più naturale valutare i tagli

fig, ax = plt.subplots(4,4)
fig.set_size_inches(40, 40)
fig.subplots_adjust(hspace = .2)


ax = ax.flatten()

for i in range(len(ax)):
    
    
    # Lo plotto
    hh = ax[i].hist2d(ph[:,i], time[:,1], cmap = "jet", bins = 100)
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("PH [ADC]", fontsize = 14)    
    ax[i].set_xlabel("Time [Tick]",fontsize = 14)
    
    fig.colorbar(hh[3], ax = ax[i])
    
    

plt.show()

Trattandosi di un calorimetro, ed avendo fatto misure a più energie, è anche possibile valutare la linearità. Inserendo le distanze rilevanti e proiettando le tracce, si possono costruire le mappe di efficienza, e in base alle situazioni in cui è installato il calorimetro, vederne un bordo, vedere varie tiles di scintillatore o altro. Basta dare spazio alla fantastia !!

FRASCATi 2023

Fortunatamente in questo caso i dati sono impacchettati in un modo più smart. Il logbook di riferimento è disponibile a questo link. I dati si trovano invece qui. Durante questo test beam è stato provato il trigger su un sottoinsieme di strip dei silici, pertanto il profilo del fascio dovrebbe corrispondere a quello su cui abbiamo triggerato, sempre descritto nel logbook. La parte relativa al test beam inizia intorno a pagina 7

In questo caso i dati sono impacchettati in degli HDF5, che possono risultare più comodi da aprire, ma come visto prima, non è assolutamente fondamentale

filepath = "/eos/project/i/insulab-como/testBeam/TB_BTF_2023_04/HDF5"
numRun = 400318

fileToLoad = os.path.join(filepath, f"run{numRun}.h5")
print(f"Caricherò il file {fileToLoad}")

Caricherò il file /eos/project/i/insulab-como/testBeam/TB_BTF_2023_04/HDF5/run400318.h5

with h5py.File(fileToLoad, 'r', libver='latest', swmr=True) as hf:
    
    # Mi faccio stampare tutte le chiavi contenute nel file di dati
    print(hf.keys())
    
    # Estraggo una certa matrice, in questo caso le coordinate
    xpos = np.array(hf["xpos"])
    ph = np.array(hf["digi_ph"])
    time = np.array(hf["digi_time"])

# Verifico che abbia 4 colonne
print(xpos.shape)
    
    
# Escludo gli eventi difettosi
"""
Un valore negativo nelle posizioni è un flag per indicare anomalie, 
In ogni caso, tali eventi vanno esclusi

I telescopi hanno un'area attiva di 1.92 cm2, quindi seleziono solo quelli con coordinate
comprese tra 0 e 2 (la seconda condizione dovrebbe essere superflua)
"""

print(xpos.min())
# In questo caso metto <2 perchè sono telescopi
condToExclude = np.all(((xpos > 0) & (xpos <2)), axis = 1)
print(condToExclude.shape)
xpos = xpos[condToExclude]
ph = ph[condToExclude]
time = time[condToExclude]


# Creo 4 vettori con le coordinate
y1 = xpos[:,0]
x1 = xpos[:,1]
y2 = xpos[:,2]
x2 = xpos[:,3]

print(xpos.shape)

<KeysViewHDF5 ['digi_ph', 'digi_time', 'iatime', 'ievent', 'itime', 'xpos']>
(48251, 4)
-4000.0
(48251,)
(44577, 4)

Vedo che ci sono una serie di matrici. In particolare a noi interessa xpos (coordiante dei silici), digi_ph, pulse height e digi_time, tempi di picco

Profilo del fascio

fig, ax = plt.subplots(2,2)
fig.set_size_inches(12,12)
fig.subplots_adjust(hspace = .2)

fig.suptitle("Profilo del fascio", fontsize = 16)

ax = ax.flatten()

lstCoord = ["y1", "x1", "y2", "x2",]

for i in range(4):
    
    # Faccio l'istogramma
    h, bins = np.histogram(xpos[:,i], bins = 384, range = (0, 0.005*384)) # Passo strip 50 um
    
    # Ottengo il centro del bin
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    # Lo plotto
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green",
               label = f"{lstCoord[i]}")
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Pos [cm]",fontsize = 14)
    
    ax[i].legend(fontsize = 14)
    
    
# Zona trigger selezionata
ax[0].axvline(x = 0.005 * 100, ls = ":", c = "k", alpha = .5)
ax[0].axvline(x = 0.005 * 190, ls = ":", c = "k", alpha = .5)
ax[2].axvline(x = 0.005 * 120, ls = ":", c = "k", alpha = .5)
ax[2].axvline(x = 0.005 * 200, ls = ":", c = "k", alpha = .5)


fig.savefig("profilo.jpg", format = "jpg")
 
plt.show()

Divergenza del fascio (con fit)

def myGauss(x, a, mu, sigma):
    return a * np.exp(- (x-mu)**2 / (2*sigma**2))

disty = 28.5 #cm
distx = 34 #cm

thetax = np.arctan( (x2-x1) / distx ) * 1e3 # mrad
thetay = np.arctan( (y2-y1) / disty ) * 1e3 # mrad

theta = [thetax, thetay]
thetaLbl = ["x", "y"]

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

for i in range(2):
    h, bins = np.histogram(theta[i], bins = 100, )
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    p0 = (np.max(h), binc[np.argmax(h)], np.std(theta[i]))
    cond = h>0
    popt, pcov = curve_fit(myGauss, binc[cond], h[cond], sigma = np.sqrt(h[cond]), absolute_sigma = True, p0 = p0)
    
    ax[i].plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = .7,
               label = fr"$\theta_{thetaLbl[i]}$")
    
    ax[i].plot(binc, myGauss(binc, *popt), ls = "--", c = "k", label = f"Fit (div = {popt[2]:.2f} mrad)")
    
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel(r"$\theta$ [mrad]",fontsize = 14)
    
    ax[i].legend(fontsize = 14)
    
    print(f"La divergenza in {thetaLbl[i]} vale {popt[2]:.2f} mrad")
        
    fig.savefig("divergenza.jpg", format = "jpg")
    
plt.show()

La divergenza in x vale 3.26 mrad
La divergenza in y vale 2.35 mrad

Istogramma PH

Chiaramente se un canale è vuoto o è indicato come rotto il suo istogramma non avrà alcun senso. Questo è un semplice esempio didattico per mostrare come caricare i dati

ph.shape

(44577, 16)

fig, ax = plt.subplots(4,4)
fig.set_size_inches(40, 40)
fig.subplots_adjust(hspace = .2)


ax = ax.flatten()

for i in range(len(ax)):
    
    # Faccio l'istogramma
    h, bins = np.histogram(ph[:,i], bins = 100,)# 384 bin perchè ho 384 strip
    
    # Ottengo il centro del bin
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    # Lo plotto
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green", )
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Energy [ADC]",fontsize = 14)
    
    ax[i].set_yscale("log")
    
    

plt.show()

Correlazione ph tempi

fig, ax = plt.subplots(4,4)
fig.set_size_inches(40, 40)
fig.subplots_adjust(hspace = .2)


ax = ax.flatten()

for i in range(len(ax)):
    
    
    # Lo plotto
    hh = ax[i].hist2d(ph[:,i], time[:,1], cmap = "jet", bins = 100)
    
    # Un po' di cose estetiche
    ax[i].grid()
    ax[i].set_ylabel("PH [ADC]", fontsize = 14)    
    ax[i].set_xlabel("Time [Tick]",fontsize = 14)
    
    fig.colorbar(hh[3], ax = ax[i])
    
    

plt.show()

Conclusioni

In questo quadernetto ho mostrato in maniera schematica principalmente come caricare i dati e come accedere alle varie quantità fisiche per poter costruire i plot di interesse. Non è da intendersi come un esempio di analisi da seguire pari pari, quanto piuttosto un punto di partenza

Riepilogo ciò che credo vada fatto

• Caratterizzare il fascio: studiarne la correlazione, il profilo, e la divergenza. Nel caso di frascati stare attenti alla "condizione di trigger scelta" in quanto il profilo del fascio dipenderà proprio da questa.
• Guardare gli istrogrammi di ph e tempi e la loro correlazione per stabilire dei tagli che consentano di decidere se una particella è stata vista o meno
• Ricostruire le traiettorie dei punti di impatto sull'oggetto di studio
• Verificare che vi sia linearità tra la risposta del calorimetro e l'energia
• Provare ad individuare le varie tiles di scintillatore. sul logbook sono descritte varie run e gli spostamenti che venivano fatti
• Provare a studiare il "Deva effect" considerando la risposta del calorimetor a multi-particelle