%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
rng = np.random.default_rng(1234)
N = 3000
xy = rng.uniform(low=0., high=1., size=[2,N])
v = np.sum(xy**2, axis=0)
w, h = plt.figaspect(1.)
plt.figure(figsize=(w,h))
plt.grid(True)
plt.xlabel('x',labelpad=0.5)
plt.ylabel('y',labelpad=0.5)
plt.scatter(xy[0,:], xy[1,:], s=2)
plt.scatter(xy[0,:][v<1], xy[1,:][v<1], s=2)
<matplotlib.collections.PathCollection at 0xffff7c1b3f10>
values = np.cumsum(np.sum(rng.uniform(0,1,size=[2,N])**2,axis=0)<1) \
/(np.arange(1,N+1))*4
plt.grid(True)
plt.xlabel('n',labelpad=0.5)
plt.ylabel('n/N*4',labelpad=0.5)
plt.plot(np.arange(1,N+1), np.ones(N)*np.pi)
plt.plot(np.arange(1,N+1), values)
plt.yticks([np.pi/2, np.pi, np.pi*1.5], [u'\u03c0/2',u'\u03c0',u'\u03c0*3/2'] )
plt.show()
you can do the plot in logscale and with $\frac{1}{\sqrt(N)}$
plt.grid(True)
plt.xlabel('n',labelpad=0.5)
plt.ylabel('n/N*4',labelpad=0.5)
plt.plot(np.arange(1,N+1), 1/np.sqrt(np.arange(1,N+1))+np.pi)
plt.plot(np.arange(1,N+1), np.ones(N)*np.pi)
plt.plot(np.arange(1,N+1), values)
plt.yticks([np.pi/2, np.pi, np.pi*1.5], [u'\u03c0/2',u'\u03c0',u'\u03c0*3/2'] )
plt.yscale('log')
plt.show()
If you setted correctly also the PYTHONPATH
environment variable
(the thisroot.sh
script should)
!echo $PYTHONPATH
import sys
for s in sys.path:
if "root" in s:
print(s)
/Users/carlo/soft/root-6.26.02-install/lib
import ROOT
Welcome to JupyROOT 6.26/02
Essentially in the same way as in C++
h = ROOT.TH1F("h2","h2",100,-20,20)
h.FillRandom('gaus',1000)
c = ROOT.TCanvas("c","c",800,600)
h.Draw()
You have to Draw
the TCanvas to see something...
c.Draw()
Old versions of ROOT do not throw an exception if there is an error in reading the file
But just a message
therefore if you read the file in a verbose script you may not see the error and it can keep running...
Let me suggest you to wrap the ROOT function to raise an exception in case the file reading has problems:
def my_read_rootfile(filename):
inputfile = ROOT.TFile(filename,'READ')
if not inputfile.IsOpen():
raise FileNotFoundError("file",filename,"not found")
if inputfile.IsZombie():
raise OSError("error opening the file",filename)
my_read_rootfile('nonexisting.root')
--------------------------------------------------------------------------- OSError Traceback (most recent call last) <ipython-input-13-0aa63365e13d> in <module> 6 raise OSError("error opening the file",filename) 7 ----> 8 my_read_rootfile('nonexisting.root') <ipython-input-13-0aa63365e13d> in my_read_rootfile(filename) 1 def my_read_rootfile(filename): ----> 2 inputfile = ROOT.TFile(filename,'READ') 3 if not inputfile.IsOpen(): 4 raise FileNotFoundError("file",filename,"not found") 5 if inputfile.IsZombie(): ~/soft/root-6.26.02-install/lib/ROOT/_pythonization/_tfile.py in _TFileConstructor(self, *args) 53 if len(args) >= 1: 54 if self.IsZombie(): ---> 55 raise OSError('Failed to open file {}'.format(args[0])) 56 57 def _TFileOpen(klass, *args): OSError: Failed to open file nonexisting.root
Error in <TFile::TFile>: file /Users/carlo/repos/pycourse/Slides/nonexisting.root does not exist
A good alternative to installing ROOT is uproot
Let's have a look on how to use it on the output of the Geant4 AnaEx01 extended example
import uproot
you can open a file as:
!wget http://www.roma1.infn.it/~mancinit/Teaching/Alghero22/AnaEx01.root
--2023-06-02 19:37:25-- http://www.roma1.infn.it/~mancinit/Teaching/Alghero22/AnaEx01.root Resolving www.roma1.infn.it (www.roma1.infn.it)... 141.108.26.150, 141.108.26.1 Connecting to www.roma1.infn.it (www.roma1.infn.it)|141.108.26.150|:80... connected. HTTP request sent, awaiting response... 302 Found Location: https://www.roma1.infn.it/~mancinit/Teaching/Alghero22/AnaEx01.root [following] --2023-06-02 19:37:25-- https://www.roma1.infn.it/~mancinit/Teaching/Alghero22/AnaEx01.root Connecting to www.roma1.infn.it (www.roma1.infn.it)|141.108.26.150|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 74210 (72K) Saving to: ‘AnaEx01.root.1’ AnaEx01.root.1 100%[===================>] 72.47K --.-KB/s in 0.03s 2023-06-02 19:37:25 (2.29 MB/s) - ‘AnaEx01.root.1’ saved [74210/74210]
file = uproot.open('AnaEx01.root')
and check the list of keys in the file:
file.keys()
['histo;1', 'histo/EAbs;1', 'histo/EGap;1', 'histo/LAbs;1', 'histo/LGap;1', 'ntuple;1', 'ntuple/Ntuple1;1', 'ntuple/Ntuple2;1']
Get one of the TTrees
tree = file['ntuple/Ntuple1']
check the keys of that TTree
tree.keys()
['Eabs', 'Egap']
and get the branches
branches = tree.arrays()
branches['Egap']
[21.4, 16.7, 21.5, 22.2, 17, 29.3, 13.4, 21.5, 19.9, 19.5, ..., 13.6, 21.4, 21.7, 22.7, 21.2, 21.7, 20.3, 14.1, 17.4] -------------------- type: 1000 * float64
type(branches['Egap'])
awkward.highlevel.Array
finally you can use them (almost) as numpy arrays
for instance, you can do an histogram
import matplotlib.pyplot as plt
plt.hist(branches['Egap'])
(array([ 1., 32., 136., 259., 247., 177., 92., 44., 9., 3.]), array([ 2.71100092, 7.56506954, 12.41913817, 17.27320679, 22.12727542, 26.98134405, 31.83541267, 36.6894813 , 41.54354992, 46.39761855, 51.25168717]), <BarContainer object of 10 artists>)
and make analysis
plt.hist(branches['Egap'])
plt.hist(branches['Egap'][branches['Eabs']<450])
(array([ 1., 5., 22., 48., 79., 70., 59., 35., 9., 3.]), array([ 2.71100092, 7.56506954, 12.41913817, 17.27320679, 22.12727542, 26.98134405, 31.83541267, 36.6894813 , 41.54354992, 46.39761855, 51.25168717]), <BarContainer object of 10 artists>)
you may have some problem with the fact that the output is not a numpy array
plt.hist2d(branches['Egap'],branches['Egap'])
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[19], line 1 ----> 1 plt.hist2d(branches['Egap'],branches['Egap']) File /usr/local/lib/python3.11/site-packages/matplotlib/pyplot.py:2669, in hist2d(x, y, bins, range, density, weights, cmin, cmax, data, **kwargs) 2665 @_copy_docstring_and_deprecators(Axes.hist2d) 2666 def hist2d( 2667 x, y, bins=10, range=None, density=False, weights=None, 2668 cmin=None, cmax=None, *, data=None, **kwargs): -> 2669 __ret = gca().hist2d( 2670 x, y, bins=bins, range=range, density=density, 2671 weights=weights, cmin=cmin, cmax=cmax, 2672 **({"data": data} if data is not None else {}), **kwargs) 2673 sci(__ret[-1]) 2674 return __ret File /usr/local/lib/python3.11/site-packages/matplotlib/__init__.py:1442, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs) 1439 @functools.wraps(func) 1440 def inner(ax, *args, data=None, **kwargs): 1441 if data is None: -> 1442 return func(ax, *map(sanitize_sequence, args), **kwargs) 1444 bound = new_sig.bind(ax, *args, **kwargs) 1445 auto_label = (bound.arguments.get(label_namer) 1446 or bound.kwargs.get(label_namer)) File /usr/local/lib/python3.11/site-packages/matplotlib/axes/_axes.py:7135, in Axes.hist2d(self, x, y, bins, range, density, weights, cmin, cmax, **kwargs) 7132 if cmax is not None: 7133 h[h > cmax] = None -> 7135 pc = self.pcolormesh(xedges, yedges, h.T, **kwargs) 7136 self.set_xlim(xedges[0], xedges[-1]) 7137 self.set_ylim(yedges[0], yedges[-1]) File /usr/local/lib/python3.11/site-packages/awkward/highlevel.py:1113, in Array.__getattr__(self, where) 1108 raise AttributeError( 1109 "while trying to get field {}, an exception " 1110 "occurred:\n{}: {}".format(repr(where), type(err), str(err)) 1111 ) from err 1112 else: -> 1113 raise AttributeError(f"no field named {where!r}") AttributeError: no field named 'T'
ROOT data can be not rectangular
numpy is designed to work with rectangular arrays
#rectangular array
[[0, 1],
[2, 3],
[4, 5],
[6, 7],
[8, 9]]
[[0, 1], [2, 3], [4, 5], [6, 7], [8, 9]]
#jagged array
[[0, 1],
[2, 3, 4],
[5],
[6, 7],
[8, 9]]
[[0, 1], [2, 3, 4], [5], [6, 7], [8, 9]]
import awkward as ak
if the data are numerical and regular can be losslessly converted to a numpy array
plt.hist2d(ak.to_numpy(branches['Egap']),ak.to_numpy(branches['Egap']));
plt.show
<function matplotlib.pyplot.show(close=None, block=None)>