Commit c44e4d44 authored by ulrich_y's avatar ulrich_y
Browse files

Added new processes section

parent 96bcd867
Pipeline #1205 passed with stage
in 16 seconds
%!TEX root=manual
/*AUXFILES*/=ee_ee2eel.f95 ee_ee2eeg.f95
/*MAIN*/=$(group)_mat_el.f95 $(group).f95
include ../makefile.conf
/*all*/: $(group).a $(group).mod .obj/tree.sha
$(OBJ): ../.obj/functions.mod
.obj/$(group)_mat_el.o .obj/$(group)_mat_el.mod: \
$(group)_mat_el.f95 $(MOD)
.obj/$(group).o .obj/$(group).mod: \
$(group).f95 .obj/$(group)_mat_el.mod $(MOD)
$(group).mod: .obj/$(group).mod
/+cp+/ $< $@
@/+echo AR+/ $@
@$(AR) $@ $^
/+rm -f .obj/*.o .obj/*.mod .obj/*.gcda .obj/*.gcno *.mod *.a+/
%!TEX root=manual
use functions
use ee_ee2eel, only: ee2eel
use ee_ee2eeg, only: ee2eeg
implicit none
FUNCTION EE2EE(p1,p2,q1,q2)
!! e-(p1) e-(p2) -> e-(q1) e-(q2)
!! for massive electrons
FUNCTION EE2EEF(p1,p2,q1,q2)
!! e-(p1) e-(p2) -> e-(q1) e-(q2)
!! massive electrons
real(kind=prec) :: p1(4),p2(4),q1(4),q2(4)
real(kind=prec) :: ee2eef, mat0, Epart
Epart = sqrt(sq(p1+p2))
mat0 = ee2ee(p1,p2,q1,q2)
ee2eef = ee2eel(p1,p2,q1,q2) + alpha / (2 * pi) * mat0 * (&
- Ieik(xieik1,Epart,p1,p2) + Ieik(xieik1,Epart,p1,q1) &
+ Ieik(xieik1,Epart,p1,q2) + Ieik(xieik1,Epart,p2,q1) &
+ Ieik(xieik1,Epart,p2,q2) - Ieik(xieik1,Epart,q1,q2))
%!TEX root=manual
use functions
use phase_space, only: ksoft, ksoftA, ksoftB
use ee_mat_el
implicit none
FUNCTION EE2EEG_S(p1,p2,q1,q2)
!! e-(p1) e-(p2) --> e-(q1) e-(q2) y(ksoft)
!! both massive and massless electrons
real(kind=prec) :: p1(4),p2(4),q1(4),q2(4)
real(kind=prec) :: ee2eeg_s
ee2eeg_s = 2 * (- eik(p1,ksoft,p2) + eik(p1,ksoft,q1) &
+ eik(p1,ksoft,q2) + eik(p2,ksoft,q1) &
+ eik(p2,ksoft,q2) - eik(q1,ksoft,q2) &
) * ee2ee(p1,p2,q1,q2)
ee2eeg_s = (4.*pi*alpha) * ee2eeg_s
%!TEX root=thesis
FUNCTION EE2EE(p1, p2, p3, p4)
!! e-(p1) e-(p2) -> e-(p3) e-(p4)
!! for massive (and massless) electrons
implicit none
real(kind=prec), intent(in) :: p1(4), p2(4), p3(p4), p4(4), ee2ee
real(kind=prec) :: den1, den2, t, scms, m2
t = sq(p1-p3) ; scms = sq(p1+p2) ; m2 = sq(p1)
den1 = sq(p1-p3) ; den2 = sq(p1-p4)
ee2ee=(8**m2**2 - 8*m2*scms + 2*s**2 + 2*scms*t + t**2)/den1**2
ee2ee=ee2ee+2*(12*m2**2 - 8*m2*scms + scms**2) / den1 / den2
ee2ee=ee2ee+(24*m2**2 + scms**2 + t**2 - 8*m2*(s + t))/den2**2
ee2ee = ee2ee * 128*pi**2*alpha**2
%!TEX root=thesis
implicit none
real (kind=prec) :: x(2),y(5)
real (kind=prec) :: single, finite, lin
real (kind=prec) :: weight
integer ido
call blockstart("ee matrix elements")
scms = 40000.
musq = me
x = (/0.75,0.5/)
call ps_x2(x,scms,p1,me,p2,me,p3,me,p4,me,weight)
call check("ee2ee" ,ee2ee (p1,p2,p3,p4), 2.273983244890001e4, threshold=2e-8)
call check("ee2eel",ee2eel(p1,p2,p3,p4), 6.964297070440638e7, threshold=2e-8)
scms = 40000.
y = (/0.3,0.6,0.8,0.4,0.9/)
call ps_x3_fks(y,scms,p1,me,p2,me,p3,me,p4,me,p5,weight)
call check("ee2eeg",ee2eeg(p1,p2,p3,p4,p5),7.864297444955537e2, threshold=2e-8)
call blockend(3)
xinormcut1 = 0.2
xinormcut2 = 0.3
call blockstart("Moller VEGAS test")
call test_INT('ee2ee0', sigma_0, 2,10,10, ??)
call test_INT('ee2eeF', sigma_0, 2,10,10, ??)
call test_INT('ee2eeR', sigma_1, 5,10,10, ??)
call blockend(3)
%!TEX root=manual
onshell = {
p.p -> M^2, q.q -> m^2, p.q -> s/2
A0 = (4GF/Sqrt[2]) "diag1"/.RunQGraf[{"mum"},{"nu","elm"},0] //. {
line[_, x_] -> x, p1->p, q1->p-q, q2->q, _(*@\mydel@*)Z | (*@\mydel@*)m -> 0
A1 = pref /. RunQGraf[{"mum"},{"nu","elm"},1] //. {
line[_, x_] -> x, p1->p, q1->p-q, q2->q, _(*@\mydel@*)Z | (*@\mydel@*)m -> 0
1/2 Z2[m] Z2[M] FermionSpinSum[
A0 /. (*@\myP@*)L -> ((*@\myO@*) - Z5 (*@\myg@*)5)/2,
A0 /. (*@\myP@*)L -> ((*@\myO@*) + Z5 (*@\myg@*)5)/2
]] /. onshell]/.{
Z2[M_] -> 1 + ((*@\mya@*)/(4(*@\myp@*))) (-3/(2(*@\mye@*))-5/2 + 3/2 Log[M^2/Mu^2]),
Z5 -> 1 - ((*@\mya@*)/(4(*@\myp@*)))
1/2 FermionSpinSum[
A1/.(*@\myg@*).k1 -> (*@\myg@*). 4[k1]+I (*@\myg@*)5 (*@\mym@*),
] /. onshell /. {
(*@\mym@*)^n_ /; EvenQ[n] -> (*@\mym@*)2^(n/2), (*@\mym@*) -> 0
4[k1]. 4[k1] -> k1.k1 + (*@\mym@*)2, 4[k1] -> k1
M1bare = Simplify[KallenExpand[LoopRefine[LoopRelease[
Coefficient[M1, (*@\mym@*)2, 0]/(16 (*@\myp@*)^2)
+ (*@\mym@*)Integrate[
Coefficient[M1, (*@\mym@*)2, 1]/(64 (*@\myp@*)^3),
]]] /. e -> Sqrt[4 (*@\myp@*)(*@\mya@*)]];
......@@ -49,6 +49,7 @@ is based on~\cite{Ulrich:2020phd} for detailed background information.}
%!TEX root=manual
\section{Implementing new processes in \mcmule{}}
In this section we will discuss how new processes can be added to
\mcmule{}. Not all of the points below might be applicable to any
particular process. Further, all points are merely guidelines that
could be deviated from if necessary as long as proper precautions are
%!TEX root=manual
\subsection{Example calculations in Mathematica}
A thorough understanding of one-loop matrix elements is crucial for
any higher-order calculation. In \mcmule{}, one-loop matrix elements
either enter as the virtual contribution to \ac{NLO} corrections or
the real-virtual contribution in \ac{NNLO} calculations. In any case,
a fast numerical routine is required that computes the matrix element.
We perform all one-loop calculations in \fdf{} as this is arguably the
simplest scheme available. For theoretical background, we refer
to~\cite{Ulrich:2020phd} and references therein.
We use {\sc Qgraf} for the diagram generation. Using the in-house
Mathematica package {\tt qgraf} we convert {\sc Qgraf}'s output for
manipulation with Package-X~\cite{Patel:2015tea}. This package is
available on request through the \ac{MMCT}
An example calculation for the one-loop calculation of
$\mu\to\nu\bar\nu e\gamma$ can be found in Listing~\ref{lst:pkgx}. Of
course this example can be made more efficient by, for example,
feeding the minimal amount of algebra to the loop integration routine.
When using {\tt qgraf} for \fdf{} some attention needs to be paid
when considering diagrams with closed fermion loops. By default, {\tt
qgraf.wl} evaluates these traces in $d$ dimensions. {\tt RunQGraf} has
an option to keep this from happening.
\caption{An example on how to calculate the renormalised one-loop
matrix element for $\mu\to\nu\bar\nu e$ in \fdf.}
There is a subtlety here that only arise for complicated matrix
elements. Because the function Package-X uses for box integrals, {\tt
ScalarD0IR6}, is so complicated, no native Fortran implementation
exists in \mcmule{}. Instead, we are defaulting to
COLLIER~\cite{Denner:2016kdg} and should directly evaluate the finite
part of the {\tt PVD} function above. The same holds true for the
more complicated triangle functions. In fact, only the simple {\tt
DiscB} and {\tt ScalarC0IR6} are natively implemented without need for
external libraries. For any other functions, a judgement call is
necessary of whether one should {\tt LoopRefine} the finite part in
the first place. In general, if an integral can be written through
logarithms and dilogs of simple arguments (resulting in real answers)
or {\tt DiscB} and {\tt ScalarC0IR6}, it makes sense to do so.
Otherwise, it is often easier to directly link to COLLIER.
%!TEX root=manual
\caption{An example implementation of $\M n0$ for M{\o}ller
scattering. Note that the electron mass and the centre-of-mass energy
are calculated locally. A global factor of $8e^4=128\pi^2\alpha^2$ is
included at the end.}
As an example, we will discuss how M{\o}ller scattering $e^-e^-\to
e^-e^-$ could be implemented.
A new process group may need to be created if the process does not
fit any of the presently implemented groups. This requires a new
folder with a makefile as well as modifications to the main
makefile as discussed in Section~\ref{sec:newpg}.
In our case, $ee\to ee$ does not fit any of the groups, so we
create a new group that we shall call {\tt ee}.
Calculate the tree-level matrix elements needed at \ac{LO} and
\ac{NLO}: $\M{n}0$ and $\M{n+1}0$. This is relatively
straightforward and -- crucially -- unambiguous as both are finite
in $d=4$. We will come back to an example calculation in
A generic matrix element file is needed to store `simple' matrix
elements as well as importing more complicated matrix elements.
Usually, this file should not contain matrix elements that are
longer than a few dozen or so lines. In most cases, this applies
to $\M n0$.
After each matrix element, the \ac{PID} needs to be denoted in a
comment. Further, all required masses as well as the
centre-of-mass energy, called {\tt scms} to avoid collisions with
the function ${\tt s(pi,pj)}=2{\tt pi}\cdot{\tt pj}$, need to be
calculated in the matrix element to be as localised as possible.
In the case of M{\o}ller scattering, a file {\tt
ee/ee\_mat\_el.f95} will contain $\M n0$. For example, $\M n0$ is
implemented there as shown in Listing~\ref{lst:mollerlo}.
Further, we need an interface file that also contains the soft
limits. In our case this is called {\tt ee/ee.f95}.
Because $\M{n+1}0$ is border-line large, we will assume that it
will be stored in an extra file, {\tt ee/ee2eeg.f95}. The required
functions are to be imported in {\tt ee/ee\_mat\_el.f95}.
Calculate the one-loop virtual matrix element $\M n1$,
renormalised in the \ac{OS} scheme. Of course, this could be done
in any regularisation scheme. However, results in \mcmule{} shall
be in the \fdh{} (or equivalently the \fdf{}) scheme. Divergent
matrix elements in \mcmule{} are implemented as $c_{-1}$, $c_0$,
and $c_1$
\M n1 = \frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\Bigg(
\frac{c_{-1}}\epsilon + c_0 + c_1\epsilon+\mathcal{O}(\epsilon^2)
For $c_{-1}$ and $c_0$ this is equivalent to the conventions
employed by Package-X~\cite{Patel:2015tea} up to a factor
$1/16\pi^2$. While not strictly necessary, it is generally
advisable to also include $c_{-1}$ in the Fortran code.
For \ac{NLO} calculations, $c_1$ does not enter. However, we wish
to include M{\o}ller scattering up to \ac{NNLO} and hence will
need it sooner rather than later anyway.
In our case, we will create a file {\tt ee/ee\_ee2eel.f95}, which
defines a function
FUNCTION EE2EEl(p1, p2, p3, p4, sing, lin)
!! e-(p1) e-(p2) -> e-(p3) e-(p4)
!! for massive electrons
implicit none
real(kind=prec), intent(in) :: p1(4), p2(4), p3(p4), p4(4)
real(kind=prec) :: ee2eel
real(kind=prec), intent(out), optional :: sing, lin
The function shall return $c_0$ in {\tt ee2eel} and, if {\tt
present} $c_{-1}$ and $c_1$ in {\tt sing} and {\tt lin}.
At this stage, a new subroutine in the program {\tt test} with
reference values for all three matrix elements should be written
to test the Fortran implementation. This is done by generating a
few points using an appropriate phase-space routine and comparing
to as many digits as possible using the routine {\tt check}.
In our case, we would construct a subroutine {\tt TESTEEMATEL} as
shown in Listing~\ref{lst:mollertest}
Define a default observable in {\tt user} for this process. This
observable must be defined for any {\tt which\_piece} that might
have been defined and test all relevant features of the
implementation such as polarisation if applicable.
Add the matrix elements to the integrands defined in {\tt
integrands.f95} as discussed above. These integrands should be
added to {\tt mcmule.f95} and a second test routine should be
written that runs short integrations against a reference value.
Because {\tt test\_INT} uses a fixed random seed, this is expected
to be possible very precisely. Unfortunately,
COLLIER~\cite{Denner:2016kdg} might produce slightly different
results on different machines. Hence, integrands involving
complicated loop functions are only required to agree up to
After some short test runs, it should be clear whether new
phase-space routines are required. Add those, if need be, to {\tt
phase\_space} as described in Section~\ref{sec:ps}.
Per default the stringent soft cut, that may be required to
stabilise the numerical integration (cf.
Section~\ref{sec:fksfor}), is set to zero. Study what the smallest
value is that still permits integration.
Perform very precise $\xc$ independence studies. Tips on how to do
this can be found in Section~\ref{sec:xicut}.
At this stage, the \ac{NLO} calculation is complete and may, after
proper integration into \mcmule{} and adherence to coding style has
been confirmed, be added to the list of \mcmule{} processes in a new
release. Should \ac{NNLO} precision be required, the following steps
should be taken
Calculate the real-virtual and double-real matrix elements
$\M{n+1}1$ and $\M{n+2}0$ and add them to the test routines as
well as integrands.
Prepare the $n$-particle contribution $\sigma_n^{(2)}$. In a
pinch, massified results can be used also for $\ieik(\xc)\M n1$
though of course one should default to the fully massive results.
Study whether the pre-defined phase-space routines are sufficient.
Even if it was possible to use an old phase-space at \ac{NLO},
this might no longer work at \ac{NNLO} due to the added
complexity. Adapt and partition further if necessary, adding more
test integrations in the process.
Perform yet more detailed $\xc$ and soft cut analyses.
\caption{Test routine for $ee\to ee$ matrix elements and integrands.
The reference values for the integration are yet to be determined.}
In the following we comment on a few aspects of this procedure such as
the $\xc$ study (Section~\ref{sec:xicut}), the calculation of matrix
elements (Section~\ref{sec:matel}), and a brief style guide for
\mcmule{} code (Section~\ref{sec:style}).
%!TEX root=manual
\subsection{Creating a new process group}
Adding M{\o}ller scattering to \mcmule{}, the example discussed above,
requires the addition of a new process group {\tt ee}. For this we
create a new folder in \mcmule{} called {\tt ee} containing a
makefile (Listing~\ref{lst:eemakefile}), a {\tt mat\_el} file
({\tt ee/ee\_mat\_el.f95}, Listing~\ref{lst:eematel}) and a module
file ({\tt ee/ee.f95}, Listing~\ref{lst:eemod}). Finally, the name of
the process group needs to be added to the {\tt GROUPS} and {\tt
WGROUPS} variables of the makefile.
\caption{The bare makefile for the new process group {\tt ee}. Large
matrix elements that are stored in extra files such as {\tt
ee/ee2eeg.f95} or {\tt ee/ee\_ee2eel.f95} need to be added to the list
of {\tt AUXFILES}}
\caption{The file {\tt ee/ee\_mat\_el.f95} imports the complicated
matrix elements {\tt ee2eel} and {\tt ee2eegl}, defines the simple
matrix element {\tt ee2ee} as per \ref{lst:mollerlo}, and provides an
interface for the $\fM{n}1$ that is called from {\tt integrands}.}
\caption{The module file {\tt ee/ee.f95} which imports all matrix
elements of {\tt ee\_mat\_el} and defines the soft limits.}
%!TEX root=manual
\subsection{Coding style and best practice}
A large-scale code base like \mcmule{} cannot live without some basic
agreements regarding coding style and operational best practice. These
range from a (recommended but not enforced) style guide over the
management of the git repository to how to best run \mcmule{} in
development scenarios. All aspects have been discussed within the
Fortran code in \mcmule{} is (mostly) written in accordance with the
following style guide. If new code is added, compliance would be
appreciated but deviation is allowed if necessary. If in doubt,
contact any member of the \ac{MMCT}.
Indentation width is two spaces. In Vim this could be implemented
by adding the following to {\tt.vimrc}
autocmd FileType fortran set tabstop=8 softtabstop=0 expandtab shiftwidth=2 smarttab
Function and subroutine names are in all-upper case.
A function body is not indented beyond its definition.
When specifying floating point literals specify the precision when
possible, i.e. \lstinline{1._prec}.
Integrands should have \lstinline{ndim} specified.
Internal functions should be used where available.
Masses and other kinematic parameters must be calculated in the
matrix elements as local variables; using the global parameters
{\tt Mm} and {\tt Me} is strictly forbidden.
These rules also hold for matrix elements.
For python code, i.e. {\tt pymule} as well as the analysis code, PEP8
compliance is strongly encouraged with the exception of {\tt E231}
(Missing whitespace after {\tt ,}, {\tt;}, and {\tt:}), {\tt E731} (Do
not assign a lambda expression, use a {\tt def}) as well, in justified
cases, i.e. if required by the visual layout, {\tt E272} (Multiple
spaces before keyword), and {\tt E131} (Continuation line unaligned for
hanging indent).
\mcmule{} uses two git repositories for version management. One
internal repository and one public-facing one. Releasing to the latter
is the responsibility of the \ac{MMCT} after sufficient vetting was
performed by squashing commits to avoid the accidental release of
embarrassing or wrong code to the public. However, even the internal
repository has certain rules attached. In general, developers are
encouraged to not commit wrong or unvetted code though this can
obviously not be completely avoided in practice. To avoid
uncontrollable growth of the git repository, large files movements are
strongly discouraged. This also means that matrix elements should not
be completely overhauled barring unanimous agreement. Instead,
developers are encouraged to add a new matrix element file and link to
that instead.
Even when running \mcmule{} for development purposes the usage of menu
files is strongly encouraged because the code will do its utmost to
automatically document the run by storing the git version as well as
any modification thereof. This allows for easy and unique
reconstruction of what was running. For production runs this is not
optional; these must be conducted with menu files after which the run
folder must be stored with an analysis script and all data on the AFS
as well as the user file library to ensure data retention.
%!TEX root=manual
\subsection{Study of \texorpdfstring{$\xc$}{xic} dependence}\label{sec:xicut}
When performing calculations with \mcmule{}, we need to check that the
dependence of the unphysical $\xc$ parameter introduced in the FKS
scheme (cf. Appendix~\ref{sec:fks}) actually drops out at \ac{NLO} and \ac{NNLO}.
In principle it is sufficient to do this once during the development
phase. However, we consider it good practice to also do this (albeit
with a reduced range of $\xc$) for production runs.
Because the $\xc$ dependence is induced through terms as
$\xc^{-2\epsilon}/\epsilon$, we know the functional dependence
of $\sigma^{(\ell)}_{n+j}$. For example, at \ac{NLO} we have
\sigma^{(1)}_{n }(\xc) &= a_{0,0} + a_{0,1}\log(\xc)\,,\\
\sigma^{(1)}_{n+1}(\xc) &= a_{1,0} + a_{1,1}\log(\xc)\,,