The automated computation of treelevel and nexttoleading order differential cross sections, and their matching to parton shower simulations
Abstract:
We discuss the theoretical bases that underpin the automation of the computations of treelevel and nexttoleading order cross sections, of their matching to parton shower simulations, and of the merging of matched samples that differ by lightparton multiplicities. We present a computer program, [email protected], capable of handling all these computations – partonlevel fixed order, showermatched, merged – in a unified framework whose defining features are flexibility, high level of parallelisation, and human intervention limited to input physics quantities. We demonstrate the potential of the program by presenting selected phenomenological applications relevant to the LHC and to a 1TeV collider. While nexttoleading order results are restricted to QCD corrections to SM processes in the first public version, we show that from the user viewpoint no changes have to be expected in the case of corrections due to any given renormalisable Lagrangian, and that the implementation of these are well under way.
CP31418, LPN14066
MCNET1409, ZUTH 14/14
1 Introduction
Quantum Chromo Dynamics is more than forty years old [1, 2], and perturbative calculations of observables beyond the leading order are almost as old, as is clearly documented in several pioneering works (see e.g. refs. [3, 4, 5, 6, 7, 8, 9]), where the implications of asymptotic freedom had been quickly understood. The primary motivation for such early works was a theoretical one, stemming from the distinctive features of QCD (in particular, the involved infrared structure, and the fact that its asymptotic states are not physical), which imply the need of several concepts (such as infrared safety, hadronparton duality, and the factorisation of universal longdistance effects) that come to rescue, and supplement, perturbation theory. On the other hand, the phenomenological necessity of taking higherorder effects into account was also acknowledged quite rapidly, in view of the structure of jet events in collisions and of the extraction of from data.
Despite this very early start, the task of computing observables beyond the Born level in QCD has remained, until very recently, a highly nontrivial affair: the complexity of the problem, due to both the calculations of the (tree and loop) matrix elements and the need of cancelling the infrared singularities arising from them, has generated a very significant theoretical activity by a numerous community. More often than not, different cases (observables and/or processes) have been tackled in different manners, with the introduction of adhoc solutions. This situation has been satisfactory for a long while, given that beyondBorn results are necessary only when precision is key (and, to a lesser extent, when large factors are relevant), and when many hard and wellseparated jets are crucial for the definition of a signature; these conditions have characterized just a handful of cases in the past, especially in hadronic collisions (e.g., the production of single vector bosons, jet pairs, or heavy quark pairs).
The advent of the LHC has radically changed the picture since, in a still relatively short running time, it has essentially turned hadronic physics into a highprecision domain, and one where events turning up in large tails are in fact not so rare, in spite of being characterised by small probabilities. Furthermore, the absence so far of clear signals of physics beyond the Standard Model implies an increased dependence of discovery strategies upon theoretical predictions for known phenomena. These two facts show that presently the phenomenological motivations are extremely strong for higherorder and multileg computations of all observables of relevance to LHC analyses.
While a general solution is not known for the problem of computing exactly the perturbative expansion for any observable up to an arbitrarily large order in , if one restricts oneself to the case of the first order beyond the Born one (nexttoleading order, NLO henceforth), then such a solution does actually exist; in other words, there is no need for adhoc strategies, regardless of the complexity of the process under study. This remarkable fact results from two equally important theoretical achievements. Namely, a universal formalism for the cancellation of infrared singularities [10, 11, 12, 13, 14], and a technique for the algorithmic evaluation of renormalised oneloop amplitudes [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], both of which must work in a process and observable independent manner. At the NLO (as opposed to the NNLO and beyond) there is the further advantage that fixedorder computations can be matched to partonshower event generators (with either the [email protected] [26] or the POWHEG [27] method – see also refs. [28, 29, 30, 31, 32, 33, 34, 35, 36, 37] for early, lessdeveloped, or newer related approaches), thus enlarging immensely the scope of the former, and increasing significantly the predictive power of the latter.
It is important to stress that while so far we have explicitly considered the case of QCD corrections, the basic theoretical ideas at the core of the subtraction of infrared singularities, of the computation of oneloop matrix elements, and of the matching to parton showers will require no, or minimal, changes in the context of other renormalisable theories, QCD being essentially a worstcase scenario. This is evidently true for treelevel multileg computations, as is proven by the flexibility and generality of tools such as MadGraph5 [38], that is able to deal with basically any userdefined Lagrangian.
In summary, there are both the phenomenological motivations and the theoretical understanding for setting up a general framework for the computation of (any number of) arbitrary observables in an arbitrary process at the tree level or at the NLO, with or without the matching to parton showers. We believe that the most effective way of achieving this goal is that of automating the whole procedure, whose technological challenges can be tackled with highlevel computer languages capable of dealing with abstract concepts, and which are readily available.
The aim of this paper is that of showing that the programme sketched above has been realised, in the form of a fully automated and public computer code, dubbed [email protected]. As the name suggests, such a code merges in a unique framework all the features of MadGraph5 and of [email protected], and thus supersedes both of them (and must be used in their place). It also includes several new capabilities that were not available in these codes, most notably those relevant to the merging of event samples with different lightparton multiplicities. We point out that [email protected] contains all ingredients (the very few external dependencies that are needed are included in the package) that are necessary to perform an NLO, possibly plus shower (with the [email protected] formalism), computation: it thus is the first public (since Dec. 16, 2013) code, and so far also the only one, with these characteristics. Particular attention has been devoted to the fact that calculations must be doable by someone who is not familiar with Quantum Field Theories, and specifically with QCD. We also show, in general as well as with explicit examples, how the construction of our framework lends itself naturally to its extension to NLO corrections in theories other than QCD, in keeping with the fact that such a flexibility is one of the standard features of the treelevel computations which were so far performed by MadGraph, and that has been inherited by [email protected].
It is perhaps superfluous to point out that approaches to automated computations constitute a field of research which has a long history, but which has had an exponential growth in the past few years, out of the necessities and possibilities outlined above. The number of codes which have been developed, either restricted to leading order (LO henceforth) predictions [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 38, 53], or including NLO capabilities [54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82] is truly staggering. The level of automation and the physics scope of such codes, not to mention other, perhaps less crucial, characteristics, is extremely diverse, and we shall make no attempt to review this matter here.
We have organized this paper as follows. In sect. 2, we review the theoretical bases of our work, and discuss new features relevant to future developments. In sect. 3 we explain how computations are performed. Section 4 presents some illustrative results, relevant to a variety of situations: total cross sections at the LHC and future colliders, differential distributions in collisions, and benchmark oneloop pointwise predictions, in the Standard Model and beyond. We finally conclude in sect. 5. Some technicalities are reported in appendices A to D.
2 Theoretical bases and recent progress
At the core of [email protected] lies the capability of computing treelevel and oneloop amplitudes for arbitrary processes. Such computations are then used to predict physical observables with different perturbative accuracies and finalstate descriptions. Since there are quite a few possibilities, we list them explicitly here, roughly in order of increasing complexity, and we give them short names that will render their identification unambiguous in what follows.

fLO: this is a tree and partonlevel computation, where the exponents of the coupling constants are the smallest for which a scattering amplitude is non zero. No shower is involved, and observables are reconstructed by means of the very particles that appear in the matrix elements.

fNLO: the same as fLO, except for the fact that the perturbative accuracy is the NLO one. As such, the computation will involve both treelevel and oneloop matrix elements.

LO+PS: uses the matrix elements of an fLO computation, but matches them to parton showers. Therefore, the observables will have to be reconstructed by using the particles that emerge from the Monte Carlo simulation.

NLO+PS: same as LO+PS, except for the fact that the underlying computation is an NLO rather than an LO one. In this paper, the matching of the NLO matrix elements with parton showers is done according to the [email protected] formalism.

MLMmerged: combines several LO+PS samples, which differ by finalstate multiplicities (at the matrixelement level). In our framework, two different approaches, called jet and shower schemes, may be employed.

FxFxmerged: combines several NLO+PS samples, which differ by finalstate multiplicities.
We would like to stress the fact that having all of these different simulation possibilities embedded in a single, processindependent framework allows one to investigate multiple scenarios while being guaranteed of their mutual consistency (including that of the physical parameters such as coupling and masses), and while keeping the technicalities to a minimum (since the majority of them are common to all types of simulations). For example, one may want to study the impact of perturbative corrections with (NLO+PS vs LO+PS) or without (fNLO vs fLO) the inclusion of a parton shower. Or to assess the effects of the showers at the LO (LO+PS vs fLO) and at the NLO (NLO+PS vs fNLO). Or to check how the inclusion of differentmultiplicity matrix elements can improve the predictions based on a fixedmultiplicity underlying computation, at the LO (MLMmerged vs LO+PS) and at the NLO (FxFxmerged vs NLO+PS).
In the remainder of this section we shall review the theoretical ideas that constitute the bases of the computations listed above in items 1–6. Since such a background is immense, we shall sketch the main characteristics in the briefest possible manner, and rather discuss recent advancements that have not yet been published.
2.1 Methodology of computation
The central idea of [email protected] is the same as that of the MadGraph family. Namely, that the structure of a cross section, regardless of the theory under consideration and of the perturbative order, is essentially independent of the process, and as such it can be written in a computer code once and for all. For example, phase spaces can be defined in full generality, leaving only the particle masses and their number as free parameters (see e.g. ref. [83]). Analogously, in order to write the infrared subtractions that render an NLO cross section finite, one just needs to cover a handful of cases, which can be done in a universal manner. Conversely, matrix elements are obviously theory and processdependent, but can be computed starting from a very limited number of formal instructions, such as Feynman rules or recursion relations. Thus, [email protected] is constructed as a metacode, that is a (Python) code that writes a (Python, C++, Fortran) code, the latter being the one specific to the desired process. In order to do so, it needs two ingredients:

a theory model;

a set of processindependent building blocks.
A theory model is equivalent to the Lagrangian of the theory plus its parameters, such as couplings and masses. Currently, the method of choice for constructing the model given a Lagrangian is that of deriving its Feynman rules, that [email protected] will eventually use to assemble the matrix elements. At the LO, such a procedure is fully automated in FeynRules [84, 85, 86, 87, 88, 89]. NLO cross sections pose some extra difficulties, because Feynman rules are not sufficient for a complete calculation – one needs at least UV counterterms, possibly plus other rules necessary to carry out the reduction of oneloop amplitudes (we shall generically denote the latter by , adopting the notation of the OssolaPapadopoulosPittau method [18]). These NLOspecific terms are presently not computed by FeynRules^{1}^{1}1We expect they will in the next public version [90], since development versions exist that are already capable of doing so – see e.g. sect. 4.3. and have to be supplied by hand^{2}^{2}2Note that these are a finite and typically small number of processindependent quantities., as was done for QCD corrections to SM processes. Therefore, while the details are unimportant here, one has to bear in mind that there are “LO” and “NLO” models to be employed in [email protected] – the former being those that could also be adopted by MadGraph5, and the latter being the only ones that permit the user to exploit the NLO capabilities of [email protected].
Given a process and a model, [email protected] will build the processspecific code (which it will then proceed to integrate, unweight, and so forth) by performing two different operations. a) The writing of the matrix elements, by computing Feynman diagrams in order to define the corresponding helicity amplitudes, using the rules specified by the model. b) Minimal editing of the processindependent building blocks. In the examples given before, this corresponds to writing definite values for particles masses and the number of particles, and to select the relevant subtraction terms, which is simply done by assigning appropriate values to particle identities. The building blocks modified in this manner will call the functions constructed in a). Needless to say, these operations are performed automatically, and the user will not play any role in them.
We conclude this section by emphasising a point which should already be clear from the previous discussion to the reader familiar with recent MadGraph developments. Namely that, in keeping with the strategy introduced in MadGraph5 [38], we do not include among the processindependent building blocks the routines associated with elementary quantities (such as vertices and currents), whose roles used to be played by the HELAS routines [91] in previous MadGraph versions [39, 51]. Indeed, the analogues of those routines are now automatically and dynamically created by the module ALOHA [92] (which is embedded in [email protected]), which does so by gathering directly the relevant information from the model, when this is written in the Universal FeynRules Output (UFO [93]) format. See sect. 2.3.1 for more details on this matter.
2.2 General features for SM and BSM physics
Since the release of MadGraph5 a significant effort, whose results are now included in [email protected], went into extending the flexibility of the code at both the input and the output levels. While the latter is mostly a technical development (see appendix B.6), which allows one to use different parts of the code as standalone libraries and to write them in various computer languages, the former extends the physics scope of [email protected] w.r.t. that of MadGraph5 in several ways, and in particular for what concerns the capability of handling quantities (e.g., form factors, particles with spin larger than one, and so forth) that are relevant to BSM theories. Such an extension, to be discussed in the remainder of this section and partly in sect. 2.3.1, goes in parallel with the analogous enhanced capabilities of FeynRules, and focuses on models and on their use. Thus, it is an overarching theme of relevance to both LO and NLO simulations, in view of the future complete automation of the latter in theories other than the SM. Some of the topics to which significant work has been lately devoted in [email protected], and which deserve highlighting, are the following:

Complex mass scheme (sect. 2.2.1).

Support of various features, of special relevance to BSM physics (sect. 2.2.2).

Improvements to the FeynRules/UFO/ALOHA chain (sect. 2.3.1).

Output formats and standalone libraries (appendix B.6).

Feynman gauge in the SM (sect. 2.4.2).

Improvements to the frontend user interface (the [email protected] shell – appendix A).

Hierarchy of couplings: models that feature more than one coupling constant order them in terms of their strengths, so that for processes with several coupling combinations at the cross section level only the assumed numericallyleading contributions will be simulated (unless the user asks otherwise) – sect. 2.4 and appendix B.1.
We would finally like to emphasise that in the case of LO computations, be they fLO, LO+PS, or merged, one can always obtain from the shortdistance cross sections a set of physical unweighted events. The same is not true at the NLO: fNLO cross sections cannot be unweighted, and unweighted [email protected] events are not physical if not showered. This difference explains why at the LO we often adopt the strategy of performing computations with small, selfcontained modules whose inputs are Les Houches event (LHE henceforth) files [94, 95], while at the NLO this is typically not worth the effort – compare e.g. appendices B.3 and B.4, where the computation of scale and PDF uncertainties at the NLO and LO, respectively, is considered. Further examples of modules relevant to LO simulations are given in sect. 2.3.3.
2.2.1 Complex mass scheme
In a significant number of cases, the presence of unstable particles in perturbative calculations can be dealt with by using the Narrow Width Approximation (NWA)^{3}^{3}3See sect. 2.5 for a general discussion of the NWA and of other related approaches.. However, when one is interested in studying either those kinematical regions that correspond to such unstable particles being very offshell, or the production of broad resonances, or very involved final states, it is often necessary to go beyond the NWA. In such cases, one needs to perform a complete calculation, in order to take fully into account offshell effects, spin correlations, and interference with nonresonant backgrounds in the presence of possibly large gauge cancellations. Apart from technical difficulties, the inclusion of all resonant and nonresonant diagrams does not provide one with a straightforward solution, since the (necessary) inclusion of the widths of the unstable particles – which amounts to a resummation of a specific subset of terms that appear at all orders in perturbation theory – leads, if done naively, to a violation of gauge invariance. While this problem can be evaded in several ways at the LO (see e.g. refs. [96, 97, 98, 99, 100, 101]), the inclusion of NLO effects complicates things further. Currently, the most pragmatic and widelyused solution is the socalled complex mass scheme [102, 103], that basically amounts to an analytic continuation in the complex plane of the parameters that enter the SM Lagrangian, and which are related to the masses of the unstable particles. Such a scheme can be shown to maintain gauge invariance and unitarity at least at the NLO, which is precisely our goal here.
In [email protected] it is possible to employ the complex mass scheme in the context of both LO and NLO simulations, by instructing the code to use a model that includes the analytical continuation mentioned above (see the Setup part of appendix B.1 for an explicit example). For example, at the LO this operation simply amounts to upgrading the model that describes SM physics in the following way [102]:

The masses of the unstable particles are replaced by .

The EW scheme is chosen where , , and (the former a real number, the latter two complex numbers defined as in the previous item) are input parameters.

All other parameters (e.g., and ) assume complex values. In particular, Yukawa couplings are defined by using the complex masses introduced in the first item.
At the NLO, the necessity of performing UV renormalisation introduces additional complications. At present, the prescriptions of ref. [103] have been explicitly included and validated. As was already mentioned in sect. 2.1, this operation will not be necessary in the future, when it will be replaced by an automatic procedure performed by FeynRules.
2.2.2 BSMspecific capabilities
One of the main motivations to have very precise SM predictions, and therefore to include higher order corrections, is that of achieving a better experimental sensitivity in the context of New Physics (NP) searches. At the same time, it is necessary to have as flexible, versatile, and accurate simulations as is possible not only for the plethora of NP models proposed so far, but for those to be proposed in the future as well. These capabilities have been one of the most useful aspects of the MadGraph5 suite; they are still available in [email protected] and, in fact, have been further extended.
As was already mentioned, the required flexibility is a direct consequence of using the UFO models generated by dedicated packages such as FeynRules or SARAH [104], and of making [email protected] compatible with them. Here, we limit ourselves to listing the several extensions recently made to the UFO format and the [email protected] code, which have a direct bearing on BSM simulations.

The possibility for the user to define the analytic expression for the propagator of a given particle in the model [105].

The possibility to define form factors (i.e., coupling constants depending on the kinematic of the vertex) directly in the UFO model. Note that, since form factors cannot be derived from a Lagrangian, they cannot be dealt with by programs like FeynRules, and have therefore to be coded by hand in the models. In order to be completely generic, the possibility is also given to go beyond the format of the current UFO syntax, and to code the form factors directly in Fortran^{4}^{4}4Which obviously implies that this option is not available should other type of outputs be chosen (see appendix B.6). Further details are given at:
https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/FormFactors.. 
Models and processes are supported that feature massive and massless particles of spin [105]. This implies that all spins are supported in the set .

The support of multifermion interactions, including the case of identical particles in the final state, and of UFO models that feature vertices with more than one fermion flow. Multifermion interactions with fermionflow violation, such as in the presence of Majorana particles, are not supported. Such interactions, however, can be implemented by the user by splitting the interaction in multiple pieces connected via heavy scalar particles, a procedure that allows one to define unambiguously the fermion flow associated with each vertex.
While not improved with respect to what was done in MadGraph5, we remind the reader that the module responsible for handling the colour algebra is capable of treating particles whose SU representation and interactions are nontrivial, such as the sextet and type vertices respectively.
2.3 LO computations
The general techniques and strategies used in [email protected] to integrate a treelevel partonic cross section, and to obtain a set of unweighted events from it, have been inherited from MadGraph5; the most recent developments associated with them have been presented in ref. [38], and will not be repeated here. After the release of MadGraph5, a few optimisations have been introduced in [email protected], in order to make it more efficient and flexible than its predecessor. Here, we limit ourselves to listing the two which have the largest overall impact.

The phasespace integration of decaychain topologies has been rewritten, in order to speed up the computations and to deal with extremely long decay chains (which can now easily extend up to sixteen particles). In addition, the code has also been optimised to better take into account invariantmass cuts, and to better handle the case where interference effects are large.

It is now possible to integrate matrix elements which are not positivedefinite^{5}^{5}5We stress that this statement is nontrivial just because it applies to LO computations. In the context of NLO simulations this is the standard situation, and [email protected] has obviously been always capable of handling it.. This is useful e.g. when one wants to study a process whose amplitude can be written as a sum of two terms, which one loosely identifies with a “signal” and a “background”. Rather than integrating , one may consider and separately, which is a more flexible approach (e.g. by giving one the possibility of scanning a parameter space, in a way that affects while leaving invariant), that also helps reduce the computation time by a significant amount. One example of this situation is the case when is numerically dominant, and thus corresponds to a very large sample of events which can however be generated only once, while smaller samples of events, that correspond to for different parameter choices, can be generated as many times as necessary. Another example is that of an effective theory where a process exists () that is also present in the SM (). In such a case, it is typically which is kept at the lowest order in (with being the cutoff scale). Finally, this functionality is needed in order to study the large expansion in multiparton amplitudes, where beyond the leading terms the positive definiteness of the integrand is not guaranteed.
In the following sections, we shall discuss various topics relevant to the calculation of LOaccurate physical quantities. Sect. 2.3.1 briefly reviews the techniques employed in the generation of treelevel amplitudes, emphasising the role of recent UFO/ALOHA developments. Sect. 2.3.2 presents the module that computes the total widths of all unstable particles featured in a given model. Sect. 2.3.3 describes reweighting techniques. Finally, in sect. 2.3.4 we review the situation of MLMmerged computations.
2.3.1 Generation of treelevel matrix elements
The computation of amplitudes at the tree level in [email protected] has a scope which is in fact broader than treelevel physics simulations, since all matrix elements used in both LO and NLO computations are effectively constructed by using treelevel techniques. While this is obvious for all amplitudes which are not oneloop ones, in the case of the latter it is a consequence of the Lcutting procedure, which was presented in detail in ref. [68] and which, roughly speaking, stems from the observation that any oneloop diagram can be turned into a treelevel one by cutting one of the propagators that enter the loop. Furthermore, as was also explained in ref. [68] and will be discussed in sect. 2.4.2, all of the companion operations of oneloop matrix element computations (namely, UV renormalisation and counterterms) can also be achieved through treelevellike calculations, which are thus very central to the whole [email protected] framework.
The construction of treelevel amplitudes in [email protected] is based on three key elements: Feynman diagrams, helicity amplitudes, and colour decomposition. Helicity amplitudes [107, 108, 109, 110, 111, 112, 113] provide a convenient and effective way to evaluate matrix elements for any process in terms of complex numbers, which is quicker and less involved than one based on the contraction of Lorentz indices. As the name implies, helicity amplitudes are computed with the polarizations of the external particles fixed. Then, by employing colour decompositions [114, 115, 116], they can be organised into gaugeinvariant subsets (often called dual amplitudes), each corresponding to an element of a colour basis. In this way, the complexity of the calculation grows linearly with the number of diagrams instead of quadratically; furthermore, the colour matrix that appears in the squared amplitude can be easily computed automatically (to any order in ) once and for all, and then stored in memory. If the number of QCD partons entering the scattering process is not too large (say, up to five or six), this procedure is manageable notwidthstanding the fact that the number of Feynman diagrams might grow factorially. Otherwise, other techniques that go beyond the Feynmandiagram expansion have to be employed [45, 117, 49, 52]. The algorithm used in [email protected] for the determination of the Feynman diagrams has been described in detail in ref. [38]. There, it has been shown that it is possible to efficiently “factorise” diagrams, such that if a particular substructure shows up in several of them, it only needs to be calculated once, thus significantly increasing the speed of the calculation. In addition, a notyetpublic version of the algorithm can determine directly dual amplitudes by generating only the relevant Feynman diagrams, thereby reducing the possible factorial growth to less than an exponential one.
The diagramgeneration algorithm of [email protected] is completely general, though it needs as an input the Feynman rules corresponding to the Lagrangian of a given theory. The information on such Feynman rules is typically provided by FeynRules, in a dedicated format (UFO). We remind the reader that FeynRules is a Mathematicabased package that, given a theory in the form of a list of fields, parameters and a Lagrangian, returns the associated Feynman rules in a form suitable for matrix element generators. It now supports renormalisable as well as nonrenormalisable theories, twocomponent fermions, spin and spin fields, superspace notation and calculations, automatic mass diagonalization and the UFO interface. In turn, a UFO model is a standalone Python module, that features selfcontained definitions for all classes which represent particles, parameters, and so forth. With the information from the UFO, the dedicated routines that will actually perform the computation of the elementary blocks that enter helicity amplitudes are built by ALOHA. Amplitudes are then constructed by initializing a set of external wavefunctions, given their helicities and momenta. The wavefunctions are next combined, according to the interactions present in the Lagrangian, to form currents attached to the internal lines. Once all of the currents are determined, they are combined to calculate the complex number that corresponds to the amplitude for the diagram under consideration. Amplitudes associated with different diagrams are then added (as complex numbers), and squared by making use of the colour matrix calculated previously, so as to give the final result. We point out that versions of MadGraph earlier than MadGraph5 used the HELAS [91, 118] library instead of ALOHA. By adopting the latter, a significant number of limitations inherent to the former could be lifted. A few examples follow here. ALOHA is not forced to deal with precoded Lorentz structures; although its current implementation of the Lorentz algebra assumes four spacetime dimensions, this could be trivially generalised to any even integer, as the algebra is symbolic and its actual representation enters only at the final stage of the output writing; its flexibility has allowed the implementation of the complex mass scheme (see sect. 2.2.1), and of generic UV and computations (see sect. 2.4.2); it includes features needed to run on GPU’s, and the analogues of the Heget [119, 120] libraries can be automatically generated for any BSM model; finally, it caters to matrixelement generators other than MadGraph, such as those now used in Herwig++ [121, 122] and Pythia8 [123]. Since its release in 2011 [92], several important improvements have been made in ALOHA. On top of the support to the complex mass scheme and to the Feynman gauge, and of specific features relevant to oneloop computations (which are discussed in sects. 2.2.1 and 2.4.2), we would like to mention here that there has been a conspicuous gain in speed, both at the routinegeneration phase as well as in the actual evaluation, thanks to the extensive use of caching. In addition, the user is now allowed to define the form of a propagator (which is not provided by FeynRules), while previously this was determined by the particle spin: nontrivial forms, such as those relevant to spin2 particles in ADD models [124], or to unparticles, can now be used.
2.3.2 Width calculator
Among the basic ingredients for Monte Carlo simulations of newphysics models are the masses and widths of unstable particles; which particle is stable and which unstable may depend on the particular parameter benchmark chosen. Masses are typically obtained by going to the mass eigenbasis and, if necessary, by evolving boundary conditions from a large scale down to the EW one. Very general codes exist that perform these operations starting from a FeynRules model, such as AsperGe [89]. The determination of the corresponding widths, on the other hand, requires the explicit calculation of all possible decay channels into lighter (SM or BSM) states. The higher the number of the latter, the more daunting it is to accomplish this task by hand. Furthermore, depending on the mass hierarchy and interactions among the particles, the computation of twobody decay rates could be insufficient, as highermultiplicity decays might be the dominant modes for some of the particles. The decay channels that are kinematically allowed are highly dependent on the mass spectrum of the model, so that the decay rates need to be reevaluated for every choice of the input parameters. The program MadWidth [106] has been introduced in order to address the above issues. In particular, MadWidth is able to compute partial widths for body decays, with arbitrary values of , at the treelevel and by working in the narrowwidth approximation^{6}^{6}6Even if those two assumptions are quite generic, there are particles for which they do not give sufficientlyaccurate results, such as the Standard Model Higgs, which has significant loopinduced decay modes.. The core of MadWidth is based on new routines for diagram generation that have been specifically designed to remove certain classes of diagrams:

Diagrams with onshell intermediate particles. If the kinematics of an particles decay allows an internal particle , that appears in an channel, to be on shell, the corresponding diagram can be seen as a cascade of two decays, particles followed by particles. It is thus already taken into account in the calculation of lowermultiplicity decay channels, and the diagram is discarded.

Radiative diagrams. Roughly speaking, if one or more zeromass particles are radiated by another particle, the diagram is considered to be a radiative correction to a lowermultiplicity decay – the interested reader can find the precise definition in ref. [106]. Such a diagram is therefore discarded, because it should be considered only in the context of a higherorder calculation. Furthermore, all diagrams with the same couplingconstant combination and the same external states are also discarded, so that gauge invariance is preserved.
MadWidth begins by generating all twobody decay diagrams, and then iteratively adds extra final state particles with the condition that any diagram belonging to either of the classes above is forbidden. This iterative procedure stops when all modes have been considered, or estimated to be numerically irrelevant^{7}^{7}7Both of these conditions can be controlled by the user.. All diagrams thus generated are integrated numerically. MadWidth uses several methods to reduce significantly the overall computation time. Firstly, it features two fast (and conservative) estimators, one for guessing the impact of adding one extra finalstate particle before the actual diagramgeneration phase, and another one for evaluating the importance of a single integration channel. Both of these estimators are used to neglect parts of the computation which are numerically irrelevant. Secondly, if the model is compatible with the recent UFO extension of ref. [106], and thus includes the analytical formulae for twobody decays, then the code automatically uses those formulae and avoids the corresponding numerical integrations.
We conclude this section by remarking that, although essential for performing BSM crosssection computations, MadWidth should be seen as a complement for the existing tools that generate models. This is because the information it provides one with must be available before the integration of the matrix elements is carried out but, at the same time, really cannot be included in a model itself, since it depends on the chosen benchmark scenario. For more details on the use of this modelcomplementing feature in [email protected], and of its massmatrix diagonalisation analogue to which we have alluded above, see the Setup part of appendix B.1.
2.3.3 Event reweighting
The generation of large samples of events for experimental analyses can be a very timeconsuming operation, especially if it involves a full simulation of the detector response. It is therefore convenient, whenever possible, to apply corrections, or to study systematics of theoretical or modelling nature, by using reweighting techniques. The reweighting of either onedimensional distributions or that performed on a eventbyevent basis are equallycommon practices in experimental physics. Although the basic idea (that a nonnull function can be used to map any other function defined in the same domain) behind these two procedures is the same, one must bear in mind that they are not identical, and in particular that the former can never be proven to be formally correct, since correlations (with other, nonreweighted variables) may be lost or modified, while the latter is correct in general, at least in the limit of a large number of events. For this reason, we consider only eventbyevent reweighting approaches in what follows.
Thanks to its flexibility and the possibility of accessing a large amount of information in a direct and straightforward manner (model and running parameters, matrix elements, PDFs, and so forth) [email protected] provides one with an ideal framework for the implementation of such approaches. The strategy is rather simple: one starts with a set of hard events, such as those contained in an LHE file, and rescales their weights:
(1) 
without modifying their kinematics. The rescaling factor is not a constant, and may change on an eventbyevent basis. This implies that, even when the original sample of events is unweighted (i.e., the ’s are all equal), the reweighted one will be in general weighted (i.e., the ’s may be different), and therefore degraded from a statistical point of view. If, however, the spread of the new weights is not too large (i.e., the ’s are close to each other, and feature a small number of outliers), the reweigthing is typically advantageous with respect to generating a new independent sample from scratch.
While eq. (1) is completely general, its practical implementation depends on the kind of problems one wants to solve. We shall consider three of them in this section, two of which constitute a direct application of the basic formula (1), and a third one which is more sophisticated. Since these procedures address different types of physics, they are conveniently associated with different modules in [email protected], but they are all fully embedded into our framework and easily accessible through it, as we shall briefly explain in what follows and in appendices B.4 and B.5.
The simplest example is that of the evaluation of the uncertainties through to the variations of the renormalisation and factorisation scales, and of the PDFs. In such a case is easily determined by using the identities and of the initialstate partons, and the power () of in the Born matrix elements^{8}^{8}8There may be cases where such matrix elements do not factorise a single factor – see e.g. sect. 2.4. We limit ourselves to discussing the simplest, and more common, situation here.:
(2) 
It should be stressed that, although scale and PDF systematics can also be computed with reweighting techniques at the NLO (see ref. [125]), in general they cannot be written in the very simple form of eq. (1), which is that of an overall rescaling. For a direct comparison of the NLO and LO techniques employed in this context by [email protected], and for fuller details about them, see appendices B.3 and B.4 respectively.
Another rather simple example is that of the case where one is interested in studying the implications of changing the modelling of a process, with the sole constraint that its initial and final states be the same. Such a situation can for example occur when the numerical values of the coupling constants are modified, or when the contributions of classes of diagrams are included or eliminated (e.g., Higgs exchange in EW vector boson scattering). The common feature of all examples of this kind is that they are associated with changes to matrix elements, all computed with the same kinematic configuration. Therefore, one can simply write:
(3) 
which for obvious reasons is dubbed matrixelement reweighting. Despite its simplicity, this method has very many different applications, from parameter scanning in newphysics models (where one starts with a single sample of events, that corresponds to a given benchmark point, and generates as many new ones as the number of parameter configurations of interest), to more advanced approaches, such as the inclusion of exact loop effects () in processes that can be also, and more simply, described by effective theories () – see refs. [126, 127] for recent results of the latter type that make use of MadGraph5 and [email protected]. Some further comments on matrixelement reweighting and its practical usage in [email protected] are given in appendix B.5.
We finally turn to discussing the matrixelement method [128, 129, 130], which can be seen as a reweighting one because the weights determined at the parton level through matrixelement computations are possibly modified by a convolution to take into account a variety of blurring effects (such as those due to a detector). On the other hand, from the practical point of view it turns out to be more convenient, rather than talking about reweighting factors, to introduce a likelihood for the observation of a given kinematic configuration (q) given a set of theoretical assumptions (). By doing so, one can just reexpress eq. (3) in a different language:
(4) 
with a suitable volume factor, and the total rate associated with the given assumptions . The advantage of eq. (4) is that it is suitable to handle cases which are much more complicated than the purelytheoretical exercise of eq. (3), which has led to its introduction here. For example, in the first approximation one may think of as the actual kinematic configuration measured by an experiment, whose accuracy is such that it can be directly used as an argument of the matrix elements, as is done in eq. (4). Note that this is a rather strong assumption, that in practice identifies hadronlevel with partonlevel quantities, and assumes that the knowledge of the final state is complete (such as that which one can ideally obtain in DrellYan production, ). It is clear that there are many ways in which this simple approximation (which is used, for example, in refs. [131, 132, 133, 134, 135]) can break down: the effects of radiation, of the imperfect knowledge of the detector, of the impossibility of a strict identification of parton with hadronlevel quantities, the uncertainties that plague the latter, the fact that all the relevant fourmomenta cannot in general be measured, are but a few examples. It is therefore necessary to generalise eq. (4), which one can do as follows:
(5) 
In eq. (5), we have denoted by p the parton kinematic configuration. All effects that may turn p into the detectorlevel quantity q (whose dimensionality therefore need not coincide with that of p) are parametrised by , called transfer function. As for any hadroncollision measurable quantity, eq. (5) features the convolution with the PDFs. The likelihood introduced in this way can be used in the context e.g. of an hypothesis test in order to determine which among various choices of is the most probable.
Although this method is both conceptually simple and very attractive, the numerical evaluation of is difficult because the transfer function behaves in a way which cannot be probed efficiently by phasespace parametrisations that work well for just the matrix elements. In order to address this issue, a dedicated program, dubbed MadWeight [136], has been introduced that includes an optimised phasespace treatment specifically designed for eq. (5). The new version of the code [137] embedded in [email protected] features several improvements w.r.t. that of ref. [136]. It includes the method for the approximate description of higherorder effects due to initialstate radiation, as proposed in ref. [138]. Furthermore, several optimizations have been achieved that render the computations much faster (sometimes by orders of magnitude); this allows one to use this approach also in the case of very involved final states, such as those relevant to Higgs production in association with a pair [139]. Further details on MadWeight and its use within [email protected] can be found in appendix B.5.
2.3.4 Treelevel merging
The goal of merging is that of combining samples associated with different parton multiplicities in a consistent manner, that avoids double counting after showering, thus allowing one to effectively define a single fullyinclusive sample. The treelevel merging algorithms implemented in [email protected] are a hybrid version of those available in Alpgen [45] and SHERPA [140]; they work for both SM and BSM hard processes, but are fully automated only when the shower phase is performed with either Pythia6 [141] or Pythia8 [123] (however, there are no reasons in principle which prevents these schemes from working with HERWIG6 [142, 143] or Herwig++ [121, 122]). They are based on the use of a measure [144] to define hardness and to separate processes of different multiplicities, and do not perform any analyticSudakov reweighting of events; rather, this operation is effectively achieved by rejecting showered events under certain conditions (see later), which implies a direct use of the welltuned showering and hadronization mechanisms of the parton shower Monte Carlos.
There are two merging schemes that can be used in conjunction with Pythia6 and Pythia8; in the case of the latter, one is also given the possibility of considering CKKWL approaches [145, 146, 147] (after having generated the samples relevant to various parton multiplicities with [email protected]); in what follows, we shall limit ourselves to briefly describe the former two methods. Firstly, one has the jet MLM scheme [148], where finalstate partons at the matrixelement level are clustered according to a jet algorithm to find the “equivalent parton shower history” of the event. In our implementation the Feynman diagram information from [email protected] is used to retain only those clusterings that correspond to actual Feynman diagrams. In order to mimic the behaviour of the parton shower, the value for each clustering vertex associated with a QCD branching is used as the renormalisation scale for in that vertex. All factorisation scales squared, and the renormalisation scale squared for the hard process (the process with the zeroextraparton multiplicity), are constructed by clustering back to the irreducible system and by using the transverse mass in the resulting frame: . The smallest value found in the jetreconstruction procedure is restricted to be larger than some minimum cutoff scale, which we denote by ; if this condition is not satisfied, the event is rejected. The hard events are then showered by Pythia: at the end of the perturbativeshower phase, finalstate partons are clustered into jets, using the very same jet algorithm as before, with the jets required to have a transverse momentum larger than a given scale , with . The resulting jets are compared to the partons at the hard subprocess level (i.e., those that result from the matrixelement computations): a jet is said to be matched to a parton if the distance between the two, defined according to ref. [144], is smaller than the minimal jet hardness: . The event is then rejected unless each jet is matched to a parton, except in the case of the largestmultiplicity sample, where extra jets are allowed if softer than the of the softest matrixelement parton in the event, . Secondly, and with the aim to give one a nonparametric way to study merging systematics, one has the shower scheme, which can be used only with Pythia’s ordered shower. In this case, events are generated by [email protected] as described above and then showered, but information is also retained on the hardest (which is also the first, in view of the ordered nature of Pythia here) emission in the shower, ; furthermore, one sets , which cannot be done in the context of the jet MLM scheme. For all samples but the largestmultiplicity one, events are rejected if , while in the case of the largestmultiplicity sample events are rejected when . This merging scheme is simpler than the jet MLM one, but it rather effectively mimics it. Furthermore, it probes the Sudakov form factors used in the shower in a more direct manner. Finally, the treatment of the largestmultiplicity sample is fairly close to that used in the CKKWinspired merging schemes. In both the jet MLM and shower methods, merging systematics are associated with variations of ; in the former case, changes to must be done by keeping a constant. For applications of the two schemes described here, see e.g. refs. [148, 149, 126, 150, 151].
2.4 NLO computations
When discussing the problem of perturbative corrections, one should bear in mind that one usually considers an expansion in terms of a single quantity (which is a coupling constant for fixedorder computations). However, this is just a particular case of the more general scenario in which that expansion is carried out simultaneously in two or more couplings, all of which are thus treated as “small” parameters – we shall refer to such a scenario as mixedcoupling expansion. Despite the fact that there is typically a clear numerical hierarchy among these couplings, a mixedcoupling situation is far from being academic; in fact, as we shall show in the following, there are cases when one is obliged to work with it. In order to study a generic mixedcoupling expansion without being too abstract, let us consider an observable which receives contributions from processes that stem from both QCD and QED interactions. The specific nature of the interactions is in fact not particularly relevant (for example, QED here may be a keyword that also understands the pureEW contributions); what matters, for the sake of the present discussion, is that may depend on more than one coupling constant. We assume that the regular function
(6) 
admits a Taylor representation:
(7) 
which is by definition the perturbative expansion of . The first few terms of the sums in eq. (7) will be equal to zero, with the number of such vanishing terms increasing with the complexity of the process under consideration – this is because and are directly related to the number of vertices that enter a given diagram. In general, it is clear that for a given pair which gives a nonvanishing contribution to eq. (7), there may exist another pair , with , , and whose contribution to eq. (7) is also non zero. It appears therefore convenient to rewrite eq. (7) with a change of variables:
(8) 
whence:
(9) 
where
(10) 
which enforces the sum over to run only over those integers whose parity is the same as that of (therefore, there are terms in each sum over for a given ). Equation (9) implies that we need to call Born the sum (over ) of all the contributions with the smallest which are nonvanishing. Hence, the NLO corrections will correspond to the sum over of all terms with . This notation is compatible with the usual one used in the context of the perturbation theory of a single coupling: the QCD or QEDonly cases are recovered by considering or respectively. In a completely general case, for any given there will exist two integers and which satisfy the following conditions:
(11) 
and such that all contributions to eq. (9) with
(12) 
vanish, and , are both nonvanishing. This implies that in the range
(13) 
there will be at least one (two if ) nonnull contribution(s) to eq. (9) (the typical situation being actually that where all the terms in eq. (13) are nonvanishing). Given eq. (13), one can rewrite eq. (9) in the following way:
(14) 
where
(15)  
(16)  
(17) 
The coefficients of eq. (14) can be expressed in terms of the quantities that appear in eq. (9), but this is unimportant here. A typical situation is where:
(18)  
(19) 
so that:
(20)  
(21)  
(22) 
whence:
(23) 
where the Born and NLO contributions correspond to and respectively. Note that eq. (23) is the most general form of the observable if one allows the possibility of having or (or both) for , since this renders eqs. (18) and (19) always true. Equation (23) has the advantage of a straightforward interpretation of the role of NLO corrections.
An example may help make the points above more explicit. Consider the contribution to dijet production due to the partonic process ; the corresponding lowestorder  and channel Feynman diagrams feature the exchange of either a gluon or a photon (or a , but we stick to the pure theory here). The Born matrix elements will therefore be the sum of terms that factorise the following coupling combinations:
(24) 
which implies , , and . Therefore, according to eq. (23), the NLO contribution will feature the following coupling combinations:
(25) 
From the procedural point of view, it is convenient to identify QCD and QED corrections according to the relationship between one coupling combination in eq. (24) and one in eq. (25), as follows:
(26)  
(27) 
which has an immediate graphic interpretation, depicted in fig. 1.
Such an interpretation has a Feynmandiagram counterpart in the case of realemission contributions, which is made explicit once one considers cutdiagrams, like those presented in fig. 2. Loosely speaking, one can indeed identify the diagram on the left of that figure as representing QED (since the photon is cut) realemission corrections to the Born contribution. On the other hand, the diagram on the right represents QCD (since the gluon is cut) realemission corrections to the Born contribution. This immediately shows that, in spite of being useful in a technical sense, QCD and QED corrections are not physically meaningful if taken separately: in general, one must consider them both in order to arrive at a sensible, NLOcorrected result. This corresponds to the fact that a given coupling combination in the bottom row of fig. 1 can be reached by means of two different arrows when starting from the top row (i.e., the Born level).
Therefore, fig. 1 also immediately shows that when one considers only the Born term associated with the highest power of (), then QCDonly (QEDonly) corrections are sensible (because only a righttoleft or lefttoright arrow is relevant, respectively): they coincide with the NLO corrections as defined above (see the paragraph after eq. (10)). It also should be clear that the above arguments have a general validity, whatever the values of , , and in eq. (23) – the former two quantities never play a role in the analogues of fig. 1, while by increasing one simply inserts more blobs (i.e., coupling combinations) in both of the rows of fig. 1. Finally, note that reading eqs. (26) and (27) in terms of diagrams, as has been done for those of fig. 2, becomes much harder when one considers virtual contributions. For example, the one whose cutdiagram is shown in fig. 3 (and its analogues) can indeed be equally well interpreted as a QED loop correction to a QCDQCD Born cutdiagram, or as a QCD loop correction to a QCDQED Born cutdiagram.
[email protected] has been constructed by having eq. (23) in mind; although the majority of the relevant features are not yet available in the public version of the code, all of them have already been thoroughly tested in the module responsible for computing oneloop matrix elements (see sects. 2.4.2 and 4.3), which is by far the hardest from this point of view, and the checks on the realemission part are also at quite an advanced stage. The basic idea is that of giving the user the choice of which coupling combinations to retain either at the Born or at the NLO level; this corresponds to choosing a set of blobs in the upper or lower row of fig. 1, respectively. [email protected] will then automatically also consider the blobs in the row not involved in the selection by the user, in order to construct a physicallymeaningful cross section, compatible with both the user’s choices, and the constraints due to a mixedcoupling expansion (the arrows in fig. 1). It should be stressed that, although the results for the coefficients can be handled separately by [email protected], such coefficients are not (all) independent from each other from a computational viewpoint, because a single Feynman diagram (an amplitudelevel quantity) may contribute to several ’s (the latter being amplitudesquared quantities). For this reason, as far as the CPU load is concerned the choice of which coupling combinations to consider can be equivalently made at the amplitude level. Indeed, this is the only option presently available in the public version of [email protected]; more detailed explanations are given in appendix B.1.
2.4.1 NLO cross sections and FKS subtraction: MadFKS
In this section, we briefly review the FKS subtraction [10, 12] procedure, and emphasise the novelties of its implementation in [email protected] w.r.t. its previous automation in MadFKS [61], the dedicated module included in [email protected].
We shall denote by the number of finalstate particles relevant to the Born contributions to a given cross section. The set of all the partonic subprocesses that correspond to these contributions will be denoted by ; each of these subprocesses can be represented by the ordered list of the identities of its partons, thus:
(28) 
The first operation performed by [email protected] is that of constructing , given the process and the theory model. For example, if one is interested in the hadroproduction of a pair in association with a light jet
(29) 
as described by the SM, [email protected] will obtain:
(30) 
Since the processes in are treelevel, [email protected] will construct them very efficiently using the dedicated algorithms (see sect. 2.3.1). Beyond the Born level, an NLO cross section receives contributions from the oneloop and realemission matrix elements. As is well known, the set of the former subprocesses coincides^{9}^{9}9This is because we are considering here only those cases where oneloop matrix elements are obtained by multiplying the oneloop amplitudes times the Born ones. Loopinduced processes, in which the LO contribution is a oneloop amplitude squared, are not to be treated as part of an NLO computation. with that of the Born, . Realemission processes are by nature treelevel, and can therefore be obtained by using the very same algorithms as those employed to generate the Born contributions. This is achieved by making the code generate all treelevel processes that have the same finalstate as the Born’s, plus one light jet – using the example of eq. (29), these would correspond to:
(31) 
Such was the strategy adopted in the original MadFKS implementation [61]. There is however an alternative procedure, which we have implemented in [email protected] because it is more efficient than the previous one in a variety of ways. Namely, for any given , one considers all possible branchings for each nonidentical with stronglyinteracting (i.e., , , and , but also , with a quark with nonzero mass). For each of these branchings, a new list is obtained by removing from , and by inserting the pair in its place. By looping over one thus constructs the set^{10}^{10}10The exceedingly rare cases of nonsingular realemission contributions can be obtained by crossing; one example is , which is the crossed process of . of realemission processes . As a byproduct, one also naturally obtains, for each , the pairs of particles which are associated with a soft and/or a collinear singularity of the corresponding matrix element (which we denote by ); by definition [61], these pairs form the set of FKS pairs, , which is central in the FKS subtraction procedure. We point out, finally, that regardless of the type of strategy adopted to construct and , it is immediate to apply it in [email protected] to theories other than QCD, such as QED.
After having obtained , [email protected] constructs the functions that are used by FKS in order to achieve what is effectively a dynamic partition of the phase space: in each sector of such a partition, the structure of the singularities of the matrix elements is basically the simplest possible, amounting (at most) to one soft and one collinear divergence. The properties of the functions are:
(32)  
(33)  
(34)  
(35) 
One exploits eq. (35) by rewriting the real matrix elements as follows:
(36) 
Thanks to eqs. (32)(34), has the very simple singularity structure mentioned above. Furthermore, the terms in the sum on the r.h.s. of eq. (36) are independent of each other, and [email protected] is thus able to handle them in parallel.
The FKS method exploits the fact that phasespace sectors associated with different functions are independent of each other by choosing different phasespace parametrisations in each of them. There is ample freedom in such a choice, bar for two integration variables: the rescaled energy of parton (denoted by ), and the cosine of the angle between partons and (denoted by ), both defined in the incomingparton c.m. frame. The idea is that these quantities are in onetoone correspondence with the soft () and collinear () singularities respectively, which renders it particularly simple to write the subtracted cross section. The body phase space is then written as follows:
(37) 
where collectively denote the independent integration variables, with
(38) 
and where
(39) 
is the set of finalstate momenta. Given eq. (38), [email protected] chooses the other integration variables and thus determines and the functional dependence according to the form of the integrand, gathered from the underlying Feynman diagrams. In general, this implies splitting the computation in several integration channels, which are independent of each other and can be dealt with in parallel. Such multichannel technique is irrelevant here, and will be understood in the following. More details can be found in refs. [152, 153] and [61]. Implicit in eq. (37) are the maps that allow one to construct soft and collinear kinematic configurations starting from a nonspecial configuration (i.e., one where no parton is soft, and no two partons are collinear). This we shall denote as follows. Given:
(40) 
[email protected] constructs its soft, collinear, and softcollinear limits with:
(41)  
(42)  
(43) 
Furthermore, all the phasespace parametrisations employed by [email protected] are such that^{11}^{11}11Equation (44) holds for all particles except the FKSpair partons; for the latter, it is the sum of their fourmomenta that is the same in the three configurations. This is sufficient owing to the underlying infraredsafety conditions.:
(44) 
which is beneficial from the point of view of the numerical stability of observables computed at the NLO, and is necessary in view of matching with the parton shower according to the [email protected] formalism. As is usual in the context of NLO computations, we call the nonspecial and the IRlimit configurations (and, by extension, the corresponding crosssection contributions) the event and the soft, collinear, and softcollinear counterevents respectively.
Given a realemission process and an contribution, the FKSsubtracted cross section consists of four terms:
(45)  
(46) 
where is the integration measure over Bjorken ’s, is the corresponding partonluminosity factor^{12}^{12}12Whose dependence on is a consequence of eq. (44) when – see ref. [26] for more details., and the shortdistance weights are reported in refs. [61, 125]. In ref. [61], in particular, an extended proof is given that all contributions to an NLO cross section which are not naturally body ones (such as the Born, virtual, and initialstate collinear remainders) can be cast in a form formally identical to that of the soft or collinear counterterms, and can thus be dealt with simultaneously with the latter. The fully differential cross section that emerges from eqs. (45) and (46) is:
(47)  
where we have understood the complete integration over the measures on the r.h.s.. [email protected] scans the phase space by generating randomly . For each of these points, an event kinematic configuration and its weight, and a counterevent kinematic configuration and its weight are given in output; with these, any number of (IRsafe) observables can be constructed. As can be seen from eq. (47), the weight associated with the single counterevent kinematics is the sum of the soft, collinear, softcollinear, Born, virtual, and initialstate collinear remainders contributions, which reduces the probability of misbinning and thus increases the numerical stability of the result.
In order for the results of eq. (47) to be physical, they must still be summed over all processes and all FKS pairs . As far as the latter sum is concerned, it is easy to exploit the symmetries due to identical finalstate particles, and thus to arrive at the set [61]:
(48) 
whose elements give nonidentical contributions to the sum over FKS pairs. Therefore:
(49) 
with a suitable symmetry factor. The sum on the r.h.s. of eq. (49) is obviously more convenient to perform than that on the l.h.s.; it is customary to include the symmetry factor in the shortdistance weights (see e.g. ref. [125]). The number of elements in can indeed be much smaller than that in . For example, when is a purely gluonic process, we have (i.e., independent of ), while . While the former figure typically increases when quarks are considered, it remains true that, for asymptotically large ’s, is a constant while scales as .
The sum on the r.h.s. of eq. (49) is what has been originally implemented in MadFKS. It emphasises the role of the realemission processes, which implies that for quantities which are naturally Bornlike (such as the Born matrix elements themselves) one needs to devise a way to map unambiguously onto . The interested reader can find the definition of such a map in sect. 6.2 of ref. [61], which we summarise here using the Born cross section as an example:
(50) 
where
(51) 
In other words, through the r.h.s. of eq. (50) one is able to define a Bornlevel quantity as a function of a realemission process. In [email protected] we have followed the opposite strategy, namely that of defining realemission level quantities as functions of a Born process. The proof that this can indeed be done is given in appendix E of ref. [61]; here, we limit ourselves to summarising it as follows, using the event contribution to the NLO cross section as an example. One has the identity:
(52) 
Here we have defined:
(53) 
where the generalised Kronecker symbol is equal to one if its arguments are equal, and to zero otherwise, and
(54) 
Although the first sum in eq. (53) might seem to involve a very large number of terms, [email protected] knows immediately which terms will give a nonzero contribution, thanks to the procedure used to construct which was outlined at the beginning of this section. On top of organising the sums over processes and FKS pairs in a different way w.r.t. the first version of MadFKS, [email protected] also performs some of these (and, specifically, those in eq. (53)) using MC techniques (whereas all sums are performed explicitly in ref. [61]). In particular, for any given , one realemission process and one FKS pair are chosen randomly among those which contribute to eq. (53); these choices are subject to importance sampling, and are thus adaptive. In summary, while the procedure adopted originally in MadFKS takes a viewpoint from the realemission level, that adopted in [email protected] emphasises the role of Born processes. The two are fully equivalent, but the latter is more efficient in the cases of complicated processes, and it offers some further advantages in the context of matching with parton showers.
2.4.2 Oneloop matrix elements: MadLoop
Both [email protected] and its predecessor [email protected] are capable of computing the virtual contribution to an NLO cross section in a completely independent manner (while still allowing one to interface to a thirdparty oneloop provider if so desired), through a module dubbed MadLoop [68]. However, there are very significant differences between the MadLoop embedded in [email protected] (i.e., the one documented in ref. [68]), and the MadLoop currently available in [email protected]; hence, in order to avoid confusion between the two, we shall call the former MadLoop4, and the latter MadLoop5. The aim of this section is that of reviewing the techniques for automated oneloop numerical computations, and of giving the first public documentation of MadLoop5.
As the above naming scheme suggests, core functionalities relevant to the handling of treelevel amplitudes were inherited from MadGraph4 in MadLoop4, while MadLoop5 uses MadGraph5. This was a necessary improvement in view of the possibility of computing virtual corrections in arbitrary renormalisable models (i.e., other than the SM). More in general, one can identify the following three items as strategic capabilities, that were lacking in MadLoop4, and that are now available in MadLoop5^{13}^{13}13Some of them are not yet public, but are fully tested.:

The adoption of the procedures introduced with MadGraph5, and in particular the UFO/ALOHA chain for constructing amplitudes starting from a given model.
It should be clear that these capabilities induce an extremely significant broadening of the scope of MadLoop (in particular, extending it beyond the SM). As a mere byproduct, they have also completely lifted the limitations affecting MadLoop4, which were described in sect. 4 of ref. [68]. It is instructive to see explicitly how this has happened. Item A. is responsible for lifting MadLoop4 limitation #1 (MadLoop4 cannot generate a process whose Born contains a fourgluon vertex, because the corresponding routines necessary in the OPP reduction were never validated, owing to the technicallyawkward procedure for handling them in MadGraph4; this step is now completely bypassed thanks to the UFO/ALOHA chain). Limitation #2 (MadLoop4 cannot compute some loops that feature massive vector bosons, which is actually a limitation of CutTools [156], in turn due to the use of the unitary gauge) is now simply absent because of the possibility of adopting the Feynman gauge, thanks again to item A. Limitation #4 (MadLoop4 cannot handle finitewidth effects in loops) is removed thanks to the implementation of the complex mass scheme, a consequence of item A. Finally, item C. lifts MadLoop4 limitation #3 (MadLoop4 cannot generate a process if different contributions to the Born amplitudes do not factorise the same powers of all the relevant coupling constants).
The advances presented in items A.–C. above are underpinned by many technical differences and improvements w.r.t. MadLoop4. Here, we limit ourselves to listing the most significant ones:

MadLoop5 is written in Python, the same language which was adopted for MadGraph5 (MadLoop4 was written in C++).

The UVrenormalisation procedure has been improved and rendered fully general (that of MadLoop4 had several hardcoded simplifying solutions, motivated by QCD).

An extensive use is made of the optimisations proposed in ref. [24] (OpenLoops).

The selfdiagnostic numericalstability tests, and the procedures for fixing numericallyunstable loopintegral reductions, have been redesigned.
More details on these points will be given in what follows. Before turning to that, we shall discuss the basic principles used by MadLoop5 for the automation of the computation of oneloop integrals.
Generalities
Given a partonic process (see eq. (28)), MadLoop computes the quantity:
(55) 
with and being the relevant treelevel and UVrenormalised oneloop amplitudes respectively; the averages over initialstate colour and spin degrees of freedom are understood. The result for is given as a set of three numbers, corresponding to the residues of the double and single IR poles, and the finite part, all in the ’t HooftVeltman scheme [157]. In the case of a mixedcoupling expansion, each of these three numbers is replaced by a set of coefficients, in the form given by eq. (23). There may be processes for which is identically equal to zero, and is finite; for these processes (called loopinduced), MadLoop computes the quantity:
(56) 
Only eq. (55) is relevant to NLO computations proper, and we shall mostly deal with it in what follows. In the current version of [email protected], the loopinduced cannot be automatically integrated (for want of an automated procedure for multichannel integration), and hence in this case the code output by MadLoop5 must be interfaced in an adhoc way to any MC integrator (including [email protected] – see e.g. ref. [127] for a recent application).
The basic quantity which MadLoop needs to compute and eventually renormalise in order to obtain that enters eq. (55) is the oneloop UVunrenormalised amplitude:
(57) 
where denotes the contribution of a single Feynman diagram after loop integration, whose possible colour, helicity, and Lorentz indices need not be specified here, and are understood; they will be reinstated later. A standard technique for the evaluation of is a socalled reduction procedure, pioneered by Passarino and Veltman [154], which can be written as follows:
(58) 
The quantities are oneloop integrals, independent of . The essence of any reduction procedure is that it is an algebraic operation that determines the coefficients and (which are functions of external momenta and of masses, and some of which may be equal to zero); the intricacies of loop integration are dealt once and for all in the computations of the ’s (which are much simpler than any ). As the leftmost equality in eq. (58) indicates, from the operator point of view is the identity; its meaning is that of replacing with the linear combination in the rightmost member of eq. (58). As the notation suggests, different reduction procedures can possibly make use of different sets of oneloop integrals.
Equation (58) is basically what one would do if one were to compute in a nonautomated manner. In automated approaches, however, additional problems arise, for example due to the necessity of relying on numerical methods, which have obvious difficulties in dealing with the noninteger dimensions needed in the context of dimensional regularisation, and with the analytical information on the integrand of , which is extensively used in nonautomated reductions. In order to discuss how these issues can be solved, let us write in the following form:
(59) 
where , and we have assumed the diagram to have propagators in the loop and have defined:
(60) 
with the mass of the particle relevant to the loop propagator, and some linear combination of external momenta. For any fourdimensional quantity , its dimensional counterpart is denoted by , and its dimensional one by . The fact that , rather than , enters eq. (60) is a consequence of the use of the ’t HooftVeltman scheme. The loop momentum is decomposed as follows:
(61) 
with similar decompositions holding for the Dirac matrices and metric tensor . One can thus define [158] the purely fourdimensional part of the numerator that appears in eq. (59):
(62) 
from whence one can obtain its dimensional counterpart:
(63) 
The quantity defined in eq. (62), not involving noninteger dimensions, can be treated by a computer with ordinary techniques. By using eq. (63) in eq. (59) one obtains:
(64) 
where
(65)  
(66) 
Both integrals in eq. (65) and eq. (66) still depend on dimensional quantities, but they do so in a way that allows one to further manipulate and cast them in a form suitable for a fully numerical treatment. In particular, one can show [158] that the computation of is equivalent to that of a treelevel amplitude, constructed with a universal set of theorydependent rules (see ref. [159], refs. [160, 161, 162], and refs. [163, 164, 165] for the QCD, QED+EW, and some BSM cases respectively), analogous to the Feynman ones and that can be derived once and for all (for each model) by just considering the oneparticleirreducible amplitudes with up to four external legs [166]. On the other hand, eq. (65) is still potentially divergent in four dimensions. The details of how this is dealt with may vary, but the common characteristic is that all of them are entirely defined by a reduction procedure. In other words, we shall use the following identity: