Tutorial calcolo efficienza

In questo quaderno ripercorro brevemente una analisi tratta dalla mia tesi. Lo scopo finale è il calcolo dell'efficienza. Essendo l'efficienza definita come il rapporto tra il numero di particelle "viste" dal mio rivelatore e il numero totale di particelle che l'hanno attraversato, non posso fare a meno di includere nella prima parte il calcolo delle soglie in ph e tempo per stabilire quando una particella sia stata vista dal mio detector e la ricostruzione delle tracce, per sapere se effettivamente l'hanno attraversato.

Allo stato attuale (13 aprile 2022, ore 11:11) ho praticamente fatto un copia-incolla brutale della mia analsi. In futuro conto di integrarlo, aggiungendo commenti in cui descrivo meglio i passaggi, e magari includendo uno schema in cui spiego come si ricostruiscono le tracce.

Dati disponibili qua

import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
from scipy.optimize import curve_fit

import copy
import matplotlib as mpl



import os
path = r"G:\Tesi\PITONE"
os.chdir(path)

filedir = "..\\CosmicData\\"

Export = False
TAGLIO = 1 # 1 x taglio quadrato, 2 x taglio diagonale
plt.close("all")

Load dei dati

I files sono sempre una linea di 38 dati divisi come segue: * 8 xy silici * 8 ph digitizer (a me interessano idx 2/3) * 8 tempi digitizer * 3 tempi generici * 4 vuoti ----------- tot = 31 * 8 righe da 10 colonne, che sono una wf Indici di Jack * PH: 10-11 * Time: 18-19

# Predispongo la lista della riga lunga. Mi interessa leggere solo quella
myData = []
"""
datiDaCaricare = ("run360074.compact", "run360080.compact", 
                  "run360086.compact")#, "run360088.compact")
"""
datiDaCaricare = ( "run360080.compact", )

# Nell'ultimo è stato cambiato il PMT, meglio non usare con gli altri tre

for myFile in datiDaCaricare:
    with open(filedir + myFile) as f:
        for i, l in enumerate(f.readlines()):
            
            # i primi 8 hanno due spazi
            ll = l.replace("  ", " ").strip()
            
    
            # Righe lunghe
            if i % 9 == 0:
                #print(l)
                
                myData.append(ll.split(" "))
            

            
            
myData = np.array(myData, dtype = np.float32)

Definizione soglie

Distribuzione dei tempi + soglia taglio

# Definisco le regioni ((start1, stop1) , (start2, stop2)) di tempi 
soglia = ((132, 150), (165, 187)) # sui bins
sogliaTempi = ((132, 150), (165, 187)) #A BREVE USERO' SOLO QUESTA. QUELLA SOPRA E' PER EVITARE PROBLEMI IMPROVVISI


hLst = []
bincLst = []

fig, ax  = plt.subplots(2,1)
fig.subplots_adjust(hspace = .6)

#fig.suptitle("Distribuzione dei tempi", fontsize = 16)

for i in range(2):
    #h, bins, _ = ax[i].hist(myData[:, 18 + i], bins = 100)
    
    # Istogrammo
    h, bins = np.histogram(myData[:, 18 + i], bins = 100)
    
    
    # Memorizzo i valori dell'istogramma
    hLst.append(h)
    binc = bins[:-1] + (bins[1] - bins[0]) / 2
    bincLst.append(binc)
    
    
    # Plotto l'istogramma
    ax[i].plot(binc, h, ds = "steps-mid", c = "tab:green", label = "Tempi di arrivo")
    
    ax[i].grid(True)
    ax[i].set_xlabel(f"Tempo (tick)", fontsize = 14)
    ax[i].set_ylabel(f"Entries", fontsize = 14)
    ax[i].set_title(f"Canale {i+1}", fontsize = 16)
    
    #ax[i].set_yscale("log")
    
    # Per la soglia
    ax[i].axvline( x = sogliaTempi[i][0], ls = ":", c = "k", label = "Soglia taglio")
    ax[i].axvline( x = sogliaTempi[i][1], ls = ":", c = "k")
    
    ax[i].legend()
    

    
if Export:
    plt.savefig('../LATEX/FIGURE/plot-cosmici/distr-tempi.eps', format='eps')
plt.show()

Distribuzione PH con e senza taglio in tempo

#sogliaPH = (1437, 810)
sogliaPH = (250, 250)



# Definiso la condizione di taglio
myLogic = []# verra tolto
logicTaglioTempi = []

fraz = [] # Metterò i due massimi da dividere


for i in range(2):
    myL = (myData[:, 18 + i] > sogliaTempi[i][0]) & (myData[:, 18 + i] < sogliaTempi[i][1])
    myLogic.append(myL)
    logicTaglioTempi.append(myL)
    
   
    
# è l'and delle due condizioni
vettTaglioTempi = logicTaglioTempi[0] & logicTaglioTempi[1]
    
    
  
fig, ax  = plt.subplots(2,1)
fig.subplots_adjust(hspace = .6)


#fig.suptitle("Distribuzione delle ph con taglio in tempo", fontsize = 16)

for i in range(2):
    h_nt, bins_nt = np.histogram(myData[:, 10 + i], bins = 300)
    h_t, bins_t = np.histogram(myData[:, 10 + i][logicTaglioTempi[i]], bins = 300)
    
    binc_nt = bins_nt[:-1] + (bins_nt[1] - bins_nt[0]) / 2
    binc_t = bins_t[:-1] + (bins_t[1] - bins_t[0]) / 2
    
    ax[i].plot(binc_nt, h_nt, ds = "steps-mid", label = "PH no taglio")
    ax[i].plot(binc_t, h_t, ds = "steps-mid", label = "PH taglio in tempi")
    
    ax[i].axvline(x = sogliaPH[i], c = "k", ls = "--", label = "Soglia di taglio")

    
    
    ax[i].grid(True)
    ax[i].set_xlabel(f"ADC", fontsize = 14)
    ax[i].set_ylabel(f"Entries", fontsize = 14)
    ax[i].set_title(f"Canale {i+1}", fontsize = 16)

    ax[i].legend()
    
    
    tmp = h_t[binc_t > sogliaPH[i]].argmax()
    fraz.append(binc_t[binc_t > sogliaPH[i]][tmp]  )
    
    
    
ax[0].set_xlim((0, 4000))
ax[1].set_xlim((0, 2000))

ax[0].set_ylim((0, 400))
ax[1].set_ylim((0, 300))


    
if Export:
    plt.savefig('../LATEX/FIGURE/plot-cosmici/distr-ph.eps', format='eps')
plt.show()
    
    
    
rapportoMassimi = fraz[0] / fraz[1]

Biplot tempo vs ph

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

g1 = ax[0].hist2d(myData[:, 10 ], myData[:, 18], bins = 100, 
                 range=((0, 6000), (120, 160)))
fig.colorbar(g1[3], ax = ax[0])
g2 = ax[1].hist2d(myData[:, 10 ], myData[:, 18], bins = 100, 
                 range=((0, 3000), (120, 160)))
fig.colorbar(g1[3], ax = ax[1])

    
for i in range(2):
    
    ax[i].set_xlabel("PH (ADC)")
    ax[i].set_ylabel("Tempi")


    

plt.show()

Per poter mettere una soglia sensata non posso basarmi solamente sugli istogrammi delle PH o dei tempi. In questo modo guardo le correlazioni tra PH e tempo. Noto che per bassi valori delle PH ho segnale per ogni tempo, mentre superata una certa soglia ho segnale solo in una certa banda: il valor medio dipende dalla lunghezza dei cavi fondamentalmente, mentre la larghezza dal timewalk del discriminatore del pretrigger. Il plot qui sotto è concettualmente uguale a quello sopra, ma può essere più comodo per capire dove ci sia effettivamente ZERO conteggi.

# Canale 1
myGrid1 = np.meshgrid(g1[1][:-1] + (g1[1][1]-g1[1][0])/2, 
                                    g1[2][:-1] + (g1[2][1]-g1[2][0])/2)

z1 = g1[0].T.flatten(order = "C")
x1 = myGrid1[0].flatten(order = "C")
y1 = myGrid1[1].flatten(order = "C")



# Canale 2
myGrid2 = np.meshgrid(g2[1][:-1] + (g2[1][1]-g2[1][0])/2, 
                                    g2[2][:-1] + (g2[2][1]-g2[2][0])/2)

z2 = g2[0].T.flatten(order = "C")
x2 = myGrid2[0].flatten(order = "C")
y2 = myGrid2[1].flatten(order = "C")






fig, ax = plt.subplots(1,2)
fig.subplots_adjust(wspace = .4)
ax[0].scatter(x1, y1, color = "k", alpha = 1, marker = ".",
            s = 1 * z1, label = "Ch 1")

ax[1].scatter(x2, y2, color = "k", alpha = 1, marker = ".",
            s = 1 * z2, label = "Ch 2")



ax[0].axvline(x = sogliaPH[0], ls = "--", c = "tab:green", lw = 3)
ax[1].axvline(x = sogliaPH[1], ls = "--", c = "tab:green", lw = 3)




for i in range(2):
    
    ax[i].set_xlabel("PH (ADC)", fontsize = 14)
    ax[i].set_ylabel("Tempi", fontsize = 14)
    ax[i].set_title(f"Canale {i+1}", fontsize = 16)


if Export:
    plt.savefig('../LATEX/FIGURE/plot-cosmici/biplot.eps', format='eps')
plt.show() 

(per verifica consistenza ) - Distribuzione PH con e senza taglio in tempo

fig, ax  = plt.subplots(1, 2)

fig.suptitle("Distribuzione delle ph con taglio in tempo", fontsize = 16)

h_t0, bins_t0 = np.histogram(myData[:, 10][logicTaglioTempi[0]], bins = 150)
h_t1, bins_t1 = np.histogram(myData[:, 11][logicTaglioTempi[1]] * rapportoMassimi, bins = 150)

binc_t0 = bins_t0[:-1] + (bins_t0[1] - bins_t0[0]) / 2
binc_t1 = bins_t1[:-1] + (bins_t1[1] - bins_t1[0]) / 2

ax[0].plot(binc_t0, h_t0, ds = "steps-mid", label = "PH taglio in tempi")
ax[1].plot(binc_t1, h_t1, ds = "steps-mid", label = "PH taglio in tempi")
    
for i in range(2):

    #ax[i].axvline(x = sogliaPH[i], c = "k", ls = "--", label = "Soglia di taglio")

    
    
    ax[i].grid(True)
    ax[i].set_xlabel(f"ADC", fontsize = 14)
    ax[i].set_ylabel(f"Entries", fontsize = 14)
    ax[i].set_title(f"Jack {i+1}", fontsize = 16)

    ax[i].legend(fontsize = 14)
    
    
    tmp = h_t[binc_t > sogliaPH[i]].argmax()
    
    
    
ax[0].set_xlim((0, 4000))
ax[1].set_xlim((0, 2000))

ax[0].set_ylim((0, 800))

    
plt.show()
    
    
    
rapportoMassimi = fraz[0] / fraz[1]

Ricostruzione tracce

Uso il silicio 1/4 (i due estremi)

Setup:

• Silici
• Spazio
• Jack
• Tavolo

Distanza:

• Tra i silici: 5,5 cm

• Tra silicio e tavolo: 36 cm

  Jack è alto 13.5 cm

Ergo:

• Distanza ("primo silicio") - Jack = 36 + 3 * 5.5 - 6.75 = 45.75 cm

Compito preliminare:

• Faccio la retta per ogni evento
• La interseco sul piano di jack, facendo finta che jack non esista

distSilici = 5.5
#distProiezione = 34.75#/33+11
#distProiezione = -16.5 #piano scinti
#sogliola = (1500, 750)

#logicTaglio = (myData[:, 10] > sogliola[0]) | (myData[:, 11] > sogliola[1])



#TAGLIO 1 quadrato: pb nel mezzo
if TAGLIO == 1:
    logicTaglio = (myData[:, 10] > sogliaPH[0]) | (myData[:, 11] > sogliaPH[1])

#TAGLIO 2 in diagonale: risolve il pb
else:
    logicTaglio = (myData[:, 10] + rapportoMassimi * myData[:, 11]) > (sogliaPH[0])


# Chiedo che almeno uno dei due canali sia sopra soglia
# Inoltre chiedo che il tempo sia sopra nel range scelto (si fa con)
logicTaglio = logicTaglio & vettTaglioTempi

Questa funzione prende tutti i dati e li interseca a distanza fissata L'origine è fissata nel primo silicio, ovvero quello più in alto Per andare in alto, si mettte una distanza negativa * Il piano dello scintillatore dista -5.5 ho supposto, come la distanza tra le barche Per andare in basso uso una distanza positiva * Il piano di jack dista

def proiettaAdistanza(distProiezione):
    mx = (myData[:,6] -  myData[:,0]) / (3 * distSilici)
    xProj = mx * distProiezione + myData[:,0]
    
    my = (myData[:,7] -  myData[:,1]) / (3 * distSilici)
    yProj = my * distProiezione + myData[:,1]
    
    return (xProj, yProj)

hist2d punti di impatto su piano scintillatore

Questo plot viene fatto per una sorta di verifica di consistenza, per verificare che sto ricostruendo le tracce nel modo giusto. Se non riuscissi a ricostruire lo scintillatore usato per il trigger, potrebbero esserci problemi nella formula (banalmente controllare se manca un segno meno, oppure se le coordinate sono incolonnate esattamente come mi aspetto...)

Infatti per definizione gli eventi salvati sono tutti e soli quelli che hanno attraversato gli scintillatori di trigger (e che sono stati visti...)

xProj, yProj = proiettaAdistanza(-8) #sht, non è giusto mail plot vien

fig, ax = plt.subplots()
fig.suptitle("Proiezione sul piano dello scintillatore", fontsize = 16)


hscinti = ax.hist2d(xProj, yProj, bins = 100, 
                            range = ((-5, 15), (-5, 15)), cmap = "plasma")

fig.colorbar(hscinti[3], ax = ax)

ax.set_xlabel("x (cm)", fontsize = 14)
ax.set_ylabel("y (cm)", fontsize = 14)


if Export:
    plt.savefig('../LATEX/FIGURE/plot-cosmici/scinti.eps', format='eps')
plt.show()

distr sul piano di BC1

fig, ax = plt.subplots()
fig.suptitle("Distribuzione punti di impatto su BC1", fontsize = 16)


gg = ax.hist2d(myData[:,0], myData[:,1], bins = 50, range = ((0,10), (0,10)),
               cmap = "plasma")

ax.set_xlabel("Posizioni x (cm)", fontsize = 14)
ax.set_ylabel("Posizioni y (cm)", fontsize = 14)

fig.colorbar(gg[3], ax = ax)

if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/profilo.eps', format='eps')
plt.show()
    



proiezx = np.nansum(gg[0], axis = 0)
proiezy = np.nansum(gg[0], axis = 1)

fig, ax = plt.subplots(2,1)
fig.subplots_adjust(hspace = .4)
fig.suptitle("Profilo dei raggi cosmici", fontsize = 16)

ax[0].plot(gg[1][:-1] + (gg[1][1]-gg[1][0])/2, proiezx,
           label = "Profilo x", c = "tab:green", ds = "steps-mid")
ax[1].plot(gg[2][:-1] + (gg[2][1]-gg[2][0])/2, proiezy,
           label = "Profilo y", c = "tab:green", ds = "steps-mid")

ax[0].set_xlabel("Asse x (cm)", fontsize = 14)
ax[1].set_xlabel("Asse y (cm)", fontsize = 14)

for i in range(2):
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)
    ax[i].legend()


if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/profilo-proiez.eps', format='eps')
plt.show()

Distribuzione angolare

L'andamento teorico è come il $\cos^2$. Tuttavia, le particelle che incidono con piccoli angoli rispetto alla verticale attraverseranno gli scintillatori di trigger indipendentemente da dove impattano, mentre le particelle con un grosso angolo rispetto alla verticale, se incidono con un certo angolo $\phi$ potranno attraversare entrambi gli scinti, mentre se incidono con angolo $\phi + \pi$ non riescono ad attraversarli entrambi (un futuro disegno chiarirà..): per questo si fitta solo la parte centrale

def cos2(x,a,b,c):
    return a * (np.cos(b*x + c) ) **2


angx = np.arctan((myData[:,6] -  myData[:,0]) / (3 * distSilici))
angy = np.arctan((myData[:,7] -  myData[:,1]) / (3 * distSilici))


hx, binsx = np.histogram(angx, bins = 100, range = (-.5, .5))
hy, binsy = np.histogram(angy, bins = 100, range = (-.5, .5))
bincx = binsx[:-1] + (binsx[1] - binsx[0]) / 2
bincy = binsy[:-1] + (binsy[1] - binsy[0]) / 2


condx = ((bincx>-.2) & (bincx<.2),)
condy = ((bincy>-.2) & (bincy<.2),)


poptx, pcovx = curve_fit(cos2, bincx[condx], hx[condx], p0 = (400, 1, 0 ))
popty, pcovy = curve_fit(cos2, bincy[condy], hy[condy], p0 = (400, 1, 0 ))




fig, ax = plt.subplots(2,1)#,figsize = (12,5))
fig.subplots_adjust(hspace = .4)

ax[0].plot(bincx, hx, ds = "steps-mid", label = "Angoli in x", c = "tab:green")
ax[0].plot(bincx[condx], cos2(bincx[condx], *poptx), ls = "--", c = "k",
           label = "Curva fittata")

ax[1].plot(bincy, hy, ds = "steps-mid", label = "Angoli in y", c = "tab:green")
ax[1].plot(bincy[condy], cos2(bincy[condy], *popty), ls = "--", c = "k",
           label = "Curva fittata")


for i in range(2):
    ax[i].grid()
    ax[i].set_ylabel("Entries", fontsize = 14)
    ax[i].set_xlabel("Angolo (rad)", fontsize = 14)
    ax[i].legend()
    


if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/distr-ang.eps', format='eps')
plt.show()


taglioX = ((angx > -.1) & (angx < .1))
taglioY = ((angy > -.1) & (angy < .1))

Proietto sul piano di jack

xProj, yProj = proiettaAdistanza(45.75)
fig, ax = plt.subplots(1,2)


ax[0].set_title("Proiezione sul piano di Jack (Taglio)", fontsize = 16)


ht = ax[0].hist2d(xProj[logicTaglio], yProj[logicTaglio], bins = 100, 
                            range = ((-15, 25), (-15, 25)), cmap = "plasma")

fig.colorbar(ht[3], ax = ax[0])


    


    
    
ax[1].set_title("Proiezione sul piano di Jack (No taglio)", fontsize = 16)


hnt = ax[1].hist2d(xProj, yProj, bins = 100, 
                            range = ((-15, 25), (-15, 25)), cmap = "plasma")

fig.colorbar(hnt[3], ax = ax[1])

plt.show()

plt.figure()
plt.hist2d(xProj[logicTaglio], yProj[logicTaglio], bins = 100, 
                            range = ((-15, 25), (-15, 25)), cmap = "plasma")
plt.colorbar()
plt.xlabel("x (cm)", fontsize = 14)
plt.ylabel("y (cm)", fontsize = 14)
plt.show()




plt.figure()
plt.hist2d(xProj, yProj, bins = 100, 
                            range = ((-15, 25), (-15, 25)), cmap = "plasma")
plt.colorbar()
plt.xlabel("x (cm)", fontsize = 14)
plt.ylabel("y (cm)", fontsize = 14)
plt.show()

Mappa di efficienza

La mappa di efficienza è definita come il rapporto tra il numero di particelle viste dal mio scinti e il numero totale di particelle che l'hanno attraversato. Lo faccio per ogni bin, per ciascun tassello della mia suddivisione spaziale. Essendo gli istogrammi 2d delle matrici di numeri, posso fare la divisione elemento per elemento

my_cmap = copy.copy(mpl.cm.plasma) 
my_cmap.set_bad(my_cmap(0))


divisione = ht[0] / hnt[0]

fig, ax = plt.subplots()

ax.set_title("Mappa di efficienza", fontsize = 16)

myPlotEff = ax.imshow(divisione.T, cmap = my_cmap, origin = "lower",
                   extent = (ht[1].min(), ht[1].max(), ht[2].min(), ht[2].max()) )
fig.colorbar(myPlotEff, ax = ax)


ax.set_xlabel("x (cm)", fontsize = 14)
ax.set_ylabel("y (cm)", fontsize = 14)




"""
ax.axvline(x = ht[1][30], c="tab:red", lw=4)
ax.axvline(x = ht[1][70], c="tab:red", lw=4)

ax.axhline(y = ht[2][40], c="tab:red", lw=4)
ax.axhline(y = ht[2][60], c="tab:red", lw=4)
"""

if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/efficienza-t{TAGLIO}.eps', format='eps')
plt.show()

C:\Users\steca\AppData\Local\Temp\ipykernel_16960\1233012966.py:5: RuntimeWarning: invalid value encountered in true_divide
  divisione = ht[0] / hnt[0]

Definisco visivamente delle regioni (una verticale ed una orizzontale) entro cui fare la media perpendicolare. Nota: non medio solo nel quadrato!!

Per esempio, considerando la regione orizzontale comrpresa tra $y=2$ e $y=8$ circa, medio ciascuna colonna, ottenendo quindi un vettore riga, al variare delle $x$, al cui centro ho il calorimentro mentre ai bordi ho il nulla

my_cmap = copy.copy(mpl.cm.plasma) 
my_cmap.set_bad(my_cmap(0))


divisione = ht[0] / hnt[0]

fig, ax = plt.subplots()

ax.set_title("Mappa di efficienza", fontsize = 16)

myPlotEff = ax.imshow(divisione.T, cmap = my_cmap, origin = "lower",
                   extent = (ht[1].min(), ht[1].max(), ht[2].min(), ht[2].max()) )
fig.colorbar(myPlotEff, ax = ax)


ax.set_xlabel("x (cm)", fontsize = 14)
ax.set_ylabel("y (cm)", fontsize = 14)





ax.axvline(x = ht[1][30], c="tab:red", lw=4)
ax.axvline(x = ht[1][70], c="tab:red", lw=4)

ax.axhline(y = ht[2][40], c="tab:red", lw=4)
ax.axhline(y = ht[2][60], c="tab:red", lw=4)


if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/efficienza-t{TAGLIO}.eps', format='eps')
plt.show()

C:\Users\steca\AppData\Local\Temp\ipykernel_16960\2542142099.py:5: RuntimeWarning: invalid value encountered in true_divide
  divisione = ht[0] / hnt[0]

print(f"La media nel quadrato 2-8 vale {np.nanmean(divisione[40:55, 50:55]):.4f}")

La media nel quadrato 2-8 vale 0.9507

Proiezioni

Effettuo le proiezioni nella regione scelta e fitto il tratto piatto con una costante per ottenere l'efficienza

def line(x,b):
    return 0*x+b


xVect = ht[1][:-1] + (ht[1][1] - ht[1][0])/2
yVect = ht[2][:-1] + (ht[2][1] - ht[2][0])/2

proiezy = np.nanmean(divisione[30:70, :], axis = 0)
proiezx = np.nanmean(divisione[:, 40:60], axis = 1)

erry = np.nanstd(divisione[30:70, :], axis = 0)/np.sqrt(40)
errx = np.nanstd(divisione[:, 40:60], axis = 1)/np.sqrt(20)




poptx, pcovx = curve_fit(line, xVect[40:60], proiezx[40:60], 
                         sigma=errx[40:60], p0 = (.9,))
popty, pcovy = curve_fit(line, yVect[40:60], proiezy[40:60], 
                         sigma=erry[40:60], p0 = (.9,))




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

ax[0].errorbar(xVect, proiezx, yerr = errx, 
               ls = "", marker = ".", c = "tab:green", zorder = 1, 
               label = "Proiezione x")
ax[1].errorbar(yVect, proiezy, yerr = erry, 
               ls = "", marker  = ".", c = "tab:green", zorder = 1, 
               label = "Proiezione y")

ax[0].plot(xVect[30:70], line(xVect[30:70], *poptx), ls = "--", c = "k", 
            zorder = 2, label = "fit costante")
ax[1].plot(yVect[40:60], line(yVect[40:60], *popty), ls = "--", c = "k", 
            zorder = 2, label = "fit costante")

for i in range(2):
    ax[i].grid()
    ax[i].legend()
    ax[i].set_ylabel("Efficienza", fontsize = 14)
    ax[i].set_xlabel("Posizione (cm)", fontsize = 14)
    


if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/efficie-proiez.eps', format='eps')
plt.show()

print(poptx, popty)
print(pcovx**.5, pcovy**.5)

C:\Users\steca\AppData\Local\Temp\ipykernel_16960\3111817874.py:8: RuntimeWarning: Mean of empty slice
  proiezy = np.nanmean(divisione[30:70, :], axis = 0)
C:\Users\steca\AppData\Local\Temp\ipykernel_16960\3111817874.py:9: RuntimeWarning: Mean of empty slice
  proiezx = np.nanmean(divisione[:, 40:60], axis = 1)
C:\Users\steca\anaconda3\lib\site-packages\numpy\lib\nanfunctions.py:1664: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,

[0.95626142] [0.95181618]
[[0.0033814]] [[0.00390742]]

Un po' di cose a caso

Correlazioni di PH

fig, ax = plt.subplots(1,2, figsize = (12,5))
fig.suptitle("Correlazioni di PH", fontsize = 16)
fig.subplots_adjust(wspace = .2)



myPlot = ax[0].hist2d(myData[:,10], myData[:,11], 
                         bins = 100, norm = mpl.colors.LogNorm(), 
                         range = ((200, 3000), (200, 3000)), cmap = my_cmap )

ax[0].set_title("Senza taglio", fontsize = 16)
fig.colorbar(myPlot[3], ax = ax[0])


myPlot = ax[1].hist2d(myData[:,10][logicTaglio], myData[:,11][logicTaglio], 
                         bins = 100, norm = mpl.colors.LogNorm(), 
                         range = ((200, 3000), (200, 3000)), cmap = my_cmap )

ax[1].set_title("Con taglio", fontsize = 16)
fig.colorbar(myPlot[3], ax = ax[1])

for i in range(2):
    ax[i].set_xlabel("Canale 1", fontsize = 14)
    ax[i].set_ylabel("Canale 2", fontsize = 14)
    
    ax[i].set_ylim((200,2000))



if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/correlazioni-ph-t{TAGLIO}.eps', format='eps')
plt.show()

Ph in funzione della posizione (Taglio)

fig, ax = plt.subplots(2,2)
ax = ax.flatten()

ax[0].hist2d(myData[:,10][logicTaglio], yProj[logicTaglio], bins = 100, range = ((0, 4000),(-10, 20)))
ax[1].hist2d(myData[:,11][logicTaglio], yProj[logicTaglio], bins = 100, range = ((0, 4000),(-10, 20)))

ax[2].hist2d(myData[:,10][logicTaglio], xProj[logicTaglio], bins = 100, range = ((0, 4000),(-10, 20)))
ax[3].hist2d(myData[:,11][logicTaglio], xProj[logicTaglio], bins = 100, range = ((0, 4000),(-10, 20)))


plt.show()

Il grafico ad x, ma come due diversi istogrammi

fig, ax = plt.subplots(1,2)
ax = ax.flatten()
fig.subplots_adjust(wspace = .4)


hist0 = ax[0].hist2d(myData[:,10][logicTaglio & taglioX], 
                     yProj[logicTaglio & taglioX], bins = 100, 
                     range = ((200, 4000),(-5, 15)))
hist1 = ax[1].hist2d(2 * myData[:,11][logicTaglio & taglioY], 
                     yProj[logicTaglio & taglioY], bins = 100, 
                     range = ((200, 4000),(-5, 15)))

for i in range(2):
    ax[i].set_ylabel("Coordinata y (cm)", fontsize = 14)
    ax[i].set_xlabel("PH normalizzata (ADC)", fontsize = 14)
    ax[i].set_title(f"Canale {i+1}", fontsize = 16)
    
if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/ph-vs-posizione-t{TAGLIO}.eps', format='eps')
plt.show()

Z plot

#Primo canale
binc0x = hist0[1][:-1] + (hist0[1][1] - hist0[1][0]) / 2
binc0y = hist0[2][:-1] + (hist0[2][1] - hist0[2][0]) / 2

myGrid = np.meshgrid(binc0x, binc0y)

z0 = hist0[0].T.flatten(order = "C")
x0 = myGrid[0].flatten(order = "C")
y0 = myGrid[1].flatten(order = "C")


#Secondo canale
binc1x = hist1[1][:-1] + (hist1[1][1] - hist1[1][0]) / 2
binc1y = hist1[2][:-1] + (hist1[2][1] - hist1[2][0]) / 2

myGrid = np.meshgrid(binc1x, binc1y)

z1 = hist1[0].T.flatten(order = "C")
x1 = myGrid[0].flatten(order = "C")
y1 = myGrid[1].flatten(order = "C")






fig, ax = plt.subplots()

ax.set_title("Pulse height vs y coord", fontsize = 16)
ax.set_ylabel("Coordinata y (cm)", fontsize = 14)
ax.set_xlabel("PH normalizzata (ADC)", fontsize = 14)


ax.scatter(x0, y0, color = "k", alpha = 1, marker = "*",
            s = 2 * z0, label = "Ch 1")

ax.scatter(x1, y1, color = "tab:green", alpha = 1, marker = ".",
            s =  2 * z1, label = "Ch2")

ax.legend(fontsize = 14)

ax.grid()


if Export:
    plt.savefig(f'../LATEX/FIGURE/plot-cosmici/zplot-t{TAGLIO}.eps', format='eps')
plt.show()