Automated one-loop calculations with GoSam

The European Physical Journal C, Mar 2012

We present the program package GoSam which is designed for the automated calculation of one-loop amplitudes for multi-particle processes in renormalisable quantum field theories. The amplitudes, which are generated in terms of Feynman diagrams, can be reduced using either D-dimensional integrand-level decomposition or tensor reduction. GoSam can be used to calculate one-loop QCD and/or electroweak corrections to Standard Model processes and offers the flexibility to link model files for theories Beyond the Standard Model. A standard interface to programs calculating real radiation is also implemented. We demonstrate the flexibility of the program by presenting examples of processes with up to six external legs attached to the loop.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1140%2Fepjc%2Fs10052-012-1889-1.pdf

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


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1140%2Fepjc%2Fs10052-012-1889-1.pdf

Gavin Cullen, Nicolas Greiner, Gudrun Heinrich. Automated one-loop calculations with GoSam, The European Physical Journal C, 2012, 1889, DOI: 10.1140/epjc/s10052-012-1889-1