#### Automated one-loop calculations with GoSam

Gavin Cullen
1
2
Nicolas Greiner
0
6
Gudrun Heinrich
0
Gionata Luisoni
5
Pierpaolo Mastrolia
0
4
Giovanni Ossola
3
8
Thomas Reiter
0
g
Francesco Tramontano
7
h
0
Max-Planck-Institut fr Physik
, Mnchen,
Germany
1
School of Physics and Astronomy, The University of Edinburgh
,
Edinburgh, UK
2
DESY, Zeuthen,
Germany
3
New York City College of Technology, City University of New York
,
New York, USA
4
Dipartimento di Fisica,
Universit di Padova
, Padova,
Italy
5
Institute for Particle Physics Phenomenology, University of Durham
,
Durham, UK
6
Department of Physics, University of Illinois at Urbana-Champaign
, Urbana,
USA
7
CERN,
Geneva, Switzerland
8
The Graduate School and University Center, City University of New York
,
New York, USA
1 Introduction
The Standard Model is currently being re-discovered at the
LHC, and new exclusion limits on Beyond the Standard
Model particlesand on the Higgs massare being
delivered by the experimental collaborations with an impressive
speed. Higher order corrections play an important role in
obtaining bounds on the Higgs boson and New Physics. In
particular, the exclusion limits for the Higgs boson would look
very different if we only had leading order tools at hand.
Further, it will be very important to have precise theory
predictions to constrain model parameters once a signal of New
Physics has been established. Therefore it is of major
importance to provide tools for next-to-leading order (NLO)
predictions which are largely automated, such that signal and
background rates for a multitude of processes can be
estimated reliably.
The need for an automation of NLO calculations has been
noticed some time ago and lead to public programs like
FeynArts [1] and QGraf [2] for diagram generation and
FormCalc/LoopTools [3] and GRACE [4] for the automated
calculation of NLO corrections, primarily in the electroweak
sector. However, the calculation of one-loop amplitudes with
more than four external legs were still tedious case-by-case
calculations. Only very recently, conceptual and technical
advances in multi-leg one-loop calculations allowed the
calculation of six-point [524] and even seven-point [25, 26]
processes at all, and opened the door to the possibility of
an automated generation and evaluation of multi-leg
oneloop amplitudes. As a consequence, already existing
excellent public tools, each containing a collection of hard-coded
individual processes, like e.g. MCFM [27, 28], VBFNLO
[29, 30], MC@NLO [31, 32], POWHEG-Box [33, 34],
POWHEL [3537], can be flanked by flexible automated
tools such that basically any process which may turn out to
be important for the comparison of LHC findings to theory
can be evaluated at NLO accuracy.
We have recently experienced major advances in the
activity of constructing packages for fully automated
oneloop calculations, see e.g. [3843]. The concepts that
lead to these advances have been recently reviewed in
[44]. Among the most important developments are the
integrand-reduction technique [45, 46] and the
generalized n-dimensional unitarity [47]. Their main outcome is
a numerical reconstruction of a representation of the tensor
structure of any one-loop integrand where the multi-particle
pole configuration is manifest. As a consequence,
decomposing one-loop amplitudes in terms of basic integrals
becomes equivalent to reconstructing the polynomial forms of
the residues to all multi-particle cuts. Within this algorithm,
the integrand of a given scattering amplitude, carrying
complete and explicit information on the chosen
dimensionalregularisation scheme, is the only input required to
accomplish the task of its evaluation. In fact, the integration is
substituted by a much simpler operation, namely by
polynomial fitting, which requires the sampling of the integrand on
the solutions of generalised on-shell conditions.
In this article, we present the program package GOSAM
which allows the automated calculation of one-loop
amplitudes for multi-particle processes. Amplitudes are expressed
in terms of Feynman diagrams, where the integrand is
generated analytically using QGRAF [2], FORM [48],
spinney [49] and haggies [50]. The individual program tasks
are steered via python scripts, while the user only needs to
edit an input card to specify the details of the process to
be calculated, and launch the generation of the source code
and its compilation, without having to worry about internal
details of the code generation.
The program offers the option to use different
reduction techniques: either the unitarity-based integrand
reduction as implemented in SAMURAI [40] or traditional tensor
reduction as implemented in Golem95C [51, 52] interfaced
through tensorial reconstruction at the integrand level [53],
or a combination of both. It can be used to calculate
oneloop corrections within both QCD and electroweak theory.
Beyond the Standard Model theories can be interfaced using
FeynRules [54] or LanHEP [55]. The Binoth-Les
Houchesinterface [56] to programs providing the real radiation
contributions is also included.
The advantage of generating analytic expressions for the
integrand of each diagram gives the user the flexibility to
organize the computation according to his own efficiency
preferences. For instance, the computing algorithm can
proceed either diagram-by-diagram or by grouping diagrams
that share a common set of denominators (suitable for a
unitarity-based reduction), and it can deal with the
evaluation of the rational terms either on the same footing as
the rest of the amplitude, or through an independent
routine which evaluates them analytically. These options and
the other features of GOSAM will be discussed in detail in
the following.
In Sect. 2, after giving an overview on the diagram
generation and on processing gauge-group and Lorentz algebra,
we discuss the code generation and the reduction strategies.
The installation requirements are given in Sect. 3, while
Sect. 4 describes the usage of GOSAM, containing all the
set-up options which can be activated by editing the input
card. In Sect. 5 we show results for processes of various
complexity. The release of GOSAM is accompanied by the
generated code for some example processes, listed in
Appendix A.
2 Overview and algorithms
2.1 Overview
GOSAM produces, in a fully automated way, all the code
required to perform the calculation of one-loop matrix
elements. There are three main steps in the process of
constructing the code: the generation of all contributing
diagrams within a process directory, the generation of the
Fortran code, and finally compiling and linking the
generated code. These steps are self-contained in the sense that
after each step all the files contained in the process directory
could be transfered to a different machine where the next
step will be carried out.
In the following sections we focus on the algorithms that
are employed for the construction of the code to produce and
evaluate matrix elements.
The first step (setting up a process directory), which
consists in the generation of some general source files and the
generation of the diagrams, is described in Sect. 2.2. The
second step (generating the fortran code) is carried out by
means of advanced algorithms for algebraic manipulation
and code optimization which are presented in Sects. 2.3
and 2.4. The third step (compilation and linking) is not
specific to our code generation, therefore will not be described
here.
The practical procedures to be followed by the user in
generating the code will be given in Sect. 4, which can be
considered a short version of the user manual.
2.2 Generation and organisation of the diagrams
For the diagram generation both at tree level and one-loop
level we employ the program QGRAF [2]. This program
already offers several ways of excluding unwanted diagrams,
for example by requesting a certain number of propagators
or vertices of a certain type or by specifying topological
properties such as the presence of tadpoles or on-shell
propagators. Although QGRAF is a very reliable and fast
generator, we extend its possibilities by adding another level of
analysing and filtering over diagrams by means of Python.
This gives several advantages: first of all, the possibilities
offered by QGRAF are not always sufficient to distinguish
certain classes of diagrams (see examples in Fig. 1); secondly,
Fig. 1 Two examples for diagrams which are difficult to isolate using
QGRAF. The diagram in Fig. 1(a) is zero in dimensional regularisation.
However, in QGRAF there is no operator to identify this type of
diagrams. In Fig. 1(b) the Z boson is emitted from a closed quark line.
These diagrams form a separate gauge invariant class and could be
treated separately from diagrams where the Z boson comes from an
external quark line
QGRAF cannot handle the sign for diagrams with Majorana
fermions in a reliable way; finally, in order to fully optimize
the reduction, we want to classify and group diagrams
according to the sets of their propagators.
Within our framework, QGRAF generates three sets of
output files: an expression for each diagram to be processed
with FORM [48], Python code for drawing all diagrams,
and Python code for computing the properties of each
diagram. The information about the model for QGRAF is either
read from the built-in Standard Model file or is generated
from a user defined LanHEP [55] or Universal FeynRules
Output (UFO) [54] file.
The Python program automatically performs several
operations:
diagrams whose color factor turns out to be zero are
dropped automatically;
the fermion flow is determined and used to compute an
overall sign for each diagram, which is relevant in the
presence of Majorana fermions;
the number of propagators containing the loop
momentum, i.e. the loop size of the diagram, the tensor rank and
the kinematic invariants of the associated loop integral are
computed;
diagrams with an associated vanishing loop integral (see
Fig. 1(a)) are detected and flagged for the diagram
selection;
all propagators and vertices are classified for the diagram
selection; diagrams containing massive quark self-energy
insertions or closed massless quark loops are specially
flagged.
Any one-loop diagram can be written in the form
D =
lN=1[(q + rl )2 ml2 + i]
where the numerator is a polynomial of tensor1 rank r .
N (q) = C0 + C11 q1 + + Cr1...r q1 qr ,
and the N N kinematic matrix is defined as
All masses can be either real or complex. Important
information about the integrals that will appear in the reduction
of each one-loop diagram is contained in the tensor rank r
of the loop integral and its kinematic matrix Sij .
We define a preorder relation on one-loop diagrams, such
that D1 D2 if their associated matrices S(D1) and S(D2)
are related by a finite (not necessarily unique) chain of
transformations
S(D2) T1 S Tm S(D1),
T2
where each transformation is one of the following:
the identity,
the simultaneous permutation of rows and columns,
the simultaneous deletion of the row and column with
the same index, which corresponds to pinching the
corresponding propagator in the diagram.
The relation can be read as appears in the
reduction of. Our algorithm groups the one-loop diagrams
D1, . . . , DD of a process into subsets V1, . . . , VG such that
V1, . . . , VG form a partition of {D1, . . . DD} and
each cell Vi contains a maximum element max Vi Vi ,
such that D max Vi , D Vi .
The partitioning procedure provides an important gain in
efficiency, because while carrying out the tensor reduction for
the diagram max Vi , all other diagrams in the same cell Vi
are reduced with virtually no additional computational cost.
The gain in efficiency can be observed when reducing the
diagram using the OPP method [45] and its implementations
in CutTools [57] and SAMURAI [40], as well as in
classical tensor reduction methods as implemented e.g. in
Golem95C [51, 52], PJFRY [58] and LoopTools [3, 59].
In order to draw the diagrams, we first compute an
ordering of the external legs which allows for a planar embedding
of the graph. Such ordering can always be found for a tree
or a one-loop graph since non-planar graphs only start to
appear in diagrams with two or more loops. After the legs have
been assigned to the vertices of a regular polygon, we use
our own implementation of the algorithms described in [60]
for fixing the coordinates of the remaining vertices; the
algorithm has been extended to determine an appealing
layout also for graphs containing tadpoles. Starting from these
1Index contractions in (2) are understood in n-dimensional space.
coordinates and using the package Axodraw [61], GOSAM
generates a LATEX file that contains graphical representations
of all diagrams.
2.3 Algebraic processing
2.3.1 Color algebra
In the models used by GOSAM, we allow one unbroken
gauge group SU(NC ) to be treated implicitly; any additional
gauge group, broken or unbroken, needs to be expanded
explicitly. Any particle of the model may be charged under the
SU(NC ) group in the trivial, (anti-)fundamental or adjoint
representation. Other representations are currently not
implemented.
For a given process we project each Feynman
diagram onto a color basis consisting of strings of generators
TiAi11 TiA1i22 TiApp1j and Kronecker deltas ij but no
contractions of adjoint indices and no structure constants f ABC .
Considering, for example, the process
u(1) + u(2) Z(3) + g(4) + g(5)
GOSAM finds the color basis
|c1 = qi(11)qj(22)g(A44)g(A55)(T A4 T A5 )j2i1 ,
|c2 = qi(11)qj(22)g(A44)g(A55)(T A5 T A4 )j2i1 ,
|c3 = qi(11)qj(22)g(A44)g(A55)j2i1 tr{T A5 T A4 },
where qi() and g( ) are the color parts of the quark and gluon
A
wave functions respectively. The dimension of this color
basis for Ng external gluons and Nqq quark-antiquark pairs is
given by [62]:
i=0
d(Ng, Nqq ) =
(1)i
(Ng + Nqq i)!.
It should be noted that the color basis constructed in this
way is not a basis in the mathematical sense, as one can find
linear relations between the vectors |ci once the number of
external partons is large enough.
Any Feynman diagram can be reduced to the form
D =
i=1
for the process specific color basis |c1 , . . . , |ck by applying
the following set of relations:
1
TiAjTkAl = TR il kj NC ij kl ,
The same set of simplifications is used to compute the
matrices ci |cj and ci |TI TJ |cj . The former is needed for
squaring the matrix element, whereas the latter is used to
provide color correlated Born matrix elements which we use
for checking the IR poles of the virtual amplitude and also
to provide the relevant information for parton showers like
POWHEG [33, 34, 63]. For the above example, GOSAM
obtains2
(NC2 1)
1
1
(NC2 1)
Similarly, the program computes the matrices ci |TI TJ |cj
for all pairs of partons I and J .
If M(0) denotes the tree-level matrix element of the
process and we have
j=1
then the square of the tree level amplitude can be written as
the n-dimensional algebra (n = 4 2) into strictly
fourdimensional objects and symbols representing the
higherdimensional remainder.
In GOSAM we have implemented the t Hooft-Veltman
scheme (HV) and dimensional reduction (DRED). In both
schemes all external vectors (momenta and polarisation
vectors) are kept in four dimensions. Internal vectors,
however, are kept in the n-dimensional vector space. We adopt
the conventions used in [49], where k denotes the
fourdimensional projection of an in general n-dimensional
vector k. The (n 4)-dimensional orthogonal projection is
denoted as k. For the integration momentum q we introduce in
addition the symbol 2 = q 2, such that
We also introduce suitable projectors by splitting the metric
tensor
In the following, we describe the t Hooft algebra in
detail. For DRED, the only differences are that the numerator
algebra is performed in four dimensions for both external
and internal vectors (i.e. q q ) and that in the very end all
appearances of q2 are replaced by q 2 2.
Wave functions and propagators GOSAM contains a
library of representations of wave functions and propagators
up to spin two.3 The exact form of the interaction vertices is
taken from the model files.
The representation of all wave functions with non-trivial
spin is based on massless spinors. Each massive external
vector pi is replaced by its light-cone projection li with
respect to a lightlike reference vector k,
For spin 1/2 particles we use the assignment of wave
functions as shown in Table 1; here, we quote the
definition of the massive spinors from [49] assuming the splitting
of (16):
= |l]
p = [l|
3Processes with particles of spin 3/2 and spin 2 have not been tested
extensively. Furthermore, these processes can lead to integrals where
the rank is higher than the loop size, which at the moment are neither
implemented in SAMURAI nor in Golem95C.
lN=1[(q + rl )2 ml2 + i]
In order to reduce the complexity at the level of the
reduction, we perform the contraction with the tree-level already
at the integrand level,
ci |cj Ci(0) Cj(1)(q),
Ci(0) Cj(0) ci |cj .
For the interference term between leading and
next-toleading order we use a slightly different philosophy. First
of all we note that it is sufficient to focus on a single group
V as defined in Sect. 2.2,
i,j=1
i,j=1
where Cj(1) is formed by the sum over the corresponding
coefficients of all diagrams D V .
2.3.2 Lorentz algebra
In this section we discuss the algorithms used by GOSAM
to transform the coefficients Ci(0) and Ci(1)(q), as defined in
the previous section, such that the result is suitable for
efficient numerical evaluation. One of the major goals is to split
2In the actual code the results are given in terms of TR and NC only.
(b) Wave functions for massless fermions
(c) Wave functions for massive fermions
u(k, +1) = v(k, 1) = |k
u(k, 1) = v(k, +1) = |k]
u(k, +1) = v(k, 1) = [k|
u(k, 1) = v(k, +1) = k|
(a) Assignment of initial and final states for
quarks and leptons
l, q
In order to preserve the condition that for any loop integral
the tensor rank does not exceed the number of loop
propagators we fix all gauge boson propagators to be in Feynman
gauge. Their wave functions are constructed as [64]
(p, +1) = q||p ] ,
2 qp
where p = p in the massless case and p = l according to
(16) in the massive case. In the latter case the third
polarisation is defined as
The wave functions and propagators for spin 3/2 and spin 2
particles correspond to those in [65].
Simplifications Once all wave functions and propagators
have been substituted by the above definitions and all
vertices have been replaced by their corresponding
expressions from the model file then all vector-like quantities
and all metric tensors are split into their four-dimensional
and their orthogonal part. As we use the t Hooft
algebra, 5 is defined as a purely four-dimensional object, 5 =
i . By applying the usual anti-commutation
relations for Dirac matrices we can always separate the
fourdimensional and (n 4)-dimensional parts of Dirac traces,
as we can use the fact that [49, 62]
= tr(1 l ) tr(l+1 l+p ).
The same logic applies to open spinor lines such as [49]
= k1|1 l |k2 tr(l+1 l+p ).
While the (n 4)-dimensional traces are reduced
completely to products of (n 4)-dimensional metric tensors
g , the four-dimensional part is treated such that the
num
ber of terms in the resulting expression is kept as small
as possible. Any spinor line or trace is broken up at any
position where a light-like vector appears. Furthermore,
Chisholm identities are used to resolve Lorentz
contractions between both Dirac traces and open spinor lines. If
any traces remain we use the built-in trace algorithm of
FORM [48].
In the final result we can always avoid the explicit
appearance of Levi-Civit tensors, noticing that any such tensor is
contracted with at least one light-like vector4 k, and we can
replace
Hence, the kinematic part of the numerator, at the end of our
simplification algorithm, is expressed entirely in terms of:
spinor products of the form ki kj , [ki kj ] or [ki | |kj
q,
dot products ki kj or ki q ,
constants of the Lagrangian such as masses, widths and
coupling constants,
the symbols 2 = q 2 q2 and = (n 4)/2.
Treatment of R2 rational terms In our representation for
the numerator of one-loop diagrams, terms containing the
symbols 2 or can lead to a so-called R2 term [66], which
contributes to the rational part of the amplitude. In general,
there are two ways of splitting the numerator function:
N (q , 2, ) = N0(q , 2) + N1(q , 2)
N (q , 2, ) = N (q ) + N (q , 2, ).
It should be noted that in (23a) the terms N1 and N2 do
not arise in DRED, where only terms containing 2
contribute to R2. Instead of relying on the construction of R2
4Any external massive vector at this point has been replaced by a pair
of light-like ones. Contractions between two Levi-Civit symbols can
be resolved to products of metric tensors.
from specialized Feynman rules [6770], we generate the
R2 part along with all other contributions without the need
to separate the different parts. For efficiency reasons,
however, we provide an implicit and an explicit construction of
the R2 terms.
The implicit construction uses the splitting of (23a) and
treats all three numerator functions Ni on equal grounds.
Each of the three terms is reduced separately in a
numerical reduction and the Laurent series of the three results are
added up taking into account the powers of .
The explicit construction of R2 is based on the
assumption that each term in N in (23b) contains at least one power
of 2 or . The expressions for those integrals are relatively
simple and known explicitly. Hence, the part of the
amplitude which originates from N is computed analytically
whereas the purely four-dimensional part N is passed to the
numerical reduction.
2.4 Code generation
2.4.1 Abbreviation system
To prepare the numerator functions of the one-loop diagrams
for their numerical evaluation, we separate the symbol 2
and dot products involving the momentum q from all other
factors. All subexpressions which do not depend on either
q or 2 are substituted by abbreviation symbols, which are
evaluated only once per phase space point. Each of the two
parts is then processed using haggies [50], which
generates optimized Fortran code for their numerical
evaluation. For each diagram we generate an interface to
SAMURAI [40], Golem95C [52] and/or PJFRY [58]. The two
latter codes are interfaced using tensorial reconstruction at the
integrand level [53].
2.4.2 Reduction strategies
In the implementation of GOSAM, great emphasis has been
put on maintaining flexibility with respect to the reduction
algorithm that the user decides to use. On the one hand, this
is important because the best choice of the reduction method
in terms of speed and numerical stability can strongly
depend on the specific process. On the other hand, we tried
to keep the code flexible to allow further extensions to new
reduction libraries, such that GOSAM can be used as a
laboratory for interfacing future methods with a realistic
environment.
Our standard choice for the reduction is SAMURAI,
which provides a very fast and stable reduction in a large
part of the phase space. Furthermore, SAMURAI reports to
the client code if the quality of the reconstruction of the
numerator suffices the numerical requirements (for details
we refer to [40]). In GOSAM we use this information to
trigger an alternative reduction with either Golem95C [52]
or PJFRY [58] whenever these reconstruction tests fail, as
shown in Fig. 2. The reduction algorithms implemented in
these libraries extend to phase space regions of small Gram
determinants and therefore cover most cases in which
onshell methods cannot operate sufficiently well. This
combination of on-shell techniques and traditional tensor
reduction is achieved using tensorial reconstruction at the
integrand level [53], which also provides the possibility of
running on-shell methods with a reconstructed numerator. In
addition to solving the problem of numerical instabilities,
in some cases this option can reduce the computational cost
of the reduction. Since the reconstructed numerator is
typically of a form where kinematics and loop momentum
dependence are already separated, the use of a reconstructed
numerator tends to be faster than the original procedure,
in particular in cases with a large number of legs and low
rank.
The flowchart in Fig. 2 summarizes all possible
reduction strategies which are currently implemented. The
strategy in use is selected by assigning the variable
reduction_interoperation in the generated Fortran
code. The availability of the branches is determined
during code generation by activating (at least one of) the
extensions (samurai, golem95, pjfry) in the input card.
Fig. 2 Reduction strategies currently implemented in GOSAM: the
reduction algorithm is chosen by setting the variable
reduction_interoperation in the generated Fortran code and can
be modified at run time. 0: SAMURAI only; 1: Golem95C only; 2:
SAMURAI with rescue option (Golem95C); 3: SAMURAI with
numerator from tensorial reconstruction; 4: same as 3 but with rescue option
(Golem95C). 11, 12 and 14 are the same as 1, 2, 3 (respectively) with
the difference that PJFRY is used instead of Golem95C.
Switching between active branches is possible at run time.
In detail, the possible choices for the variable
reduction_interoperation are the following:
0 the numerators of the one-loop diagrams are reduced by
SAMURAI, no rescue system is used in case the
reconstruction test fails;
1 the tensor coefficients of the numerators are
reconstructed using the tensorial reconstruction at the
integrand level, the numerator is expressed in terms of
tensor integral form factors which are evaluated using
Golem95C;
2 the numerators are reduced by SAMURAI; whenever the
reconstruction test fails, numerators are reduced using
the option 1 as a backup method;
3 tensorial reconstruction is used to compute the tensor
coefficients; SAMURAI is employed for the reduction of
the reconstructed numerator, no rescue system is used;
4 as in option 3, SAMURAI is used to reduce the
reconstructed numerator, Golem95C is used as backup
option;
11 same as 1 but PJFRY is used instead of Golem95C;
12 same as 2 but PJFRY is used instead of Golem95C;
14 same as 4 but PJFRY is used instead of Golem95C.
It is difficult to make a statement about the optimal
reduction method because this depends on the process under
consideration. For multi-leg processes, e.g. bbbb
production, we found that SAMURAI is clearly superior to tensor
reduction in what concerns timings and size of the code.
Concerning points which need a special treatment, we did
not make extensive studies using traditional tensor reduction
only, but one can certainly say that the combination of
SAMURAI and tensorial reconstruction seems to be optimal in
what concerns the avoidance of numerical instabilities due
to inverse Gram determinants.
2.5 Conventions of the amplitudes
In this section we briefly discuss the conventions chosen for
the results returned by GOSAM. Depending on the actual
setup for a given process, in particular if an imported model
file is used, conventions may be slightly different. Here we
restrict the discussion to the case where the user wants to
compute QCD corrections to a process and in the setup files
he has put gs = 1. In this case, the tree-level matrix element
squared can be written as
2
|M|tree = A0A0 = (gs )2b a0.
The fully renormalised matrix element at one-loop level, i.e.
the interference term between tree-level and one-loop, can
be written as
= A1A0 + A0A1 = 2 (A0A1)
2 2 2 2 2
= |M|bare + |M|ct,mQ + |M|ct,s + |M|wf,g + |M|wf,Q
A call to the subroutine samplitude returns an array
consisting of the four numbers (a0, c0, c1, c2) in this order.
The average over initial state colours and helicities is
included in the default setup. In cases where the process is
loop induced, i.e. the tree level amplitude is absent, the
pro
gram returns the values for A1A1 where a factor
has been pulled out.
After all UV-renormalisation contributions have been
taken into account correctly, only IR-singularities remain,
which can be computed using the routine
ir_subtractions. This routine returns a vector of length two,
containing the coefficients of the single and the double pole, which
should be equal to (c1, c2) and therefore can be used as a
check of the result.
Ultraviolet renormalisation in QCD For
UV-renormalisation we use the MS scheme for the gluon and all massless
quarks, whereas a subtraction at zero momentum is chosen
for massive quarks [71]. Currently, counterterms are only
provided for QCD corrections. In the case of electroweak
corrections only unrenormalised results can be produced
automatically.
For computations involving loop propagators for
massive fermions, we introduced the automatic generation of
a mass counter term needed for the on-shell
renormalisation of the massive particle. Here, we exploit the fact that
such a counter term is strictly related to the massive fermion
self energy bubble diagrams (see Fig. 3). As described in
Sect. 2.2, the program GOSAM analyzes all generated
diagrams. In that step also self-energy insertions of massive
Fig. 3 Feynman diagram of a
massive quark self energy in
QCD. For this type of diagram
GOSAM automatically
generates UV-counterterms
quarks are detected, where we make the replacement
(q/ + /r + m) g
[(q + r)2 m2]q2
(q/ + /r + m) g m 6q r + 3(r2 m2)
[(q + r)2 m2]q2 + 4 m2
[(q + r)2 m2]q2
The symbol 1HV is one in the t Hooft Veltman scheme and
zero in DRED.
Performing the integral, contracting the expression with
the QCD vertices at both sides and multiplying the
missing factor of (2 )1 we retrieve the expression for the mass
counter-term,
Furthermore, the renormalisation of s leads to a term of
the form
2 s (4 ) 2 0
|M|ct,s = b 2 (1 ) |M|tree
q=Nf +1
with 0 = (11CA 4TRNf )/6, Nf being the number of
light quark flavours, Nf,h the number of heavy flavours, and
b is the power of the coupling in the Born amplitude as
defined in (24). The last term of (28) provides the finite
renormalisation needed to compensate the scheme dependence
of s ,
A further contribution consists of the wave-function
renormalisation of massive external quark lines. If we
denote the set of external massive quark lines by Qh =
{Q1(m1), . . . , Qp(mp)} we obtain
2 s (4 ) CF
|M|wf,Q = 2 (1 ) 2
Q(m)Qh
loops. If Ng is the number of external gluon lines, this
contribution can be written as
2 s (4 ) 2TR
|M|wf,g = 2 (1 ) Ng 3
q=Nf +1
At the level of the generated Fortran code the presence
of these contributions can be controlled by a set of variables
defined in the module config.f90. The variable
renormalisation can be set to 0, 1, or 2. If
renormalisation=0, none of the counterterms are present. If
renormalisation=2 only |M|c2t,mQ is included, which is the
counterterm stemming from all terms of the type of (27)
contributing to the amplitude.
In the case where renormalisation=1 a more
finegrained control over the counterterms is possible.
renorm_logs: if set to false, in all counterterms the
generation of logarithms is disabled, i.e. factors of the form
() in (27) to (31) are replaced by one.
2
renorm_beta: if set to false, the counterterm |M|ct,s is
set to zero.
2
renorm_mqwf: if set to false, the counterterm |M|wf,Q is
set to zero.
2
renorm_mqse: if set to false, the counterterm |M|ct,mQ
is set to zero.
renorm_decoupling if set to false, the counterterm
2
|M|wf,g is set to zero.
The default settings for renormalisation=1 are true
for all the renorm options listed above.
Finite renormalisation of 5 in QCD In the t Hooft
Veltman scheme, a finite renormalisation term for 5 is required
beyond tree level. The relevant terms are generated only if
fr5 is added in the input card to the list of extensions
before code generation. Currently, the automatic generation of
this finite contribution is not performed if model files
different from the built-in model files are used. In agreement
with [72] and [73] we replace the axial component at each
vertex,
5 21 Zaxial 5 5 ,
Zaxial = 1 2 2s CF 1HV.
Finally, also the wave function of the gluon receives a
contribution from the presence of heavy quarks in closed fermion
Once it is generated, this contribution can be switched on
and off at run-time through the variable renorm_gamma5,
which is defined in the module config.f90.
Conversion between the schemes In GOSAM we have
implemented two different schemes, the t Hooft Veltman
scheme and dimensional reduction. By default, the former
is used, while the latter can be activated by adding the
extension dred. If a QCD computation has been done in
dimensional reduction the result can be converted back to the
t Hooft Veltman scheme by adding a contribution for each
external massless parton,
I =1
with qDR = qDR = CF /2 and gDR = CA/6. This
conversion can be switched on by setting convert_to_cdr to
true in the module config.f90. At one-loop level, the
t Hooft Veltman scheme and conventional dimensional
regularisation (CDR) are equivalent in the sense that It HV = 0
for all partons.
3 Requirements and installation
3.1 Requirements
The program GOSAM is designed to run in any modern
Linux/Unix environment; we expect that Python ( 2.6),
Java (1.5) and Make are installed on the system. Further
more, a Fortran 95 compiler is required in order to
compile the generated code. Some Fortran 2003 features are
used if one wants to make use of the Les Houches
interface [56]. We have tried all examples using gfortran
versions 4.1 and 4.5.
On top of a standard Linux environment, the programs
FORM [48], version 3.3 (newer than Aug. 11, 2010) and
QGRAF [2] need to be installed on the system. Whereas
spinney [49] and haggies [50] are part of GOSAM and
are not required to be installed separately, at least one of the
libraries SAMURAI [40] and Golem95C [52] needs to be
present at compile time of the generated code. Optionally,
PJFRY [58] can be used on top of Golem95C.
3.2 Download and installation
QGRAF The program can be downloaded as Fortran
source code from
http://cfif.ist.utl.pt/~paulo/qgraf.html.
After unpacking the tar-ball, a single Fortran 77 file
needs to be compiled.
FORM The program is available at http://www.nikhef.nl/~form/ both as a compiled binary for many platforms and as a tar
ball. The build process, if built from the source files, is
controlled by Autotools.
SAMURAI and Golem95C These libraries are available as
tar-balls and from subversion repositories at
http://projects.hepforge.org/samurai/
http://projects.hepforge.org/golem/95/
respectively. For the users convenience we have
prepared a package containing SAMURAI and Golem95C
together with the integral libraries OneLOop [74],
QCDLoop [75] and FF [59]. The package
gosam-contrib1.0.tar.gz containing all these libraries is available for
download from:
http://projects.hepforge.org/gosam/
GOSAM The user can download the code either as a
tarball or from the subversion repository at
http://projects.hepforge.org/gosam/.
The build process and installation of GOSAM is controlled
by Python Distutils, while the build process for the
libraries SAMURAI and Golem95C is controlled by
Autotools.
Therefore the installation proceeds in two steps:
1. For all components which use Autotools, the
following sequence of commands installs them under the user
defined directory MYPATH.
./configure --prefix=MYPATH
make FC=gfortran F77=gfortran
make install # or sudo make install
If the configure script is not present, the user needs
to run sh ./autogen.sh first.
2. For GOSAM which is built using Distutils, the user
needs to run
python setup.py install \
--prefix MYPATH
If MYPATH is different from the system default (e.g.
/usr/bin), the environment variables PATH,
LD_LIBRARY_PATH and PYTHONPATH might have to be set
accordingly. For more details we direct the user to the
GOSAM reference manual and to the documentation of the
beforementioned programs.
4 Using GOSAM
4.1 Setting up a simple process
GOSAM is a very flexible program and comes with a wide
range of configuration options. Not all of these options are
relevant for simple processes and often the user can leave
most of the settings at their default values. In order to
generate the code for a process, one needs to prepare an input file,
which will be called process card in the following, which
contains
process specific information, such as a list of initial and
final state particles, their helicities (optional) and the order
of the coupling constants;
scheme specific information and approximations, such as
the regularisation and renormalisation schemes, the
underlying model, masses and widths which are set to zero,
the selection of subsets of diagrams; the latter might be
process dependent;
system specific information, such as paths to programs
and libraries or compiler options;
optional information for optimisations which control the
code generation.
In the following we explain how to set up the required files
for the process qq gZ0 g ee+. The example
computes the QCD corrections for the uu initial state, where
me = 0 and Nf = 5 massless quarks are assumed. For our
example, we follow an approach where we keep the
different types of information in separate filesprocess.rc,
scheme.rc and system.rcand use GOSAM to
produce a process card for this process based on these files. This
is not requiredone could also produce and edit the process
card directlyit is however more convenient to store
system specific information into a separate, re-usable file, and
it makes the code generation more transparent.
Process specific information The following listing
contains the information which is specific to the process. The
syntax of process cards requires that no blank character is
left between the equals sign and the property name.
Commentary can be added to any line, marked by the #
character. Line continuation is achieved using a backslash at the
end of a line.5
5The line numbers are just for reference and should not be included in
the actual files.
The first line defines the (relative) path to the directory
where the process files will be generated. GOSAM expects
that this directory has already been created. Lines 2 and 3
define the initial and final state of the process in terms of
field names, which are defined in the model file. Besides the
field names one can also use PDG codes [76, 77] instead.
Hence, the following lines would be equivalent to lines 2
and 3 in Listing 1:
in=2,-2
out=21,11,-11
Line 4 describes the helicity amplitudes which should
be generated. If no helicities are specified, the program
defaults to the generation of all possible helicity
configurations, some of which may turn out to be zero. The
different helicity amplitudes are separated by commas; within one
helicity amplitude there is one character (usually +, -
and 0) per external particle from the left to the right. In the
above example for the reaction
we have the following assignments:
With the above value for helicities we generate all
non-vanishing helicities for the partons but keep the
lepton helicities fixed. In more complicated examples this way
of listing all helicities explicitly can be very tedious.
Therefore, we introduced the option to generate sets of helicities
using square brackets. For example, if the gluon helicity is
replaced by [+-], the bracket is expanded automatically to
take the values +,-.
A further syntactical reduction can be achieved for the
quarks. The current expansion of a square bracket and its
opposite value can be assigned to a pair of variables as in
[xy=+-]. If the bracket expands to + then x is assigned
+ and y is assigned the opposite sign, i.e. -. If the bracket
expands to - the assignments are x=- and y=+. Hence, the
helicity states of a massless quark anti-quark pair are
generated by [qQ=+-]Q, and the selection of helicities in our
example can be abbreviated to
which is equivalent to the version of this line in Listing 1.
Listing2 Filescheme.rc
1 extensions=dred
2 qgraf.options=onshell
3 zero=mU,mD,mC,mS,mB,me,wT
4 one=gs
gosam.py --template gosam.in \
--merge process.rc \
--merge scheme.rc --merge system.rc
6In this example we assume that the user has defined an environment
variablePREFIX.
The generated file can be processed with gosam.py
directly but requires the process directory to be present.
All further steps are controlled by the generated make files;
in order to generate and compile all files relevant for the
matrix element one needs to invoke
make compile
The generated code can be tested with the program
matrix/test.f90. The following sequence of commands
will compile and run the example program.
cd matrix
make test.exe
./test.exe
The last lines of the program output should look as follows7
# LO:
# NLO, finite part
# NLO, single pole
# NLO, double pole
# IR, single pole
# IR, double pole
0.3450350717601E-06
-10.77604823456547
-19.98478948141949
-5.666666665861926
-19.98478948439310
-5.666666666666666
The printed numbers are, in this order, a0, c0/a0, c1/a0,
c2/a0 and the pole parts calculated from the infrared
insertion operator [78, 79].
One can generate a visual representation of all generated
diagrams using the command
make doc
which generates the file doc/process.ps using a
Python implementation of the algorithm described in [60] and
the LATEX package AXODRAW [61].
4.1.1 Further options
GOSAM provides a range of options which influence the
code generation, the compilation and the numerical
evaluation of the amplitude. Giving an exhaustive list of all options
would be far beyond the scope of this article and the
interested user is referred to the reference manual. Nonetheless,
we would like to point out some of GOSAMs capabilities
by presenting the corresponding options.
7The actual numbers depend on the random number generator of the
system because the phase space point is generated randomly; however,
the pole parts should agree between the matrix element and the infrared
insertion operator given that the matrix element is fully renormalised.
Generating the R2 term When setting up a process the user
can specify if and how the R2 term of the amplitude should
be generated by setting the variable r2 in the setup file.
Possible options for r2 are implicit, which is the
default, explicit, off and only. The keyword
implicit means that the R2 term is generated along with
the four-dimensional numerator as a function in terms of q ,
2 and and is reduced at runtime by sampling different
values for 2. This is the slowest but also the most general
option. Using the keyword explicit carries out the
reduction of terms containing 2 or during code generation
(see Appendix B). The keyword off puts the R2 term to
zero which is useful if the user wants to provide his own
calculation for these terms. Conversely, using r2=only
discards everything but the R2 term (reducing it as in the
case explicit) and puts GOSAM in the position of
providing R2 terms for external codes which work entirely in
four dimensions.
Diagram selection GOSAM offers a two-fold way of
selecting and discarding diagrams. One can either influence
the way QGRAF generates diagrams or apply filters to the
diagrams after they have been generated by QGRAF or
combine the two methods. Let us assume that in the above
example we want to remove the third generation of quarks
completely. Hence, all closed quark loops would be massless and
therefore the second generation is just an exact copy of the
first one. We can therefore restrict the generation of closed
quark loops to up and down quarks. GOSAM has a filter
precisely for this purpose, which takes the field names of the
flavours to be generated as arguments.
This filter can be combined with the already existing filter
selecting only diagrams containing a Z-propagator using the
AND function:
filter.nlo=AND( NFGEN(U,D), \
IPROP([Z]) == 1 )
A further feature of the code generated by GOSAM is the
possibility of selecting diagrams at runtime. For example,
we would like to distinguish at runtime three different gauge
invariant sets of diagrams at one-loop level:
1. diagrams with a closed quark loop where the Z is
at
tached to the loop;
2. diagrams with a closed quark loop where the Z is emitted
from the external quark line;
3. diagrams without a closed quark loop.
In order to provide the code for a diagram selection at run- 4.2 Interfacing the code
time one simply replaces the above filter by a list of filters
as follows
The first mandatory arguments of this routine are the
external momenta vecs, where vecs(i,:) contains the
momentum of the i-th particle as a vector [Ei , pix , piy , piz],
and we use in-out kinematics, i.e. p1 + p2 p3 + +
pN . Maximal numerical stability is achieved if the beam
axis is chosen along the z-axis. The second argument,
scale2 = 2R , is the square of the renormalisation scale.
As a third argument the routine expects a vector which
accepts the result in the format [a0, c0, c1, c2] with the
coefficients being defined in (24) and (25). The optional
argument ok may be used in order to report the outcome of
the reconstruction tests in samurai if no rescue method has
been chosen (see Sect. 2.4.2). The last argument allows one
to select a single helicity subamplitude; the index h runs
from zero to the number of helicities minus one. The
labeling of the helicities is documented for each process in the
file doc/process.ps.
exitgolem This routine should be called once after the
last amplitude evaluation in the program. It closes all open
log files and gracefully terminates the reduction and loop
libraries.
subroutine exitgolem(exit_libs)
use config, only: ki
logical, optional,
& intent(in) :: exit_libs
end subroutine
end interface
The optional argument exit_libs should only be set if
multiple calls to this routine (e.g. for different matrix
elements) are necessary and the dependent libraries should be
terminated only once.
A small program which computes the amplitude for a set
of phase space points is automatically generated with the
amplitude code in the file test.f90 in the subdirectory
matrix. The script config.sh in the process directory
returns suitable compilation and linking options for the
generated matrix element code.
4.3 Using the BLHA interface
The so-called Binoth Les Houches Accord (BLHA) [56]
defines an interface for a standardized communication between
one-loop programs (OLP) and Monte Carlo (MC) tools. The
communication between the two sides is split into two main
phases: an initialisation phase and a runtime phase. During
initialisation the two programs establish an agreement by
exchanging a set of files and typically initiate the code
generation. The OLP runtime code is then linked to the MC
program and, during the runtime phase, called through a
welldefined set of routines providing NLO results for the phase
space points generated by the MC. According to this
standard, it is the responsibility of the MC program to provide
results for the Born matrix element, for the real emission and
for a suitable set of infrared subtraction terms. A schematic
overview on this procedure is shown in Fig. 4.
GOSAM can act as an OLP in the framework of the
BLHA. In the simplest case, the MC writes an order filein
this example it is called olp_order.lhand invokes the
script gosam.py as follows:
gosam.py --olp olp_order.lh
Further, GOSAM specific options can be passed either in a
file or directly at the command line. One can, for example,
use autotools for the compilation by modifying the above
line as follows.
gosam.py --olp olp_order.lh \
extensions=autotools
The contract file is given the extension .olc by default
and would be olp_order.olc in this example.
Alternatively, the name can be altered using the -o option.
If successful, the invocation of gosam.py generates
a set of files which can be compiled as before with a
generated make file. The BLHA routines are defined in
the Fortran module olp_module but can also be
accessed from C programs.9 The routines OLP_Start and
OLP_EvalSubProcess are defined exactly as in the
BLHA proposal [56]. For convenience, we extended the
interface by the functions OLP_Finalize(), which
terminates all reduction libraries, and OLP_Option(char*,
int*), which can be used to pass non-standard options at
runtime. For example, a valid call in C to adjust the Higgs
mass would be
int ierr;
OLP_Options("mH=146.78", &ierr);
A value of one in ierr indicates that the setting was
successful. A value of zero indicates an error.
4.4 Using external model files
With a few modifications in the process description files,
GOSAM can immediately make use of model files generated
by either FeynRules [80] in the UFO format [54] or by
LanHEP [55]. In both cases, the following limitations and
differences with respect to the default model files, sm and
smdiag, apply:
9A header file is provided in olp.h.
Fig. 4 Schematic overview
over the interaction between
Monte Carlo tool and one-loop
program in the Binoth Les
Houches Accord
As usual, particles can be specified by their PDG code.
The field names, as used by QGRAF, are parti and
antii for the particles with the PDG code i and i
respectively. For example, the W + and the W boson
would be called part24 and anti24.
All model parameters are prefixed by the letters mdl in
order to avoid name clashes with existing variable names
in the matrix element code.
The variable model.options and the extension fr5
are not guaranteed to work with models other than the
built-in models.
Importing models in the UFO format A model description
in the UFO format consists of a Python package stored in
a directory. In order to import the model into GOSAM one
needs to set the model variable in the input card to specify
the keyword FeynRules in front of the directory name,
where we assume that the model description is in the
directory $HOME/models/MSSM_UFO.
Importing models in the LanHEP format LanHEP model
descriptions consist of a set of plain text files in the same
directory with a common numbering (such as func4.mdl,
lgrng4.mdl, prtcls4.mdl, vars4.mdl). A
LanHEP model can be loaded by specifying the path and
the common number in the model variable. Assuming
the files are situated in the directory $HOME/models/
MSSM_LHEP one would set the variable as follows.
Details about the allowed names for the table columns are
described in the GOSAM reference manual. Precompiled
MSSM_UFO and MSSM_LHEP files can also be found in the
subdirectory examples/model.
5 Sample calculations and benchmarks
The codes produced by GOSAM have been tested on several
processes. In this section we describe some examples of
applications. Additional results, whose corresponding code is
also included in the official distribution of the program, will
be reported in Appendix B.
5.1 pp W + j with SHERPA
In Sect. 4.3 the BLHA interface of GOSAM was presented.
This interface allows one to link the program to a Monte
Carlo event generator, which is, in general, responsible for
supplying the missing ingredients for a complete NLO
calculation of a physical cross section. Among the different
general purpose Monte Carlo event generators, SHERPA
[81] is one of those which offers these tools: computing the
LO cross section, the real corrections with both the
subtraction terms and the corresponding integrated counterparts
[8284]. Furthermore, SHERPA offers the possibility to
match a NLO calculation with a parton shower [85, 86].
Using the BLHA interface, we linked GOSAM with SHERPA
to compute the physical cross section for W + 1 jet at NLO.
The first steps to perform this linking is to write a
SHERPA input card for the desired process. Instructions and many
examples on how to write this can be found in the on-line
manual [87]. Running the code for the first time will
produce an order file OLE_order.lh which contains all the
necessary information for GOSAM, to produce the desired code
for the loop part of the process. This includes a list of all
partonic subprocesses needed. In parallel to the production of
the needed SHERPA libraries with the provided script, one
can at this point run the gosam.py command with the flag
-olp and the correct path to the order file as explained in
Sect. 4.3. Further options may be specified. Among them it
is useful to have a second, GOSAM-specific, input card with
all the important GOSAM options. Since, at the end,
SHERPA needs to be linked to a dynamic library, it is convenient
to run GOSAM with the autotools extension, which allows
the direct creation of both static and dynamic libraries,
together with the test routine test. The gosam.py script
creates all the files needed for interfacing GOSAM with
the Monte Carlo event generator together with the code for
the one-loop computation of all needed subprocesses, and
a makefile to run them. The different parton-level
subprocesses are contained in different subdirectories. At this point
the user simply has to run the makefile to generate and
compile the code. Once the one-loop part of the code is ready, the
produced shared library must be added to the list of needed
libraries in the SHERPA input card as follows.
With this operation the generation of the code is completed.
The evaluation of the process and the physical analysis can
then be performed at the users discretion following the
advice given in the SHERPA on-line documentation [87].
We tested the BLHA interface by computing W + 1 jet
and producing distributions for several typical observables.
result ud Wg
Fig. 5 NLO calculation of W + 1 jet production at LHC using GOSAM interfaced with SHERPA via the BHLA interface. The comparison to
MCFM is also shown
In Figs. 5(a) and 5(b) the inclusive transverse momentum
and rapidity of the jets is shown. These distributions were
compared with similar ones produced using the program
MCFM [27, 28], and perfect agreement was found.
5.2 pp W + j , EW corrections
As a first example of an electroweak calculation, we
computed the virtual one-loop corrections to ud W g. A
complete analytical calculation for this process was presented in
Ref. [88].
91.1876
0.88156596117995394232
MW
For the kinematic point given in Table 2 and the above
parameters we obtain the following result:
2.812364835883295
94.52525523327047
17.84240236996827
0.5555555555555560
(67, 70) of Ref. [88]
4.7438251678146885
0.5555555555555555
The poles have been renormalized using (49)(64) in
Sects. 3.3 and 3.4 of [88]. Our result is agreement with (67),
(70) of Ref. [88] and with Ref. [89] for the infrared
divergences that remain after renormalisation.
Table 2 Kinematic point used in pp W + j , EW
The process in the Standard Model first arises at
the one-loop order, and proceeds through a closed loop of
fermions and W bosons. Of the 16 helicity amplitudes
contributing to it, only three are independent and their analytic
expressions can be found in [90]. The pure QED
contribution, involving a fermion loop, is contained in
samurai1.0 [40] and will not be repeated here. Instead, we show the
results of the W -loop contribution to the independent
helicity amplitudes, as an example of EW corrections that can be
handled with GOSAM.
With the above parameters and the kinematics of Table 3 we
obtain the following results.
12.025419045962
7.3804060437434
982.78049397093
GOSAM (dred)
12.02541904626610
7.380406043429961
982.7804939723322
As an example for the usage of GoSam with a model file
different from the Standard Model we calculated the QCD
corrections to neutralino pair production in the MSSM. The
model file has been imported via the interface UFO
(Universal FeynRules Output) [54] which facilitates the
import of Feynman rules generated by FeynRules [80] to
programs generating one-loop amplitudes. To import such
files within the GoSam setup, all the user has to do is to
give the path to the corresponding model file in the input
card.
For this example, we combined the one-loop amplitude
with the real radiation corrections to obtain results for
differential cross sections. A calculation of neutralino pair
production for the LHC presenting total cross sections at NLO
is given in [91].
For the infrared subtraction terms the program
MadDipole [92, 93] is used, the real emission part is
calculated using MadGraph/MadEvent [94]. The virtual matrix
element is renormalized in the MS scheme, while massive
particles are treated in the on-shell scheme. The
renormalisation terms specific to the massive MSSM particles have
been added manually.
In Fig. 6 we show the differential cross section for the
m1010 invariant mass, where we employed a jet veto to
sup0 0
press large contributions from the channel qg 1 1 q
which opens up at order 2s , but for large pjet belongs to
T
the distinct process of neutralino pair plus one hard jet
production at leading order. We used Nf = 5 massless quark
flavours and the MSTW08 [95] parton distribution
functions. For the SUSY parameters we use the modified
benchmarks point SPS1amod suggested in [96], and we use s =
7 TeV.
For reference, we also give the result for the
unrenormalised amplitude at one specific phase space point for
0 10 in the DRED scheme, using the following
pauu 1
rameters and momenta:
All widths have been set to zero; for further settings we
refer to the model parameter files contained in the
subdirectory examples/model/MSSM_UFO. We have checked
that the pole terms of the renormalised amplitude cancel
with the infrared poles from MadDipole. For the phase
0 10 with a jet veto on jets with pjet > 20 GeV and
process pp 1 T
< 4.5. The band gives the dependence of the result on = F = R
between 0/2 and 20. We choose 0 = MZ. The black line gives the
bin error for the value at the central scale
space point given in Table 4 we obtain the following
numbers.
0.8680577964243597103
31.9136615197871
13.4374663711899
2.6666666666667
As an example of a QED calculation, we compared the
virtual QED corrections for the process e+e e+e with
the results provided in [97]. The results compared in the
table are the bare unrenormalised amplitudes in the t Hooft
Veltman scheme. No counterterms or subtraction terms have
been added to the result.
7.2973525376 103
0.51099891 103
Using the parameters given above and the kinematics of
Table 5 we obtain the following results.
0.7586101468103622
0.5005827938274887
0.0474506407008029
0
0.7586101468103619
0.5005828268263969
0.0474506427003504
0
0.4999997388800458
0.4999997388800458
0.05006809884093004
0.1832142729949070
0.1331461741539769
s
result uu ttH
t
t
Table 6 Kinematic point used in pp ttH
Table 7 Kinematic point used for gg ttZ
136.35582793693018
181.47665951104506
182.16751255202476
6270.1855170414337
6925.5258180925930
804.28866486597315
5.6 pp t t H
This process has been compared with the results given
in [39]. The partonic subprocesses uu t tH and gg
t tH where computed both in the t Hooft Veltman scheme
and in dimensional reduction and the fully renormalised
results were successfully compared as an internal
consistency check. Apart from wave function renormalisation and
mass counterterms, Yukawa coupling renormalisation is also
needed here. Yukawa coupling counterterms are in this case
equal to the wave function counterterms. The Yukawa top
mass is set equal to its pole mass.
The kinematics used to obtain the results below is given
in Table 6. The results are given in the t Hooft Veltman scheme, and are fully renormalised.
15.133871809486299
20.889486679044587
36.023358488530903
27.986733991031045
50.105625289561424
22.118891298530357
4977.7694025303863
5306.3374282745517
328.56802574416463
806.93726196887712
1281.8763412410237
474.93907927214622
250.0
26.088703626953386
14.002628607367491
40.091332234320859
7000.0
3725.2619580634337
4258.3185872039012
533.05662914046729
6.127399805961155
9.006680638719660
2.986347664537282
6.000000000000004
6.127400074872043
9.006680836410272
2.9863477301662056
6.000000131659877
On an Intel Core i7 950 at 3 GHz the evaluation of
a single phase space point took 44 ms in the uu channel
and 223 ms in the gg channel. The code was compiled with
gfortran without optimisations.
This amplitude, fully renormalised, has been compared with
the results given in [37].
result gg ttH
result gg ttZ
2.200490364806190
15.29615178164782
1.640361500121837
2.666666666666666
2.2004904613782828
15.29615211731521
1.640361536072381
2.666666725182165
0.1531395190212139
204.9208290898557
50.62939646427283
5.999999999999997
0.1531395190212831
204.920829867328
50.6293965717156
6.00000000000003
The evaluation of a single phase space point took 1433 ms
on a 2 GHz processor. The code was compiled with
gfortran -O2.
5.8 pp bbbb + X
A detailed discussion of this process can be found in
[98, 99]. In this section we focus on the parts that are
relevant in the context of the virtual corrections. In particular
we compared our result to the one given in [38], which is
the fully renormalised amplitude including the mass
counterterms for the top-quark contribution.
The results below are obtained for the phase space point of
Table 8 using the above parameters.
1.022839601391936
36.97653243659754
34.01491655155776
11.33333333333512
Table 8 Kinematic point used in pp bbbb
Table 9 Kinematic point used in pp ttbb
250.0
250.0
147.5321146846735
108.7035966213640
194.0630765341365
49.70121215982584
250.0
250.0
190.1845561691092
182.9642163285034
100.9874727883170
25.86375471407044
1.022839601391910
36.97653243473214
34.01491655142099
11.33333333333343
5.753293428094349
22.19223384585620
20.89828996870689
8.000000000000199
5.753293428094391
22.19223384564902
20.89828996857439
8.000000000000037
On an Intel Xeon E7340 the running time for the
calculation of a single phase space point was 19.6 s for the gluon
initiated channel and 440 ms for the quark channel.
5.9 pp t t bb + X
This process has been compared with the results given
in [38]. We have set up the process both in the t Hooft
Veltman scheme and in dimensional reduction and successfully
compared the fully renormalised results as an internal
consistency check. The results below are given in the t Hooft
2
Veltman scheme, and only the counterterms for |M|ct,mt
are included.
Using the above parameters and the phase space point of
Table 9 we obtain the following results.
result uu ttbb
result gg ttbb
2.201164677187755
8.880263116574282
4.730495922109534
5.333333333333468
8.279470201927135
21.83922035777929
12.59181277770347
8.666666666666764
On an Intel Core i7 950 at 3 GHz the evaluation of a
single phase space point took 393 ms in the uu channel
and 12.3 s in the gg channel. The code was compiled with
gfortran without optimisations.
5.10 pp W +W bb
have been calculated both in [38] and [39]. Accordingly,
the results below are given in the t Hooft Veltman scheme,
5
1
0
1
2.44140351
0
s
W
With the above parameters and the kinematics defined in
Table 10 we obtain the following results.
Table 10 Kinematic point used in pp W +W bb
2.201164677187727
8.880263117410131
4.730495921691266
5.333333333333190
8.279470201927128
21.83922035648926
12.59181277853837
8.666666666666549
result gg W +W bb
result uu W +W bb
1.549796787502985
17.80558461276584
19.61125131175888
8.666666666666668
2.338048681706755
5.936151367348438
10.44868110371249
5.333333333333312
1.549795815702494
17.80558440908488
19.611251301307803
8.66666666666661
2.338048676370483
5.936151368788066
10.44868110378090
5.333333333333336
71.68369328357026
43.36246487297426
16.53704884550758
0.4920930146078141
5.11 ud W +ggg
The amplitude ud W +ggg is an important channel in the
calculation of the process pp W + + 3 jets. The QCD
corrections to this process have been presented in Refs. [69].
The subprocess with one quark pair and three gluons
consists of more than 1500 Feynman diagrams. We have
computed the amplitude including the leptonic decay of the W
boson and compared our result to [38].
s
Furthermore, the values for the dependent parameters are
cos2 W = MW2 /MZ2 and = 2GF MW2 sin2 W / . For the
phase space point of Table 11 we obtain the numbers below.
result ud W +ggg
8.552735739069321
36.45372625230239
34.70010131004584
11.66666666666747
36.4536949986367
34.70007155977844
11.666656664302845
250.0
158.5329205583824
1.357120667754454
Table 11 Kinematic point used in ud W +ggg
On an Intel Core 2 i5 Laptop at 2.0 GHz the evaluation of a
single phase space point took about 2.5 s for ud e+eggg
and about 7.5 s for on-shell Ws without decay. The code
was compiled with gfortran -02.
The process ud W +bb, with an on-shell W -boson, has
been studied in [100], while the effects of the W -decay
have been recently accounted for in [101], and implemented
within MCFM. We consider the latter process, and
compare the renormalised amplitude evaluated by MCFM. The
b-quark is treated as massive in all diagrams except in the
vacuum-polarisation like contributions.
Using the above parameters and the kinematics given in
Table 12 we obtain the following results.
Table 12 Kinematic point used in ud W +bb
250.0
250.0
162.5391101447744
104.0753327455388
185.8004692730082
47.58508783667868
76.084349979114506
1998.0331337409114
953.55303294091811
190.20402007017753
417.39085287123652
360.80087787946474
76.084349979114506
1998.0331337409114
17.060211586132972
39.063065018596490
40.625785184424117
96.749061789153586
1.88443466774536441
41.217129894410029
I R1
I R2
1.884434667673654
41.21712989438873
26.60367070701196
2.666666666666624
26.60367070701218
2.666666666666667
The evaluation of a single phase space point took 9.12 ms
on a 2 GHz processor. The code was compiled with
gfortran -O2.
6 Conclusions
We have presented the program package GOSAM which
produces, in a fully automated way, the code required to
perform the evaluation of one-loop matrix elements for
multiparticle processes. The program is publicly available at
http://projects.hepforge.org/gosam/ and can be used to
calculate one-loop amplitudes within QCD, electroweak
theory, or other models which can be imported via an interface
to LanHEP and UFO, also included in the release. Monte
Carlo programs for the real radiation can be easily linked
through the BLHA interface.
GOSAM is extremely flexible, allowing for both
unitarity-based reduction at integrand level and traditional
tensor reduction, or even for a combination of the two
approaches when required. The amplitudes are generated in
terms of Feynman diagrams within the dimensional
regularisation scheme, and optionally the calculation can be
carried out either in the t Hooft Veltman or in the dimensional
reduction variant. The user can choose among different
libraries for the master integrals, and the setup is such that
other libraries can be linked easily.
The calculation of the rational terms is very modular and
can proceed either along with the same numerical
reduction as the rest of the amplitude, or independently, before
any reduction, by using analytic information on the integrals
which can potentially give rise to a rational part. In the
current version of the code, UV-renormalisation counterterms
are provided for QCD corrections only. Further
improvements concerning the full automatisation of electroweak
corrections are planned.
Different systems to detect and rescue numerical
instabilities are implemented, and the user can switch between them
without having to re-generate the source code. Due to a
careful organisation of the calculation both at the code
generation stage and at the reduction stage, the runtimes for
multiparticle amplitudes are very satisfactory. Moreover, the
GOSAM generator can also produce codes for processes that
include intermediate states with complex masses.
Within the context of the automated matching of Monte
Carlo programs to NLO virtual amplitudes, GOSAM can
be used as a module to produce differential cross sections
for multi-particle processes which can be compared directly
to experiment. Therefore we believe that GOSAM can
contribute to the goal of using NLO tools as a standard
framework for the LHC data analysis at the TeV scale.
Acknowledgements We would like to thank the SHERPA
collaboration for the support, in particular Jennifer Archibald, Frank Krauss
and Marek Schnherr. We also would like to thank Rikkert
Frederix, Adam Kardos, Stefano Pozzorini, Zoltan Trocsanyi, Thomas
Schutzmeier and Christian Sturm for their input to various
comparisons and clarifying discussions, and Edoardo Mirabella for important
feedback on Samurai. G.C. and G.L. were supported by the British
Science and Technology Facilities Council (STFC). The work of G.C.
was supported by DFG Sonderforschungsbereich Transregio 9,
Computergesttzte Theoretische Teilchenphysik. N.G. was supported in
part by the U.S. Department of Energy under contract No.
DE-FG0291ER40677. P.M. and T.R. were supported by the Alexander von
Humboldt Foundation, in the framework of the Sofja Kovaleskaja Award
Project Advanced Mathematical Methods for Particle Physics,
endowed by the German Federal Ministry of Education and Research.
The work of G.O. was supported in part by the National Science
Foundation under Grant PHY-0855489 and PHY-1068550. The research
of F.T. is supported by Marie-Curie-IEF, project: SAMURAI-Apps.
We also acknowledge the support of the Research Executive Agency
(REA) of the European Union under the Grant Agreement number
PITN-GA-2010-264564 (LHCPhenoNet).
Open Access This article is distributed under the terms of the
Creative Commons Attribution License which permits any use,
distribution, and reproduction in any medium, provided the original author(s)
and the source are credited.
Appendix A: Examples included in the release
In the following we give results for the processes listed in
the examples directory. Unless stated otherwise, we
assume that the coupling constants (e and gs in the standard
model) have been set to one in the input card. The
conventions for the returned numbers (a0, c0, c1, c2) are as
stated in Sect. 2.5. Dimensionful parameters are understood
to be in powers of GeV.
As an illustration of the potential of GOSAM, we display
in Table 13 the timings required by a wide list of benchmark
processes. The first value provided in the table is the time
required for the code generation (Generation, given in
seconds): we remind the reader that this operation only needs
to be performed once per process. The second value is the
timing for the full calculation of the amplitude at one
phasespace point (Evaluation, in milliseconds). Results are
obtained with an Intel(R) Core(TM) i7 CPU 950 @ 3.07 GHz.
Table 13 Time required for code generation and calculation of
one phase-space point. The results were obtained with an Intel(R)
Core(TM) i7 CPU 950 @ 3.07 GHz. The time for the evaluation of
a phase space point is taken as the average of the time obtained from
the evaluation of 100 random points generated using RAMBO [102],
where the code was compiled using gfortran without any
optimisation options. The generation of the R2 term was set to explicit
dd t t (DRED)
dg dg (DRED)
e+e t t
e+e t t (LanHEP)
e+e uu
gg gg (DRED)
gg gg (LanHep)
gg t t (DRED)
gg t t (UFO)
ud W +W +cs
A.1 How to run the examples
The example directories only define the system
independent part of the setup. All settings which are defined in the
file system.rc (see Sect. 4) must be put either in a file
called $HOME/.gosam or in the file setup.in in the
GOSAM examples/ directory. A script runtests.sh
is provided to generate, compile and run the test programs.
The names of the directories respectively examples to be run
should be specified at the command line, e.g.
./runtests.sh eeuu bghb
If the script is invoked without arguments it will loop over all
subdirectories. A second script, summarize.sh, can be
used in order to collect the test results and print a summary
to the screen. The command
./summarize.sh
will produce an output like the following one.
+ bghb (succeeded)
+ eeuu (succeeded)
grep: ./ddtt/...: No such file ...
The examples e+e t t have an explicit dependence on
the Golem95C library and will therefore fail if the
extension golem95 is not added.
A.2 e+e uu
The following parameters and momenta have been used to
produce the numerical result:
Table 14 Kinematic point used in e+e tt
result e+e uu
A.3 e+e t t
GOSAM
3.7878306213027528
1.86960440108932 CF
3.0000000000000 CF
2.0000000000000 CF
result e+e tt
A.4 uu dd
This example has been produced twice: once with the
default model file and once with a model file imported from
LanHEP [55]. Thus it also can serve as an example of how to
import model parameters from LanHEP. The result is given
in dimensional reduction, and no renormalisation terms are
included.
The following results are obtained with the above
parameters and the kinematic point of Table 14.
6.3620691850584166
13.182472828297422
12.211527682024421
0
6.3620691850631061
13.182472828302023
12.211527682032367
0
This example has been produced twice: once in the t Hooft
Veltman (HV) scheme and once with dimensional
reduction (DRED). Only the result in the HV scheme will be
74.7646520969852
6067.88254935176
5867.13826404309
275.508937405653
74.7646520969852
6067.88254935176
5862.12966020487
130.988237049907
listed below, for the DRED calculation see the directory
uudd_dred.
Using the above parameters and the phase space point of
Table 15 we obtain the following numbers.
0.28535063700913421
2.7940629929270155
6.4881359148866604
5.3333333333333
0.28535063700913416
2.7940629929268876
6.4881359148866391
5.3333333333333
A.5 gg gg
This example has been produced both with the default model
file and with a model file imported from LanHEP.
Further, it has been calculated in the t Hooft Veltman scheme
and in the dimensional reduction scheme. Only the results
in the t Hooft Veltman scheme are listed below, for
further details please see the subdirectories gggg_dred and
gggg_lhep. The result is for the helicity configuration
g(+)g(+) g()g(), and pure Yang-Mills theory, i.e.
fermion loops are not included.
Table 15 Kinematic point used in uu dd
Table 16 Kinematic point used in gg gg
102.6289752320661
102.6289752320661
102.6289752320661
102.6289752320661
220.9501779577791
220.9501779577791
220.9501779577791
220.9501779577791
Evaluating the amplitude for above parameters and the
phase space point given in Table 16 we obtain the
following results.
result gg gg
GOSAM(HV)
14.120983050796795
124.0247557942351
55.003597347101078
12.00000000000000
A.6 gg gZ
14.120983050796804
124.02475579423495
55.003597347101035
12
As this process has no tree level amplitude, the result is for
the one-loop amplitude squared.
With the above parameters and the kinematics given in
Table 17 we obtain the following result.
0.1075742599502829 0.10757425995048300
This example has been calculated in the t Hooft Veltman
scheme and in the dimensional reduction scheme. Only the
102.6289752320661
102.6289752320661
54.70017191625945
54.70017191625945
220.9501779577791
220.9501779577791
30.55485589367430
30.55485589367430
Table 17 Kinematic point used in gg gZ
Table 19 Kinematic point used in gg tt
results in the t Hooft Veltman scheme are listed below, for
the renormalised amplitude with Nf = 5 and the top mass
renormalised on-shell.
For further details please see the subdirectories ddtt and
ddtt_dred.
a0 0.43024349783870747
c0/a0 22.526901042662193
c1/a0 10.579577611830414
c2/a0 2.6666666666666599
With the above parameters and the kinematics given in
Table 18 we obtain the following results.
Ref. [27, 106] (MCFM)
0.43024349783867882
22.526901042658068
10.579577611830567
2.666666666666721
With the above parameters and the kinematics given in
Table 19 we obtain the following results.
a0 4.5576116986983433
c0/a0 15.352143751168184
c1/a0 27.235240992743407
c2/a0 6.0
A.9 bg H b
Ref. [27, 106] (MCFM)
4.5576116986983424
15.352143750919995
27.235240936279297
6.0
For this process the mass of the b-quark is set to zero.
However, in order to have a coupling between the b-quark and
the Higgs boson, the following Yukawa coupling is
implemented in the model file:
The result is for the renormalised amplitude in the HV
scheme.
Table 18 Kinematic point used in dd tt
74.7646520969852
6067.88254935176
5867.13826404309
275.508937405653
137.84795086008967
3161.1731634194916
3058.6441209877348
240.37699329184659
74.7646520969852
6067.88254935176
5862.12966020487
130.988237049907
137.84795086008967
3161.1731634194916
3049.2945357402382
25.969323180836145
With the above parameters and the kinematics given in
Table 20 we obtain the following results.
Refs. [39, 107]
2.09926265849001642 2.09926265848997195
24.131948141318752 24.131948141995107
11.957924609547224 11.957924605423791
5.6666666666666643 5.6666666666666670
The decay width H of this loop induced process is
known analytically at lowest order. For comparison we used
the equations including the top loop and the bosonic
contribution given in [108, 109]. The decay width can be
expressed as
where i = m2H /(4mi2) for i = W, t .
Table 20 Kinematic point used in bg H b
Refs. [108, 109]
This example has been calculated in the t Hooft Veltman
scheme and in the dimensional reduction scheme. Only the
results in the t Hooft Veltman scheme are listed below, for
the renormalised amplitude. In addition to a calculation with
the default model file, calculations using LanHEP [55] and
UFO [54] are also contained in the examples directory.
s
With the above parameters and the kinematics given in
Table 21 we obtain the following results.
1.4138127601912656
5.4861229357937624
0.18879169932851950
2.666666666666667
We list the renormalised amplitude in the HV scheme.
1.4138127601912673
5.4861229357937660
0.18879169932852413
2.6666666666666665
With the above parameters and the kinematics given in
Table 22 we obtain the following results.
GOSAM(HV)
2.8398509625435832
8.6052919370147745
18.722010655600936
5.6666666666666
2.8398509625435922
8.6052919368774248
18.722010655557121
5.66666666666667
We list the renormalised result in the dimensional reduction
scheme.
With the above parameters and the kinematics given in
Table 23 we obtain the following results.
e
e
1187.7086110647201
2293.0435558834492
509.48956356743611
1282.3236278741238
Table 23 Kinematic point used in g b e e t
86.3112218694181
629.81047833131981
144.72113807954338
774.53161641086319
Ref. [27, 106] (MCFM)
8.52301540708130106
79.879718569273024
26.570185487963364
4.3333333331689596
a0 102 8.52301540675800134
c0/a0 79.879718568538991
c1/a0 26.570185488790770
c2/a0 4.3333333333333401
A.14 u d W +W + c s
Results are given for the unrenormalised amplitude in the
dimensional reduction scheme.
With the above parameters and the kinematics given in
Table 24 we obtain the following results.
result u d W +W + c s
23.3596455167118
13.6255429251954
5.333333333333
5.58083951102529
142.048678636208
258.58120146220904
278.46456389968404
Ref. [20, v3]
500
451.975082051212
1187.7086110647201
2897.148136260289
2189.6399870328105
488.098411670514
968.29887350775562
Table 24 Kinematic point used in u d W +W + c s
Appendix B: Explicit reduction of R2 rational terms
In this Appendix we list all integrals which give rise to R2
terms as we use these expressions in their explicit
construction. We use the definition
2 2 2
Sij = (ri rj ) mi mj .
500
31.1330162081798
27.0607980217775
8.22193223977868
57.8599829481937
65.9095476235891
7.92796656791140
98.5198083786150
171.863734086635
49.8952157196287
l=1
77.0725048002414
5.61185898481311
155.113300136598