Process heavy enubet data in chunk¶
In [1]:
import uproot
import numpy as np
import os
import matplotlib.pyplot as plt
import matplotlib as mpl
from numba import njit, prange
In [2]:
filepath = "/eos/experiment/neutplatform/enubet/testbeam2025/picosec_data/sampic_runs/rootSampicData/processed"
numRun = 28
fileToLoad = os.path.join(filepath, f"sampic_run{numRun}_merged.root")
print(f"File to be loaded --> {fileToLoad}")
File to be loaded --> /eos/experiment/neutplatform/enubet/testbeam2025/picosec_data/sampic_runs/rootSampicData/processed/sampic_run28_merged.root
In [3]:
# Quick look to file
with uproot.open(fileToLoad)["picoTree"] as f:
f.show()
name | typename | interpretation
---------------------+--------------------------+-------------------------------
Board | int32_t[140] | AsDtype("('>i4', (140,))")
Channel | int32_t[140] | AsDtype("('>i4', (140,))")
Cell0TimeStamp | double[140] | AsDtype("('>f8', (140,))")
TimeInstant | double[140] | AsDtype("('>f8', (140,))")
TOTValue | float[140] | AsDtype("('>f4', (140,))")
PeakValue | float[140] | AsDtype("('>f4', (140,))")
Baseline | float[140] | AsDtype("('>f4', (140,))")
Waveform | float[140][64] | AsDtype("('>f4', (140, 64))")
Amplitude | float[140] | AsDtype("('>f4', (140,))")
HitFeb | int32_t[3] | AsDtype("('>i4', (3,))")
ArraySize | int32_t | AsDtype('>i4')
Quick load and pack into 2d histos¶
In [4]:
%%time
# --- Parameters ---
filename = fileToLoad
treename = "picoTree"
chunk_size = 200
n_pixels = 140
n_real_pixels = 135
n_samples = 64
adc_min, adc_max = 0, 1 # That's how it seemed to me
n_bins_adc = 100
adc_edges = np.linspace(adc_min, adc_max, n_bins_adc + 1)
adc_bin_width = adc_edges[1] - adc_edges[0]
# Final histogram: pixel, time, amplitude
hist_all = np.zeros((n_real_pixels, n_samples, n_bins_adc), dtype=np.uint32)
# --- Quick method to accumulate ---
@njit(parallel=True, fastmath=True)
def accumulate_histograms_fast(wf_chunk, hist_all, adc_min, adc_bin_width, n_bins_adc):
n_events, n_pix, n_samp = wf_chunk.shape
for pix in prange(n_pix):
if pix >= hist_all.shape[0]:
continue
for ev in range(n_events):
if wf_chunk[ev, pix, 0] == -1:
continue
for s in range(n_samp):
val = wf_chunk[ev, pix, s]
if val >= adc_min:
b = int((val - adc_min) / adc_bin_width)
if 0 <= b < n_bins_adc:
hist_all[pix, s, b] += 1
# --- Read root file in chunks ---
with uproot.open(filename) as f:
tree = f[treename]
for i, (data, report) in enumerate(tree.iterate(
["Waveform"],
step_size=chunk_size,
library="np",
report=True,
)):
wf = data["Waveform"] # (chunk, 140, 64)
accumulate_histograms_fast(wf, hist_all, adc_min, adc_bin_width, n_bins_adc)
if i%25 == 0:
print(f"Processed chunk {i}")
Processed chunk 0 Processed chunk 25 Processed chunk 50 Processed chunk 75 Processed chunk 100 Processed chunk 125 Processed chunk 150 Processed chunk 175 Processed chunk 200 Processed chunk 225 Processed chunk 250 Processed chunk 275 Processed chunk 300 Processed chunk 325 Processed chunk 350 Processed chunk 375 Processed chunk 400 Processed chunk 425 Processed chunk 450 Processed chunk 475 Processed chunk 500 Processed chunk 525 Processed chunk 550 Processed chunk 575 Processed chunk 600 Processed chunk 625 Processed chunk 650 CPU times: user 1min 46s, sys: 5 s, total: 1min 51s Wall time: 43.2 s
Persistence plot¶
In [5]:
%%time
n_rows, n_cols = 8, 17
fig, ax = plt.subplots(n_rows, n_cols, figsize=(60, 40))
ax = ax.flatten()
for iChan in range(n_real_pixels):
H = hist_all[iChan].T
hh = ax[iChan].imshow(
np.log1p(H),
origin="lower",
aspect="auto",
extent=[0, n_samples, adc_min, adc_max],
cmap="jet",
norm=mpl.colors.LogNorm()
)
ax[iChan].set_title(f"Pix {iChan}", fontsize=20)
fig.colorbar(hh, ax=ax[iChan], fraction=0.046, pad=0.04)
for j in range(n_real_pixels, len(ax)):
ax[j].axis("off")
plt.tight_layout()
plt.show()
CPU times: user 32.8 s, sys: 758 ms, total: 33.6 s Wall time: 33.6 s