Commit 75e0f7e6 authored by janlukas.bosse's avatar janlukas.bosse
Browse files

asdfasdf asdf <

parents daa5ceb7 aa9bfbef
This diff is collapsed.
This diff is collapsed.
# Advanced Computational physics
Libraries and reference data needed for the Advanced Computational Physics Lab Course 2020
\ No newline at end of file
doi = {10.1016/0021-9991(78)90004-9},
url = {},
year = {1978},
month = may,
publisher = {Elsevier {BV}},
volume = {27},
number = {2},
pages = {192--203},
author = {G Peter Lepage},
title = {A new algorithm for adaptive multidimensional integration},
journal = {Journal of Computational Physics}
doi = {10.5281/ZENODO.3897199},
url = {},
author = {Lepage, Peter},
title = {gplepage/vegas: vegas version 3.4.5},
publisher = {Zenodo},
year = {2020},
copyright = {Open Access}
\ No newline at end of file
%!TEX root = report.tex
% ----------------------------------------------------------------------------
% Use the KomaScript styling on A4 paper
% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
% Use UTF-8 encoding
% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
% pretty fonts
% ----------------------------------------------------------------------------
\usepackage[sc]{mathpazo} % Use the Palatino font
\usepackage[T1]{fontenc} % Use 8-bit encoding that has 256 glyphs
\linespread{1.05} % Line spacing - Palatino needs more space between lines
\usepackage{microtype} % Slightly tweak font spacing for aesthetics
% ----------------------------------------------------------------------------
% Headings and spacing
% ----------------------------------------------------------------------------
% ----------------------------------------------------------------------------
% physics maths and code related things
% ----------------------------------------------------------------------------
\usepackage{amsmath, amsthm, amssymb, amsfonts}
\renewcommand{\vec}[1]{\boldsymbol{#1}} % bold vectors
% ----------------------------------------------------------------------------
% Figures etc
% ----------------------------------------------------------------------------
% whack shit to get the pngs in the pgf to be included correctly
% ----------------------------------------------------------------------------
% Citing and Referencing
% ----------------------------------------------------------------------------
\usepackage[backend=biber, style=phys, biblabel=brackets,
citestyle=numeric-comp, doi=false, sorting=none]{biblatex}
\usepackage[pdfstartview=FitH, % Oeffnen mit fit width
breaklinks=true, % Umbrueche in Links, nur bei pdflatex default
bookmarksopen=true, % aufgeklappte Bookmarks
bookmarksnumbered=true, % Kapitelnummerierung in bookmarks
colorlinks=true, urlcolor=blue, pdfborder={0 0 0},
% ----------------------------------------------------------------------------
% The title page
% ----------------------------------------------------------------------------
\title{Project Report --- Monte-Carlo Simulations for Particle Physics}
\author{Jan Lukas Bosse}
\Large \textbf{Advanced Computational Physics Lab}
%{\LARGE \bfseries
% \titlefont{Computational Many-Body Methods} \\[0.5cm]}
{\LARGE \bfseries
Project Report --- Monte-Carlo Simulations for Particle Physics \\[1.5cm]}
\textbf{\large Jan Lukas Bosse}\\[3mm]
% ----------------------------------------------------------------------------
% And the main part
% ----------------------------------------------------------------------------
\section{Simulation of the \(e^+ e^- \to q \overline{q}\)
process at fixed-order QED}
In this first part we integrate the leading order matrix element for
the \(e^+ e^- \to q \overline{q}\) process using Monte-Carlo integration. The
probability for this process to happen at a given scattering angle \(\theta\),
center of mass energy \(s\) and azimuthal angle \(\phi\) is proportional
to the differential cross-section
f_{q\overline{q}}(s, \cos \theta, \phi) :=
\frac{\dd{\sigma}_{q\overline{q}}}{\dd{s} \dd{\cos \theta} \dd{\phi}}
= p(s) \frac{1}{8\pi} \frac{1}{4 \pi} \frac{1}{2s}
|\mathcal{M}_{q\overline{q}}(s, \cos \theta, \phi)|^2
where \(p(s)\) is the distribution of the center of mass spectrum of the beam
and \(\mathcal{M}_{q\overline{q}}\) is the leading order matrix element
as given in the exercise sheet. To find the total cross-section
\(\sigma_{q\overline{q}}\) of the process
one has to integrate \(f_{q\overline{q}}\) over \(s, \theta\)
and \(\phi\). A closer look at \(\mathcal{M}_{q\overline{q}}\) shows that it
is not dependent on \(\varphi\), so we can skip the \(\varphi\)-integration
and replace it by multiplying the end-result with \(2 \pi\). Furthermore,
\(\cos \theta\) is distributed uniformly when integrating in spherical coordinates,
while \(\theta\) is not. So we will integrate over \(\cos \theta\) from
-1 to 1 instead of over \(\theta\). Besides \(s\) and \(\cos \theta\) the
integrand also depends on the quark flavour of the outgoing quark. Since
all outgoing flavours are equally likely and we integrate using Monte-Carlo
integration, we can simply pick a random flavor for each sample to achieve the
averaging over different quark flavors.
\caption{The integrands \(f_{q \overline{q}}(s, \cos \theta, 0)\) for all
possible outgoing quark flavours. Sadly, it is not easy to tell from this plot
which of the quark flavours corresponds to which of the surfaces. But the
general features are clear: Depending on the outgoing flavour \(q\) the
integrand takes one of two very similar forms. In either case
\(f_{q\overline{q}}\) has its highes values along the \(s=M_Z^2\) line.}
\todo{Plot the integrand here?}
\subsection{Integrating \(\dd{\sigma}_{q\overline{q}}\) at fixed \(s\)}
First, we calculate \(\sigma_{q\overline{q}}\) for \(p(s) = \delta(s - M_Z^2)\),
where \(M_Z\) is the mass of the Z-boson. Now the the \(s\)-dependence
of the integrating also drops out and we only have to perform a 1-dimensional
integral. The code for the Monte-Carlo integrator is found in
\texttt{integrator.MCIntegrator} and as a sampler we use a
\texttt{integrator.UniformSampler} in the range
\((\cos \theta, y) \in [-1,1] \times [0, 30]\) with \(N_{try} = 100,000\)
samples. After multiplication with
\(2\pi\) for the \(\phi\) integration and other conversion factors we find that
\sigma_{q\overline{q}}(s=M_Z^2) = 42510 \pm 200 \, \mathrm{pb}.
\subsection{Integrating over \(\cos \theta\) and \(s\)}
Although the integrand has its largest values at the \(s=M_Z^2\) line, one
should also integrate over \(s\) to get the correct result for the
cross-section. So now we integrate of \(\cos \theta\) and \(s\), the latter
in the interval \([(M_Z - 3 \Gamma_Z)^2, (M_Z + 3\Gamma_Z)^2]\). Since
the integrand is fastly decaying in \(s\) this range should suffice to
approximate the integral for \(s \in [0, \infty]\).
We use mutatis mutandis the same integrator and sampler as in the previous part.
Because we now perform a 2D integral instead of a 1D integral we use now
\(N_{try} = 1,000,000\) samples. As a result we now obtain
\sigma_{q\overline{q}} = 9950 \pm 100 \, \mathrm{pb}.
Little surprisingly this result is smaller than in the previous section. In
the previous section we chose \(p(s) = \frac{1}{2} \delta(s-M_Z^2)\), which is
supported exactly along the ridge of highest values of the integrand,
while in this section we had
\(p(s) = N \chi_{[(M_Z - 3 \Gamma_Z)^2, (M_Z + 3\Gamma_Z)^2]}(s)\), where
\(\chi_{[a,b]}(x)\) denotes the characteristic function on the interval \([a,b]\)
and \(N\) is a normalization constant.
Hence we sampled alot more smaller function values than in the previous section
and one expects the integral to be smaller.
\caption{Comparison of the acceptance probability as a function of \(s\)
with the marginal integral
\(\int f_{q\overline{q}}(s, \cos \theta, 0) \dd{\cos \theta}\). For comparison
we also show a Breit-Wigner distribution centenred at \(M_Z^2\) with scale
\(M_Z^2 \Gamma_Z^2\). All three
distributions are normalized such that the interal over them is 1.
To check consistency with the 1D integration from the previous section one can
also compute the marginal 1D integrals
\(\int f_{q\overline{q}}(s, \cos \theta, 0) \, \mathrm{d} \cos \theta\) for
sufficiently closely spaced \(s\)-values and then compare their values with the
acceptance probabilities \(p_{acc}(s)\) of sample points as a function of their \(s\)-value
in the 2D integration. This is done in \Cref{fig:histograms}.
For comparison we also show the Breit-Wigner distribution that we will
use later for importance sampling. We see that both functions agree
qualitatively in their overall shape. But \(p_{acc}(s)\) is more peaked than
the marginal integrals; More than could be explained by statistical fluctuations.
This hints, that there still must be some mistake in the code. See
\cref{sec:discussion} for more details.
\subsection{Monte-Carlo errors and number of samples}
One of the advantages of Monte-Carlo integration over e.g. quadrature rules is
that the error of the results scales as \(O\left(N_{try}^{-\frac{1}{2}}\right)\)
of the dimension of the integration domain while for said quadrature rule the
neccesary number of function evaluations scales exponentially in the dimension
of the integration domain\footnote{I know, apples and oranges... Still, this is
what people claim as the advntage of Monte Carlo integration.}.
We verify this claim by plotting the Monte-Carlo error estimate as a function of
samples \(N_{try}\) \Cref{fig:monte_carlo_errors}. Little surprisingly, it
follows the \(\frac{1}{\sqrt{N}}\) prediction almost perfectly.
\caption{Monte-Carlo error estimate for the 1D integral in
\cref{sec:1d_integration} (left) and for the 2D integral (right) as a function
of samples \(N_try\). For better comparison we also plot a
\(\frac{1}{\sqrt{N_{try}}}\) line in both cases.}
\subsection{Vegas integration}
One of the shortcomings of the naive Monte-Carlo integration method used in
\cref{sec:2d_uniform_integration} is that all \(x, y\) pairs get sampled with the
same probability. This means that most of the samples in regions where the
integrand is small get rejected thus decreasing the sampling efficiency and
increasing the Monte-Carlo error. One way around this is the Vegas algorithm
\cite{VegasAlgorithm}, which adapts its sampling grid based on where samples
got previously accepted. We use the \texttt{vegas} python implementation
\cite{vegas_python} here with 10,000 samples split into 10 iterations.
The result is
\sigma_{q\overline{q}} = 9930 \pm 20 \, \mathrm{pb}.
Again, this result agrees with the literature value given. More importantly,
the standard error is smaller than in \cref{sec:2d_uniform_integration}
despite us taking only one tenth of the samples.
This clearly shows the advantage of adapting the sampling algorithm to the integrand.
\caption{Sampling grid of the \texttt{vegas} algorithm after 10 iterations with
1,000 samples. The orange line is along \(s=M_Z^2\), the region where
the integrand is largest and over which we integrated in
\cref{sec:1d_integration}. Each rectangle in the grid has the same sampling
\begin{lstlisting}[frame=none, float,
caption={Summary of the \texttt{vegas} algorithm after
10 iterations},
captionpos=b, label=lst:vegas_log, xleftmargin=.15\textwidth]
itn integral wgt average chi2/dof Q
1 0.0002562(40) 0.0002562(40) 0.00 1.00
2 0.0002582(20) 0.0002578(18) 0.21 0.64
3 0.0002573(17) 0.0002575(13) 0.13 0.88
4 0.0002562(17) 0.0002570(10) 0.23 0.88
5 0.0002564(15) 0.00025682(83) 0.20 0.94
6 0.0002548(15) 0.00025635(73) 0.44 0.82
7 0.0002560(16) 0.00025629(66) 0.37 0.90
8 0.0002564(14) 0.00025631(60) 0.32 0.95
9 0.0002560(15) 0.00025627(56) 0.28 0.97
10 0.0002571(14) 0.00025639(52) 0.29 0.98
In \cref{fig:vegas_grid,lst:vegas_log} we show the Vegas integration grid
and summary after the 10 iterations. One can clearly see, that the grid
is denser around the maximal function values at the \(s=M_Z^2\) line. There
is also some non-trivial binning in the \(\cos \theta\) direction, namely
the values \(\cos \theta \approx \pm 1\) get sampled more often than
\(\cos \theta \approx 0\). Both make sense, when considering
\cref{fig:integrand}. In the first iteration, we get non-sensical \(Q\) and
\(\chi^2\)-values of 1.0 and 0.0 respectively, but the value of the integral is
already close to the final value. In the subsequent iterations value of the
integral does not change by much anymore. This hints, that the integrand was
"too simple" for vegas and it didn't need many samples to create an almost
perfectly adapted grid.
\subsection{Importance sampling}
If one already has a priori knowledge about the rough shape of the integrand
one can use importance sampling instead of the Vegas algorithm to increase
the sampling efficiency. The rough idea is to take more \(x\)-samples in the
region where one knows the integrand to be large and restrict the \(y\)-sampling
range according to the value of the integrand (plus some leeway to make sure
that one could also obtain sample points above the integrand). The
\(x\) and \(y\) samples have to be compatible that they fill the region around
the integrand uniformly.
Since in our case the integrand is dominated by the Breit-Wigner distribution
around \(s=M_Z^2\) we sample the \(s\) values from this distribution and
the \(\cos \theta\) values again uniformly in \([-1, 1]\). Such a sampler
is implemented in \texttt{integrator.BreitWignerSampler}. Running the importance
sampling algorithm with 10,000 samples (the same number as in
\cref{sec:2d_uniform_integration}) we obtain
\sigma_{q\overline{q}} = 11310 \pm 50 \, \mathrm{pb}.
Sadly, this is not the same result as in the previous parts. Also after lengthy
debugging sessions we couldn't find the reason why it isn't and eventually gave
up. More details are in \cref{sec:discussion}. Nevertheless, it is important
to note that the importance sampling algorithm increased the sampling
efficiency significantly, see \cref{tab:comparison}.
\subsection{Comparison of Results}
& & & & & & & \\
& & & & & & & \\
& & & & & & & \\
& & & & & & & \\
& & & & & & &
\caption{An empty table}
\subsection{Ideas for the next years}
\caption{Oh no!}
from vector import Vec4
from particle import Particle
from histogram import Histo1D, Scatter2D
import math as m
from utils.vector import Vec4
from utils.particle import Particle
from utils.histogram import Histo1D, Scatter2D
from itertools import combinations
from typing import List
# Why the fuck turn this into a class? Simply to share self.ecm2??
class Algorithm:
def Yij(self, p, q):
def Yij(self, p: Vec4, q: Vec4):
pq = p.px * q.px + * + p.pz * q.pz
return (2. * min(p.E, q.E)**2
* (1.0 - min(max(pq / m.sqrt(p.P2() * q.P2()), -1.0), 1.0))
/ self.ecm2)
def Cluster(self, event):
def Cluster(self, event: List[Particle]):
# TODO: implement the clustering algorithm here, and return a list of
# splitting scales Yij, ordered from smallest Yij to largest Yij
......@@ -25,7 +26,29 @@ class Algorithm:
# version of what we'd like plot). Instead, keep clustering until only
# two jets are left.
return []
self.ecm2 = sum( for e in event).M2()
ys = []
while (N := len(event)) > 2:
yij = {}
for i, j in combinations(range(N), 2):
yij[(i, j)] = self.Yij(event[i].mom, event[j].mom)
# find the minimum distance among the remaining particles
# and corresponding particles
dist = 1000000000.
i, j = 0, 0
for p1, p2 in combinations(range(N), 2):
if (newdist := yij[(p1, p2)]) < dist:
dist = newdist
i, j = p1, p2
# and combine them:
event[i].mom += event[j].mom
self.ecm2 = sum( for e in event).M2()
return ys
class Analysis:
import constants as const
import numpy as np
def chi1(s):
out = const.kappa * s * (s - const.M_Z**2)
out /= (s - const.M_Z**2)**2 + const.Gamma_Z**2 * const.M_Z**2
return out
def chi2(s):
out = const.kappa**2 * s**2
out /= (s - const.M_Z**2)**2 + const.Gamma_Z**2 * const.M_Z**2
return out
def eeqq_matrix_element(s : float, costheta : float, phi : float, q):
"""The squared ee -> qq matrix elemenent.
s : The centre of mass energy
costheta : The angle between the incoming e+ and the outgoing q
phi : The azimuthal angle of outgoing quark in the plane perpendicular
to the beam
q : one of ["u", "c", "d", "s", "b"], the flavour of the outgoing quarks
prefactor = (4 * np.pi * const.alpha)**2 * const.N_C
term1 = 1 + costheta**2
term1 *= (const.Q["e"]**2 * const.Q[q]**2
+ 2 * const.Q["e"] * const.Q[q] * const.V["e"] * const.V[q]
* chi1(s)
+ (const.A["e"]**2 + const.V["e"]**2)
* (const.A[q]**2 + const.V[q]**2)
* chi2(s)
term2 = np.copy(costheta)
term2 *= (4. * const.Q["e"] * const.Q[q] * const.A["e"] * const.A[q]
* chi1(s)
+ 8. * const.A["e"] * const.V["e"] * const.A[q] * const.V[q]
* chi2(s)
return prefactor * (term1 + term2)
def eeqq_mel_random_q(s, costheta, phi):
"Wrapper around `eeqq_matrix_element` that picks a random flavour"
q = np.random.choice(const.out_flavours)
return eeqq_matrix_element(s, costheta, phi, q)
\ No newline at end of file
import sys
import matplotlib.pyplot as plt
import numpy as np
class Bin1D:
def __init__(self, xmin, xmax):
self.xmin = xmin
self.xmax = xmax
self.w = 0.
self.w2 = 0.
self.wx = 0.
self.wx2 = 0.
self.n = 0.
def __repr__(self):
return str(self)
def __str__(self):
return "{0:10.6e}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5:10.6e}\t{6}".format(
self.xmin, self.xmax, self.w, self.w2, self.wx, self.wx2, int(self.n))
def Format(self, tag):
return "{0}\t{0}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5}".format(
tag, self.w, self.w2, self.wx, self.wx2, int(self.n))
def Fill(self, x, w):
self.w += w
self.w2 += w * w
self.wx += w * x
self.wx2 += w * w * x
self.n += 1.
def ScaleW(self, scale):
self.w *= scale
self.w2 *= scale * scale
self.wx *= scale
self.wx2 *= scale * scale
class Point2D:
def __init__(self, xmin, xmax):
self.x = (xmax + xmin) / 2
self.xerr = [(xmax - xmin) / 2] * 2
self.y = 0
self.yerr = [0] * 2
def __repr__(self):
return str(self)
def __str__(self):
return "{0:10.6e}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5:10.6e}".format(
self.x, self.xerr[0], self.xerr[1], self.y, self.yerr[0], self.yerr[1])
def ScaleY(self, scale):
self.y *= scale
self.yerr = [err * scale for err in self.yerr]
class Histo1D:
def __init__(self, nbin, xmin, xmax, name="/MC/untitled"): = name
self.bins = []
self.uflow = Bin1D(-sys.float_info.max, xmin)
width = (xmax - xmin) / nbin
for i in range(nbin):
self.bins.append(Bin1D(xmin + i * width, xmin + (i + 1) * width))
self.oflow = Bin1D(xmax, sys.float_info.max)