
import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl
import scipy.optimize as opt
import copy
import os
import h5py
import json
from matplotlib import ticker, cm
from matplotlib.ticker import LogFormatter


config_file = r"./config.json"


with open(config_file, "r") as f:
            dizi = json.load(f)
        
dizi

{'RunAllignment': 580223,
 'd_12': 925,
 'd_1c': 953,
 'd_13': 2009,
 'offset_x2': 0.02773369382613908,
 'offset_y2': -0.10399194207684116,
 'offset_x3': -0.03472523001319236,
 'offset_y3': 0.061180428232952125,
 'div_12x': 3.459271249720224e-05,
 'div_12y': 4.1801429042218794e-05,
 'div_13x': 3.446746723655453e-05,
 'div_13y': 4.1447333251267285e-05}


# Carico i dati

nRun = 580226

with h5py.File(f"/eos/project/i/insulab-como/testBeam/TB_2023_06_H8_GALORE/HDF5/run{nRun}.h5", 
               'r', libver='latest', swmr=True) as hf:
    print(hf.keys())
    for k in hf.keys():
        print(hf[k].shape)
    xpos = np.array(hf['xpos'])
    xinfo = np.array(hf['xinfo'])

# Tolgo gli eventi "strani"
logico = (xpos > 0) & (xpos < 5)
logico2 = logico.all(axis = 1)

xpos = xpos[logico2]
xinfo = xinfo[logico2]


# Allineo


xpos[:,2] -= dizi['offset_x2']
xpos[:,3] -= dizi['offset_y2']
xpos[:,4] -= dizi['offset_x3']
xpos[:,5] -= dizi['offset_y3']

<KeysViewHDF5 ['ievent', 'info_plus', 'nstrip', 'xinfo', 'xpos']>
(2358010,)
(2358010, 2)
(2358010, 6)
(2358010, 5)
(2358010, 6)


# Mi preparo i vettori

x1 = xpos[:,0]
y1 = xpos[:,1]
x2 = xpos[:,2]
y2 = xpos[:,3]
x3 = xpos[:,4]
y3 = xpos[:,5]

rot = xinfo[:,0]


fig, axes = plt.subplots(3,2)
fig.subplots_adjust(hspace=.4)
fig.set_size_inches(20,20)

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

for i in range(3):
    ax = axes[i,0]
    ax.set_title(f"Tele {i+1} - x")
    
    h, c = np.histogram(xpos[:,2*i], bins = np.linspace(0, 384*0.005, 384))

    c = c[:-1] + (c[1]-c[0])/2
    ax.plot(c, h, ds = 'steps-mid', c = 'hotpink')
    
    ax.set_xlabel("x pos [cm]")
    ax.set_ylabel("Counts")
    ax.grid()
    
    ax = axes[i,1]
    ax.set_title(f"Tele {i+1} - y")
    
    h, c = np.histogram(xpos[:,2*i+1], bins = np.linspace(0, 384*0.005, 384))
    c = c[:-1] + (c[1]-c[0])/2
    
    ax.plot(c, h, ds = 'steps-mid', c = 'hotpink')    
    
    ax.set_xlabel("y pos [cm]")
    ax.set_ylabel("Counts")
    ax.grid()
plt.show()


def proiettaDistZ(z):
    
    mx = (x2-x1)/dizi['d_12']
    xProj = x1 + mx * z
    
    my = (y2-y1)/dizi['d_12']
    yProj = y1 + my * z
    
    return (xProj, yProj)

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


# Proietto le tracce sul piano del cristallo

xc, yc = proiettaDistZ(dizi['d_1c'])

fig, axes = plt.subplots(1,2)
fig.subplots_adjust(hspace=.4)
fig.set_size_inches(12,6)


ax = axes[0]
ax.set_title(f"Proiezioni cristallo - x")
    
h, c = np.histogram(xc, bins = np.linspace(0, 384*0.005, 384))

c = c[:-1] + (c[1]-c[0])/2
ax.plot(c, h, ds = 'steps-mid', c = 'hotpink')
    
ax.set_xlabel("x pos [cm]")
ax.set_ylabel("Counts")
ax.grid()
    
ax = axes[1]
ax.set_title(f"Proiezioni cristallo - y")
    
h, c = np.histogram(yc, bins = np.linspace(0, 384*0.005, 384))
c = c[:-1] + (c[1]-c[0])/2
    
ax.plot(c, h, ds = 'steps-mid', c = 'hotpink')    

ax.set_xlabel("y pos [cm]")
ax.set_ylabel("Counts")
ax.grid()
plt.show()


# Distribuzione angolare tele in => tele 1-2

thetaY_12 = np.arctan((y2-y1)/dizi['d_12'])*1e6 #murad
thetaX_12 = np.arctan((x2-x1)/dizi['d_12'])*1e6 #murad

fig, ax = plt.subplots(1,2)
fig.subplots_adjust(hspace=.4)
fig.set_size_inches(18,8)

fig.suptitle("Distribuzione angolare tele1-2", fontsize = 20)

# Distribuizione in x

h_thetaX,b_thetaX = np.histogram(thetaX_12,bins = 400)
b_thetaX = b_thetaX[0:-1] + (b_thetaX[1]-b_thetaX[0])/2
ax[0].plot(b_thetaX, h_thetaX, ds = 'steps-mid', color = 'hotpink')

sigma = np.std(thetaX_12)
mu = np.mean(thetaX_12)
a = np.max(h_thetaX)

logi = (b_thetaX > (mu - 2*sigma)) & (b_thetaX < (mu + 2*sigma))
p0 = [a, mu, sigma]

x_fit = b_thetaX[logi]
y_fit = h_thetaX[logi]
oP, pC = opt.curve_fit(gaus,xdata = x_fit, ydata = y_fit, sigma=np.sqrt(y_fit), p0 = p0, absolute_sigma=True)
ax[0].plot(x_fit, gaus(x_fit,  *oP), c ='k',ls = '--')

y_mu = gaus(oP[1],*oP)
mu = oP[1]

testo = f'$\mu = $ {mu:.2f} $\mu$rad \n $\sigma$ = {oP[2]:.2f}'
ax[0].plot(mu, y_mu,  '*', color = 'k', ms = 11, label = testo)

ax[0].set_xlabel('$\\theta_X$ [$\mu$rad]')
ax[0].set_ylabel('Entries')
ax[0].grid()
ax[0].legend()

#=====================================================================================================================#


# Distribuizione in y

h_thetaY,b_thetaY = np.histogram(thetaY_12,400)
b_thetaY = b_thetaY[0:-1] + (b_thetaY[1]-b_thetaY[0])/2
ax[1].plot(b_thetaY, h_thetaY, ds = 'steps-mid', color = 'hotpink')

sigma = np.std(thetaY_12)
mu = np.mean(thetaY_12)
a = np.max(h_thetaY)

logi = (b_thetaY > (mu - 1*sigma)) & (b_thetaY < (mu + 1*sigma))
p0 = [a, mu, sigma]
x_fit = b_thetaY[logi]
y_fit = h_thetaY[logi]

oP, pC = opt.curve_fit(gaus, xdata = x_fit, ydata = y_fit, sigma = np.sqrt(y_fit), p0 = p0)

ax[1].plot(x_fit, gaus(x_fit, *oP), c ='k', ls = '--')

y_mu = gaus(oP[1],*oP)
mu = oP[1]

testo = f'$\mu = $ {mu:.2f} $\mu$rad \n $\sigma$ = {oP[2]:.2f}'
ax[1].plot(mu, y_mu,  '*', color = 'k', ms = 11, label = testo)
# 

ax[1].set_xlabel('$\\theta_Y$ [$\mu$rad]')
ax[1].set_ylabel('Entries')
ax[1].grid()
ax[1].legend()

plt.show()


# Distribuzione angolare tele in => Cristallo - tele3
d_c3 = dizi['d_13']- dizi['d_1c']

thetaY_c3 = np.arctan((y3-yc)/d_c3)*1e6 #murad
thetaX_c3 = np.arctan((x3-xc)/d_c3)*1e6 #murad

fig, ax = plt.subplots(1,2)
fig.subplots_adjust(hspace=.4)
fig.set_size_inches(18,8)

fig.suptitle("Distribuzione angolare cristallo - tele 3", fontsize = 20)


# Distribuizione in x

h_thetaX,b_thetaX = np.histogram(thetaX_c3,bins = 400)
b_thetaX = b_thetaX[0:-1] + (b_thetaX[1]-b_thetaX[0])/2
ax[0].plot(b_thetaX, h_thetaX, ds = 'steps-mid', color = 'hotpink')

sigma = np.std(thetaX_c3)
mu = np.mean(thetaX_c3)
a = np.max(h_thetaX)

logi = (b_thetaX > (mu - 2*sigma)) & (b_thetaX < (mu + 2*sigma))
p0 = [a, mu, sigma]

x_fit = b_thetaX[logi]
y_fit = h_thetaX[logi]
oP, pC = opt.curve_fit(gaus,xdata = x_fit, ydata = y_fit, sigma=np.sqrt(y_fit), p0 = p0, absolute_sigma=True)
ax[0].plot(x_fit, gaus(x_fit,  *oP), c ='k',ls = '--')

y_mu = gaus(oP[1],*oP)
mu = oP[1]

testo = f'$\mu = $ {mu:.2f} $\mu$rad \n $\sigma$ = {oP[2]:.2f}'
ax[0].plot(mu, y_mu,  '*', color = 'k', ms = 11, label = testo)

ax[0].set_xlabel('$\\theta_X$ [$\mu$rad]')
ax[0].set_ylabel('Entries')
ax[0].grid()
ax[0].legend()

#=====================================================================================================================#


# Distribuizione in y

h_thetaY,b_thetaY = np.histogram(thetaY_c3,400)
b_thetaY = b_thetaY[0:-1] + (b_thetaY[1]-b_thetaY[0])/2
ax[1].plot(b_thetaY, h_thetaY, ds = 'steps-mid', color = 'hotpink')

sigma = np.std(thetaY_c3)
mu = np.mean(thetaY_c3)
a = np.max(h_thetaY)

logi = (b_thetaY > (mu - 2*sigma)) & (b_thetaY < (mu + 2*sigma))
p0 = [a, mu, sigma]
x_fit = b_thetaY[logi]
y_fit = h_thetaY[logi]

oP, pC = opt.curve_fit(gaus, xdata = x_fit, ydata = y_fit, sigma = np.sqrt(y_fit), p0 = p0)

ax[1].plot(x_fit, gaus(x_fit, *oP), c ='k', ls = '--')

y_mu = gaus(oP[1],*oP)
mu = oP[1]

testo = f'$\mu = $ {mu:.2f} $\mu$rad \n $\sigma$ = {oP[2]:.2f}'
ax[1].plot(mu, y_mu,  '*', color = 'k', ms = 11, label = testo)
# 

ax[1].set_xlabel('$\\theta_Y$ [$\mu$rad]')
ax[1].set_ylabel('Entries')
ax[1].grid()
ax[1].legend()

plt.show()


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

fig.suptitle("Divergenza angolare VS posizione cristallo in x", fontsize = 20)

# Seleziono in x dove si trova il cristallo

x_sx = 0.911
x_dx = x_sx + 0.005

# Creo la mia mappa 

my_cmap = copy.copy(mpl.cm.jet) # copy the default cmap
my_cmap.set_bad(my_cmap(0))

#Istogrammo i dati

hh = ax.hist2d(xc, (thetaX_c3-thetaX_12) , bins = 400, norm = mpl.colors.LogNorm(), cmap = my_cmap)

ax.axvline(x_sx, c = 'white')
ax.axvline(x_dx, c = 'white')

ax.set_xlim(0.83,1.05)
ax.set_ylim(-250, 200)

ax.set_xlabel('x crys [cm]')
ax.set_ylabel('x divergence [$\mu$rad]')
fig.colorbar(hh[3], ax = ax)

plt.show()


logi = (xc > x_sx) & (xc < x_dx)
div_X = thetaX_c3-(thetaX_12)
fig, ax = plt.subplots()
fig.set_size_inches(8, 8)

fig.suptitle("Bruco Plot", fontsize = 20)

my_cmap = copy.copy(mpl.cm.jet) # copy the default cmap
my_cmap.set_bad(my_cmap(0))

bins=[500, 500]
hh = ax.hist2d((rot+thetaX_12)[logi], (div_X )[logi], bins = bins,  cmap = my_cmap, norm = mpl.colors.LogNorm()) 

ax.set_ylim(-100, 100)
ax.set_xlim(-50900, -50750)
ax.xaxis.set_major_locator(ticker.MultipleLocator(30))  # Imposta il tick ogni 100 valori

ax.set_xlabel(r'$\theta_X$ in [$\mu$rad]')
ax.set_ylabel('x divergence [$\mu$rad]')

fig.colorbar(hh[3], ax = ax)
plt.show()


formatter = LogFormatter(10, labelOnlyBase=False) 


logi = (xc > x_sx) & (xc < x_dx)
divX = thetaX_c3-thetaX_12

# Istogrammo

hh = np.histogram2d((rot+thetaX_12)[logi], (divX)[logi], bins = [500,700]) 

#Creo la mia griglia

myGrid = np.meshgrid((hh[1][:-1] + (hh[1][1]-hh[1][0])/2), 
                      (hh[2][:-1] + (hh[2][1]-hh[2][0])/2))
a = hh[1][:-1] + (hh[1][1]-hh[1][0])/2

z1 = hh[0].T

zP = hh[0].flatten(order = "A")
x1 = myGrid[0]
y1 = myGrid[1]

fig, ax = plt.subplots()
fig.set_size_inches(8, 8)

fig.suptitle("Bruco Plot Smooth", fontsize = 20)
z1[z1<1] = 1
h = ax.contourf(x1, y1, z1,  cmap = my_cmap,levels = np.logspace(0, np.log(max(zP)), 100, base = np.e),  norm = mpl.colors.LogNorm())

ax.set_ylim(-100, 100)
ax.set_xlim(-50900, -50750)

ax.xaxis.set_major_locator(ticker.MultipleLocator(100))  # Imposta il tick ogni 100 valori

ax.set_xlabel(r'$\theta_X$ in [$\mu$rad]')
ax.set_ylabel('x divergence [$\mu$rad]')

fig.colorbar(h, ax = ax, format=formatter, )

plt.show()

Distribuzioni angolari¶

Divergenza angolare in funzione della posizione in x del cristallo¶

Bruco plot¶