<jupyter_start><jupyter_text>Files, data, event processing **Quick intro to the following packages**- `uproot` - ROOT I/O in pure Python and NumPy.- `awkward-array` - manipulate arrays of complex data structures as easily as NumPy. &nbsp;        ROOT I/O in pure Python and NumPy  **What is ``uproot``?****Effectively what connects HEP data (ROOT format) with the Python scientific ecosystem around NumPy!**`uproot` provides very fast, efficient, and convenient access to ROOT trees.- Pure Python + NumPy implementation of ROOT I/O.- An array-centric view of ROOT TTree data:  - Branches of simple types are simple arrays.  - Branches of complex types are “jagged arrays”.-  High performance for large baskets, despite Python’s slowness (because all per-entry operations are performed in NumPy).- Greatest benefits: simplicity, minimal installation, set-up, and affinity with machine learning interfaces.  **Why does it exist?**1. To extract columnar data (branches) from a ROOT file without invoking theevent-handling infrastructure of the ROOT framework.3. To express the semantics and conventions of the ROOT file format independentlyof ROOT, in lieu of a formal specification.  **Why Python + NumPy?**- As stressed several times, the scientific Python ecosystem, including much of ML, is designed around a fundamental abstraction called the NumPy array.- Working with computer scientists is easier when you can say, "*pip install uproot*".- Implemented correctly, Python + NumPy doesn't have to be slow.  - Finding the columnar data in a ROOT file may be done in slow Python, as long asdecompression and array manipulations are done by compiled code, see the now-old-ish performance study below.   **1. Getting data from a simple ROOT TTree**By "simple" we mean a file without *jagged structures*, or nested structures with branch sizes depending on an event-by-event basis.<jupyter_code>import uproot
f = uproot.open('data/Zmumu.root')

f<jupyter_output><empty_output><jupyter_text>ROOT files, directories, and trees are like Python dicts with keys() and values().<jupyter_code>f.keys()
t = f[b'events']
t.keys()
t['M']<jupyter_output><empty_output><jupyter_text>Uproot's main purpose is to read branches from ROOT files as NumPy arrays:<jupyter_code>t['M'].array()<jupyter_output><empty_output><jupyter_text>All branches can be looked at with `t.arrays()`. A subset is specified e.g. as `t.arrays(['Run', 'Event'])`:<jupyter_code>t.arrays()<jupyter_output><empty_output><jupyter_text>One can now start performing calculations such as the below. But this is not necessary the best way, think NumPy type of operations ;-).<jupyter_code>import numpy

for px,py,pz in t.iterate(["px1","py1","pz1"], outputtype=tuple):
    pt = numpy.sqrt(px**2 + py**2)
    eta = numpy.arctanh(pz / numpy.sqrt(px**2 + py**2 + pz**2))
    phi = numpy.arctan2(py, px)
    print(pt)
    print(eta)
    print(phi)<jupyter_output>[44.7322 38.8311 38.8311 ... 32.3997 32.3997 32.5076]
[-1.21769 -1.05139 -1.05139 ... -1.57044 -1.57044 -1.57078]
[ 2.74126   -0.440873  -0.440873  ...  0.0370275  0.0370275  0.0369644]<jupyter_text>**Create a Pandas DataFrame**<jupyter_code>import pandas

df = t.pandas.df()
df
t.pandas.df(['Run', 'Event', 'pt1', 'pt2']).head()<jupyter_output><empty_output><jupyter_text>**2. Data with jagged structure**<jupyter_code>f = uproot.open('data/uproot-tutorial-file.root')

f.keys()
branches = f[b'Events'].arrays(namedecode='utf-8')
                               
branches<jupyter_output><empty_output><jupyter_text>The jagged structure here comes from the number of muons per event, which is variable:<jupyter_code>branches['nMuon']<jupyter_output><empty_output><jupyter_text>This becomes evident when checking for example the $p_T$ of all muons:<jupyter_code>branches['Muon_pt']<jupyter_output><empty_output><jupyter_text>Print the $p_T$ for the muons in the first 10 events to trivially see the jagged structure:<jupyter_code>print(' \n'.join([str(elm) for elm in branches['Muon_pt'][:10]]))<jupyter_output>[10.763697 15.736523] 
[10.53849  16.327097] 
[3.2753265] 
[11.429154  17.634033   9.624728   3.5022252] 
[ 3.2834418  3.6440058 32.911224  23.721754 ] 
[3.566528 4.572504 4.371863] 
[57.6067  53.04508] 
[11.319675 23.906353] 
[10.193569 14.204061] 
[11.470704   3.4690065]<jupyter_text>We will get back to jagged arrays in a sec. Let's first show that `uproot` also has (limited) writing functionality.  **3. Interoperability - writing files with uproot** `uproot` version 3 started support for writing files!<jupyter_code>f = uproot.recreate("tmp.root")
f["name"] = numpy.histogram(numpy.random.normal(0, 1, 100000), bins=20)
f.closed
f.keys()
ar = f["name"]

ar.show()<jupyter_output>0                                               18257
                           +---------------------------------------------------+
[-inf, -4.8553)      0     |                                                   |
[-4.8553, -4.411)    1     |                                                   |
[-4.411, -3.9667)    6     |                                                   |
[-3.9667, -3.5224)   18    |                                                   |
[-3.5224, -3.0782)   94    |                                                   |
[-3.0782, -2.6339)   294   |*                                                  |
[-2.6339, -2.1896)   997   |***                                                |
[-2.1896, -1.7453)   2607  |*******                                            |
[-1.7453, -1.301)    5634  |****************                                   |
[-1.301, -0.85676)   9811  |***************************                        |
[-0.85676, -0.41247) 14276 |[...]<jupyter_text>But ... can the file actually be read back in ROOT ...? (Yes it can!)If you have ROOT installed, try```pythonimport ROOTf = ROOT.TFile("tmp.root")h = f.Get("name")c = ROOT.TCanvas("myCanvasName","The Canvas Title",800,600)h.Draw('hist')c.Draw()```If not, and you aren't working in Windows, try to install ROOT via Conda,```python!conda install root -c conda-forge -y```and then try the above. &nbsp;Manipulate arrays of complex data structures as easily as NumPy<jupyter_code>import awkward
branches
table = awkward.Table(branches)

table
table.nMuon[0]
branches['nMuon'][0]
import matplotlib.pyplot as plt

plt.hist(table.nMuon, bins=10, range=(0, 10))
plt.xlabel('Number of muons in event')
plt.ylabel('Number of events');<jupyter_output><empty_output><jupyter_text>How many muon entries are there in total?<jupyter_code>branches['nMuon'].sum()
len(branches['Muon_pt']), len(branches['Muon_pt'].flatten())  # 235286 muons in 1e5 events<jupyter_output><empty_output><jupyter_text>Plot the $p_T$ and $\eta$ of all muons:<jupyter_code>plt.hist(table.Muon_pt.flatten(), bins=100, range=(0, 100))
plt.xlabel('Muon pT')
plt.ylabel('Number of candidates')
plt.yscale('log');
plt.hist(table.Muon_eta.flatten(), bins=100, range=(-2.5, 2.5))
plt.xlabel('Muon $\eta$')
plt.ylabel('Number of candidates');
%%timeit
len(branches['Muon_pt'].flatten())
%%timeit
branches['nMuon'].sum()<jupyter_output>44.9 µs ± 2.56 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>Selections are done via masks. Let's create one that singles out events with a single muon:<jupyter_code>branches['nMuon'] == 1
single_muon_mask = branches['nMuon'] == 1
single_muon_mask.sum()<jupyter_output><empty_output><jupyter_text>Just checking:<jupyter_code>len(branches['Muon_pt'][single_muon_mask]), branches['Muon_pt'][single_muon_mask]
plt.hist(branches['Muon_pt'][single_muon_mask].flatten(), bins=100, range=(0, 100))
plt.xlabel('Muon $p_{\mathrm{T}}$ [MeV]')
plt.ylabel('Number of single muons / 1 MeV')
plt.yscale('log')
plt.show()<jupyter_output><empty_output><jupyter_text>Mask to selection muons within $|\eta| <2$:<jupyter_code>eta_mask = abs(branches['Muon_eta']) < 2
eta_mask
eta_mask.flatten().sum()
eta_mask.sum()
eta_mask.sum().sum()
plt.hist(branches['Muon_eta'].flatten(), bins=50, range=(-2.5, 2.5))
plt.title('No selection')
plt.xlabel('Muon $\eta$')
plt.ylabel('Number of muons')
plt.show()

plt.hist(branches['Muon_eta'][eta_mask].flatten(), bins=50, range=(-2.5, 2.5))
plt.title('With $|\eta| < 2$ selection')
plt.xlabel('Muon $\eta$')
plt.ylabel('Number of muons')
plt.show()
len(single_muon_mask & eta_mask)
plt.hist([branches['Muon_pt'][single_muon_mask & eta_mask].flatten(),
          branches['Muon_pt'][single_muon_mask & ~eta_mask].flatten()],
         label=['$|\eta| < 2$', '$|\eta| \geq 2$'],
         density=True,
         bins=25, range=(0, 50))
plt.xlabel('Muon $p_{\mathrm{T}}$ [GeV]')
plt.ylabel('Number of single muons / 2 GeV')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Concentrate on 2-muon events:<jupyter_code>two_muons_mask = branches['nMuon'] == 2
two_muons_table = table[two_muons_mask]
# Better methods to deal with Lorentz vectors are soon to be available

import uproot_methods

two_muons_p4 = uproot_methods.TLorentzVectorArray.from_ptetaphim(two_muons_table['Muon_pt'],
                                                                 two_muons_table['Muon_eta'],
                                                                 two_muons_table['Muon_phi'],
                                                                 two_muons_table['Muon_mass'])
two_muons_p4
len(two_muons_p4)
first_muon_p4 = two_muons_p4[:, 0]
second_muon_p4 = two_muons_p4[:, 1]
first_muon_p4.delta_r(second_muon_p4)
plt.hist(first_muon_p4.delta_r(second_muon_p4), bins=100)
plt.xlabel('$\Delta R$ between muons')
plt.ylabel('Number of two-muon events')
plt.show()<jupyter_output><empty_output><jupyter_text>Further refine to $\mu^+\mu^-$ pairs - you see where're getting ;-):<jupyter_code>sum_p4 = first_muon_p4 + second_muon_p4
opposite_sign_muons_mask = two_muons_table['Muon_charge'][:, 0] != two_muons_table['Muon_charge'][:, 1]
dimuon_p4 = sum_p4[opposite_sign_muons_mask]
dimuon_p4
import numpy as np

figsize_l, figsize_h = plt.rcParams["figure.figsize"]
plt.figure(figsize=(figsize_l*2.5, figsize_h*3.))

(yvals, binedges, patches) = plt.hist(dimuon_p4.mass, bins=np.logspace(np.log10(0.1), np.log10(1000), 200))

plt.xlabel('Dimuon invariant mass [GeV]')
plt.ylabel('Number of dimuon events')
plt.xscale('log')
plt.yscale('log')

import particle.particle.literals as lpart
from hepunits import GeV
    
list_particles = [getattr(lpart,name) for name in ('eta', 'rho_770_0', 'omega_782','phi_1020','Jpsi_1S', 'Z_0')]

for p in list_particles:
    binnumber = np.searchsorted(binedges, p.mass/GeV)
    plt.text(p.mass/GeV, yvals[binnumber-1]*1.02, '${}$'.format(p.latex_name), horizontalalignment='center', fontsize=16)

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14836961599551526, 'fraction_numerical': 0.05736709333831636, 'mean_word_length': 3.639791937581274, 'pattern_counts': {'":': 0, '<': 77, '<?xml version=': 0, '>': 73, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 25, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '433825', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4140, 'n_tokens_neox': 3707, 'n_words': 898, 'path': '05_Scikit-HEP_files-data-eventprocessing.ipynb', 'repo_name': 'eduardo-rodrigues/2020-10-13_LIV.DAT_school', 'size': '10703'}