Tutorial analisi waveform CNAO - Soluzione¶

Questo quaderno è la seconda parte (cfr. per la prima) in cui vengono descritti una serie di problemi

I dati sono costituiti dagli HBOOK convertiti in ROOT usando il comando h2root e quindi concatenati.
Il lavoro è stato fatto col seguente script

#!/bin/bash

for f in $(ls ./hbook/*.hbook); do

    # Solo nome (e.g. run016491_001.hbook)
    bb=$(basename $f)

    # Prima parte del nome (e.g. run016491)
    nn=$(echo $bb | cut -d'_' -f1)

    # Nome completo no estensione (e.g. run016491_001)
    nome=${bb/.hbook/}


    # Root singolo (e.g. ./root/run016491_001.root)
    fileParziale=./root/${nome}.root

    # Root unito (e.g. ./root_merged/run016491.root)
    fileTotale=./root_merged/${nn}.root

    h2root $f $fileParziale
    hadd -a $fileTotale $fileParziale


done

La guida di riferimento è questo mio altro esempio realizzato per il test beam di NanoCal, le spiegazioni dettagliate si trovano lì

In [1]:
import sys, os, re, glob
import uproot
import numpy as np
from matplotlib import pyplot as plt
In [2]:
dataPath = r"/eos/project/i/insulab-como/testBeam/CNAO_2023_06/root_merged"
dataPath = r"/eos/project/i/insulab-como/testBeam/CNAO_2023_06/root"
numRun = 16496

fileToLoad = os.path.join(dataPath, f"run{numRun:06d}_001.root")

Ispeziono la struttura del root¶

In [3]:
with uproot.open(fileToLoad) as f:
    for k in f:
        print(k)
        try:
            for kk in f[k]:
                print(kk)
        except:
            pass

h11;1
<TBranch 'I15h1' at 0x7fa59a81bca0>
<TBranch 'Itime' at 0x7fa598f5f520>
<TBranch 'Iatime' at 0x7fa598f5fcd0>
<TBranch 'Iutil' at 0x7fa598f664c0>
<TBranch 'Nword0_742' at 0x7fa598f66c70>
<TBranch 'Idata0_742' at 0x7fa598f6b460>
<TBranch 'Info_onl' at 0x7fa598f6bc10>
h1;1
<TBranch 'I15g1' at 0x7fa598f71670>
<TBranch 'Itime' at 0x7fa598f78280>
<TBranch 'Iatime' at 0x7fa598f78a30>
<TBranch 'Iutil' at 0x7fa598f7e220>
<TBranch 'Nword0_730' at 0x7fa598f7e9d0>
<TBranch 'Idata0_730' at 0x7fa598f831c0>
<TBranch 'Info_onl' at 0x7fa598f83970>
h2;1
<TBranch 'I15g1' at 0x7fa598f898b0>
<TBranch 'Itime' at 0x7fa598f8e130>
<TBranch 'Iatime' at 0x7fa598f8e8e0>
<TBranch 'Iutil' at 0x7fa598f930d0>
<TBranch 'Nword0_730' at 0x7fa598f93880>
<TBranch 'Idata0_730' at 0x7fa598f93fa0>
<TBranch 'Info_onl' at 0x7fa598f98820>
h3;1
<TBranch 'I15g1' at 0x7fa598f20220>
<TBranch 'Itime' at 0x7fa598f20cd0>
<TBranch 'Iatime' at 0x7fa598f274c0>
<TBranch 'Iutil' at 0x7fa598f27c70>
<TBranch 'Nword0_730' at 0x7fa598f2b460>
<TBranch 'Idata0_730' at 0x7fa598f2bc10>
<TBranch 'Info_onl' at 0x7fa598f32400>
h12;1
<TBranch 'I15h1' at 0x7fa598f37190>
<TBranch 'Itime' at 0x7fa598f379d0>
<TBranch 'Iatime' at 0x7fa598f3e1c0>
<TBranch 'Iutil' at 0x7fa598f3e970>
<TBranch 'Nword0_742' at 0x7fa598f42160>
<TBranch 'Idata0_742' at 0x7fa598f42910>
<TBranch 'Info_onl' at 0x7fa598f49100>
h13;1
<TBranch 'I15h1' at 0x7fa598f499d0>
<TBranch 'Itime' at 0x7fa598f4e700>
<TBranch 'Iatime' at 0x7fa598f4eeb0>
<TBranch 'Iutil' at 0x7fa598f536a0>
<TBranch 'Nword0_742' at 0x7fa598f53e50>
<TBranch 'Idata0_742' at 0x7fa598f5a640>
<TBranch 'Info_onl' at 0x7fa598f5adf0>
h4;1
<TBranch 'I15g1' at 0x7fa598ee0730>
<TBranch 'Itime' at 0x7fa598ee5430>
<TBranch 'Iatime' at 0x7fa598ee5be0>
<TBranch 'Iutil' at 0x7fa598eeb3d0>
<TBranch 'Nword0_730' at 0x7fa598eebb80>
<TBranch 'Idata0_730' at 0x7fa598ef0370>
<TBranch 'Info_onl' at 0x7fa598ef0b20>
h110000;1
h120000;1
h5101;1
h5251;1
h5151;1
h5102;1
h5252;1
h5152;1
h5103;1
h5253;1
h5153;1
h5104;1
h5254;1
h5154;1
h5105;1
h5255;1
h5155;1
h5106;1
h5256;1
h5156;1
h5107;1
h5257;1
h5157;1
h5108;1
h5258;1
h5158;1

Carico la branch con le waveform¶

Io ipotizzo che sia h11, ma lo dimostriamo

In [4]:
with uproot.open(fileToLoad) as f:
    f["h11"].show()
    
    wf = f["h11"]["Idata0_742"].array(library = "np")
    nw = f["h11"]["Nword0_742"].array(library = "np")
    
    print(wf.shape, nw.shape, )
    print(wf[0].shape, nw[0].shape, )
    

name                 | typename                 | interpretation                
---------------------+--------------------------+-------------------------------
I15h1                | int32_t                  | AsDtype('>i4')
Itime                | int32_t                  | AsDtype('>i4')
Iatime               | int32_t                  | AsDtype('>i4')
Iutil                | int32_t[16]              | AsDtype("('>i4', (16,))")
Nword0_742           | int32_t                  | AsDtype('>i4')
Idata0_742           | int32_t[]                | AsJagged(AsDtype('>i4'))
Info_onl             | uint16_t[100]            | AsDtype("('>u2', (100,))")
(500,) (500,)
(32992,) ()

Descrizione formato dati¶

Nword0_742¶

Contiene il numero di parole, ovvero il "numero di colonne" che compongono la corrispondente. Questo valore è sempre uguale a $32992$, come si evince dal grafico sopra

Struttura di una singola waveform¶

Una singola waveform è costituita da 1024 campionamenti, preceduti da 5 valori da saltare, e seguiti da altri 2 valori da saltare, per un totale di 1031 punti. Il digitizer V1742 ha 32 canali, quindi si avrà una ripetizione di 32 volte questa unità elementare $1031 * 32 = 32992 $.
In realtà vengono letti due digitizer, un odoscopio ha infatti 64 canali. Quindi ci sono $1031*64$ parole, dove parole da 16 bit vengono a due a due compattate in una parola da 32 bit.

Idata0_742¶

Ora la matrice Idata0_742 contiene $500$ righe, una per evento e $32992 = 32 * 1031$ colonne, una volta che ho trasformato l'array di array in una matrice. Le colonne sono da interpretare come la ripetizione della seguente unità elementare

  • 5 colonne da saltare
  • 1024 punti di waveform
  • 2 colonne da saltare

Questa struttura di $1031$ parole si ripete per ciascun canale.

LA VERA COSA IMPORTANTE¶

Come già osservato, i dati provenienti dal digitizer V1742 hanno risoluzione a 12 bit, per la maggior parte dei digitizer invece a 14 bit. Si potrebbe sicuramente osservare che il tipo di dato int, ovvero intero a 32 bit sia eccessivo. La cosa che più balza all'occhio è l'asse verticale completamente sballato del quaderno precedente.
La chiave, trovata da Valerio, è che ciascun intero a 32 bit è da interpretare come due interi a 16 bit giustapposti. Quella che vedavamo non era altro che la rappresentazione decimale di un numero a 32 bit, ma che avrebbe dovuto essere interpretato come due interi a 16 bit!

Spacchetto le waveform in un modo intelligente¶

In [5]:
# Trasformo da array di array a matrice
waveform = np.stack(wf, axis = 0)
print(waveform.shape)

(500, 32992)
In [6]:
# Creo una terza dimensione per avere indice di evento e indice di canale
waveform2 = waveform[:, np.newaxis, :]
print(waveform2.shape)

(500, 1, 32992)
In [7]:
# Effettuo il reshaping (evento, canale, tempo)
waveformMatrix = waveform2.reshape(waveform.shape[0], 32, 1031)
print(waveformMatrix.shape)

(500, 32, 1031)

Estraggo dunque i due blocchi da 16 bit¶

  • Facendo l'AND logico con 0xFFFF ottengo i 16 bit meno significativi
  • Righ-shiftando a destra di 16 bit, mi rimangono solo i 16 bit più significativi, che ora si trovano in posizione dei 16 meno significativi
In [8]:
least = np.bitwise_and(waveformMatrix[:,:,:], 0xFFFF)
most =  np.right_shift(waveformMatrix[:,:,:], 16)
least.shape, most.shape
Out[8]:
((500, 32, 1031), (500, 32, 1031))

Attualmente mi trovo con due matrici, aventi il primo indice come quasi numero di evento, il secondo indice rappresenta il canale e il terzo indice il campionamento, dove al solito

  • 5 header da saltare
  • 1024 campionamenti
  • 2 trailer da saltare

Una riga di Idata0_742 contiene le informazioni relative a 64 canali, ovvero due digitizer

A questo punto dovrei pescare un elemento da most e uno da least. Alla fine avrò una matrice col doppio dei canali

In realtà inizio a mettere le due parole su una riga, una riga avrà dunque due eventi

In [9]:
# Inizializzo la matrice temporanea, con due waveform per riga
wfMatrix = np.empty((least.shape[0], least.shape[1], least.shape[2]*2), dtype = "int16")
print(wfMatrix.shape, wfMatrix.dtype)

(500, 32, 2062) int16
In [10]:
# Lista per inserire i campionamenti pari
evenIdx = 2 * np.arange(0, least.shape[2], 1, dtype = "int")
print(evenIdx)

[   0    2    4 ... 2056 2058 2060]
In [11]:
# Lista per inserire i campionamenti dispari
oddIdx = 2 * np.arange(0, least.shape[2], 1, dtype = "int") + 1
print(oddIdx)

[   1    3    5 ... 2057 2059 2061]
In [12]:
# Check che tutto sia corretto
print(wfMatrix[:,:,evenIdx].shape) # Una delle due matrici in cui vado ad inserire
print(most.shape)                  # Matrice che vado ad inserire

(500, 32, 1031)
(500, 32, 1031)

Questo punto meriterebbe una discussione dettagliata. Io, ignorantemente mi aspettavo che si riempisse most-least, esattamente come leggo un numero scritto su un foglio. Empiricamente, ho verificato che pur tagliando nel modo giusto, tutte le waveform presentavano un primo punto a zero, e qualcosa di strano alla fine. Invertendo l'ordine, invece le waveform hanno proprio la faccia che devono avere.

Il lettore curioso può invertire e scoprire cosa succede

In [13]:
# Alterno most0-least0-most1-least1....
# wfMatrix[:,:,evenIdx] = most
# wfMatrix[:,:,oddIdx] = least

wfMatrix[:,:,evenIdx] = least
wfMatrix[:,:,oddIdx] = most
In [14]:
# La struttura ora è: 5 - 1024 - 2 - 5 - 1024 - 2
# Mi serve dunque 5:1024+5 U 5+1024+2+5:5+1024+2+5+1024
a = np.arange(5, 1024+5)
print(a.shape)
b = np.arange(5+1024+2+5, 5+1024+2+5+1024)
print(b.shape)

idxToKeep = np.concatenate((a,b))
print(idxToKeep.shape)

# Matrice a cui rimuovo header and trailer
wfMatrixCrop = wfMatrix[:,:,idxToKeep]

print(wfMatrix.shape)
print(wfMatrixCrop.shape)

(1024,)
(1024,)
(2048,)
(500, 32, 2062)
(500, 32, 2048)

Impacchetto gli eventi nei canali uno sotto l'altro¶

Costruisco dunque una matrice col doppio di eventi e metà colonne, e con lo stesso trucchetto alterno eventi pari ed eventi dispari

In [15]:
# Inizializzo la matrice finale
wfMatrixFinal = np.empty((wfMatrixCrop.shape[0], wfMatrixCrop.shape[1] * 2, int(wfMatrixCrop.shape[2]/2)), dtype = "int16")
print(wfMatrixFinal.shape, wfMatrixFinal.dtype)

(500, 64, 1024) int16
In [16]:
# Lista per popolare gli eventi pari
evenIdxRow = 2 * np.arange(0, wfMatrixCrop.shape[1], 1, dtype = "int")
print(evenIdxRow)

[ 0  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46
 48 50 52 54 56 58 60 62]
In [17]:
# Lista per popolare gli eventi dispari
oddIdxRow = 2 * np.arange(0, wfMatrixCrop.shape[1], 1, dtype = "int") + 1
print(oddIdxRow)

[ 1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47
 49 51 53 55 57 59 61 63]
In [18]:
# Notare attentamente la posizione dei : nella terza dimensione..
wfMatrixFinal[:,evenIdxRow,:] = wfMatrixCrop[:,:,:wfMatrixFinal.shape[2]]
wfMatrixFinal[:,oddIdxRow,:] = wfMatrixCrop[:,:,wfMatrixFinal.shape[2]:]

print(wfMatrixFinal.shape)

(500, 64, 1024)

Faccio alcuni plot e verifico che tutto funzioni¶

In [19]:
# Verifico che siano coerenti con un digitizer a 12 bit...
print(wfMatrixFinal.min())
print(wfMatrixFinal.max())
print(wfMatrixFinal.mean())
print(wfMatrixFinal.std())

1457
3693
3501.9082320861817
327.5669701957176

Plot su matrice non croppata e con 2062 numeri per riga, ovvero 2 wf¶

Piuttosto inutile, l'ho usato per stabilire il corretto ordine least-most

In [20]:
print(wfMatrix.shape)

idxEv = 11

fig, ax = plt.subplots(8,4)
ax = ax.flatten()
fig.set_size_inches(60, 40)

for idxChan in range(len(ax)):
    ax[idxChan].plot(wfMatrix[idxEv, idxChan, :])
    
    ax[idxChan].set_xlabel("Time [tick]")
    ax[idxChan].set_ylabel("Voltage [ADC]")
    ax[idxChan].set_title(f"Chan. {idxChan}", fontsize = 20)
    
plt.show()

(500, 32, 2062)

Plot di quella che mi aspetto sia una wf completa¶

In [21]:
idxEv = 10

fig, ax = plt.subplots(16,4)
ax = ax.flatten()
fig.set_size_inches(80, 40)

for idxChan in range(len(ax)):
    ax[idxChan].plot(wfMatrixFinal[idxEv, idxChan, :], c = "tab:green")
    
    ax[idxChan].set_xlabel("Time [tick]")
    ax[idxChan].set_ylabel("Voltage [ADC]")
    ax[idxChan].set_title(f"Chan. {idxChan}", fontsize = 20)
    
plt.show()
In [22]:
# Anche questa cella è stata usata per debuggare l'alternanza least-most
print(wfMatrixFinal[10, 20, :] < 50)
print(wfMatrixFinal[10, 20, :])
print(wfMatrixFinal[10, 20, :].shape)

print()

print(wfMatrix[5, 20, :] < 50)
print(wfMatrix[5, 20, :])
print(wfMatrix[5, 20, :].shape)

[False False False ... False False False]
[3531 3531 3533 ... 3540 3539 3538]
(1024,)

[False False False ... False False False]
[2056 1024  697 ... 3538 2639 2057]
(2062,)

Plot di una waveform, che bellezza!!!¶

In [23]:
idxEv = 10
idxChan = 8

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

ax.plot(np.arange(0,wfMatrixFinal.shape[2],1), wfMatrixFinal[idxEv, idxChan, :], c = "tab:green")

ax.set_xlabel("Time [tick]")
ax.set_ylabel("Voltage [ADC]")
ax.set_title(f"Chan. {idxChan}", fontsize = 20)
    
plt.show()

Plot a posteriori, istogrammi delle componenti¶

In [24]:
# Ispezione le prime 6 componenti di tutte le waveform
# 0-1-2-3-4 sono l'header

fig, ax = plt.subplots(1,6, figsize=(21,3))

for i in range(6):
    ax[i].hist(wfMatrix[:,:,i].flatten(), 100)
    ax[i].set_title(f"pos. {i}")

plt.show()
In [25]:
# 6 componenti successive: sono tutti "dati"

fig, ax = plt.subplots(1,6, figsize=(21,3))

for i in range(6):
    pos = i+6
    ax[i].hist(wfMatrix[:,:,pos].flatten(), 100)
    ax[i].set_title(f"pos. {pos}")
plt.show()
In [26]:
# 1029-1030 sono la fine del primo digitizer
# 1031-1032-1033-1034-1035 sono l'header del secondo digitizer

fig, ax = plt.subplots(1,6, figsize=(21,3))

for i in range(6):
    pos = i+1028
    ax[i].hist(wfMatrix[:,:,pos].flatten(), 100)
    ax[i].set_title(f"pos. {pos}")
plt.show()
In [27]:
# 1031-1032-1033-1034-1035 sono l'header del secondo digitizer

fig, ax = plt.subplots(1,6, figsize=(21,3))

for i in range(6):
    pos = i+1028+6
    ax[i].hist(wfMatrix[:,:,pos].flatten(), 100)
    ax[i].set_title(f"pos. {pos}")
plt.show()
In [28]:
# Gli ultimi due sono l'end del secondo digitizer

fig, ax = plt.subplots(1,6, figsize=(21,3))

for i in range(6):
    pos = -i
    if pos == 0: continue
    ax[i].hist(wfMatrix[:,:,pos].flatten(), 100)
    ax[i].set_title(f"pos. {pos}")
plt.show()
In [29]:
wfMatrixFinal.dtype
Out[29]:
dtype('int16')