In [1]:
import numpy as np
import sys, os, re, glob
import uproot
from matplotlib import pyplot as plt

Definisco il file da caricare¶

In [2]:
path = "/eos/project/i/insulab-como/testBeam/TB_2023_06_T9_NANOCAL/H2ROOT"
numRun = 650320

nameFile = f"run{numRun}_00010.root"

fileToLoad = os.path.join(path, nameFile)

Visualizzo le varie ntuple e ciò che contengono¶

Fate riferimento al logbook

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

h1;1
<TBranch 'Ievent' at 0x7f3ab1b02670>
<TBranch 'Itime' at 0x7f3ab1b02eb0>
<TBranch 'Iatime' at 0x7f3ab16266a0>
<TBranch 'N1' at 0x7f3ab1626e50>
<TBranch 'N2' at 0x7f3ab162d640>
<TBranch 'N3' at 0x7f3ab162ddf0>
<TBranch 'N4' at 0x7f3ab16325e0>
<TBranch 'N5' at 0x7f3ab1632d90>
<TBranch 'N6' at 0x7f3ab1637580>
<TBranch 'N7' at 0x7f3ab1637d30>
<TBranch 'N8' at 0x7f3ab163e520>
<TBranch 'N9' at 0x7f3ab163ecd0>
<TBranch 'N10' at 0x7f3ab16434c0>
<TBranch 'N11' at 0x7f3ab1643ca0>
<TBranch 'Ivfas_strip1' at 0x7f3ab1648490>
<TBranch 'Ivfas_strip2' at 0x7f3ab1648c40>
<TBranch 'Ivfas_strip3' at 0x7f3ab164f430>
<TBranch 'Ivfas_strip4' at 0x7f3ab164fbe0>
<TBranch 'Ivfas_strip5' at 0x7f3ab16553d0>
<TBranch 'Ivfas_strip6' at 0x7f3ab1655b80>
<TBranch 'Ivfas_strip7' at 0x7f3ab165b370>
<TBranch 'Ivfas_strip8' at 0x7f3ab165bb20>
<TBranch 'Ivfas_strip9' at 0x7f3ab15e0310>
<TBranch 'Ivfas_strip10' at 0x7f3ab15e0ac0>
<TBranch 'Ivfas_data1' at 0x7f3ab15e62b0>
<TBranch 'Ivfas_data2' at 0x7f3ab15e6a60>
<TBranch 'Ivfas_data3' at 0x7f3ab15ec250>
<TBranch 'Ivfas_data4' at 0x7f3ab15eca00>
<TBranch 'Ivfas_data5' at 0x7f3ab15f21f0>
<TBranch 'Ivfas_data6' at 0x7f3ab15f29a0>
<TBranch 'Ivfas_data7' at 0x7f3ab15f7190>
<TBranch 'Ivfas_data8' at 0x7f3ab15f7940>
<TBranch 'Ivfas_data9' at 0x7f3ab15fd130>
<TBranch 'Ivfas_data10' at 0x7f3ab15fd8e0>
<TBranch 'Ivfas_data11' at 0x7f3ab16030d0>
<TBranch 'Xinfo' at 0x7f3ab1603880>
<TBranch 'Info_plus' at 0x7f3ab1603fa0>
<TBranch 'Info_digi' at 0x7f3ab1608820>
<TBranch 'Icrys' at 0x7f3ab1608fd0>
h2;1
<TBranch 'Idigi2' at 0x7f3ab160fd90>
<TBranch 'Idtime2' at 0x7f3ab1615610>
<TBranch 'Idatime2' at 0x7f3ab1615dc0>
<TBranch 'Nwdigi_730' at 0x7f3ab161b5b0>
<TBranch 'Idigi_730' at 0x7f3ab161bd60>
h3;1
<TBranch 'Idigi3' at 0x7f3ab161f6a0>
<TBranch 'Idtime3' at 0x7f3ab15a63a0>
<TBranch 'Idatime3' at 0x7f3ab15a6b50>
<TBranch 'Nwdigi_742a' at 0x7f3ab15ac340>
<TBranch 'Idigi_742a' at 0x7f3ab15acaf0>
<TBranch 'Nwfast_742a' at 0x7f3ab15b12e0>
<TBranch 'Idigifast_742a' at 0x7f3ab15b1a90>
<TBranch 'Ifast_base742a' at 0x7f3ab15b6280>
<TBranch 'Ifast_time742a' at 0x7f3ab15b6a30>
<TBranch 'Ifast_time50_742a' at 0x7f3ab15bd220>
<TBranch 'Ifast_ph742a' at 0x7f3ab15bd9d0>
h100;1
h200;1
h101;1
h201;1
h102;1
h202;1
h103;1
h203;1
h104;1
h204;1
h105;1
h205;1
h106;1
h206;1
h107;1
h207;1
h108;1
h208;1
h109;1
h209;1
h110;1
h210;1
h111;1
h211;1
h112;1
h212;1
h113;1
h213;1
h114;1
h214;1
h115;1
h215;1
h116;1
h216;1
h117;1
h217;1
h1000;1
h1001;1
h1002;1
h1003;1
h1004;1
h1005;1
h1006;1
h1007;1
h1008;1
h1009;1
h4;1
<TBranch 'Idigi4' at 0x7f3ab14625b0>
<TBranch 'Idtime4' at 0x7f3ab1462df0>
<TBranch 'Idatime4' at 0x7f3ab14675e0>
<TBranch 'Nwdigi_742b' at 0x7f3ab1467d90>
<TBranch 'Idigi_742b' at 0x7f3ab146d580>
<TBranch 'Nwfast_742b' at 0x7f3ab146dd30>
<TBranch 'Idigifast_742b' at 0x7f3ab1473520>
<TBranch 'Ifast_base742b' at 0x7f3ab1473c10>
<TBranch 'Ifast_time742b' at 0x7f3ab1478400>
<TBranch 'Ifast_time50_742b' at 0x7f3ab1478bb0>
<TBranch 'Ifast_ph742b' at 0x7f3ab147f3a0>

Carico concretamente i dati¶

In [4]:
with uproot.open(fileToLoad)["h3"] as f:
    f.show()
    for k in f:
        print(k)
    Nwdigi_742a = f["Nwdigi_742a"].arrays(library = "np")["Nwdigi_742a"]
    Idigi_742a = f["Idigi_742a"].arrays(library = "np")["Idigi_742a"]
    
    print(Nwdigi_742a.shape, Idigi_742a.shape)
    
# Condizione da soddisfare, (direi che è sempre vera)
conda = (Nwdigi_742a == 8248)#32992)
print(f"Good events: {conda.sum()}\nTotal: {conda.shape}")

#Nwdigi_742a = Nwdigi_742a[conda]
Idigi_742a = Idigi_742a[conda]

name                 | typename                 | interpretation                
---------------------+--------------------------+-------------------------------
Idigi3               | int32_t                  | AsDtype('>i4')
Idtime3              | int32_t                  | AsDtype('>i4')
Idatime3             | int32_t                  | AsDtype('>i4')
Nwdigi_742a          | uint16_t                 | AsDtype('>u2')
Idigi_742a           | uint16_t[40000]          | AsDtype("('>u2', (40000,))")
Nwfast_742a          | uint16_t                 | AsDtype('>u2')
Idigifast_742a       | uint16_t[4112]           | AsDtype("('>u2', (4112,))")
Ifast_base742a       | int32_t[32]              | AsDtype("('>i4', (32,))")
Ifast_time742a       | int32_t[32]              | AsDtype("('>i4', (32,))")
Ifast_time50_742a    | int32_t[32]              | AsDtype("('>i4', (32,))")
Ifast_ph742a         | int32_t[32]              | AsDtype("('>i4', (32,))")
<TBranch 'Idigi3' at 0x7f3ab1494940>
<TBranch 'Idtime3' at 0x7f3ab14991c0>
<TBranch 'Idatime3' at 0x7f3ab1499970>
<TBranch 'Nwdigi_742a' at 0x7f3ab149e160>
<TBranch 'Idigi_742a' at 0x7f3ab149e910>
<TBranch 'Nwfast_742a' at 0x7f3ab1423100>
<TBranch 'Idigifast_742a' at 0x7f3ab14238b0>
<TBranch 'Ifast_base742a' at 0x7f3ab1423fd0>
<TBranch 'Ifast_time742a' at 0x7f3ab1429850>
<TBranch 'Ifast_time50_742a' at 0x7f3ab1429f70>
<TBranch 'Ifast_ph742a' at 0x7f3ab142e7f0>
(299,) (299, 40000)
Good events: 299
Total: (299,)

Ora la matrice Idigi_742a contiene $191$ righe, una per evento e $40000$ colonne. 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. Nella modalità in questione abbiamo solo 8 canali letti: infatti $1031*8=8248$, che è quello che trovate in Nwdigi_742a. Se leggeste tutti e $32$ i canali, trovereste $1031*32=32992$

Snippet in caso voleste leggere anche l'altro digitizer,

with uproot.open(r"rootdigi/run400857.root")["h4"] as f:
    for k in f:
        print(k)
    Nwdigi_742b = f["Nwdigi_742b"].arrays(library = "np")["Nwdigi_742b"]
    Idigi_742b = f["Idigi_742b"].arrays(library = "np")["Idigi_742b"]

    print(Nwdigi_742a.shape, Idigi_742a.shape)

# Condizione da soddisfare, non ho ben capito perché
condb = (Nwdigi_742b == 32992)
print(f"Good events: {condb.sum()}\nTotal: {condb.shape}")

#Nwdigi_742a = Nwdigi_742a[condb]
Idigi_742b = Idigi_742b[condb]

Ispeziono come sono gli oggetti¶

Chiaramente se vuoi caricarti più files, puoi fare un np.vstack()

In [5]:
print(type(Idigi_742a))
print(Idigi_742a.shape, Idigi_742a.dtype)

<class 'numpy.ndarray'>
(299, 40000) uint16
In [6]:
Idigi_742a
Out[6]:
array([[ 512, 1024, 2202, ...,    0,    0,    0],
       [ 512, 1024, 2394, ...,    0,    0,    0],
       [ 512, 1024, 1805, ...,    0,    0,    0],
       ...,
       [ 512, 1024, 3344, ...,    0,    0,    0],
       [ 512, 1024, 1670, ...,    0,    0,    0],
       [ 512, 1024, 2319, ...,    0,    0,    0]], dtype=uint16)

A titolo esemplificativo, plotto le prime 10 wf del canale 1 (python-like)¶

In [7]:
# Parametri caratteristici che ho appena descritto sopra

nPtsDigi = 1024
nWordSingleChannel = 5+nPtsDigi+2
nChannels = 8 # 32


def getWaveform(matrice, canale, evento):
    """
    matrice: matrice dei dati
    canale: numero di canale, contando da 0
    evento: numero di evento (riga), contando da 0
    """
    i = evento
    j = canale
    return matrice[i, (5 + nWordSingleChannel*j) : (5 + nWordSingleChannel*j + nPtsDigi)]
In [8]:
xVect = np.arange(0, nPtsDigi) * 0.2 #ns, sampling 5 GHz

chanToPlot = 0

for i in range(10):
    fig, ax = plt.subplots()
    fig.set_size_inches(12, 5)
    wf = getWaveform(Idigi_742a, chanToPlot, i)
    
    ax.plot(xVect, wf, ls = ":", c = "tab:green")
    
    ax.grid()
    ax.set_xlabel("Time [ns]", fontsize = 14)
    ax.set_ylabel("[ADC]", fontsize = 14)
    
    ax.set_title(f"Chan. {chanToPlot} -- Ev. {i}", fontsize = 16)
    
    plt.show()
In [ ]: