#!/usr/bin/env python
# coding: utf-8

# <style>
# @import url(https://www.numfys.net/static/css/nbstyle.css);
# </style>
# <a href="https://www.numfys.net"><img class="logo" /></a>
# 
# # Electron Positron Annihilation
# 
# ### Examples - Particle and Nuclear Physics
# <section class="post-meta">
# By Jonas Tjemsland, Eilif Sommer Øyre and Jon Andreas Støvneng
# </section>
# Last edited: April 15th 2018 
# ___
# 
# In this notebook we will discuss electron positron scattering at the $Z$-resonance. In particular, we will use an event generator to run Monte Carlo simulations for the $e^++e^-\to Z\to ?$ annihilation process and then visualize the results as energy spectra. The event generator will take care of the cascading and hadronization. Physical experiments have been performed on the subject. Electron-positron scattering at the Z-resonance was studied by the ALEPH collaboration at the Large Electron-Positron collider (LEP) at CERN [1]. 
# 
# The event generator Pythia 8.2 [2, 3, 4] will be used. It is a standard tool for the generation of high-energy collisions, and is considered accurate for $>10\;\mathrm{GeV}$. The program works for hadron-hadron and lepton-lepton collitions. Pythia 8 is written in C++, but there exists a wrapper to Python.
# 
# In the following we always refer to the Standard Model. We start by importing necesarry packages. The installation of the Python wrapper for Pythia is described at the end of this notebook. The Pythia installation contains several examples you can play around with.

# In[1]:


# Import Pythia
import sys
cfg = open("Makefile.inc")
lib = ""
for line in cfg:
    if line.startswith("PREFIX_LIB="): lib = line[11:-1]; break
sys.path.insert(0, lib)
import pythia8

# Import other packages
import numpy as np
import matplotlib.pyplot as plt
import progressbar
get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


# Set common figure parameters
fontSize = 14
newparams = {'figure.figsize': (15, 6),
             'font.size': fontSize, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral',
             'lines.linewidth': 2.0}
plt.rcParams.update(newparams)

# Set constants
me = 0.000511 # GeV. Electron mass
mp = 0.93827  # GeV. Proton mass
mZ = 91.1876  # GeV. Z mass


# ## Setting up the process, $e^+ + e^- \rightarrow Z \rightarrow ?$
# 
# #### Settings
# 
# We will be simulating the scattering of a positron and an electron. That is, we will collide a positron and an electron at some center of mass energy $E_\mathrm{CM}$ and see what kind of particles is produced after hadronization and cascades.
# 
# We need to define a `Pythia`-object, which will handle the simulation of the process.

# In[3]:


# Initialise a pythia object
pythia = pythia8.Pythia()


# The settings for the event generator can be stored in a file and loaded using `pythia.readFile()`. However, here we will use `pythia.readString()` to read single settings.
# 
# The setting for the center of mass energy is called `Beams:eCM` and is given in units of $\mathrm{GeV}$. We choose that the particles collide at the $Z$-resonance, $E_\mathrm{CM}=91.1876\;\mathrm{GeV}$.
# 
# The inital particles are defined using the settings `Beams:idA` and `Beams:idB`, and is given by the [Monte Carlo Particle Numbering Scheme](http://pdg.lbl.gov/2017/reviews/rpp2017-rev-monte-carlo-numbering.pdf) [5]. In this numbering scheme, each particle is given an integer. Particles are given positive numbers and anti-particles are given negative numbers. For example, electrons is $11$ and positrons are $-11$. Some values are shown in the table below (see [5] for a complete list and the details of the scheme).
# 
# <center>**Table 1. ** *List of some elementary particles and their Monte Carlo Particle Numbering [5].*
# 
# | Quarks |    | Leptons    |    | Gauge Bosons |    | Hadrons |      |
# |:------:|:--:|:----------:|:--:|:------------:|:--:|:-------:|:----:|
# | $d$    | 1  | $e^-$      | 11 | $\gamma$     | 22 | $p$     | 2212 |
# | $u$    | 2  | $\nu_e$    | 12 | $Z^0$        | 23 | $n$     | 2112 |
# | $s$    | 3  | $\mu^-$    | 13 | $W^+$        | 24 | $\pi^0$ | 111  |
# | $c$    | 4  | $\nu_\mu$  | 14 | $h^0$        | 25 | $\pi^+$ | 211  |
# | $b$    | 5  | $\tau^-$   | 15 |              |    | $K^0_L$ | 130  |
# | $t$    | 6  | $\nu_\tau$ | 16 |              |    | $K^+$   | 321  |

# In[4]:


# Set collision properties
pythia.readString("Beams:idA = 11")
pythia.readString("Beams:idB = -11")
pythia.readString("Beams:eCM = 91.1876");


# In Pythia, particles with a (nomial) lifetime $\tau_0 < 10^3\;\mathrm{mm/}c$ decay by default. However, we will only consider stable particles. Thus, we let even the long-lived particles $\mu^{\pm}$, $\pi^\pm$, $K^\pm$, $K^0_L$ and $n$ decay. This is achieved using the `???:mayDecay` setting. The neutron $n$, for example, has an lifetime of $\tau= 880\;\mathrm{s}\approx 3\times10^8\;\mathrm{km}/c$, which in most applications and detectors can be treated as stable. However, when the neutrons have an astroparticle origin, they will decay before reaching Earth.

# In[5]:


pythia.readString("13:mayDecay   = true")
pythia.readString("211:mayDecay  = true")
pythia.readString("321:mayDecay  = true")
pythia.readString("130:mayDecay  = true")
pythia.readString("2112:mayDecay = true");


# #### Possible processes and interactions
# 
# Finally, we need to choose which processes to turn on. There are four fundamental forces in nature, shown in table 2.
# 
# <br></br>
# <center> **Table 2. ** *The four fundamental forces (taken from [6]).*</center>
#     
# | Force | Strength | Theory | Mediator |
# |:-----:|:--------:|:------:|:--------:|
# | Strong | 10 | Chromodynamics | Gluon, $g$ |
# | Electromagnetic | $10^{-2}$ | Electrodynamics | Photon, $\gamma$ |
# | Weak | $10^{-13}$ | Flavourdynamics | $W^\pm$ and $Z$ |
# | Gravitational | $10^{-42}$ | Geometrodynamics | Graviton |
# <br>
# 
# The initial state $e^+ + e^-$ is charge neutral. The electron and positron may thus decay into an intermediate photon (electromagnetism) or an intermediate $Z$ boson (weak interaction) (we neglect the coupling to the higgs). The corresponding fundamental vertices is shown is figure 2. A more complete list of vertices in Quantum Electrodynamics (QED), Quantum Chromodynamics (QCD) and Electroweak theory (GWS) is shown at the end of this notebook.
# 
# <br>
# ![Possible fundamental vertices](images/vertices_z_gamma.png)
# <center> **Figure 1. ** *Two of the fundemental vertices for a charge neutral electroweak interaction.*</center>
# 
# Note that the processes shown in figure 1 cannot be physical due to conservation of momentum. In order to get a physical process, two or more such vertices must be combined following a set of rules (e.g. charge and baryon number conservation for all interactions and color conservation in electromagnetic and weak interactions; see [6] for more information and deeper insight). This is known as a *Feynman diagram*. All elementary particle interactions can be described in this way. Note that the intermediate $\gamma$ and $Z$ are virtual particles, which means that they cannot be observed directly.
# 
# Each interaction vertex contribute with a factor $\propto \lambda$ to the final amplitude. This factor called the coupling constant is assumed to be small, such that the problem can be treated as a perturbation. That is, the diagrams with few vertices contribute the most. As an example, consider the process $e^++e^-\to e^++e^-$. This is known as *Babha scattering*. The first order Feynman diagrams for this process is shown in figure 2 for the photon. 
# 
# ![Babha scattering](images/bhabha.png)
# 
# <center> **Figure 2** *The first order Feynman diagrams to the Babha scattering. The upper diagram is called $t$-channel and the $s$-channel is shown on the bottom. The $u$-channel is not shown. How does the $u$-channel look like?*</center>
# 
# We are now ready to choose which processes to turn on. A complete list of processes and which settings they corresponds to is shown at the [Pythia 8.2 website](http://home.thep.lu.se/Pythia/pythia82html/ProcessSelection.html). As we have seen, the electron and positron can annihilate to either a $Z$ boson or a photon. That is, we must turn on some electroweak processes that takes two fermions to either $Z$ or $\gamma$. This is exactly what `WeakSingleBoson:ffbar2gmZ = on` does. Note that the photon branch will we suppressed at the $Z$-resonance, which is what we consider. We can check later that this is the case.

# In[6]:


pythia.readString("WeakSingleBoson:ffbar2ffbar(s:gmZ) = on");
pythia.readString("WeakSingleBoson:all = on");


# At last, we initialize the event generator.

# In[7]:


pythia.init();


# ## Running the simulations
# 
# We are now ready to run the simulations. We will use $10^5$ events in the Monte Carlo. This might be a bit much for our purposes, but it gives smooth energy spectra later on. The pythia class contains everything that is known about the current event (the $e^++e^-$ scattering). This includes e.g. initial, intermediate and final particles, the id of the particles and the four-momentum of the particles. A new event is generated using `pythia.next()`. For each event, we iterate through all the particles, search for a given set of particles ($\gamma$, $e^+$, $e^-$, $p$ and $\bar p$) and store the energy of the particle in an array, which later will be used to compute the energy spectrum.
# 
# **NOTE: ** When running the simulations, Pythia prints messages, used settings etc. in the terminal in which Jupyter was run. These messages can be turned off and additional can be turned on. Check them out!

# In[8]:


def run_simulations(pythia, iEvent):
    # First element is mass
    eGamma = [0]
    eE = [0.000511]
    eN = [.93957 ]
    eP = [.93828]
    eNu = [0] # We treat the neutrinos as massless
    eRest = []

    w = [progressbar.Bar('=', '[', ']'), ' ', progressbar.Percentage(), ' ', progressbar.ETA()]
    bar = progressbar.ProgressBar(widgets=w)

    for i in bar(range(iEvent)):
        # Generate next event. Skip if fail
        if not pythia.next(): continue

        # Iterate through the particles in the event. For each final particle,
        # store the energy in an array
        for i in range(pythia.event.size()):
            if pythia.event[i].isFinal():
                idAbs = pythia.event[i].idAbs() # Absolute id
                eI = pythia.event[i].e()        # Energy
                if idAbs == 22:
                    # PHOTON
                    eGamma.append(eI)
                elif idAbs == 2212:
                    # PROTON
                    eP.append(eI)
                elif idAbs == 2112:
                    # NEUTRON
                    eN.append(eI)
                elif idAbs == 11:
                    # e+e-
                    eE.append(eI)
                elif idAbs == 12 or idAbs == 14 or idAbs == 16:
                    # NEUTRINOS
                    eNu.append(eI)
                else:
                    eRest.append(eI)
    return eGamma, eE, eN, eP, eNu, eRest


# In[9]:


iEvent = 100000


# In[10]:


eGamma, eE, eN, eP, eNu, eRest = run_simulations(pythia, iEvent)


# ## Energy Spectrum
# 
# The final results from the simulations above are a set of arrays containing energies in GeV from the final particles in the events. The results must be transformed such that they become independent of the number of events used. The usual way to present the result is some kind of *energy spectra* such as
# $$\frac{\mathrm{d}N}{\mathrm{d}E}\; [1/\mathrm{GeV}],\quad \text{or}\quad T\times \frac{\mathrm{d}N}{\mathrm{d}T} \; [1/\mathrm{GeV}],$$
# where $N$ is the number density, $E$ is the energy and $T$ is the kinetic energy. We will use the latter. This is achieved by creating a histogram of the data and dividing by the number of events and the bin width of the histogram. A histogram can be created using the `histrogram()` function from `numpy`.

# In[11]:


def get_hist(E, bins):
    """ Creates energy spectrum given an array of energies and the bins. """
    if len(E) < 2: return np.array([])
    bin_width = bins[1:] - bins[:-1]
    T = np.array(E[1:]) - E[0]
    y, bin_edges = np.histogram(T, bins=bins)
    return y/(iEvent*bin_width)


# Let's create the histograms for the results. We should use logarithmically distributed bins.

# In[12]:


histNum = 100
bins = np.logspace(-5, 3, histNum)
bin_centers = 0.5*(bins[1:] + bins[:-1])

histGamma = get_hist(eGamma, bins)
histP = get_hist(eP, bins)
histNu = get_hist(eNu, bins)
histE = get_hist(eE, bins)
histN = get_hist(eN, bins)
if len(eRest) > 0:
    print("Additional particles were detected!")
else:
    print("All the different particles have been accounted for!")


# We can now plot the results.

# In[13]:


def plot_spectrum(bin_centers, histGamma, histP, histE, histNu):
    # Plot histograms loglog
    fig, ax = plt.subplots()

    if len(histGamma) > 0: ax.plot(bin_centers, bin_centers*histGamma, label=r"$\gamma$")
    if len(histP) > 0:     ax.plot(bin_centers, bin_centers*histP, label=r"$p\bar p$")
    if len(histE) > 0:     ax.plot(bin_centers, bin_centers*histE, label=r"$e^-e^+$")
    if len(histNu) > 0:    ax.plot(bin_centers, bin_centers*histNu, label=r"$\nu\bar\nu$")
    if len(histN) > 0:     ax.plot(bin_centers, bin_centers*histN, label=r"$n\bar n$")

    ax.set_xscale("log", nonposx='clip')
    ax.set_yscale("log", nonposy='clip')
    ax.legend(loc=2)
    plt.xlabel(r"$T$ [GeV]")
    plt.ylabel(r"Td$N$/d$T$")
    plt.show()


# In[14]:


plot_spectrum(bin_centers, histGamma, histP, histE, histNu)


# Note that all of the graphs goes, except the massless photons and the light neutrinos, goes to zero as the energy decreases. This change occurs rapidly. Moreover, when the energy increases, the energy spectrum decreases as well. There is also a sudden peak at $~5\times 10^1\;\mathrm{GeV}$ for the electron and neutrino spectrum.
# 
# Let's discuss some properties of the shapes of the graphs, in particular why the energy spectra are 0 for $T\gtrsim 50\;\mathrm{GeV}$.
# 
# ### Proton and antiproton $p\bar p$
# 
# The energy spectrum for the protons and anti-protons becomes zero at

# In[15]:


print("T = %.2f GeV." % (np.max(eP) - mp))


# The initial particles has an energy
# 
# $$E_i = 2m_e + E_\mathrm{CM} = 2\cdot 0.511\;\mathrm{MeV} + 91.1876\;\mathrm{GeV} \approx 91.2\;\mathrm{GeV}.$$
# 
# A proton and an anti-proton at rest has the energy
# 
# $$E_f = 2m_p = 2\cdot 0.93827231 \;\mathrm{GeV} \approx 1.88 \;\mathrm{GeV}.$$
# 
# If the process is $e^+ + e^- \to p + \bar p$, the kinetic energy of the final particles must be
# 
# $$T_p = T_\bar{p} = \frac{E_i - E_f}{2}=44.7 \;\mathrm{GeV}.$$
# 
# Thus, if we increase the number of events, we observe the energy spectrum for the protons and the anti-protons to become zero at $E_{0, 2}\to 44.7\;\mathrm{GeV}.$

# ### Electron and positron $e^-e^+$
# We repeate the same calculations as above:

# In[16]:


print("T = %.5f GeV" % (np.max(eE)))
print("T = %.5f GeV (analytical)" % ((2*me + mZ - 2*me)/2))


# Note also that the graph for $e^-e^+$ diverges for $T~45.6\;\mathrm{GeV}$. This corresponds to the process $e^+ + e^- \to Z \to e^+ + e^-$. The details are beyond the scope of this notebook.

# ### Photon and neutrino, $\gamma$ and $\nu$

# In[17]:


print("T = %.5f GeV" % (np.max(eGamma)))
print("T = %.5f GeV (analytical)" % ((2*me + mZ)/2))


# ### Further work and exercises
# 
# Play around with the code and concepts introduced. For instance, let the long lived particles be considered stable, try with other initial particles (remember to use the right processes) or try other initial energies.
# 
# - Where do other characteristics of the graphs such as the peak at $\sim 45\;\mathrm{GeV}$ and the maximum values come from?
# - Let the neutron be stable. Why are the energy spectrum for the neutron and proton similar? *Hint. Conservation law.*
# - Study the the Babha scattering in detail. For example, how are the angular distribution from Pythia? What do you obtain at tree-level in QED (similar calculations are found in e.g. [5])?
# - Electron-positron scattering at the Z-resonance was studied by the ALEPH collaboration at the Large Electron-Positron collider (LEP) at CERN. 
# 
# ### Resources and Further Reading
# 
# [1] Check out the ALEPH website (http://aleph.web.cern.ch/) or the CERN website (https://home.cern/about/experiments/aleph)
# 
# [2] Sjostrand, Torbjorn et al.: A Brief Introduction to PYTHIA 8.1. Comput.Phys.Commun. 178 (2008) 852-867 [arXiv:0710.3820](https://arxiv.org/abs/0710.3820) [hep-ph] CERN-LCGAPP-2007-04, LU-TP-07-28, FERMILAB-PUB-07-512-CD-T  
# 
# [3] Sjostrand, Torbjorn et al.: PYTHIA 6.4 Physics and Manual. JHEP 0605 (2006) 026 hep-ph/0603175 FERMILAB-PUB-06-052-CD-T, LU-TP-06-13
# 
# [4] Pythia 8.2 online manual: http://home.thep.lu.se/Pythia/pythia82php/Welcome.php (accessed March 2018)
# 
# [5] C. Patrignani et al. (Particle Data Group), Chin. Phys. C, 40, 100001 (2016) and 2017 update.
# 
# [6] Griffiths, David J. 2008. Introduction to elementary particles. Second, revised edition. Wiley-vch.  

# ---
# ## Appendix
# 
# ### Installation of Pythia and Setup for Python
# 
# Download [Pythia](http://home.thep.lu.se/~torbjorn/pythia8/pythia8230.tgz) and extract content. In Linux, you can use the following terminal commands to extract Pythia in the home folder:
# ```
# cd ~
# wget http://home.thep.lu.se/~torbjorn/pythia8/pythia8230.tgz
# tar -xvf pythia8230.tgz
# cd pythia8230
# ```
# Configure Python before configuring pythia.
# ```
# ./configure --with-python-include=/usr/include/python2.7 --with-python-bin=/usr/bin
# ```
# Now, configure Pythia (see http://home.thep.lu.se/Pythia/pythia82html/PythonInterface.html and http://home.thep.lu.se/Pythia/ for more settings) by running
# ```
# make
# ```
# Navigate to the folder containing the python file, `cd /folder/containing/python/file`. Then, copy the `Makefile.inc` file using
# ```
# cp ~/pythia8230/examples/Makefile.inc .
# ```
# For more information, see the information page for the [Python Interface](http://home.thep.lu.se/Pythia/pythia82html/PythonInterface.html).
# 
# ### Some Additional Fundamental Vertices
# Some of the fundamental vertices in the electroweak theory is shown in the figure below. The photon is charge neutral and couples to two of any electrically charged particle, such as the quarks or the electron. The $Z$ boson is also charge neutral, and couples to any lepton. The $W^\pm$ bosons are charged, and must therefore couple to two particles with a net charge. This is a lepton and corresponding neutrino, $(l^-, \nu_l)$, or a quark and a quark of oposite flavour such as the up and down quark. The flavour of the quarks is thus not conserved in electroweak theory. Quantum chromodynamics is not discussed in this notebook, but one can note that the gluon couples to any quark and that gluons can couple with other gluon in a three- or a four-vertex creating a so-called gluon ball. Gluons can also couple to photons.
# 
# ![vertices](images/vertices.png)