ROOT Histograms in Python

ROOT is a data analysis framework developed by CERN that is well-suited for the analysis of certain large scientific data sets such as particle collision events and astronomical data.

One useful ROOT feature are the histograms. A simple program that histograms 100,000 Gaussian random numbers is shown here:

from __future__ import division, print_function
from ROOT import gRandom, TCanvas, TH1F

c1 = TCanvas('c1', 'Example', 200, 10, 700, 500)
hpx = TH1F('hpx', 'px', 100, -4, 4)

for i in xrange(100000):
px = gRandom.Gaus()
hpx.Fill(px)

hpx.Draw()
c1.Update()

# Save the plot as an image
c1.Print("myhistogram.eps")

It’s output is shown here: A second program is given on page~\pageref{root-histograms} that shows how the ROOT TCanvas can be split into regions to show multiple plots on the same canvas. It’s output is shown here:

We can also split ROOT’s TCanvas into regions to show multiple plots on a single canvas. A simple program that does that is shown here:

from __future__ import division, print_function
from ROOT import gRandom, TCanvas, TH1F

# Create a canvas
c1 = TCanvas('c1', 'My Basic Histograms', 200, 10, 700, 500)

# Divide the canvas into two plot areas
c1.Divide(1,2)

# Create a histograms
hist1 = TH1F('hist1', 'First', 100, -4, 4)
hist2 = TH1F('hist2', 'Second', 100, -4, 4)

for i in xrange(100000):
px = gRandom.Gaus()
hist1.Fill(px)

for i in xrange(100000):
px = gRandom.Gaus()
hist2.Fill(px)

# Go to plot area 1
c1.cd(1)

# Plot the first histogram
hist1.Draw()

# Go to plot area 2
c1.cd(2)

# Plot the second histogram
hist2.Draw()

c1.Update()

# Save the plot as an image
c1.Print("myhistogram.eps")

ROOT’s TLorentzVector Class in Python

ROOT is a data analysis framework developed by CERN that is well-suited for the analysis of certain large scientific data sets such as particle collision events and astronomical data.

One useful part of ROOT is their TLorentzVector class—objects that store and can work with relativistic four-vectors such as the space-time 4-vector or the energy-momentum 4-vector.

Below is a basic Python (2.7) program exploring the TLorentzVector class.

#
# This file explores the use of ROOT's TLorentzVector class, which is useful
# for working with energy-momentum 4-vectors.
#

from __future__ import division, print_function
from ROOT import TLorentzVector

# Create the energy-momentum vector for a photon with energy 5 GeV
photon = TLorentzVector(0, 0, 5.0, 5.0)

# Create the 4-vector for a proton at rest
proton = TLorentzVector(0, 0, 0, 0.938)

# Create vectors for two pi+ pions and one pi- pion
pip1, pip2, pim = TLorentzVector(), TLorentzVector(), TLorentzVector()

# Set the vector components for the pions
pip1.SetPxPyPzE(-0.226178, -0.198456, 1.946048, 1.974144)
pip2.SetPxPyPzE(0.554803, -0.301158, 1.301439, 1.453219)
pim.SetPxPyPzE(-0.07765, 0.072333, 1.372624, 1.38382)

# Print the magnitudes of each 4-vector
print("The magnitude of the photon 4-vector is", photon.Mag())
print("The magnitude of the proton 4-vector is", proton.Mag())
print("The magnitude of the pip1 4-vector is", pip1.Mag())
print("The magnitude of the pip2 4-vector is", pip2.Mag())
print("The magnitude of the pim 4-vector is", pim.Mag())

# Create a new 4-vector by adding/subtracting the others
diff = photon + proton - ( pip1 + pip2 + pim )

diffMass = diff.Mag()
print("The invariant mass of the missing particle is", diffMass)

Barnsley Fern in Python

The Barnsley Fern is a beautiful fractal that can easily be generated in Python. If we zoom in on one branch, we see that the pattern is repeated: The python code follows.

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt

# Main body of program
points = 10000  # The number of points to use.

X = []  # A list of x-coordinates
Y = []  # A list of y-coordinates

# Set the starting point
point = [0.5, 0.0]
X.append(point)
Y.append(point)

def new_point(p):
"""
This function takes in a point (x, y) and generates a new point according
to the given equations.
"""
r = np.random.uniform(0, 1)
if r < 0.02:
p = [0.5, 0.27*p]
elif 0.02 <= r <= 0.17:
p = [-0.139*p + 0.263*p + 0.57, 0.246*p + 0.224*p - 0.036]
elif 0.17 < r <= 0.3:
p = [0.17*p - 0.215*p + 0.408, 0.222*p + 0.176*p + 0.0893]
elif 0.3 < r < 1.0:
p = [0.781*p + 0.034*p + 0.1075, -0.032*p + 0.739*p + 0.27]

return p

# Generate a large number of points
for i in range(points):
point = new_point(point)
X.append(point)
Y.append(point)

# Plot the results
plt.scatter(X, Y, c='g', s=.05)
plt.axis('Off')
plt.axes().set_aspect('equal')
plt.title("Barnsley Fern")
plt.savefig("barnsley_fern.png")
plt.show()

Prime Number Sieve in Python

This program computes all the prime numbers up to 10,000 using an efficient algorithm.

Instead of checking if a number n is prime by dividing by all previous numbers, or even all previous primes, this program only divides by all previous primes less than or equal to the square root of n.

from __future__ import division, print_function

primes = 
print(2)

for n in range(3, 10001):
i, k = 0, 0
while primes[i] <= n**(1/2):    # No need to check against ALL smaller primes
if n % primes[i] == 0:      # Check if divisible by previous primes
k = 1
i += 1
if k == 0:                      # If good, at it to the list and print it
primes.append(n)
print(n)

Newton’s Cradle in Visual Python

This simple program shows an animation of a two-pendulum version of Newton’s cradle. One pendulum begins at some starting angle and the other pendulum is stationary. When the moving pendulum slams into the stationary one, its momentum is transferred to the other pendulum.

from visual import *
import numpy as np

# Constants
g = 9.80            # (m/s^2)
L = 10              # Length of the pendulums (m)
initialAngle = 1.2  # In radians

# Create the pendulum bob and rod
pend = sphere(pos=(L*np.sin(initialAngle), -L*np.cos(initialAngle), 0), radius=1, color=color.yellow)
rod = cylinder(pos=(0, 0, 0), axis=(pend.pos, pend.pos, 0), radius=0.1)
pend2 = sphere(pos=(-2, -L, 0), radius=1, color=color.red)
rod2 = cylinder(pos=(-2, 0, 0), axis=(pend2.pos+2, pend2.pos, 0), radius=0.1)

def position(right, t):
"""
Only one of the pendulums is in motion at a given time. This function
moves the moving pendulum to its new position. We use the equation:
theta(t) = theta_0*cos(sqrt(g/L)*t)
"""
theta = initialAngle*np.cos((g/L)**(1/2)*t)

if right:
pend.pos = [L * np.sin(theta), -L * np.cos(theta), 0]  # Update position of bob
rod.axis = [pend.pos, pend.pos, 0]  # Update rod's position
else:
pend2.pos = [L * np.sin(theta) - 2, -L * np.cos(theta), 0]  # Update position of bob
rod2.axis = [pend2.pos + 2, pend2.pos, 0]  # Update rod's position

# Once the moving pendulum reaches theta = 0, switch to the other one
if theta <= 0:
return False  # Return
else:
return True

# Increment time
i = 0
right = True  # The right pendulum is the first in motion
while True:
rate(200)
right = position(right, i)
i += 0.01

If you want to include damping, you can change the line

theta = initialAngle*np.cos((g/L)**(1/2)*t)

to

theta = initialAngle*exp(-mu*t)*np.cos((g/L)**(1/2)*t)

with some value for mu.