GoSam-2.0: a tool for automated one-loop calculations within the Standard Model and beyond

We present the version 2.0 of the program package GoSam for the automated calculation of one-loop amplitudes. GoSam is devised to compute one-loop QCD and/or electroweak corrections to multi-particle processes within and beyond the Standard Model. The new code contains improvements in the generation and in the reduction of the amplitudes, performs better in computing time and numerical accuracy, and has an extended range of applicability. The extended version of the"Binoth-Les-Houches-Accord"interface to Monte Carlo programs is also implemented. We give a detailed description of installation and usage of the code, and illustrate the new features in dedicated examples.


Introduction
After the great achievement of discovering a new boson at the LHC [1,2], the primary goal is now to study its properties in detail, and to detect the slightest hints for possible extensions of the Standard Model. Certainly, precise theory predictions are indispensable to achieve this aim, which calls for calculations at next-to-leading order (NLO) accuracy and beyond.
NLO predictions nowadays should be considered as the standard for experimental data analysis. Ideally, matching NLO results to a parton shower and merging different jet multiplicities should be aimed for. However, this also requires fast and highly automated NLO tools to be available, to be compared to a vast amount of measurements, most of them dealing with multi-particle final states.
The integrand-reduction method [14][15][16] has changed our way of addressing the decomposition of amplitudes in terms of master integrals, whose coefficients can be determined by applying algebraic projections to polynomial functions.
The principle of an integrand-reduction method, which is valid at any order in perturbation theory [17][18][19][20][21], is the underlying multi-particle pole expansion for the integrand of any scattering amplitude, or, equivalently, a representation where the numerator of each Feynman integral is expressed as a combination of products of the corresponding denominators, with polynomial coefficients. These coefficients correspond to the residue of the integrand at the multiple-cut. Each residue is a multivariate polynomial in the irreducible scalar products formed by the loop momenta and either external momenta or polarization vectors constructed out of them.
GoSam is a code which was designed to maximally exploit both the integrand reduction for dimensionally regulated one-loop amplitudes [14,22], as implemented in Samurai [23], as well as improved tensor reduction methods as developed in [24,25]. The algebraic expression of the integrands are automatically generated by means of the Golem technology [24,[26][27][28].
The polynomial structure of the multi-particle residues is a qualitative information that turns into a quantitative algorithm for decomposing arbitrary amplitudes in terms of master integrals at the integrand level. In fact, in the context of an integrand-reduction, any explicit integration procedure is replaced by a simpler operation like polynomial fitting, which in Samurai is implemented via Discrete Fourier Transform [29][30][31].
GoSam produces analytic expressions for the integrands. Because of this feature, it is suitable to be interfaced with a new library, called Ninja [32,33], implementing an ameliorated integrand-reduction method, where the decomposition in terms of master integrals is achieved by Laurent expansion through semi-analytic polynomial divisions [34]. With the new reduction algorithm, GoSam-2.0 can produce results for NLO virtual corrections that are more accurate and less time consuming than the ones provided by version 1.0.
In this paper we present the new version 2.0 of the program GoSam [10], which has been used already to produce a multitude of NLO predictions both within [33,[35][36][37][38][39][40][41][42][43][44][45] and beyond [46,47] the Standard Model. The new version contains important improvements in speed, numerical robustness, range of applicability and user-friendliness. GoSam can be linked to different Monte Carlo programs via the Binoth-Les-Houches-Accord BLHA [48], where the extended version BLHA2 [49] is also implemented. The program can be downloaded from [50] http://gosam.hepforge.org The structure of the paper is the following. In Section 2 we give a brief overview of the program structure. The new features of the program are presented in Section 3. Section 4 describes the installation and usage of GoSam, while in Section 5 we give examples illustrating some of the new features, before we conclude in Section 6. The appendices contain a commented example of an input card for convenience of the user, and some details about higher rank integrals.

Overview of the program
GoSam can be used either as a standalone code producing one-loop (and tree level) amplitudes, or it can be used as a One Loop Provider (OLP) in combination with a Monte Carlo (MC) program, where the interface is automated, based on the standards defined in [48,49]. The main workflow of GoSam is shown in Fig. 1 for the standalone version and in Fig. 2 for GoSam as an OLP within a Monte Carlo setup.
In the standalone version, the user will fill out a process run card which we call process.in, where the process is defined, together with some options. Then the code for the virtual amplitudes is generated by invoking gosam.py process.in . After running the above command with an appropriate run card, all the files which are relevant for code generation will be created. The command make source will invoke QGRAF [51] and FORM [52,53] to generate the diagrams and algebaric expressions for the amplitudes, using also spinney [54] for the spinor algebra within FORM and haggies [55] for code generation. In version 2.0 of GoSam, the production of optimized code however is largely relying on the new features of FORM version ≥ 4. The command make compile will finally compile the produced Fortran90 code.
In the OLP version, the information for the code generation is taken from the order file generated by the Monte Carlo program. Depending on the MC, the whole generation can be invoked automatically and steered by its setup. This is shown schematically in Fig. 2 and explained in more detail in Section 4.3. The amplitudes are evaluated using D-dimensional reduction at integrand level [14,29,56], which is available through two different reduction procedures and libraries: Samurai [23,31] and Ninja [33,34]. Alternatively, tensorial reconstruction [25] is also available, based on the libraries Golem95C [28,57,58] and OneLOop [59].
It should be emphasized that all the reduction libraries used in GoSam-2.0 are included in the program package, and the installation script described in Section 4.1 will take care of compilation and linking, such that the user does not have to worry about installing them separately. Interfacing other tensor integral libraries, such as LoopTools [4,60], PJFRY [61,62] or Collier [63], should be straightforward, due to the modular structure of our setup.
More details about the reduction procedures implemented in GoSam will be given in Section 3.

New features
The version 2.0 of GoSam comes with several new features, which lead to an improvement in speed for both the generation and the evaluation of the amplitudes, more compact code, and more stable numerical evaluation. Further, the range of applicability of the code is extended, in particular to deal with effective theories and physics Beyond the Standard Model. We will describe some of the new features in more detail below. While in version 1.0 of GoSam the Fortran code for the amplitudes was written using haggies [55], we now largely use the features provided by FORM version 4.x [53] to produce optimized code. This leads to more compact code and a speed-up in amplitude evaluation of about a factor of ten. The option to use haggies is still available by setting the extension noformopt.

Grouping/summing of diagrams which share common subdiagrams
Already in the first release of GoSam, the diagrams were analyzed according to their kinematic matrix S ij and grouped together before reduction. This lead to an important gain in efficiency, both with reduction based on integrand reduction methods, as well as with classical tensor reduction techniques. Details about the way diagrams are grouped can be found in [10]. This feature is still present when Samurai or Golem95C are used for computing the amplitudes.
In the new release an option called diagsum combines diagrams which differ only by a subdiagram into one "meta-diagram" to be processed as an entity. This allows one to further reduce the number of calls to the reduction program and therefore to increase the computational speed.

General Information
The maximum effective rank in this group is 5. Fig. 3 Example of diagrams sharing a common tree part, which are summed when the diagsum option is set to diagsum=true. Fig. 4 Example of diagrams sharing a common loop propagator, but with different particle content in the loop, which are summed when the diagsum option is set to diagsum=true.
When the option diagsum is active, diagrams which differ only by a propagator external to the loop, as is the case e.g. for the Z/γ propagator in QCD corrections to the production of Z+jets, are summed together before being processed by FORM. Similarly, diagrams which differ only by an external tree part, but which share exactly the same set of loop propagators, are summed together prior the algebraic manipulation. An example is shown Figure 3. Finally, diagrams which share the same set of propagators, but have different particles circulating in the loop, as shown in Figure 4, are also summed into one "meta-diagram". The default setting for this option is diagsum=true.

Numerical polarisation vectors
The use of numerical polarisation vectors for massless gauge bosons (gluons, photons) is activated by default. This means that the various helicity configurations for the massless bosons will be evaluated numerically, based on a unique code containing generic polarisation vectors, rather than producing separate code for each helicity configuration. To switch off this default setting, for example if the user would like to optimize the choice of reference vectors for each helicity configuration, the option polvec=explicit should be given in the input file process.in.

Improvements in the reduction
The algebraic generation of the integrands in GoSam is tailored to the maximal exploitation of the D-dimensional integrand reduction algorithm.
In the previous version, GoSam-1.0, Samurai has been the default library for the amplitude decomposition in terms of master integrals. Within the original integrand reduction algorithm, implemented in Samurai, the determination of the unknown coefficients multiplying the master integrals requires: i) to sample the numerator on a finite subset of the on-shell solutions; ii) to subtract from the integrand the non-vanishing contributions coming from higher-point residues; and iii) to solve the resulting linear system of equations. Gauss substitutions and the integrand subtractions enforce a triangular system solving strategy, which proceeds top-down, from the pentuple-cut to the single-cut. In this fashion, because of the integrand subtractions, the integrand which has to be evaluated numerically gets updated at any level, cut-by-cut, by the subtraction of the polynomial residues determined at the previous step. The sampling and the determination of the coefficients in Samurai proceeds with a projection technique based on the Discrete Fourier Transform [29,31].
In the new version GoSam-2.0, the amplitude decomposition is obtained by a new integrand-reduction method [34], implemented in the C++ code Ninja [32,33], which is the default reduction library.
In Ref. [34] an improved version of the integrand reduction method for one-loop amplitudes was presented. This method allows, whenever the analytic dependence of the integrand on the loop momentum is known, to extract the unknown coefficients of the residues by performing a Laurent expansion of the integrand with respect to one of the free loop components which are not constrained by the corresponding on-shell conditions.
Within the Laurent expansion approach, the system of equations for the coefficients becomes diagonal. In fact, in the asymptotic limit, both the integrand and the higher-point subtractions exhibit the same polynomial behaviour as the residue. Therefore one can identify the unknown coefficients with the ones of the expansion of the integrand, corrected by the contributions coming from higherpoint residues. In other words, the subtractions of higher-point contributions are replaced by corrections at the coefficient level. Because of the universal structure of the residues, the parametric form of these corrections can be computed once and for all, in terms of a subset of the higher-point coefficients.
This novel D-dimensional unitarity-based algorithm is lighter than the original integrand reduction method, because less coefficients need to be determined, and turns out to be faster and numerically more accurate.
The integrand reduction via Laurent expansion has been implemented in the library Ninja, where the Laurent expansions of the integrands are performed by a polynomial division between some parametric expansions of the numerator and the uncut denominators. The expansions of the numerator, required by Ninja as input, are efficiently generated by GoSam using FORM, after collecting the terms that do not depend on the loop momentum into global abbreviations.
Ninja and the new version of Samurai, as well as Golem95C, all distributed with the GoSam-2.0 package, can deal with processes where the masses of the internal particles are complex, and where the rank r of the numerator of the integrands can exceed the number N of denominators by one unit, i.e. r ≤ N + 1, as it may happen e.g. in effective theories (see also Section 3.5.1).

The extension derive
The derive feature generates code to access the tensor coefficients of each diagram or group of diagrams individually. While it has been among the possible keywords for the extensions options in GoSam-1.0 already, it now has been promoted to be used by default in the context of tensorial reconstruction [25]. It improves both the speed and the precision of tensorial reconstruction and makes connection to other reduction methods.
The idea behind it is to compute the numerator N (q) from a Taylor series In this form one can read off a one-to-one correspondence between derivatives at q = 0 and the coefficients of the tensor integrals. At a technical level, the files helicity*/d*h*l1d.f90 contain the routines derivative(µ 2 , [i 1 ], [i 2 ],...) and reconstruct d*(coeffs), where the latter is only generated in conjunction with the extension golem95, and coeffs is a type which comprises all coefficients of a diagram of a certain rank. The number of optional indices i 1 , i 2 , . . . determine which derivative should be returned. The subroutine reconstruct d* also takes into account the proper symmetrisation.

Electroweak scheme choice
Regularisation and renormalisation within the Standard Model can be performed using various schemes, which also may differ in the set of electroweak parameters considered as input parameters, while other electroweak parameters are derived ones. Within GoSam Table 1 Possible choices to select the electroweak scheme. To simplify the notation we write the sine of the weak mixing angle as sw. The lists of derived parameters contain only the parameters which are computed from the input parameters and used in the expressions for the amplitudes.
by setting appropriately the flag model.options, depending on whether the scheme might be changed after the generation of the code or not. By default, when the flag is not set in the input card, GoSam generates a code which uses m W , m Z and α as input parameters, allowing however to change this in the generated code, by setting the variable ewchoice in the configuration file common/config.f90 to the desired value. The user can choose among 8 different possibilities, which are listed in Table 1. When the electric charge e is set algebraically to one, the schemes 6 − 8 cannot be used.
The flag model.options in the input card allows one also to directly set the values of the different parameters appearing in the model. If the values of exactly three electroweak parameters are specified, GoSam automatically takes them as input parameters. In that case, in order to be able to switch among different schemes after code generation, the variable ewchoose also must be added to the model.options flag.

Stability tests and rescue system
Within the context of numerical and semi-numerical techniques, we should be able to assess in real time, for each phase space point, the level of precision of the corresponding one-loop matrix element. Whenever a phase space point is found in which the quality of the result falls below a certain threshold, either the point is discarded or the evaluation of the amplitude is repeated by means of a safer, albeit less efficient procedure. This procedure is traditionally called "rescue system".
Apart from improvements in the stability of the reduction itself, which are provided by the new versions of Samurai and Golem95C, and in particular by the new reduction algorithm Ninja, the new version of GoSam also has a more refined rescue system as compared to version 1.0.
Looking at the literature, we observe that various techniques for detecting points with low precision have been implemented within the different automated tools for the evaluation of one-loop virtual corrections.
A first commonly used approach relies on the comparison between the numerical values of the infrared pole coefficients computed by the automated tool with their known analytic results dictated by the universal behaviour of the infrared singularities [64]. We refer to this as the pole test.
The main advantages of this method are its broad applicability to all amplitudes and the fact that it requires a negligible additional computation time. However, since not all integrals which appear in the reconstruction of the amplitude give a contribution to the double and single poles, this method often provides an overestimate of the precision, which might result in keeping phase space points whose finite part is less precise than what is predicted by the poles.
Different techniques have been proposed that target directly the precision of the finite part. Using the symmetry properties of scattering amplitudes under scaling of all physical scales, or alternatively the invariance under rotation of the momenta, we can build pairs of points that should provide identical results, both for the finite parts and for the poles, and use the difference between them as an estimator of the precision.
The scaling test [65], is based on the properties of scaling of scattering amplitudes when all physical scales (momenta, renormalisation scale, masses) are rescaled by a common multiplicative factor x. As shown in [65], this method provides a very good correlation between the estimated precision and the actual precision of the finite parts.
The rotation test [33] exploits the invariance of the scattering amplitudes under an azimuthal rotation about the beam axis, namely the direction of the initial colliding particles. Whenever the initial particles are not directed along the beam axis, one can perform a rotation of all particles by an arbitrary angle in the space of momenta. A validation of this technique, and the corresponding correlation plots, has been presented in [33].
While the scaling and the rotation test provide a more reliable estimate of the precision of the finite parts that enter in the phase space integration, their downside is that they require two evaluations of the same matrix element, therefore leading to a doubling in the computational time.
Additional methods have been proposed, within the context of integrandreduction approaches, which target the relations between the coefficients before integration, namely the reconstructed algebraic expressions for the numerator function before integration, known as N = N tests [16,23]. This kind of tests can be applied to the full amplitude (global N = N test) or individually within each residue of individual cuts (local N = N test). The drawback of this technique comes from the fact that the test is applied at the level of individual diagrams, rather than on the final result summed over all diagrams, making the construction of a rescue system quite cumbersome.
For the precision analysis contained in GoSam-2.0, and to set the trigger for the rescue system, we decided to employ a hybrid method, that takes advantage of the computational speed of the pole test, combined with the higher reliability of the rotation test. This hybrid method requires setting three different thresholds. After computing the matrix elements, GoSam-2.0 checks the precision δ pole of the single pole with the pole test. Comparing the single pole S IR that can be obtained from the general structure of infrared singularities and the one provided by GoSam-2.0, which we label S, we define (2) 11 The corresponding estimate of the number of correct digits in the result is provided by P pole = − log 10 (δ pole ). This step does not require any increase in computational time. The value of P pole is then compared with two thresholds P high and P low . If P pole > P high the point is automatically accepted. Given the high quality of the computed pole, the finite part is very unlikely to be so poor that the point should be discarded.
If P pole < P low the point is automatically discarded, or sent to the rescue system. If already the pole has a low precision, we can expect the finite part to be of the same level or worse.
In the intermediate region where P high > P pole > P low , it is more difficult to determine the quality of the result solely based on the pole coefficients. Only in this case the point is recalculated using the rotation test, which requires additional computational time.
If we call the finite part of the amplitudes evaluated before and after the rotation A fin and A fin rot respectively, we can define the error δ rot estimated with the rotation as and the corresponding estimate on the number of correct digits as P rot = − log 10 (δ rot ). P rot provides a reliable estimate of the precision of the finite part [33], and can be compared with a threshold P set to decide whether the point should be accepted or discarded.
The values of the three thresholds P high , P low and P set can be chosen by the user, to adjust the selection mechanism to the fluctuations in precision which occur between different processes. In the input card, P high , P low and P set correspond to PSP chk th1, PSP chk th2 and PSP chk th3, respectively, see appendix A. It is worth to notice that the rotation test can be bypassed simply by setting the initial thresholds P high = P low . In this case the selection is performed solely on the basis of the pole test.

Higher rank integrals
The libraries Ninja, Golem95C and Samurai all support integrals with tensor ranks r exceeding the number of propagators N . Such integrals occur for example in effective theories (a prominent example is the effective coupling of gluons to the Higgs boson), or in calculations involving spin-two particles beyond the leading order. These extensions are described in detail in Refs. [31,34,58] and are contained in the distribution of GoSam-2.0. The additional integrals will be called automatically by GoSam if they occur in an amplitude, such that the user can calculate amplitudes involving higher rank integrals without additional effort. Ninja and Samurai provide higher rank integrals for rank r = N + 1, version 1.3 of Golem95C [58] provides higher rank integrals and the tensorial reconstruction routines up to r = N + 1 for N ≤ 6, as well as form factors up to rank ten for N ≤ 4. More details about the higher rank integrals are given in Appendix B.

Production of colour-/spin correlated trees
GoSam can also generate tree level amplitudes in a spin-and colour-correlated form. Colour correlated matrix elements are defined as and we define spin-correlated matrix elements as The spin-correlated matrix element (as well as the colour correlated matrix element) contains implicitly the sum over all non-specified helicities, while only the helicities with the indices i and j are fixed, i.e.
These matrix elements are particularly useful in combination with Monte Carlo programs which use these trees to build the dipole subtraction terms for the infrared divergent real radiation part. With these modified tree level matrix elements GoSam is able to generate all necessary building blocks for a complete NLO calculation. Such a setup has been used successfully in combination with the framework of Herwig++/Matchbox [66][67][68].

Support of complex masses
The integral libraries contained in the GoSam package as well as the GoSam code itself fully support complex masses. The latter are needed for the treatment of unstable fermions and gauge bosons via the introduction of the corresponding decay width. A fully consistent treatment of complex W -and Z-boson masses requires the use of the complex mass scheme [69]. According to this scheme the boson masses become Gauge invariance requires that the definition of the weak mixing angle has to be modified accordingly: To make use of the complex mass scheme, we introduce two new model files, sm complex and smdiag complex, which contain the Standard Model with complex mass scheme, the first with a full CKM matrix, the latter with a diagonal unit matrix for the CKM matrix. An example dealing with a complex top quark mass is given in Section 5.

Installation
The user can download the code either as a tar-ball or from a subversion repository at http://gosam.hepforge.org The installation of GoSam-2.0 is very simple when using the installation script. The latter can be downloaded by Upon installation, the installer will ask some questions, which are described in more detail in the manual [50], which also can be downloaded from the webpage given above.
To use the default installation all the questions can be "answered" by pressing the ENTER key.
In particular, the installer will check if QGRAF [51] and FORM [52,53] already exist on the system. If they are not found, one can either press ENTER to have them installed by the script, or provide a path to the binary (tab-completion can be used). If they are found, their version is checked, and if needed the installation of a version which has been tested to run with GoSam is suggested.
As soon as all questions are answered, the main installation process will start. The components will be downloaded, built and installed. The whole procedure can take about 10-30 minutes.
At the end, a script gosam setup env.sh will be created in the bin/ subdirectory of the install location, which will set (temporarily) all environment variables as soon as the script is sourced into a shell (with the command source [path]/gosam setup env.sh). The installer also gives a recommendation how these environment variables can be set permanently. The script can be used in all tcshand bash-compatible shells.
All files which have been installed are tracked in the installer-log.ini file. It is important to keep this file and the install script. They are needed to update and uninstall GoSam. For the default installation, internet access is required.
The program GoSam is designed to run in any modern Unix-like environment (Linux, Mac). The system requirements are Python (≥ 2.6), a Fortran compiler (gfortran or ifort), a C/C++ compiler (gcc/icc), and (GNU) make. Compatibility with gcc versions 4.2.-4.9 as well as clang has been tested. By default, GoSam uses the gfortran/gcc compilers from the GNU Compiler Suite. To use an Intel compiler (ifort/icc), the --intel option can be used. Specific paths to the compilers can be provided using the --fc, --cc, --cxx options. All further options can be listed by invoking the installation script with the flag help: gosam installer.py --help.

Using GoSam
We first start describing the use of GoSam in the standalone version. For the use in combination with a Monte Carlo program, based on the BLHA interface, we refer to Section 4.3.
In order to generate the matrix element for a given process the user should create a process specific setup file, process.in, which we call process card. An example process card for the process e + e − → tt is given in Appendix A, where we explain each entry in detail.
It is recommended to generate and modify a template file for the process card. This can be done by invoking the shell command gosam.py --template process.in which generates the file process.in with some documentation for all defined options. The options are filled with the default values, and some paths are set by the installation script. User-defined options changing the default values can also be set in a global configuration file. The script will search in the GoSam-2.0 directory, in the user's home directory and in the current working directory for a file named '.gosam' or 'gosam.in'.
In order to generate the Fortran code for the process specified in the input card one needs to invoke gosam.py process.in

Structure of the generated code
The generated process directory will have the following sub-directory structure: codegen: This directory contains files which are only relevant for code generation. These files will therefore not be included in a tar-ball created with make dist.
common: Contains Fortran files which are common to all helicity amplitudes and to the constructed matrix element code. The file config.f90 contains some global settings, the file model.f90 contains the definitions and settings for the model parameters. This directory is always compiled first. doc: Contains files which are necessary for creating doc/process.pdf, which displays all Feynman diagrams of the Born level and one-loop amplitude, together with colour and helicity info.
helicity*: These directories contain all files for a specific helicity amplitude.

Interfacing to Monte Carlo programs
The interface of GoSam with a Monte Carlo event generator program is based on the Binoth-Les Houches Accord (BLHA) standards. GoSam-2.0 supports both BLHA1 [48] and BLHA2 [49]. Certainly, a dedicated interface without using the BLHA is also possible, and such an interface with MadGraph/MadDipole/MadEvent [70][71][72][73] has been built and applied successfully in various phenomenological applications [38,41,43,46,47]. If GoSam is used as a One Loop Provider (OLP), the Monte Carlo program is steering the different stages of the calculation, in particular the phase space integration and the event generation, as illustrated in Fig. 2. Therefore, the user frontend will depend on the user interface of the Monte Carlo program. The latter will call GoSam at runtime to provide the corresponding value of the one-loop amplitude at the given phase space points.
A number of phenomenological results produced by combining GoSam with various Monte Carlo programs can be found in the literature, e.g in combination with Sherpa [37,[40][41][42]45], PowHeg [39], Herwig++/Matchbox [66]. Examples how to run GoSam with Sherpa can also be found on the Sherpa manual webpage [74] and on the GoSam process packages webpage [75]. For the interface with PowHeg, a detailed description can be found in the appendix of Ref. [39]. The interface with Herwig++/Matchbox is described in [66].

Using external model files
The GoSam-2.0 package comes with the built-in model files sm, smdiag, smehc, sm complex, smdiag complex, where the latter two are needed in the case of complex masses and couplings, see Section 3.5.3. The model files smehc contain the effective Higgs-gluon couplings.
Other models can be imported easily in the UFO (Universal FeynRules Output) [76] format. The model import in the UFO format can be used in the standalone as well as the OLP mode of GoSam, where both the BLHA1 and BLHA2 standards are supported for the syntax of the model import.
A model description in the UFO format consists of a python package which the user can either generate using FeynRules [77,78] or write himself and store in any directory. In order to import the model into GoSam one needs to set the model variable in the process card (line 5 in the example process card of Appendix A), specifying the keyword FeynRules in front of the path pointing to the python files defining the model. For example, if we assume that the model description is in the directory $HOME/models/MSSM UFO, the process card should contain the line model= FeynRules,$HOME/models/MSSM UFO . The import of model files generated by LanHEP [79] is also supported. More details about the import from LanHEP are given in the

gg → H+1 jet in the heavy top mass limit
Recently GoSam was used to compute the virtual corrections for the production of a Higgs boson in association with 2 and 3 jets [37,41] in the infinite top-mass limit. As an example for this type of processes, where a special model file is needed, containing the Feynman rules for the effective vertices, which furthermore give rise to higher rank loop integrals, we consider here the process gg → H g.
An example process card for the generation of this process, and a test routine comparing a phase space point with results from analytical amplitude representations is provided among the examples of the GoSam-2.0 distribution. In the following we will refer to that example to describe some feature of this process.
In order to compute amplitudes using the effective gluon-gluon-Higgs vertices, the model smehc has to be used. This model contains also the effective vertex for the Higgs boson decaying to two photons. When setting the powers of the strong coupling using the order flag, one has to remember that the effective vertex counts as two powers of the strong coupling. To compute the virtual corrections for H + 1 jet we therefore have to set order=QCD, 3, 5.
The inclusion of the effective gluon-gluon-Higgs coupling at NLO also requires corrections of the Wilson coefficient. At NLO the Wilson coefficient is given by [80]  where the effective Lagrangian is given by The smehc model file also contains the effective vertex for the Higgs decay into a pair of photons via top-and W -loops. For the vertex factor we use the formula given by FeynRules [77]: 2 987 where x t = m H mt , x W = m H m W and the corresponding effective Lagrangian is given by Table 3 contains numerical results for gg → Hg at the phase space point shown in Table 2.

Single top production
An example containing complex masses in the loop propagators is the so called s-channel single top quark production, where the top quark has a width w t = 1.5 GeV. In Table 5 we give numerical results for the subprocess ud → νee + bb at the phase space point given in Table 4. The b-quarks are taken to be massless and the comparison has been performed against the HELAC-NLO code [8].  Table 4 Kinematic point used for ud → νee + bb. The W -boson and top-quark mass and width are set to m W = 80.25 GeV, w W = 0, mt = 170.9 GeV and wt = 1.5 GeV.
Results for ud → νee + bb with the kinematic point from Table 4.   Table 6 Kinematic point used for uū → G → γγ.

Graviton production within models of large extra dimensions
As an example for the usage of GoSam with a model file different from the Standard Model we consider the QCD corrections to graviton production in ADD models [81,82] with large extra dimensions (LED). The corresponding model files in UFO [76] format, which we generated using FeynRules [77,78], are located in the subdirectory examples/model/LED UFO. To import new model files within the GoSam setup, the user should specify the path to the model file in the process card. In the given example, this already has been done, i.e. the process card contains the line model=FeynRules,[gosampath]/examples/model/LED UFO. The example process we included in the GoSam-2.0 distribution is uū → G → γγ, where G denotes a graviton, and the program calculates the virtual QCD corrections. Note that this example also involves integrals where the rank exceeds the number of propagators, due to the spin-2 nature of the graviton. Running make test in the subdirectory examples/uu graviton yy should produce the result shown in Table 7, using the phase space point given in Table 6. The full process, including an additional jet, has been calculated in [47], where we refer to for details about the parameter settings.
Results for uū → G → γγ with the kinematic point from Table 6.  Table 7 Result for the virtual amplitude including the QCD corrections to uū → G → γγ within ADD models of large extra dimensions.

Conclusions
We have presented the program package GoSam-2.0, which is a highly automated tool to calculate one-loop multi-particle amplitudes. As the amplitudes at a first stage are produced in an algebraic form, the program offers a lot of flexibility concerning the particle content and the couplings, the choice of the reduction method and the treatment of the rational parts.
GoSam-2.0 can be used to calculate NLO QCD corrections both within and Beyond the Standard Model, as well as electroweak corrections, in combination with a Monte Carlo program providing the tree-level and NLO real radiation parts. The latter can be interfaced using the Binoth-Les-Houches-Accord, where both BLHA1 and BLHA2 standards are supported. The automated interface to various Monte Carlo programs also offers the possibility to produce parton showered events and to compare different shower Monte Carlo event generators at NLO level.
We also note that the structure of the code is favourable to be used as a building block for the one-loop virtual times singly unresolved real radiation part entering NNLO calculations.
GoSam-2.0 contains many important new features. The installation procedure is extremely simple: all dependencies are provided in one package, and an install script is building the whole package in a completely automated way. Setting up a process is also very user-friendly: the user only has to fill out a well documented text file, the process card, where the program automatically chooses appropriate default values for unspecified options.
Improvements in the code generation compared to version 1.0 lead to more compact and faster code. GoSam-2.0 also contains a new integrand reduction method, the integrand decomposition via Laurent expansion, implemented in the library Ninja, which leads to a considerable gain in stability and speed, in particular for amplitudes containing internal masses.
The range of applicability of GoSam also has been extended considerably. In particular, integrals where the rank exceeds the number of propagators (needed e.g. in effective theories) are fully supported, and propagators for spin-2 particles are implemented. The complex mass scheme is supported, including the complexification of the couplings, and several electroweak schemes can be chosen. Moreover, a new system for stability tests and the rescue of 'unstable' phase space points has been implemented. In addition, the program offers the possibility to produce spin-and colour correlated tree-level matrix elements. As a consequence, GoSam-2.0 can provide all the building blocks needed by modern Monte Carlo programs to construct a full NLO event generator, for QCD corrections both within and beyond the Standard Model, as well as electroweak corrections.
Therefore, to follow the strive for precision in the next phases of LHC data taking as well as at a future Linear Collider, not only regarding QCD corrections, GoSam-2.0 can serve as a highly valuable tool.
Here we give a commented example for the process e + e − → tt.
In the following it is assumed that the process e + e − → tt should be calculated to order O(ααs) (QCD corrections). We neglect the exchange of a Z or a Higgs boson and treat the electron as massless. The output directory is assumed to be in the relative path eett. A template file for a generic process card (called eett.in here) can be generated by invoking the shell command gosam.py --template eett.in The template file eett.in then should be edited by the user to define the process specifications. All lines starting with # are comments.
At this point we would like to emphasize that almost all specifications in the process card are options, which will take default values if they are not filled in by the user. The paths to the libraries will be inserted automatically by the install script. The only mandatory fields are the in and out particles, the perturbative order and the path where to store the process files. Therefore, a minimal process card can look like this: In order to populate the process subdirectory specified under process path with files for code generation one invokes gosam.py eett.in This will create the subdirectory structure described in Section 4.2.
In the following, we will give detailed comments to all the fields and options available in the process card for the example eett.in. (Please note that the line numbers on the left are only included for better readability and should not be included in your input file).  The comments to the file eett.in are as follows.
1 Setting a process name is optional but recommended. All module names will be prefixed with the process name (e.g. precision → eett precision). This will avoid name conflicts if at a later stage more than one matrix elements are linked into one executable. 2 The item process path specifies the directory to which all generated files and directories are written. Specification of a process path is mandatory. 3-4 The items in and out specify the particles of the initial and final state. The particle names must be defined in the selected model file. As the model files usually define mnemonics for the particle names there might be several ways of specifying the same process. Instead of 'e+' one could have written 'ep' or 'positron'. For a complete list of alternative particle names please refer to the documentation of the corresponding model file. Specifying in and out particles is mandatory. 5 The option model specifies which model files should be used in order to generate and evaluate the diagrams. How to import models in UFO or LanHep format is described in Section 4.4. The default for this field is smdiag, i.e. the built-in Standard Model file with a diagonal CKM matrix. 6 The option model.options can be used to pass options which are specific to a certain model. The default is ewchoose, which means that the electro-weak scheme is selected automatically according to the given input parameters. 7 The item order is a comma separated list with three entries. The first entry specifies a symbol that denotes a coupling constant. In the Standard Model file sm the only two possibilities are 'gs' for the strong coupling constant gs and 'e' for the electroweak coupling. The second number is the power of the chosen coupling constant for the tree-level diagrams and the third number specifies the power of that coupling constant for the one-loop diagrams. Note that the numbers refer to the powers in the diagrams of the amplitude rather than the squared amplitude. In the above example the string 'gs, 0, 2' specifies that the tree-level diagrams should be of order g 0 s and the one-loop diagrams should be of order g 2 s and an unspecified power of e in both cases. If there is no tree level, i.e. the process is loop induced, the keyword NONE should be put as second item in the list, instead of the tree level power of the coupling. The values of order are translated into a vsum constraint in the file qgraf.dat. This field is mandatory. 8-9 The keywords zero and one specify a set of symbols that should be treated as zero (resp. one). These simplifications are applied at the symbolical level.
Only symbols that appear in the FORM interface of the model file should be specified here (masses, couplings, CKM-matrix elements, etc). In the example we specify the electron mass 'me' to be zero and we do not keep the coupling constants in the calculation explicitly (gs = e = 1). These options can be omitted. 10 The option regularisation scheme allows to choose the dimensional regularisation scheme, in our example dred for dimensional reduction, which is the default. cdr for "conventional dimensional regularisation" is also possible. 11 helicites: a comma separated list of helicities to be calculated. An empty list means that all possible helicities should be generated. The characters correspond to particles 1, 2, ... from left to right. Example: e + e − → γγ: Only three helicity configurations are required; the other ones are either zero or can be obtained by symmetry transformations. This corresponds to helicities=+-++,+-+-,+---Multiple helicities can be encoded in patterns, which are expanded at the time of code generation. For more details we refer to the manual. 12 qgraf.options=onshell,notadpole,nosnail: a list of options which is passed to QGRAF via the 'options' line. Possible values (as of qgraf.3.1.1) are the following keywords: onepi, onshell, nosigma, nosnail, notadpole, floop, topol. In our example, it means that external lines are on-shell, i.e. do not contain selfenergy corrections, and that tadpole and snail diagrams are discarded. We refer to the QGRAF documentation for more details. 13-16 The value of the option qgraf.verbatim is passed verbatim to the file qgraf.dat.
In our example we suppress the generation of diagrams containing Higgs and Z bosons. As these commands are passed verbatim to QGRAF, no mnemonic names are allowed here, e.g. the Higgs particle has to be denoted by 'H' and cannot be replaced by 'h'. For a complete list of available options, please consult the QGRAF manual. For a complete list of particle names we refer to the documentation of the corresponding model file. These options can be omitted. 17 polvec: by default (polvec=numerical), numerical polarisation vectors are used for the massless gauge bosons, rather than producing separate code for each helicity (see Section 3.1.3). To switch off the use of numerical polarisation vectors, use polvec=explicit. 18 diagsum: if True, one-loop diagrams sharing some propagators are combined before the algebraic reduction. The default is diagsum = True. 19 The option reduction programs allows to choose the amplitude reduction method.
If several choices are given, the code is produced such that the reduction methods can be switched at runtime. The default is ninja, golem95. 20 extensions: this option contains a list of useful extensions to the core of the program, which operate at the code generation stage. The currently available extensions are In this example, the gluon (particle 1) takes the momentum k 2 as reference momentum for the polarisation vector. The massive top quark (particle 3) uses the light-cone projection l 4 of the W-boson as reference direction for its own momentum splitting. Similarly, the momentum of the W-boson is split into a direction l 4 and one along l 3 . 45 abbrev.limit: maximum number of instructions per subroutine when calculating abbreviations. The default is 0, which means that no maximum is set. 46 templates: path pointing to the directory containing the template files for the process. If not set, the program uses the directory gosam path /templates. The directory must contain a file called template.xml. 47 qgraf.bin: path to the QGraf executable. The default path will be set by the installation script. 48 form.bin: path to the FORM executable. The default path will be set by the installation script. 49 form.threads: the number of FORM threads when using tform, the parallel version of FORM. The default is 2. 50 form.tempdir: the temporary directory where FORM can store (large) intermediate files. the default is /tmp. 51 haggies.bin: path to the haggies executable. The default path will be set by the installation script. 52 fc.bin: path to the Fortran compiler. The default path will be set by the installation script. 53 python.bin: path to the python executable. The default path will be set by the installation script. 54 ninja.fcflags: compiler flags to compile with Ninja. The default will be set by the installation script. 55 ninja.ldflags: ldflags required to link the Ninja library. The default will be set by the installation script. 56 samurai.fcflags: compiler flags to compile with Samurai. The default will be set by the installation script. 57 samurai.ldflags: ldflags required to link the Samurai library. The default will be set by the installation script. 58 golem95.fcflags: compiler flags to compile with Golem95C. The default will be set by the installation script. 59 golem95.ldflags: ldflags required to link the Golem95C library. The default will be set by the installation script. 60 r2: treatment of the rational part R 2 . The possibilities are: implicit: µ 2 terms are kept in the numerator and reduced at runtime, explicit: µ 2 terms are reduced analytically, off: all µ 2 terms are set to zero. The default is r2=explicit. 61 symmetries: this information is used when the list of helicity configurations is generated. An empty list means that all helicity configurations will be generated, even if some of them could be mapped onto each other. Possible values are: flavour: assumes that no flavour changing interactions are present. When calculating the list of helicities, fermions with PDG codess 1-6 are assumed not to mix. n = h : restriction of particle helicities, e.g. 1=-, 2=+ specifies helicities of particles 1 and 2. -% n = h restriction by PDG code, e.g. %23=+-specifies the helicity of all Z-bosons to be '+' and '-' only (no '0' polarisation), % n refers to both +n and −n, %+ n refers to +n only, %n refers to −n only. 62 crossings: a list of crossed processes derived from this process. For each process in the list a module similar to matrix.f90 is generated. Example: process name=ddx uux in=1,-1 out=2,-2 crossings = dxd uux: -1 1 → 2 -2, ud ud: 2 1 → 2 1

B: Higher rank integrals
Higher rank integrals are implemented in all reduction libraries included in GoSam. Ninja and Samurai are based on integrand reduction, as described in Section 3.2, Golem95C provides tensor integrals, using a tensor reduction method and a basis of scalar integrals which has been designed to provide numerical stability in problematic phase space regions, for example in the limit of small Gram determinants.
In the following we briefly sketch the main features of the higher rank extensions for both approaches, more details can be found in [31,34,58].

B.1: Integrand reduction approach
If the rank r of a one-loop integrand is not larger than the number of propagators N , the respective integral can be written as the following combination of known master integrals In the case where r = N + 1 the integral is generalized as The three integrals in Eq. (B.4) whose numerator is proportional to µ 2 are finite and contribute to the rational part of the amplitude. They have been computed in Ref. [27,34] and they read where s ij ≡ (p i − p j ) 2 . The tadpole of rank 2 appearing in Eq. (B.4) can be written in terms of the scalar tadpole I i1 as The latter can be computed using the formulas of Ref. [84], as a function of form factors of scalar integrals B 0 . In the special case with s i2i1 = 0 we use Eq. (A.6.2) and (A.6.3) of that reference. For the general case s i2i1 = 0 we use instead the following formula [32] B 111 (s i2i1 , m 2 i1 , m 2 (B.10)

B.2: Tensor reduction approach
In the tensor reduction approach, the tensor integrals are written in terms of linear combinations of scalar form factors and all possible combinations of external momenta and metric tensors carrying the Lorentz structure. The form factors themselves are then reduced to a convenient set of basis integrals. It is well known that, due to the 4-dimensionality (resp. D = 4 − 2 dimensionality in dimensional regularisation) of space-time, integrals with N ≥ 6 can be reduced iteratively to 5-point integrals. Therefore form factors for N ≥ 6 are never needed. The general form factor decomposition of an arbitrary tensor integral can be written as where ∆ µ ij = r µ i − r µ j are differences of external momenta r, and qa = k + ra. The notation [· · · ] {µ1···µr} {a1···ar} stands for the distribution of the r Lorentz indices µ i , and the momentum labels a i , to the vectors ∆ µi j ai and metric tensors in all distinguishable ways. Note that the choice r N = 0, a i = N ∀ i leads to the well known representation in terms of external momenta where the labels a i are not necessary, but we prefer a completely shift invariant notation here.
S denotes an ordered set of propagator labels, corresponding to the momenta forming the kinematic matrix S, defined by We should point out that the form factors of type D N,r j1...jr−6 and beyond, i.e. form factors associated with three or more metric tensors, are not needed for integrals where the rank r does not exceed the number N of propagators, no matter what the value of N is, because integrals with N ≥ 6 can be reduced algebraically to pentagons.
The program Golem95C reduces the form factors A, . . . , D internally to a set of basis integrals, i.e. the endpoints of the reduction (they do not form a basis in the mathematical sense, as some of them are not independent). The choice of the basis integrals can have important effects on the numerical stability in certain kinematic regions. Our reduction endpoints are 4-point functions in 6 dimensions I 6 4 , which are IR and UV finite, 4-point functions in D + 4 dimensions, and various 3-point, 2-point and 1-point functions. A special feature of Golem95C is that the algebraic reduction to scalar basis integrals is automatically replaced by a stable and fast one-dimensional numerical integration of parametric integrals corresponding to tensor rather than scalar integrals in kinematic situations where a further reduction would lead to spurious inverse Gram determinants tending to zero. This leads to improved numerical stability.
The extension of Golem95C to higher rank integrals [58] follows the reduction formalism as outlined in [24]. However, the extension of the formalism to rank six pentagons required some care, as the latter develop an UV divergence, and therefore O( ) terms occurring in the reduction need to be taken into account at intermediate stages.
The rational parts of all the integrals contained in Golem95C can be extracted separately, and analytic formulae for r ≤ N are provided in [85]. The results for those integrals which are relevant for the higher rank extension can be extracted from [27], where formulae for all possible rational parts are given in a general form. The ones which are relevant for the higher rank extension which have not been given above already are listed explicitly here, where the notation conventions are [g µ1µ2 g µ3µ4 + g µ1µ3 g µ2µ4 + g µ1µ4 g µ2µ3 ] + O( ) ,