Stereogramma

Questo quaderno mostra come realizzare uno stereogramma, ovvero uno scatter plot rot vs crad, dove ogni punto rappresenta la Pulse Height media di un qualche rivelatore (Calorimetro, APC...) per ciascuna posizione (rot, crad) del goniometro.

• È sufficiente fornire i numeri di run da usare e la corrispondente modalità di scansione
• Per punti singoli, è indifferente se mettere rot o crad, ma deve comunque essere specificato
• Questo quaderno è predisposto per costruire due stereogrammi usando le PH di due diversi rivelatori. Inoltre, verrà effettuata la media sopra una certa soglia sogliola, che deve essere fornita
• TODO: allo stato attuale delle cose, non supporta scan diagonali o meshgrid (ma neanche la DAQ...)

Questo quaderno è stato realizzato in collaborazone con Giulia

import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl

import h5py
import json

import os

from matplotlib.patches import Circle, Ellipse

Definisco il percorso dei file

filePath = "/eos/project-i/insulab-como/testBeam/TB_T9_2022/HDF5"

Definisco e visualizzo il file JSON con le varie run e scan type

In realtà, sarebbe sufficiente avere o un file json esterno, o scriversi il dizionario qui dentro. Diciamo che è facilmente adattabile ad altre necessità, soprattutto se diventasse necessario passare più argomenti per ciascuna run.

dizi = {
    '550299' : 'rot',
    '550300' : 'crad',
    '550303' : 'rot',
    '550304' : 'crad'
}

JSONfile = r"./stereogramma.json"

with open(JSONfile, "w") as f:
    json.dump(dizi, f, indent=4)

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

print(dizi.keys())

dict_keys(['550299', '550300', '550303', '550304'])

Scelta dei rivelatori da usare

Definisco:

• Gli indici nella matrice delle pulse height
• Le soglie da usare
• Le label che compariranno

# Indici digitizer: 1 -> LG, 4 -> APC2

indici = [1, 4]
sogliola0 = 100
sogliola1 = 100
label0 = 'LG'
label1 = 'APC2'

Itero su ciascuna run

Calcolo la PH media per ogni step, la appendo alla lista generale, e visualizzo gli istogrammi

# Preparo le liste per i dati che voglio tirarmi fuori
mean0_lst = []
mean1_lst = []
rot_lst = []
crad_lst = []

# Ciclo sui run
for k in dizi.keys():

    # Percorso completo del file da caricare
    fileName = os.path.join(filePath, f"run{k}.h5")
    
    # Importo i dati
    with h5py.File(fileName,'r', libver='latest', swmr=True) as hf:
        xpos = np.array(hf['xpos'])
        digi_ph = np.array(hf['digi_ph'])
        xinfo = np.array(hf['xinfo'])
        info_plus = np.array(hf['info_plus'])

    # Tolgo gli eventi problematici
    # Tengo solo gli eventi validi, ed estraggo le sotto componenti utili per il resto dell'analisi
    logico = (xpos > 0) & (xpos < 10)
    logico2 = logico.all(axis = 1)

    # PH dei due rivelatori da usare
    ph0 = digi_ph[logico2][:, indici[0]]
    ph1 = digi_ph[logico2][:, indici[1]]

    # Rot e Crad del goniometro
    gonio = [xinfo[logico2][:,0], xinfo[logico2][:,1]]
    
    # Numero di step, e valori unici
    step = info_plus[logico2][:, 1]
    unique_step = np.unique(step)

    # Valori di rot e crad ad ogni step
    rot_step  = [gonio[0][step == i][0] for i in unique_step]
    crad_step = [gonio[1][step == i][0] for i in unique_step]

    # Guardo che tipo di scan è (rot o crad) => DA SISTEMARE
    scan_type = dizi[k]
    if scan_type == 'rot':
        scan_step = rot_step
    elif scan_type == 'crad':
        scan_step = crad_step
    else:
        print("MEGA ERRORE") 

    # PH per ogni step
    phxstep0 = [ph0[step == i] for i in unique_step]
    phxstep1 = [ph1[step == i] for i in unique_step]

    # Media delle PH per ogni step
    mean0 = [np.mean(phxstep0[i][phxstep0[i]>sogliola0]) for i in unique_step]
    mean1 = [np.mean(phxstep1[i][phxstep1[i]>sogliola1]) for i in unique_step]

    ######################## Figura ##############################

    
    fig, ax = plt.subplots(2,2)
    fig.set_size_inches(12,12)
    ax = ax.flatten()
    
    
    fig.suptitle(f"nrun = {k} - scan {scan_type}")


    # Plotto gli istogrammi per ciascun step
    for i in unique_step:
        ax[0].hist(phxstep0[i], 100, label = f"{scan_step[i]:.2f} -> {mean0[i]:.2f}", histtype = 'step', lw = 3)
        ax[1].hist(phxstep1[i], 100, label = f"{scan_step[i]:.2f} -> {mean1[i]:.2f}", histtype = 'step', lw = 3)

    opts = {"ls" : "--", "c" : "k"}
    ax[0].axvline(x = sogliola0, **opts)
    ax[1].axvline(x = sogliola1, **opts)
    
    ax[0].set_title(label0)
    ax[1].set_title(label1)
    

    for i in [0, 1]:
        a = ax[i]
        a.set_yscale('log')
        a.set_xlabel('PH [ADC]')
        a.set_ylabel('Counts')
        a.legend()
    


    # Plotto le medie
    ax[2].plot(scan_step, mean0, '--')
    ax[3].plot(scan_step, mean1, '--')

    for i in [2, 3]:
        a = ax[i]
        a.set_xlabel(scan_type + ' [murad]')
        a.set_ylabel("PH mean [ADC]")
    
    
    plt.show()

    # Aggiungo le medie della run corrente alle liste generali

    mean0_lst += mean0
    mean1_lst += mean1
    rot_lst   += rot_step
    crad_lst  += crad_step


# Converto in array di numpy
mean0 = np.array(mean0_lst)
mean1 = np.array(mean1_lst)
rot = np.array(rot_lst)
crad = np.array(crad_lst)

Stereogramma

Plotto lo stereogramma e cerchio il punto più alto

# Opzioni per gli ellissi
opts = {"facecolor":'none', "lw":3, "edgecolor":"k"}

def plottaStereo(currMean, currLabel):

    fig, ax = plt.subplots()
    fig.set_size_inches(12, 10)
    
    # Plotto lo stereogramma
    gonioimg = ax.scatter(rot,crad,c=currMean, cmap="jet")
    
    
    # ax.set_title(currLabel)
    ax.set_xlabel(r"Rot [$\mu$rad]", fontsize = 14)
    ax.set_ylabel(r"Crad [$\mu$rad]", fontsize = 14)
    fig.colorbar(gonioimg, ax = ax)
    
    
    # Ottengo il massimo dello scan e lo cerchio
    myXcenter = rot[np.argmax(currMean)]
    myYcenter = crad[np.argmax(currMean)]
    myXmin, myXmax = ax.get_xlim()
    myYmin, myYmax = ax.get_ylim()
    myPercent = .4
    print(myXcenter, myYcenter)
    
    # Per avere il cerchio che è veramente rotondo, indipendentemente dall'aspect ratio della figura
    bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    width, height = bbox.width, bbox.height
    
    
    ax.add_patch(Ellipse(xy = (myXcenter, myYcenter), width = (myXmax - myXmin) * myPercent / width, height =  (myYmax - myYmin) * myPercent / height, angle=0,  **opts))
    
    ax.grid()
    fig.savefig(f"stereo-{currLabel}.pdf", format = "pdf",bbox_inches='tight', pad_inches=0)
    plt.show()
    


plottaStereo(mean0, label0)
plottaStereo(mean1, label1)

31142.121 -10425.97

32249.232000000004 -11305.403999999999

Bisognerebbe verificare che in questo modo sia davvero rotondo....

# Opzioni per gli ellissi
opts = {"facecolor":'none', "lw":3, "edgecolor":"k"}


# ++++ PLOT 1 +++++
currMean = mean0
currLabel = label0

fig, ax = plt.subplots()
fig.set_size_inches(24, 6)

# Plotto lo stereogramma
gonioimg = ax.scatter(rot,crad,c=currMean, cmap="jet")


ax.set_title(currLabel)
ax.set_xlabel("Rot [murad]")
ax.set_ylabel("Crad [murad]")
fig.colorbar(gonioimg, ax = ax)


bbox = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
width, height = bbox.width, bbox.height



# Ottengo il massimo dello scan e lo cerchio
myXcenter = rot[np.argmax(currMean)]
myYcenter = crad[np.argmax(currMean)]
myXmin, myXmax = ax.get_xlim()
myYmin, myYmax = ax.get_ylim()
myPercent = 1
print(myXcenter, myYcenter)

ax.add_patch(Ellipse(xy = (myXcenter, myYcenter), width = (myXmax - myXmin) * myPercent/width, height =  (myYmax - myYmin) * myPercent/height, angle=0,  **opts))

ax.grid()

plt.show()

31142.121 -10425.97