Esempio analisi dati HDF5¶

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

import sys, os, time, re
import h5py

from scipy.optimize import curve_fit

Load dati¶

Definisco il file da caricare (lavorando su SWAN posso caricarlo direttamente da eos)

In [2]:
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
In [3]:
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"])

# Verifico che abbia 4 colonne
print(xpos.shape)
    
    
# Escludo gli eventi difettosi
"""
Un valore negativo nelle posizioni è un flag per indicare anomalie, dovrebbe essere -999
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())
condToExclude = np.all(((xpos > 0) & (xpos <2)), 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]

print(xpos.shape)

<KeysViewHDF5 ['digi_ph', 'digi_time', 'iatime', 'ievent', 'itime', 'xpos']>
(48251, 4)
-4000.0
(48251,)
(44577, 4)
In [4]:
hf
Out[4]:
<Closed HDF5 file>

Profilo del fascio¶

In [5]:
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 = "darkgreen", lw = .7,
               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()
In [6]:
fig, ax = plt.subplots(2,2)
fig.set_size_inches(12,12)
fig.subplots_adjust(hspace = .2)

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

ax = ax.flatten()

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

for i in range(4):
    h, bins = np.histogram(xpos[:,i], bins = 384, range = (0, 0.005*384)) # Passo strip 50 um
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    
    ax[i].plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = .7,
               label = f"{lstCoord[i]}")
    
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)    
    ax[i].set_xlabel("Pos [cm]",fontsize = 14)
    
    ax[i].legend(fontsize = 14)
    
    ax[i].set_yscale("log")
    
    
# 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("profiloLog.jpg", format = "jpg")
 
plt.show()

Divergenza del fascio¶

In [7]:
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
In [ ]:
In [ ]: