Altro esempio di analisi dati, spero a breve arriveranno dei commenti più sensati

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

import matplotlib as mpl
from scipy.optimize import curve_fit
import copy
In [2]:
my_cmap = copy.copy(mpl.cm.jet) # copy the default cmap
my_cmap.set_bad(my_cmap(0))

%matplotlib inline

Load dei dati¶

Carico i dati dal file compact, ricordando che c'è una riga di dati interessanti

  • 4 coord silici (x1 inv, y1, x2, y2)
  • 6 ph (canada, vuoto, piani)
  • 6 tempi
  • 3 valori tipo tempo run, tempo assoluto e numero evento

e poi ci sono 6 righe di waveform

In [3]:
myData=[]
dati=("run360388.compact",)
for myFile in dati:
    with open(myFile) as f:
        for i, I in enumerate(f.readlines()):
            II=I.replace("  ", " ").strip()
            if i % 7==0:
                myData.append(II.split(" "))
            
            
myData=np.array(myData, dtype=np.float32)
In [4]:
silici = myData[:,0:4]
ph = myData[:,4:10]   # 6 colonne di pulse height
tempi = myData[:,10:16] # 6 colonne di tempi
absTime = myData[:,17] # Tempo unix

# Elimino possibili punti oltre le coordinate dei silici
var = (silici>0) & (silici<10)
condizione=np.all(var,axis=1)

silici = silici[condizione]
ph = ph[condizione]
tempi = tempi[condizione]
absTime = absTime[condizione]


x1 = silici[:,0]
y1 = silici[:,1]
x2 = silici[:,2]
y2 = silici[:,3]

Correlazioni coordinate¶

Verifico se ci sia una coordinata da invertire

In [5]:
fig, ax = plt.subplots(2,1)

ax[0].hist2d(x1,x2, cmap = "jet", bins = 100)
ax[1].hist2d(y1,y2, cmap = "jet", bins = 100)

plt.show()

Siccome noto che le x sono anticorrelate, ne inverto una per ottenerle correlate. In realtà non importa quale, basta che siano d'accordo

visto che dobbiamo fare un lavoro di fino, la larghezza esatta delle camere non è 10 cm, ma 0.0242384 (larghezza delle strip n° di strip) quindi per "girare" più correttamente la prima x, bisogna fare x1 --> 0.0242*384 - x1

In [6]:
x1 = 0.0242*384 - x1

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

ax[0].hist2d(x1,x2, cmap = "jet", bins = 100)
ax[1].hist2d(y1,y2, cmap = "jet", bins = 100)

plt.show()

Ph vs tempo run¶

Plotto l'istogramma della Ph vs il tempo di run: in caso ci fossero stati cali di tensione, spegnimenti o altre amenità, questo plot può essere un buon indicatore

In [7]:
fig, ax = plt.subplots()

hh = ax.hist2d(absTime, ph[:,2], cmap = "jet", bins = 100)
ax.set_ylim((0,400))
fig.colorbar(hh[3], ax = ax)

plt.show()

Distr. PH canale 1 (che è disconnesso)¶

Posso notare gli spike: mi permette di fissare una ulteriore soglia su questo canale mantenendo al tempo stesso più bassa la soglia sulle ph dei miei piani

In [8]:
sogliaSpike = 44

fig, ax = plt.subplots()

h, bins = np.histogram(ph[:,1], bins = 200, range = (0, 100))
binc = bins[:-1] + (bins[1] - bins[0])/2


ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Spettro ch 1")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)

ax.axvline(x = sogliaSpike, ls = "--", c = "k")

ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_yscale("log")
ax.legend(fontsize = 14)
ax.set_title(f"Canale digi disconnesso", fontsize = 16)

    
plt.show()

Tagli in ph e tempo¶

I primi due canali sono uno scintillatore e vuoto. A me interessano gli indici 2-5, che sono i nostri piani.

Al solito, le soglie vengono pattuite sugli istogrammi 2d, per comodità riporto anche quelli monodimensionali

In [9]:
%matplotlib inline
#soglie = (65, 63, 64, 66) # Prima versione
soglie = (16, 18, 29, 28) # Più basse, se taglio gli spike dall'uno

bincVectPh = []
hVectPh = []

fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,8)
ax = ax.flatten()

fig.suptitle("Distribuzioni ph", fontsize = 16)
fig.subplots_adjust(hspace = .4)


for i, a in enumerate(ax):
    h, bins = np.histogram(ph[:,i+2], bins = 200, range = (0, 1000))
    binc = bins[:-1] + (bins[1] - bins[0])/2
    
    bincVectPh.append(binc)
    hVectPh.append(h)
    
    a.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Pulse height")
    a.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)
    
    a.axvline(x = soglie[i], ls = "--", c = "k")
    
    a.grid()
    a.set_xlabel("Energia [ADC]", fontsize = 14)
    a.set_ylabel("Entries", fontsize = 14)
    a.set_yscale("log")
    a.legend(fontsize = 14)
    a.set_title(f"Piano {i+1}", fontsize = 16)
    
    
plt.show()

bincVectPh = np.array(bincVectPh)
hVectPh = np.array(hVectPh)
In [10]:
%matplotlib inline
tempiMin = (150, 139, 135, 139)
tempiMax = (242, 237, 217, 241)

bincVectTime = []
hVectTime = []


fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,8)
ax = ax.flatten()

fig.suptitle("Distribuzioni tempi", fontsize = 16)
fig.subplots_adjust(hspace = .4)


for i, a in enumerate(ax):
    h, bins = np.histogram(tempi[:,i+2], bins = 200, range = (0, 1000))
    binc = bins[:-1] + (bins[1] - bins[0])/2
    
    bincVectTime.append(binc)
    hVectTime.append(h)
    
    a.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Tempi")
    a.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)
    
    a.axvline(x = tempiMin[i], ls = "--", c = "k")
    a.axvline(x = tempiMax[i], ls = "--", c = "k")
    
    a.grid()
    a.set_xlabel("Tempo [tick]", fontsize = 14)
    a.set_ylabel("Entries", fontsize = 14)
    a.set_yscale("linear")
    a.legend(fontsize = 14)
    a.set_title(f"Piano {i+1}", fontsize = 16)
    
    
plt.show()

bincVectTime = np.array(bincVectTime)
hVectTime = np.array(hVectTime)
In [11]:
%matplotlib qt
%matplotlib inline


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

ax = ax.flatten()

for i, a in enumerate(ax):
    a.hist2d(ph[:,i+2], tempi[:,i+2], bins = (200, 200), range = ((0, 600), (0, 1000)),
            norm = mpl.colors.LogNorm(), cmap = my_cmap)
    
    opts = {"ls":"--", "c":"r", "lw":3}
    
    a.axvline(x = soglie[i], **opts)
    
    a.axhline(y = tempiMin[i], **opts)
    a.axhline(y = tempiMax[i], **opts)
    
    
    a.set_xlabel("Energia [ADC]", fontsize = 14)
    a.set_ylabel("Tempo [tick]", fontsize = 14)
    a.set_title(f"Piano {i+1}", fontsize = 16)
    
plt.show()

Ricostruzione tracce¶

  • La seguente funzione ritorna le proiezioni ad una data quota. Lo zero è fissato in BC1
  • La distanza tra i silici è 5.6 cm
  • La distanza tra il silicio 1 (che sta sopra) e i piani di scinti è 74.6 cm
In [12]:
z2 = 15.6 #cm


def proiettaDistZ(z):
    mx = (x2-x1)/z2
    xProj = x1 + mx * z
    
    my = (y2-y1)/z2
    yProj = y1 + my * z
    
    return (xProj, yProj)

Implemento i tagli¶

I tagli scelti sono

  • Ph maggiore della soglia (quella più bassa)
  • Tempo nell'intervallo scelto
  • Canale 1 disconnesso, metto una soglia per evitare gli spike
In [13]:
logicTagli = []

for i in range(4):
    p = ph[:,i+2]
    t = tempi[:,i+2]
    
    myLogic = (p > soglie[i]) & (t > tempiMin[i]) & (t < tempiMax[i]) & (ph[:,1] < sogliaSpike)
    
    logicTagli.append(myLogic)

Mappe dei punti di impatto¶

In [14]:
xProj, yProj = proiettaDistZ(-5)

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 = "jet")

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

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

plt.show()
In [15]:
# Wide
xProj, yProj = proiettaDistZ(84.6)
#xProj, yProj = proiettaDistZ(86)

fig, ax = plt.subplots()
fig.suptitle("Proiezione sul piano dei piani (wide)", fontsize = 16)


hnt = ax.hist2d(xProj, yProj, bins = 100, 
                            cmap = "jet")

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

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

plt.show()
In [16]:
myRange = ((-15, 25), (-15, 25))
myBins = 50

xProj, yProj = proiettaDistZ(84.6)

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


hnt = ax.hist2d(xProj, yProj, bins = myBins, 
                            range = myRange, cmap = "jet")

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

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

plt.show()
In [17]:
vectPuntiImpattoTagli = []

fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,8)
fig.subplots_adjust(hspace = .4)
ax = ax.flatten()


for i in range(4):
    
    taglio = logicTagli[i]
    
    hscinti = ax[i].hist2d(xProj[taglio], yProj[taglio], bins = myBins, 
                                range = myRange, cmap = "jet")
    

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

    ax[i].set_xlabel("x (cm)", fontsize = 14)
    ax[i].set_ylabel("y (cm)", fontsize = 14)
    ax[i].set_title(f"Piano {i+1}", fontsize = 16)

    
    vectPuntiImpattoTagli.append(copy.copy(hscinti))
    
plt.show()
    
vectPuntiImpattoTagli = np.array(vectPuntiImpattoTagli)

C:\Users\steca\AppData\Local\Temp\ipykernel_16900\9005217.py:28: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  vectPuntiImpattoTagli = np.array(vectPuntiImpattoTagli)

Mappe di efficienza¶

In [18]:
%matplotlib qt
%matplotlib inline

ht = vectPuntiImpattoTagli[0] # Per esempio, uso lo zero, tanto hanno i range utuali

vectDivisioni = []

# Definisco gli indici entro cui sommare
xcoordmin = -3
xcoordmax = 13
ycoordmin = 5
ycoordmax = 17

minx = np.nonzero(ht[1] > xcoordmin)[0][0]
maxx = np.nonzero(ht[1] > xcoordmax)[0][0]
miny = np.nonzero(ht[2] > ycoordmin)[0][0]
maxy = np.nonzero(ht[2] > ycoordmax)[0][0]

print(minx,maxx,miny,maxy)


opts = {"ls":"--", "c":"lime", "lw":3}

fig, ax = plt.subplots(2,2)
fig.set_size_inches(16,8)
fig.subplots_adjust(hspace = .4, wspace = .2)

ax = ax.flatten()

for i in range(4):
    ht = vectPuntiImpattoTagli[i]
    divisione = ht[0] / hnt[0]
    
    vectDivisioni.append(divisione)
    
    myPlotEff = ax[i].imshow(divisione.T, cmap = my_cmap, origin = "lower",
                   extent = (ht[1].min(), ht[1].max(), ht[2].min(), ht[2].max()))
    
    
    ax[i].axvline(x = xcoordmin, **opts)
    ax[i].axvline(x = xcoordmax, **opts)
    
    ax[i].axhline(y = ycoordmin, **opts)
    ax[i].axhline(y = ycoordmax, **opts)
    
    """
    ax[i].axvline(x = hnt[1][15], **opts)
    ax[i].axvline(x = hnt[1][35], **opts)
    
    ax[i].axhline(y = hnt[2][40], **opts)
    ax[i].axhline(y = hnt[2][25], **opts)
    """
    
    fig.colorbar(myPlotEff, ax = ax[i])
    
    ax[i].set_xlabel("x (cm)", fontsize = 14)
    ax[i].set_ylabel("y (cm)", fontsize = 14)
    ax[i].set_title(f"Piano {i+1}", fontsize = 16)
    

                 
                 
plt.show()

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

16 36 26 41

Scritta in modo orribile, urge rivederla

In [19]:
def line(x,b):
    return 0*x+b

ht = vectPuntiImpattoTagli[1] #Per esempio, uso lo zero, tanto hanno i range utuali

for j, divisione in enumerate(vectDivisioni):

    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[minx:maxx, :], axis = 0)
    proiezx = np.nanmean(divisione[:, miny:maxy], axis = 1)

    erry = np.nanstd(divisione[minx:maxx, :], axis = 0)/np.sqrt(20)
    errx = np.nanstd(divisione[:, miny:maxy], axis = 1)/np.sqrt(15)




    poptx, pcovx = curve_fit(line, xVect[minx:maxx], proiezx[minx:maxx], 
                             sigma=errx[minx:maxx], p0 = (.9,))
    popty, pcovy = curve_fit(line, yVect[miny:maxy], proiezy[miny:maxy], 
                             sigma=erry[miny:maxy], p0 = (.9,))




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

    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[15:35], line(xVect[15:35], *poptx), ls = "--", c = "k", 
                zorder = 2, label = "fit costante")
    ax[1].plot(yVect[25:40], line(yVect[25:40], *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)
    fig.suptitle(f"Piano {j+1}", fontsize = 16)


    plt.show()

    print(poptx, popty)

[0.77456844] [0.76780788]

[0.84215955] [0.84252488]

[0.76719587] [0.75365857]

[0.77877104] [0.77514257]

I valori sono leggermente più bassi di quelli trovati da luca, in my opinion, perchè le distanze non sono precisissime, ma leggermente approssimative

WORK IN PROGRESS

Suggerimento di Luca: anziché proiettare la mappa di efficienza, faccio direttamente gli spettri¶

Istogramma dei punti di impatto¶

Va sistemato...

In [20]:
xProj, yProj = proiettaDistZ(84.6)
myBins = 50

# Fette in cui calcolare l'efficienza
myxMin = 5#-5
myxMax = 10#15
myyMin = 5#0
myyMax = 15#20

ulterioreTagliox = (xProj > myxMin) & (xProj < myxMax)
ulterioreTaglioy = (yProj > myyMin) & (yProj < myyMax)

xProj = xProj[ulterioreTaglioy]
yProj = yProj[ulterioreTagliox]



# Istogrammo i punti di impatto
hx, binsx = np.histogram(xProj, bins = myBins)
bincx = binsx[:-1] + (binsx[1] - binsx[0])/2

hy, binsy = np.histogram(yProj, bins = myBins)
bincy = binsy[:-1] + (binsy[1] - binsy[0])/2



# Plotto gli spettri dei punti di impatto
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,5)

ax[0].plot(bincx, hx, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Spettro")
ax[0].fill_between(bincx, hx, step = "mid", color = "lime", alpha = 1)
ax[1].plot(bincy, hy, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Spettro")
ax[1].fill_between(bincy, hy, step = "mid", color = "lime", alpha = 1)

ax[0].set_title("Coordinata x", fontsize = 16)
ax[1].set_title("Coordinata y", fontsize = 16)


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


plt.show()
In [21]:
fig, ax = plt.subplots(4,4)
fig.set_size_inches(20,20)
ax = ax.flatten()

for i in range(4):
    
    # Definisco il taglio. Oltre a quello in ph serve anche una condizione sulla fetta
    cutx = logicTagli[i][ulterioreTaglioy]
    cuty = logicTagli[i][ulterioreTagliox]
    
    # Istogrammo gli spettri col taglio
    hxt, binsxt = np.histogram(xProj[cutx], bins = myBins)
    bincxt = binsxt[:-1] + (binsxt[1] - binsxt[0])/2

    hyt, binsyt = np.histogram(yProj[cuty], bins = myBins)
    bincyt = binsyt[:-1] + (binsyt[1] - binsyt[0])/2
    
    
    # Plotto gli spettri tagliati
    ax[4*i + 0].plot(bincxt, hxt, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Hit taglio x")
    ax[4*i + 0].fill_between(bincxt, hxt, step = "mid", color = "lime", alpha = 1)
    ax[4*i + 1].plot(bincyt, hyt, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Hit Taglio y")
    ax[4*i + 1].fill_between(bincyt, hyt, step = "mid", color = "lime", alpha = 1)

    ax[4*i + 0].set_title("Coordinata x - taglio", fontsize = 16)
    ax[4*i + 1].set_title("Coordinata y - taglio", fontsize = 16)
    
    
    # Calcolo il rapporto degli spettri
    rappx = hxt / hx
    rappy = hyt / hy
    
    
    # Plotto i rapporti
    ax[4*i + 2].plot(bincxt, rappx, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Efficienza")
    ax[4*i + 3].plot(bincyt, rappy, ds = "steps-mid", c = "darkgreen", lw = 2, label = f"Efficienza")
    
    ax[4*i + 2].set_ylim((-.1, 1.1))
    ax[4*i + 3].set_ylim((-.1, 1.1))
    
    
    
for i in ax:
    i.grid()
    i.legend()
    i.set_ylabel("Entries", fontsize = 14)
    i.set_xlabel("Posizione (cm)", fontsize = 14)
    
plt.show()

C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: divide by zero encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: invalid value encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: divide by zero encountered in true_divide
  rappy = hyt / hy
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: invalid value encountered in true_divide
  rappy = hyt / hy
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: divide by zero encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: invalid value encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: divide by zero encountered in true_divide
  rappy = hyt / hy
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: invalid value encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: divide by zero encountered in true_divide
  rappy = hyt / hy
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: invalid value encountered in true_divide
  rappy = hyt / hy
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: divide by zero encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:30: RuntimeWarning: invalid value encountered in true_divide
  rappx = hxt / hx
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: divide by zero encountered in true_divide
  rappy = hyt / hy
C:\Users\steca\AppData\Local\Temp\ipykernel_16900\1051747373.py:31: RuntimeWarning: invalid value encountered in true_divide
  rappy = hyt / hy
In [ ]:
In [ ]: