Load dei dati¶

I dati sono disponibili a questo link

Data format:¶

1.  bit legato allo smart readout. E' una funzione che non siamo ancora riusciti ad utilizzare in maniera efficace, lo dovresti poter ignorare. E' un singolo bit, quindi o 0 o 1;
2.  hitmap del pacchetto. Sono 8 pixel, quindi 8 byte. E' one-hot encoded, un "1" in una data posizione indica che quel pixel è stato colpito, uno "0" che non è stato colpito. Ogni pacchetto ha almeno un 1. In teoria i pixel sono bottom left a bottom right poi top left e top right partendo dal LSB. Lo slave è più vicino al LSB, poi c'è il master. Lo slave è fisicamente più in basso. Romualdo, comunque, sa benissimo come sono, in caso date un'occhiata al suo SW. Comunque da 0 a 0xFF;
3.  corepr è il numero del core che si è acceso. La matrice è alta 512 pixel, e un core è 2x4, quindi ci sono in totale 512/4 = 128 pixel regions. Questo numero va da 0 a 0x7F;
4.  col è il numero della colonna. La matrice è larga 512 pixel, divisa in 16 sezioni, ognuna divisa in 16 colonne da 2 pixel. Quindi la colonna va da 0 a 0xF;
5.  sec è la sezione, come detto prima va da 0 a 0xF anche questa;
6.  ser è il serializzatore da cui proviene il dato. Se tutto funziona come dovrebbe, il dato della n-esima colonna dovrebbe arrivare dall'n-esimo serializzatore. Questo parametro di controllo ci ha portato a trovare diversi errori in realtà. Va comunque da 0 a 0xF;
7.  falling è un booleano. Serve per implementazioni future, in questa versione dovrebbe, se ricordo bene, essere sempre false;
8.  ts è il timestamp dato dal chip vero e proprio. È 1 byte, quindi va da 0 a 0xFF;
9.  ts_fpga è il timestamp dato dall'FPGA. Viene associato al dato quando questo esce dal chip. un tick del timestamp dell'FPGA sono 20 cicli di clock dell'FPGA, a 80 MHz, vuol dire che sono 250 ns. Va da 0 a 0xFFFFFF ( 6 byte);
10. ts_sw è il timestamp associato dal software al dato. E' legato ai secondi che passano durante l'acquisizione. Dovrebbe essere un normale intero python;
11. ts_ext è l'"extended timestamp". E' il reale timestamp del dato, tiene conto del timestamp del chip, di quello dell'FPGA e delle colte che l'FPGA è andata in overflow e ha ricominciato da zero. Dovrebbe essere un normale intero python;
12. E' il numero del trigger - incrementale. Dovrebbe essere un normale intero python.
In [1]:
import numpy as np
from matplotlib import pyplot as plt

import copy
import os
import time
from sys import getsizeof

from itertools import chain

from io import StringIO
import pandas as pd

from itertools import chain
In [2]:
listaFiles = ("run__57__18_03_52.csv",
              "run__16__13_15_14.csv", 
              "run__18__13_27_50.csv")

fileToLoad = rf"./Data/{listaFiles[2]}"
print(f"File to load: {fileToLoad}")

File to load: ./Data/run__18__13_27_50.csv

Parsing del file¶

"(1, 60, 11, 14, 11, 11, False, 98, 2210921, 0, 2210914, 6)","(0, 224, 11, 13, 11, 11, False, 102, 2210925, 0, 2210918, 6)" 1, 60, 11, 14, 11, 11, False, 98, 2210921, 0, 2210914, 6 0, 224, 11, 13, 11, 11, False, 102, 2210925, 0, 2210918, 6

Replace

  • ," -> \n
  • "( -> (
  • )" ->
  • ( ->
In [3]:
t = time.time()

# Lista delle righe parsate
newtext = []

with open(fileToLoad) as f:
    for line in f:
        newline = line.replace(",\"","\n").replace("\"(","(").replace(")\"", "").replace("(","")#.replace("False", "0").replace("True","1")
        #print(newline)
        newtext.append(newline)

# Sarebbe il file scritto bene, una riga per evento, separata da ,
testo_nuovo_file = "".join(newtext)

print(f"Time elapsed: {time.time()-t:.2f}")

# Creo un dizionario contenente l'header
header = eval(testo_nuovo_file.split("\n")[0].replace("\"",""))
# header.keys() # test per vedere se il dizionario è ok


# Creo un oggetto file del testo in questione
#csvStringIO = StringIO("".join(newtext))
csvStringIO = StringIO(testo_nuovo_file)


# Creo un dataframe
nomiCol = ("smartReadout", "hitmap", "corepr", "col", 
           "sec", "ser", "falling", "ts", 
           "ts_fpga", "ts_sw", "ts_ext", "nev")

tipiCol = { 0: np.uint8,    # Smart readout
            1: np.uint8,    # hitmap
            2: np.uint8,    # corepr = pixel region (0, 128)
            3: np.uint8,    # col (0,16)
            4: np.uint8,    # sec (0, 16)
            5: np.uint16,   # ser = serializzatore (dovrebbe essere = col)
            6: "category",  # falling (booleano per implementazioni future)
            7: np.uint8,    # ts = timestam dato dal chip vero e proprio (1 byte)
            8: np.uint64,   # ts_fpga (6 bytes)
            9: np.uint32,   # ts_sw 
            10: np.uint32,  # ts_ext
            11: np.uint32}  # numero di trigger

print(f"Time elapsed: {time.time()-t:.2f}")
           
df = pd.read_csv(csvStringIO, sep=",", header=None, skiprows=1, 
                 names = nomiCol, dtype = tipiCol,)


print(f"Time elapsed: {time.time()-t:.2f}")
print(f"Spazio occupato: {df.memory_usage(deep=True).sum()/1024**2:.2f} Mb")

Time elapsed: 4.52
Time elapsed: 7.25
Time elapsed: 10.94
Spazio occupato: 123.54 Mb

Passo a numpy...¶

In [4]:
sec = df.sec.values
col = df.col.values
corepr = df.corepr.values
hitmap = df.hitmap.values
nev = df.nev.values

eventi = np.unique(nev)

#hitmapStr = hitmapstr.values

xStart = 32*sec + 2*col
yStart = 4 * corepr

maxLen = len(sec)

Ricostruzione pixels colpiti¶

Vado a costruire una matrice che contiene le coordinate: vale nan se in un certo evento un pixel non si è acceso, mentre contiene le coordinate se si è acceso

In [5]:
t = time.time()

# x0, y0, x1 ... x7, y7
# Preinizializzo la matrice che conterrà le coordinate di hit per i pixel
# accesi, e nan per i pixel spenti
coordHit = np.empty((maxLen, 16), dtype = np.int16) * np.nan


# Devo farlo solo se il bit è 1
myBool = np.empty((maxLen, 8), dtype = bool)
for i in range(8):
    myBool[:,i] = hitmap & 1 << i != 0
    
"""
# OPPURE SE VOGLIO LA MATRICE DA 16
myBool = np.empty((maxLen, 16), dtype = bool)
for i in range(8):
    pippo = hitmap & 1 << i != 0
    myBool[:,2*i] = pippo
    myBool[:,2*i+1] = pippo
"""  


# Inizializzo le x
# coordHit[:,::2][myBool] = np.repeat(xStart[:,np.newaxis], 8, axis = 1)[myBool]
# coordHit[:,1::2][myBool] = np.repeat(yStart[:,np.newaxis], 8, axis = 1)[myBool]


for i in range(8):
    # Memorizzo le coordinate iniziali solo per i pixel accesi
    coordHit[:,2*i][myBool[:,i]] = xStart[myBool[:,i]]
    coordHit[:,2*i+1][myBool[:,i]] = yStart[myBool[:,i]]

    
    # Aggiungo la correzione sulla posizione
    coordHit[:,2*i] += i%2
    coordHit[:,2*i+1] += int(np.floor(i/2))
    
    
print(f"Time elapsed: {time.time()-t:.2f}")

Time elapsed: 3.85

Mostro un frame¶

In [6]:
mat = np.zeros((512, 512))

aaa = coordHit[nev == eventi[100]] 
bbb = aaa[~np.isnan(aaa)]
ss = np.reshape(bbb, (int(bbb.size/2), 2)).astype(int)


# Costruisco il frame da mostrare
for i in range(ss.shape[0]): 
    print(ss[i,0], ss[i,1])
    mat[ss[i,0], ss[i,1]] = 1
    
plt.close("all")
fig, ax = plt.subplots(1,2)
for a in ax:
    a.imshow(mat.T, origin = "lower", cmap = "viridis")



# Per ora è un approccio barbaro
ax[1].set_xlim((int(ss[:,0].mean() - 10), 
                int(ss[:,0].mean() + 10)))
ax[1].set_ylim((int(ss[:,1].mean() - 10), 
                int(ss[:,1].mean() + 10)))



plt.show()

193 127
195 126
194 127
197 124
196 125
198 124
199 123
200 122
203 118
204 116
205 116
204 117
205 115
In [ ]:

Approccio datato e poooco efficiente¶

In [7]:
print(f"Andrò a ciclare su {maxLen} elementi: wish me good luck")


currIdx = 0
listaHit = [] # Un elemento per evento
tmpList = [] # Contiene le coordinate di un certo evento


# Ciclo su ciascun evento salvato
t = time.time()
for i in range(maxLen):
    
    if nev[i] != currIdx:
        currIdx = nev[i]
        
        if i != 0:
            
            pippo = np.array(tmpList)
            listaHit.append(pippo)
            #listaHit.append(list(chain.from_iterable(tmpList)))
            tmpList = []
            
    
    # Controllo quale degli 8 pixels si sia acceso
    for j in range(8):
        if hitmap[i] & 1 << j != 0:

            myX = xStart[i] + j%2
            myY = yStart[i] + int(np.floor(j/2))
            
            tmpList.append((myX, myY))
            
    if i==1000: break
            
print(f"Time elapsed: {time.time()-t:.2f}")

Andrò a ciclare su 4466758 elementi: wish me good luck
Time elapsed: 0.05
In [8]:
mat = np.zeros((512, 512))
arrayInQuestione = listaHit[100]

for i in arrayInQuestione: # per esempio voglio vedere il 3
    print(i[0], i[1])
    mat[i[0],i[1]] = 1
    
plt.close("all")
fig, ax = plt.subplots(1,2)
for a in ax:
    a.imshow(mat.T, origin = "lower", cmap = "viridis")


# Per ora è un approccio barbaro
ax[1].set_xlim((int(arrayInQuestione[:,0].mean() - 10), 
                int(arrayInQuestione[:,0].mean() + 10)))
ax[1].set_ylim((int(arrayInQuestione[:,1].mean() - 10), 
                int(arrayInQuestione[:,1].mean() + 10)))



plt.show()

193 127
195 126
194 127
197 124
196 125
198 124
199 123
200 122
203 118
204 116
205 116
204 117
205 115

Altro approccio per parsing dati¶

Al momento salto la prima linea, ma in linea di principio, non è obbligatorio

In [9]:
superLista = []



t = time.time()
with open(fileToLoad) as f:
    for i, line in enumerate(f):
        
        # Salto la prima riga di HEADER.
        # A dire la verità ci sono anche un po' di dati...
        if i == 0: 
            #print("===== HEADER =====")
            #print(line)
            continue

        # Sostituisco le stringhe vero e falso con una rappr. numerica
        line = line.replace("False", "0")
        line = line.replace("True", "1")
        
        # Splitto ciascuna linea in corrispondenza di )","(
        # Stando attento a saltare apici e parentesi ad inizio e fine riga
        superLista.append( line[2:-3].split(')","(') )

    
# Let's flatten the list
superLista = np.array([np.array(i.split(","), dtype = np.uint64) for sublist in superLista for i in sublist])
#superLista = np.array([np.fromstring(i, sep = ",", dtype = np.int64) for sublist in superLista for i in sublist])
#superLista = np.array(list(chain.from_iterable(superLista)))

print(f"Time elapsed: {time.time()-t:.2f}")

print(superLista.shape)
print(type(superLista))

print(f"Dimensione su disco: {getsizeof(superLista)/1024**2:.2f} Mb")

Time elapsed: 23.95
(4466747, 12)
<class 'numpy.ndarray'>
Dimensione su disco: 408.94 Mb