13. A Brief Introduction to Monte Carlo Event Generators in Particle Physics

13.1. Integrating for the Cross Section

Here, we will extend our Monte Carlo integration method to perform an integral to get the cross section for a specific process.

Remember that the cross section is defined as (See Mark Thomson, page 70):

\(\sigma=\) (number of particles per unit time per target particle)/(incident flux),

and an expression for \(2 \rightarrow 2\) scattering was derived, in the center-of-mass (C.O.M.) frame:

\(\sigma = \frac{1}{64 \pi^2 s}\frac{p_f^*}{p_i^*} \int |\mathcal{M}_{fi}|^2 \mathrm{d} \Omega^*\),

where the \(*\) indicates that the quantities are defined in the C.O.M. frame, and \(s\) is the Mandelstam variable.

All the physics of the process is encapsulated in tha matrix element \(\mathcal{M}_{fi}\). The rest of the expression is related to the “kinematics” of the process, i.e. energy and momentum conservation.

\(\frac{ \mathrm{d}{\sigma} } { \mathrm{d} \Omega^*}\) is known as a differential cross section, so in order to find the total cross section \(\sigma\), we need to perform an integral over the solid angle \(\Omega^*\). Recall that \(\mathrm{d} \Omega^* = \mathrm{d} \phi^* \mathrm{d} (\cos \theta^*)\) (i.e. the solid angle element is the angular part of the spherical polar coordinate volume element).

Generally speaking, integrals over differential cross sections will be complicated multi-differential and multi-dimensional functions of several variables, and therefore the Monte Carlo technique is the ideal method for solving them. In the first example we will consider here, the cross section can also be calculated analytically. However, we are going to use this as a first step towards understanding Monte Carlo event generators.

13.2. Monte Carlo Event Generators

This part of the tutorial is primarily based on the article https://arxiv.org/pdf/1412.4677.

In a similar way to the Monte Carlo integration that we have previously discussed, we can scan over a function \(f(x)\), and “collect” phase-space points, along with their associated “probabilities”, that are proportional to the weight of each in the calculation for the integral. These points effectively correspond to possible “events”, with their weights corresponding to their probability of occurring.

However, if we want to use these events, e.g. to perform an experimental analysis, then we must always carry the associated weight around for use in histograms, averages and so on. This can be inconvenient, but also very inefficient, as time may be wasted in some latter part of the simulation (e.g. a simulation of the effects of the detector) to events that possess only a very small weight.

The “hit-or-miss” method aims to equalize the weights of different events as far as possible.

Since the weight of each event is proportional to the probability of occurring, we can unweighted the events by keeping only a fraction of them, according to their weights. We do this by finding the maximum weight which occurs in the integration region, while we are performing the Monte Carlo integration.

Let’s do this for \(f(x) = \sqrt{1-x^2}\):

import matplotlib.pyplot as plt # import matplotlib, a conventional module name is plt
import numpy as np
from math import pi

x = np.linspace(-1, 1, 200) # creates a NumPy array from -1 to 1, 200 equallys-paced points  
y = np.sqrt(1-x**2) # take the NumPy array and create another one, where each term is now the sine of each of the elements of the above NumPy array

fig, ax = plt.subplots() # create the elements required for matplotlib. This creates a figure containing a single axes.

# set the labels and titles:
ax.set_xlabel(r'$x$', fontsize=20) # set the x label
ax.set_ylabel(r'$f(x)$', fontsize=20) # set the y label. Note that the 'r' is necessary to remove the need for double slashes. You can use LaTeX! 

# make a one-dimensional plot using the above arrays, add a custom label
ax.plot(x, y, label=r'$\sqrt{1-x^2}$') 

# construct the legend:
ax.legend(loc='upper right')  # Add a legend

plt.show() # show the plot here

First let’s write the Monte Carlo integration function:

import math
import random

# Let's define a function that performs one-dimensional MC integration of an arbitrary for N points
# We'll also make it capable of performing the integration in an interval a,b
def mcint(func, a, b, N):
    """Calculates the one-dimensional Monte Carlo integral of func in [a,b] for N points"""
    sumw = 0 # we will use this variable for the sum of f(x_i)
    sumwsq = 0 # and this one for the sum of f(x_i)^2, used in the error calculation
    for i in range(int(N)):
        xi = (b-a) * random.random() + a
        wi = (b-a)*func(xi)
        sumw += wi
        sumwsq += wi**2 
    # now calculate the average value of f (i.e. the integral):
    I = sumw/N
    # and the error: 
    sigmaIsq = (1/N) * ( (1/N) * sumwsq - I**2 ) # this is the variance (i.e. the error squared)
    sigmaI = math.sqrt(sigmaIsq) # this is the actual error
    return I, sigmaI # return the integral and its error

And apply it:

# Now let's also define the function that we wish to integrate:
def f(x):
    return np.sqrt(1-x**2)

# and use our mcint() function to integrate it in [0,1]:
N = 100000 # with 10000 points we expect O(1/sqrt(N)) ~ O(1%) error
Int, Err = mcint(f,-1,1,N) # this way you can access both the integral and its error

# print the results:
print("Analytic value of the integral=", math.pi/2)
print("The MC integral of f(x) = sqrt(1-x^2) in [-1,1] is", Int, "+-", Err)
print("The fractional uncertainty is", Err/Int)
print("The fractional difference from the analytical value is", abs(Int-math.pi/2)/(math.pi/2))

Analytic value of the integral= 1.5707963267948966
The MC integral of f(x) = sqrt(1-x^2) in [-1,1] is 1.569254005501887 +- 0.0014145090696274995
The fractional uncertainty is 0.0009013894912284157
The fractional difference from the analytical value is 0.0009818722304733582

# Let's define a function that performs one-dimensional MC integration of an arbitrary for N points
# We'll also make it capable of performing the integration in an interval a,b
def mcint(func, a, b, N):
    """Calculates the one-dimensional Monte Carlo integral of func in [a,b] for N points"""
    sumw = 0 # we will use this variable for the sum of f(x_i)
    sumwsq = 0 # and this one for the sum of f(x_i)^2, used in the error calculation
    maxw = -1E99 # the maximum weight found during the integration
    xmaxw = -1E99 # the location of the maximum weight 
    for i in range(int(N)):
        xi = (b-a) * random.random() + a
        wi = (b-a)*func(xi)
        sumw += wi
        sumwsq += wi**2 
        if wi > maxw:
            maxw = wi
            xmaxw = xi
    # now calculate the average value of f (i.e. the integral):
    I = sumw/N
    # and the error: 
    sigmaIsq = (1/N) * ( (1/N) * sumwsq - I**2 ) # this is the variance (i.e. the error squared)
    sigmaI = math.sqrt(sigmaIsq) # this is the actual error
    return I, sigmaI, maxw, xmaxw # return the integral and its error and the maximum weight and its location

# Use the function again to get the integral and the maximum weight: 
N = 100000 # with 10000 points we expect O(1/sqrt(N)) ~ O(1%) error
Int, Err, MaxW, xMaxW = mcint(f,-1,1,N) # this way you can access both the integral and its error

print('Maximum weight found=', MaxW, 'at', xMaxW)

Maximum weight found= 1.999999999840376 at 1.263423735009539e-05

Let’s now see how to generate a sample of \(x\)’s distributed according to \(f(x) = \sqrt{1-x^2}\) (appropriately normalized).

To proceed, we choose to keep (“accept”) each \(x\) (i.e. an “event”), with probability \(w(x)/w_\mathrm{max})\), where \(w(x)\) is the weight of the point and \(w_\mathrm{max})\). The rest are thrown away (“rejected”). All accepted events are given a weight equal to the average weight \(<w>\), used to calculate the integral.

Let’s try this with our given function, and compare to the expected (normalized) function.

# event generation function according to func, given the number of events to be generated Ngen in the interval a, b,
# and number of integration events Nint
def evtgen(func, a, b, Ngen, Nint):
    # first find the integral using the integration function: 
    Int, Err, wmax, xwmax = mcint(func,a,b,Nint) 
    # define a counter to check when to terminate the event generation:
    counter = 0
    # an empty list to hold the Ngen events
    eventslist = []
    while counter < Ngen:
        # get a random point and calculate its weight: 
        xi = (b-a) * random.random() + a
        wi = (b-a)*func(xi)
        if wi/wmax > random.random():
            # accept the event, increment counter
            eventslist.append(xi)    
            counter += 1
        else:
            # reject the event and check another one
            continue 
    # return the integral, the eror and the events list
    return Int, Err, eventslist
    

Let’s try it out!

# integral number of points:
Nintegral = 100000
# number of events to be generated: 
Nevents = 10000
sigma, sigmaerr, events = evtgen(f,-1,1,Nevents,Nintegral)

# check that the number of events generated is as requested:
print('Number of events generated=', len(events))

Number of events generated= 10000

And now let’s compare the generated \(x\)’s to the normalized distribution:

import matplotlib.pyplot as plt # import matplotlib, a conventional module name is plt
import numpy as np

fig, ax = plt.subplots() # create the elements required for matplotlib. This creates a figure containing a single axes.

# set the labels and titles:
ax.set_xlabel(r'$x$', fontsize=20) # set the x label
ax.set_ylabel(r'frequency', fontsize=20) # set the y label
ax.set_title(r'$f(x)=\sqrt{1-x^2}$', fontsize=10) # set the title 

# make one-dimensional plots using the above arrays, add a custom label, linestyles and colors:
ax.hist(events, color='blue', bins=25, density=True, fill=False, histtype='step', label="hit-or-miss") 

# plot the distribution itself:
xv = np.linspace(-1, 1,1000)
gaussianpoints = f(xv) / (np.pi/2) # pi/2 is the analytic integral 

ax.plot(xv, gaussianpoints, ls='--', color='red', label='analytic')

# construct the legend:
ax.legend(loc='upper right')  # Add a legend

plt.show() # show the plot here

We are now ready to proceed to the particle physics case!

13.3. A Monte Carlo Event Generator for \(e^+ e^- \rightarrow \mu^+ \mu^-\)

The differential cross section for the process \(e^+ e^- \rightarrow \gamma \rightarrow \mu^+ \mu^-\) is given by:

\(\frac{\mathrm{d}\sigma}{\mathrm{d} \Omega} = \frac{ \alpha^2 } { 4s} (1+ \cos^2\theta)\),

where \(\alpha = e^2/(4\pi)\) is the QED running coupling.

Since the expression does not depend on the angle \(\phi\), we may integrate over it, introducing a multiplicative factor of \(2\pi\) on the right-hand side.

The integration to obtain the cross section is in fact trivial, since we know how to integrate cosine functions analytically, and the \(e^+ e^-\) center-of-mass energy squared, \(s\) is fixed. The total cross section in this case is given by:

\(\sigma = \frac{ 4 \pi \alpha^2 } { 3 s}\).

As a reminder, the conversion factor from GeV to pb is given by: \(3.894 \times 10^8 ~\mathrm{pb/GeV}^{−2}\).

Despite the fact that this exaple is rather easy to handle analytically, we wish begin our venture into event generation with this simple exercise, as it provides an insight to the basic building blocks of an event generator.

Let’s first calculate the total cross section analytically, using the integration routines that we have developed thus far.

Let’s start by copying the MC integration routine, which also finds the maximum.

Let’s define the function we wish to integrate, i.e. the differential cross section:

# conversion factor from GeV^-2 to pb:
pbconvert=3.894E8
# this is the QED running coupling at 90 GeV:
alpha = 1/132.507 

# Now let's also define the function that we wish to integrate:
def f(x,s=30**2): # x = costheta, and set the COM energy to 90 GeV
    return 2*math.pi*((alpha**2)/(4 * s))*(1+x**2)

And integrate!

# fix the COM energy globally
ECM = 30 # in GeV
s=ECM**2

# Use the function again to get the integral and the maximum weight: 
N = 100000 # with 10000 points we expect O(1/sqrt(N)) ~ O(1%) error
Int, Err, MaxW, xMaxW = mcint(f,-1,1,N) # this way you can access both the integral and its error
print('The cross section in pb is:', Int*pbconvert, '+-', Err*pbconvert)
print('Maximum weight found=', MaxW, 'at', xMaxW)

The cross section in pb is: 103.20920489046665 +- 0.07318641769311837
Maximum weight found= 3.976083870415857e-07 at 0.9999903771364542

Let’s also calculate the analytical expectation for the total cross section:

# analytical expression for the cross section for e+e- -> gamma -> mu+mu-
def sigma_analytical():
    return 4 * math.pi * alpha**2 / 3 / s 

sigma_analyt = sigma_analytical() 
sigma_analytical_pb = sigma_analyt*pbconvert
print("The analytical cross sectin for e+e- -> gamma -> mu+mu- at 30 GeV is", sigma_analytical_pb)

The analytical cross sectin for e+e- -> gamma -> mu+mu- at 30 GeV is 103.22013054444798

Evidently these two are compatible!

We are now ready to proceed to the event generation stage. This will consist simply of generating the \(\cos \theta\) values and then “reconstructing” the momenta according to those values, for each event. The angle \(\phi\) can be chosen uniformly and randomly in the interval \((0,2 \pi\)). Let’s start by generating a list of \(\cos \theta\), distributed according to our differential cross section. This list is essentially the list of events.

Let’s use the evtgen function from the introduction to this section!

# integral number of points:
Nintegral = 100000
# number of events to be generated: 
Nevents = 10000
sigma, sigmaerr, events = evtgen(f,-1,1,Nevents,Nintegral)

# check that the number of events generated is as requested:
print('Number of events generated=', len(events))

Number of events generated= 10000

Let’s check whether it corresponds to what we expect! In this case we actually can also plot the analytical form, as we did earlier.

import matplotlib.pyplot as plt # import matplotlib, a conventional module name is plt
import numpy as np

fig, ax = plt.subplots() # create the elements required for matplotlib. This creates a figure containing a single axes.

# set the labels and titles:
ax.set_xlabel(r'$\cos\theta$', fontsize=20) # set the x label
ax.set_ylabel(r'frequency', fontsize=20) # set the y label
ax.set_title(r'Differential cross section for $e^+ e^- \rightarrow \gamma \rightarrow \mu^+ \mu^-$', fontsize=10) # set the title 

# make one-dimensional plots using the above arrays, add a custom label, linestyles and colors:
ax.hist(events, color='blue', bins=25, density=True, fill=False, histtype='step', label="hit-or-miss") 

# plot the distribution itself:
xv = np.linspace(-1, 1,1000)
analytic = f(xv) / sigma_analyt# pi/2 is the analytic integral 

ax.plot(xv, analytic, ls='--', color='red', label='analytic')

# construct the legend:
ax.legend(loc='upper right')  # Add a legend

plt.show() # show the plot here

Therefore, our events have indeed the right distribution: they follow the differential cross section, as expected!

Let’s now proceed to set up the momenta of the particles in each event. This is relatively easy in this case, since we are actually working directly in the center of mass frame.

If we assume that the incoming beans are traveling in the positive and negative \(z\)-directions, e.g. \(e^-\) in the positive and \(e^+\) in the negative, then the incoming particle four-momenta can be written as:

\(p_{e^-} = \frac{\sqrt{s}}{2} ( 1, 0, 0, 1)\),

\(p_{e^+} = \frac{\sqrt{s}}{2} ( 1, 0, 0, -1)\),

and the outgoing four-momenta according to the generated \((\theta, \phi)\):

\(p_{\mu^-} = \frac{\sqrt{s}}{2} ( 1, \sin \theta \cos \phi, \sin \theta \sin\phi, \cos\theta)\),

\(p_{\mu^+} = \frac{\sqrt{s}}{2} ( 1, -\sin \theta \cos \phi, -\sin \theta \sin\phi, -\cos\theta)\),

where one can clearly observe that energy and momentum are conserved.

Let’s now “reconstruct” (i.e. build) the four-vectors in each event according to the list of \(\cos \theta\) that we have generated.

# takes a list of cos thetas and reconstructs the four-momenta
# for the process e- e+ -> mu- mu+ in the COM
# returns a list of four-vectors, where a four-vector is a list of four numbers
def reconstruct(costhetas):
    recoevents = []
    sqrtso2 = math.sqrt(s)/2 # appears often so no need to recalc.
    for costh in costhetas:
        # generate a random phi:
        phi = random.random()*2*math.pi
        # get sintheta as well:
        sinth = math.sqrt(1-costh**2)
        pem = [sqrtso2, 0, 0, sqrtso2]
        pep = [sqrtso2, 0, 0, -sqrtso2]
        pmm = [sqrtso2, sqrtso2*sinth * math.cos(phi), sqrtso2*sinth * math.sin(phi), sqrtso2*costh]
        pmp = [sqrtso2, -sqrtso2*sinth * math.cos(phi), -sqrtso2*sinth * math.sin(phi), sqrtso2*costh]
        recoevents.append([pem, pep, pmm, pmp])
    return recoevents

# get the reconstructed events: 
reconstructed_events = reconstruct(events)
    

Let’s see an example of an event, e.g. the first one:

for p in reconstructed_events[0]:
    print('p=',p)

p= [15.0, 0, 0, 15.0]
p= [15.0, 0, 0, -15.0]
p= [15.0, 0.843942317719295, 14.957281930476524, -0.7533117658732513]
p= [15.0, -0.843942317719295, -14.957281930476524, -0.7533117658732513]

Seems pretty reasonable!

Since there is effectively only one degree of freedom, there’s not many variables that are interesting in this simple case. We can imagine reconstructing the ``transverse momentum’’, the magnitude of the momentum in the \(x\) and \(y\) directions, defined as:

\(p_T = \sqrt{p_x^2 + p_y^2}\)

Let’s calculate this per event for the outgoing \(\mu^-\) and plot it!

pT_list = [] # the list that will contain all the transverse momenta
for ev in reconstructed_events:
    pmm = ev[2] # the mu- is the third entry in the event in this case
    pT = math.sqrt(pmm[1]**2 + pmm[2]**2)
    pT_list.append(pT)

# PLOT! 
fig, ax = plt.subplots() # create the elements required for matplotlib. This creates a figure containing a single axes.

# set the labels and titles:
ax.set_xlabel(r'$p_T$', fontsize=20) # set the x label
ax.set_ylabel(r'frequency', fontsize=20) # set the y label
ax.set_title(r'Transverse momentum of the muon in $e^+ e^- \rightarrow \gamma \rightarrow \mu^+ \mu^-$', fontsize=10) # set the title 

# make one-dimensional plots using the above arrays, add a custom label, linestyles and colors:
ax.hist(pT_list, color='blue', bins=25, density=True, fill=False, histtype='step', label="hit-or-miss") 

# construct the legend:
ax.legend(loc='upper right')  # Add a legend

plt.show() # show the plot here

An interesting next step is to consider the full process for \(e^+ e^- \rightarrow \mu^+ \mu^-\), which includes a contribution from the \(Z\) boson. The main difference in the way that the \(Z\) boson interacts with the leptons is that it couples with different strengths to the left- and right-handed fermions.

See the details discussed in secton 3.1.2 of https://arxiv.org/pdf/1412.4677. Following these, you should be able to address section 3.2, Exercise 1.