NLO matrix elements and truncated showers

In this publication, an algorithm is presented that combines the ME+PS approach to merge sequences of tree-level matrix elements into inclusive event samples with the POWHEG method, which combines exact next-to-leading order matrix element results with the parton shower. It was developed in parallel to the MENLOPS technique and has been implemented in the event generator Sherpa. The benefits of this approach are exemplified by some first predictions for a number of processes, namely the production of jets in e+ e- annihilation, in deep-inelastic ep scattering, in association with single W, Z or Higgs bosons, and with vector boson pairs at hadron colliders.


Introduction
In the past decade, the incorporation of higher-order corrections into parton-shower simulations has been in the centre of formal improvements of existing general-purpose Monte-Carlo event generators like HERWIG [5], PYTHIA [6], and SHERPA [4]. As parton showers approximate higher-order matrix elements only in the soft and collinear limits of the real-emission phase space, they do inherently not yield a good prediction in the hard domain. This deficiency can lead to rather large discrepancies between Monte-Carlo predictions and experimental data and should therefore be corrected. In this endeavour, two somewhat orthogonal approaches have been pursued. To improve the description of hard QCD radiation, "merging algorithms" (ME+PS), have been proposed [7], which were shown to be correct up to next-to-leading logarithmic accuracy in e + e − -annihilation into hadrons. These methods were improved and extended in a series of publications [8,9]. A reformulation which generically maintains the logarithmic accuracy of the parton shower was achieved in [1] 1 , providing also the means to classify the various other methods and implementations according to their formal accuracy. It makes use of "truncated" parton showers, which are a concept initially introduced in [2] in the context of the POWHEG method. An alternative merging technique, known as "MLM merging", was suggested in [12] and described in more detail and compared with the other algorithms in [13]. While the ME+PS methods succeed at improving the simulation of multijet events, they do not address the apparent problem that the total cross sections of the simulated inclusive samples are still of leading-order The individual terms in the sum are given by where |M B | 2 ({ a}) denotes the partonic matrix element squared, dΦ B ({ a}) is the corresponding differential n-particle partonic phase-space element, S({ f }) is the symmetry factor, F ({ a}) is the flux factor, and L is the parton luminosity. In the case of leptonic initial states, and ignoring QED initial-state radiation, the parton distribution functions f (x, µ 2 ) are replaced by δ(1 − x).
In a similar fashion, the real-emission part of the QCD next-to-leading order cross section can be written as a sum, depending on parton configurations {a 1 , . . . , a n+1 }, by replacing the Born-level matrix elements B with the real-emission matrix elements R and the Born-level phase space dΦ B with the real-emission phase-space dΦ R . Furthermore, it is useful to introduce a notation for mappings from real-emission parton configurations to Born-level parton configurations and vice versa. They are given by (cf. [19]) 3) The map b ij,k (a) combines the partons a i and a j into a common "mother" parton a ı , in the presence of the spectator a k by defining a new flavour f ı and by redefining the particle momenta. The inverse map, r ı,k (a) determines the parton configuration of a real-emission subprocess from a Born parton configuration and a related branching process ı,k → ij, k. The radiative variables Φ R|B are thereby employed to turn the n-parton momentum configuration into an (n + 1)-parton momentum configuration. Following [19], the real-emission matrix elements, R({ a}), can be decomposed into a number of terms R ij,k as , (2.4) where D ij,k are the dipole terms defined ibidem. Therefore, the real-emission differential cross section can be rewritten as a sum of trivially factorised contributions . (2.6) and where R ij,k = L R ij,k , cf. [19]. The no-branching probability of a parton-shower algorithm is derived, based on Eq. (2.6) and the decomposition of the dipole terms D ij,k in the parton-shower approximation. It reads where the dipole-dependent no-branching probabilities are given by 2  is a constrained no-branching probability, as x 1/2 ≤ 1 before and after the branching process.
The quantities K ij,k are the parton-shower evolution kernels, which depend on the parton flavours f i , f j and f k and on the radiative phase space Φ ij,k R|B = {t, z, φ}, cf. [19]. The variable t, with t ∝ 2 p i p j , is the evolution parameter while z is called the splitting variable of the parton-shower model and J ij,k is the Jacobian factor associated with the transformation of integration variables.

The POWHEG approach
Using a simple corrective weight, it is possible to modify the parton-shower such that it produces the O(α s ) radiation pattern of the matrix element [19]. The corresponding dipole-dependent no-branching probability reads where the actual matrix elements R and B, which also include the respective parton luminosity factors, replace the parton shower kernel K and the parton luminosities of the equation above, Eq. (2.8).
The key point of the POWHEG method is, to supplement Monte-Carlo event samples from such matrix-element corrected parton showers with an approximate next-to-leading order weight to arrive at NLO accuracy. This is achieved by multiplying Eq. (2.11) Note that the second term in the square bracket describes resolved emissions, simulated by the matrix-element corrected parton shower, while the first term incorporates unresolved emissions and virtual corrections.
To reveal the fixed order properties of (2.11) it is useful to inspect its expansion keeping terms up to O(α s ). (2.12) It is imperative to note that neither the presence nor the precise form of ∆ (ME) in the resolved emission term influence the fixed-order accuracy of the method at O(α s ). Similarily, different choices of scales at which α s is evaluated in R, theB-function, and the matrix-element corrected parton shower emission terms contribute at O(α 2 s ) because α s (µ 1 ) = α s (µ 2 ) (1 + O(α s )) [2]. In order to obtain the correct leading logarithmic behaviour of the real-emission term, α s and the parton luminosities in ∆ (ME) and the resolved emission term of Eq. (2.11) must be evaluated at scale k 2 T . NLL accuracy can be restored for processes with no more than three coloured partons by means of the replacement [21] α s → α s 1 + α s 2π where the MS expression of α s should be used. In our Monte-Carlo implementation of the POWHEG method, we follow this approach.

The ME+PS approach
The ME+PS approach for the inclusion of matrix-element corrections into the parton-shower relies on a twofold generation of radiative corrections: Through an ordinary parton shower on the one hand and through real-radiation tree-level matrix elements on the other hand. In contrast to the POWHEG method, which only corrects the first emission off the core interaction, the ME+PS technique can be employed for arbitrary higher-order tree-level configurations.
A first algorithm to achieve this was presented in [7]. The solution there is based on separating the radiative phase space into a region of soft/collinear emissions, the parton-shower (PS) region, and a region of hard emissions, the matrix-element (ME) region. By demanding each region to be filled by the respective way of generating radiation and some reweighting double counting and other problems can be avoided. The original formulation has been tremendously improved in [1], by realising that, formally, separate splitting kernels in the ME and PS regions can be defined, which add up to the full splitting kernel: The functional form of the separation criterion Q ij,k is in principle arbitrary as long as it identifies soft and collinear divergences in the real-radiation matrix elements. The approach, fully outlined in [1], then replaces the splitting kernels in the ME region by the ratio of the real-emission and Born matrix elements, just like this is done in a matrix-element corrected parton shower. However, in contrast to a reweighting technique, only emission terms are modified and no correction is applied to the no-emission probabilities. The ME+PS technique can be implemented in a master formula for the first emission, describing the expectation value of an arbitrary infrared safe observable O, similar to the POWHEG case: (2. 16) There are three components to the differential cross section: The term describing unresolved emissions, which is generated in the standard parton-shower approach, and the resolved part, which is now split between the PS and the ME domain. Within the ME domain, the matrix-element generator is directly invoked to define the real-emission configuration. This is possible due to the restricted phase space, removing all infrared divergent regions by applying the cut in Q ij,k and rendering the matrix element finite. In this case, the Sudakov form factor ∆ (PS) , which makes the matrix element exclusive, must be added explicitly. It can either be calculated analytically, like in the original formulation of [7], or by utilising the shower itself to generate the correct probabilities. This latter option is commonly referred to as the pseudoshower approach [8,1]. A complication arises if the phase-space separation criterion Q ij,k is different from the parton-shower evolution variable t. This can imply the possibility of a shower emission Q < Q cut being allowed "between" two branchings at Q > Q cut in the parton-shower history of the matrix element. In such cases, in order not to spoil the logarithmic accuracy of the parton shower, the existing branchings need to be embedded into the subsequent parton-shower evolution. This leads to a truncated shower algorithm [2,1].
where (cf. [19]) If w({ a}) equals one, i.e. if the parton shower approximation is equal to the real-radiation matrix element, the t-integral can be performed easily and the square bracket in Eq. (2.17) equals one. The total cross section is therefore identical to the leading order result, as a direct consequence of the unitarity constraint for the parton shower. The more common configuration will however be, that w({ a}) = 1. In this case, the square bracket in Eq. (2.17) is different from one and the total cross section is not equal to the leading order result. The origin of the mismatch is a difference in the emission rate of the parton shower compared to the ratio R/B of matrix elements. While the former is exponentiated into the Sudakov form factor, the latter appears in the differential real-radiation probability only. We refer to this effect as an "emission-rate difference" in the following.
To investigate its consequences on arbitrary observables, Eq. (2.16) can be expanded in powers of α s , resulting in While the first and second terms represent the usual parton-shower prediction, the third term encodes the emission-rate differences effected by matrix element corrections in the region of well separated partons. In the ME+PS approach they enter at O(α s ). Emission-rate differences seem to be an undesirable side-effect of the ME+PS method at first. However, given Eq. (2.19), they can serve as an indicator for the relevance of higher-order real-emission corrections. We will elaborate on this fact in some more detail in Sec. 4.

Merging POWHEG and ME+PS -The MENLOPS approach
In this section, the two master equations for the POWHEG (Eq. (2.11)) and ME+PS (Eq. (2.16)) approaches are combined into one single expression, defining the MENLOPS approach. The aim of this combination algorithm is to simultaneously have NLO accuracy in the cross section, leading logarithmic accuracy as implemented in the parton shower and hard higher-order emissions corrected using tree-level matrix elements. Our method of choice is to simply replace the unresolved and the PS resolved part in Eq. (2.16) with the respective POWHEG expression. This essentially amounts to the replacement of the parton-shower no-emission probability with the corresponding POWHEG result, ∆ (PS) → ∆ (ME) and a substitution of the leading-order weight B byB, like in the POWHEG method itself. The ME part of the cross section is then generated separately, starting from real-emission matrix elements, as described in Sec. 2.2. This immediately implies that it will not automatically benefit from a POWHEG implementation regarding the local K-factorB/B, and it is therefore necessary to supply this K-factor explicitly. There is no a-priori definition of a Born-level parton configuration in this context, because the ME event is defined in terms of a real-emission configuration. One rather has to identify a branching history { a} → b ı,k ({ a}) such thatB/B can be computed depending on b ı,k ({ a}). The definition is achieved by clustering the real-emission configuration using an algorithm which is similar to a sequential recombination jet scheme and which determines the node to be clustered according to the related branching probability in the parton shower. For more details on this technique we refer the reader to [1].
Implementing these ideas, the master formula for the first emission in MENLOPS is obtained as (3.1) Note that the arguments of Q ij,k have been suppressed for ease of notation. The resolved-ME part of the expression in square brackets exhibits an additional factor compared to the POWHEG master formula. This makes the emission-rate difference in the MENLOPS method explicit. However, the expectation value of O is still determined correct to O(α s ), as can be seen by explicitly expanding Eq. (3.1) in powers of α s : Comparing this with Eq. (2.12), we find that the O(α s ) accuracy of the MENLOPS method is identical to what is obtained from POWHEG. The potential mismatch between exact higher-order tree-level matrix elements and their respective parton-shower approximation, leading to a difference between ∆ (ME) and ∆ (PS) , contributes terms of O(α 2 s ) or higher as long as Θ (Q ij,k − Q cut ) enforces t > t 0 . The precise value of Q cut must therefore be chosen such that this constraint is satisfied. Apart from being of order α s , the term in the square brackets of Eq. (3.2) should be rather small in practice, as potential differences between matrix-element and parton-shower expressions merely lie in subleading logarithmic and power corrections. This corresponds to saying that the mismatch between POWHEG and MENLOPS results is at most of order α 2 s log µ 2 /Q 2 cut if the parton shower has LL accuracy, and of order α 2 s if it has NLL accuracy. Generating a second emission using the ME+PS method supplemented with the above metioned local Kfactor of course introduces additional emission-rate differences, as described by Eq. (2.19). However, because such terms are of O(α 2 s ), they do not spoil the next-to-leading order accuracy of the method.
At this point we would like to stress that in their publication Hamilton and Nason arrived at the same ideas [3].

Results
This section collects results obtained with an implementation of the algorithm described in the previous sections in the SHERPA framework. It aims at detailing the improved description of data collected in various collider experiments and at quantifying some of the systematic uncertainties inherent to the MENLOPS method, in particular those related to the merging of the multijet tree-level contributions. Note again, that the MENLOPS approach is designed to merge the next-to-leading order accurate description of a given core interaction (like for example e + e − → qq) through the POWHEG method with higher-order tree-level contributions (like e + e − → qqgg) described in the ME+PS approach. Since the total cross section is essentially defined by the POWHEG expression of the core process in question, uncertainties like those related to the choice of scales are encoded mostly there. They have been discussed in our parallel publication [19], while uncertainties related to the ME+PS method were discussed for example in [1,11]. However, a comparison with results of the ME+PS and POWHEG techniques alone is extremely useful to assess the quality of the approach and the improvements related to it. The precise setup of SHERPA for this comparison, including in particular a parton shower based on Catani-Seymour subtraction terms [22] and an automated implementation of the Catani-Seymour subtraction method [23] in the matrix-element generator AMEGIC++ [24] was described in detail in our parallel publication [19]. Throughout the paper, we use the CTEQ6.6 parton distribution functions [25], (and, correspondingly, the MS subtraction scheme) with α S (m Z ) = 0.118 and running at two-loop. If not stated otherwise, hadronisation is not accounted for. Multiple parton interactions are not included in the simulation. We use the Rivet program package [26] and the HZTool library [27] for analyses and comparison with data.

Merging Systematics
As pointed out in the previous section, the ME+PS approach violates the unitarity of the parton-shower simulation. This discrepancy is directly inherited by the MENLOPS method. The extent of this effect depends entirely on the quality of the parton-shower algorithm, as can be seen in Eq. (3.2): If the parton-shower approximation to the real-emission matrix element is good, the correction factor, Eq. (3.2), is close to one.
We test the quality of our algorithms in the reactions e + e − → hadrons, deep-inelastic lepton-nucleon scattering, Drell-Yan lepton-pair production, W -and Higgs-boson production and W + W − -production by varying the phase-space separation cut, Q cut , and the maximum number of partons, N max , which is simulated with matrix elements in the MENLOPS approach. This is in close correspondence to the sanity checks of the ME+PS method which have been presented in [1]. The respective results are summarised in Tabs. 1-6 and depicted in Fig. 1. It is interesting to note that differences in the total cross-section are smallest for the MENLOPS samples with N max = 1, and that they increase steadily for larger N max . This indicates that the parton shower tends to underestimate the cross section of higher-order tree-level contributions. We observe two important effects, which allow us to judge the quality of the MENLOPS approach with respect to NLO accuracy: Firstly, for N max = 1 the emission-rate differences never exceed the size of the NLO corrections. Secondly, for any given value of Q cut and N jet , the relative difference between the cross sections from MENLOPS and POWHEG is always smaller than the one between the cross sections from ME+PS and LO+PS. This is best seen in Fig. 1, and it gives some confidence that the MENLOPS technique can help to improve perturbative QCD predictions from parton-shower Monte Carlo. The above analysis can be seen from a different perspective as follows: Usually, the biggest intrinsic uncertainty of the ME+PS approach stems from the freedom to choose the phase-space separation cut, Q cut , as explained and exemplified in a number of processes in [1]. Since the MENLOPS method relies on identical ideas to separate the real-emission phase space, it naturally inherits this source of uncertainty. Deviations of MENLOPS results from results with different values of Q cut are to be expected. However, their small size in a reasonable range of Q cut is a sign of the algorithm working well. The following rule of thumb can be applied: If the value of Q cut is chosen too large, too much extra emission phase space is left to the POWHEG simulation, typically leading to an underestimation of jet rates, since POWHEG only simulates the first emission through matrix elements. If, on the other hand, this value is too small, too much phase space is filled by matrix elements with large final-state multiplicity, which may lead to noticeable emission-rate differences.        (1) 1.305(4) 1.369 (5) 1.391 (6) 1.469 (7) 1.405 (7) 1.470 (9)     The value of Q cut should therefore lie well between the parton-shower cutoff and the factorisation scale of the core process, with some margin on either side of this interval. We exemplify the stability of our MENLOPS implementation with respect to variations of Q cut in Fig. 2. Due to their similarity to Q ij,k , the differential jet rates shown there are extremely sensitive to the details of the radiation pattern and thus to the accuracy of the ME+PS implementation. They tend to expose even the slightest mismatch between PS and ME subsamples, which then shows up as a kink in the distribution. However, when varying Q cut in a rather wide range, we observe no sizable discrepancies between the respective MENLOPS predictions, which is a very encouraging result regarding the quality of the algorithm and its implementation in SHERPA.

e + e − → jets
In this section we focus on electron-positron annihilation into hadrons at LEP energies ( √ s =91.25 GeV).
The core process of the simulation is therefore the reaction e + e − → qq. A full wealth of experimental data has been provided by the LEP experiments, which allows to assess the quality of the MENLOPS approach in this simplest realistic scenario. Although the improvements discussed in this paper concern only the perturbative QCD part of the Monte-Carlo simulation, our results account for hadronisation effects using the Lund model [30,6] to make them comparable to experimental data. Virtual matrix elements needed for the simulation were supplied by code provided by the BlackHat collaboration [31]. Figure 3 highlights the improvement in the description of jet data. In the hard-emission region the MENLOPS results for the 2 → 3-, the 3 → 4-and the 4 → 5-jet rate are generally closer to the data than the POWHEG ones, which hints at the success of the simulation. Deviations in the 5 → 6-jet rate are most likely due to the fact that matrix elements for six-jet production are not included. Note that these distributions are normalised to the total cross section, such that no rate difference between the ME+PS and the MENLOPS samples can be observed. Figures 4 and 5 show examples of event-shape variables, which are all very well described in the hard-emission region by the MENLOPS simulation. Several distributions for jet angular correlations in 4-jet production, that have been important for the analysis of QCD and searches for physics beyond the Standard Model are investigated in Fig. 6. The good fit to those data proves that correlations amongst the final-state partons are correctly implemented by the higher-order matrix elements.

Deep-inelastic lepton-nucleon scattering
Deep-inelastic scattering (DIS) is one of the best understood processes in perturbative QCD. However, it has been an obstacle for a very long time to properly simulate hadronic final states in DIS using general-purpose Monte Carlo based on collinear factorisation. Only recently, a consistent approach was presented [11], that allows to describe jet data throughout the experimentally accessible range of Q 2 , the negative virtuality of the exchanged virtual γ * /Z-boson. It is absolutely mandatory for this method that a large number of final-state partons can be described by hard matrix elements in order to lift the severe restrictions on the real-emission phase space of the parton shower, which are imposed by the factorisation theorem.
In all our simulations we use the core process e + q → e + q. We present results for two analyses. The first is the measurement of inclusive jet production in [32], which covers different ranges of jet-pseudorapidity in the laboratory frame, η lab , in the low-Q 2 domain 5 < Q 2 < 100 GeV 2 . Jets are defined using the inclusive k T -algorithm [34] and are constrained to E T,B > 5 GeV and the pseudorapidity range −1 < η lab < 2.8, where E T,B is the jet transverse energy in the Breit frame. The second analysis corresponds to the measurement of dijet production in [33], which covered a wider range of Q 2 and produced many doubly differential jet spectra. The acceptance region is 5 < Q 2 < 15000 GeV 2 and −1 < η lab < 2.5. Jet transverse energies are subject to the cuts E T,B 1,2 > 5 GeV and E T,B 1 + E T,B 2 > 17 GeV. The latter requirement is introduced to avoid E T,B 1 ≈ E T,B 2 , which is the region of the phase space where next-to-leading order corrections are unstable due to implicit restrictions on soft emissions [35]. As outlined in [11], a crucial observable is given by the inclusive jet cross section, differential with respect to E 2 T,B /Q 2 . For E 2 T,B /Q 2 > 1 it probes a part of the phase space where leading order Monte-Carlo models without the inclusion of low-x effects are bound to fail in their description of jet spectra. Another very good observable to validate the proper Monte-Carlo simulation is the dijet cross section as a function of Q 2 . While still a relatively inclusive quantity, it is an important indicator for the correct simultaneous implementation  Figure 4: Total jet broadening and jet mass difference at LEP compared to data taken by the ALEPH experiment [28]. T MC/data Figure 5: C parameter and thrust distribution at LEP compared to data taken by the ALEPH experiment [28].  Figure 6: Angles between the leading (in energy) four jets defined using the Durham algorithm with y cut = 0.008. Results at the parton level are compared to data from the OPAL experiment [29].  Figure 7: Left: The inclusive jet cross section as a function of E 2 T,B /Q 2 in bins of η lab , compared to data from the H1 collaboration [32]. E 2 T,B is the jet transverse energy in the Breit frame, while η lab denotes the jet rapidity in the laboratory frame. Right: The dijet cross section as a function of Q 2 in bins of E T,1 + E T,2 , compared to data from the H1 collaboration [33].  Figure 9: The dijet cross section as a function of the dijet mass m jj , compared to data from the H1 collaboration [33].
of inclusive DIS and the additional production of hard QCD radiation. The high quality of the MENLOPS prediction for the two above observables is confirmed in Fig. 7. Discrepancies in the description of the E 2 T,B /Q 2 -spectrum in the forward region can be attributed to the fact that the simulation is limited to three additional partons in the hard matrix elements. This restriction is imposed by the usage of the matrixelement generator AMEGIC++ [24]. Figures 8 and 9 exemplify again that the MENLOPS simulation correctly predicts multijet differential distributions in all regions of the phase space, while the POWHEG approach fails in the low-Q 2 domain.

Drell-Yan lepton-pair production
Results for lepton-pair production through the Drell-Yan process are compared to data from the Tevatron at √ s = 1.96 TeV in Figs. 10-13, using the core process qq → , where = e, µ. The invariant mass of the lepton pair was restricted to be within 66 < m /GeV < 116 in the simulation. The MENLOPS and ME+PS samples use tree-level matrix elements up to Z + 3 jets with a merging cut of Q cut = 20 GeV. Virtual matrix elements are provided by BlackHat [31]. The Z → + − decay is corrected for QED next-to-leading order and soft-resummation effects using the Yennie-Frautschi-Suura (YFS) approach [36].
The Tevatron experiments provide a wealth of measurements sensitive to QCD corrections in Drell-Yan production. Fig. 10 shows the transverse momentum distribution of the lepton pair in two different analyses from the DØ experiment. The left hand plot displays a very recent analysis using the Z → µµ channel [37] to measure the Z-p ⊥ distribution normalised to the inclusive cross section. It requires muons with p ⊥ > 15 GeV in a mass window of 65 < m µµ /GeV < 115 and with |η| < 1.7. The muon signal is corrected to the particle level including photons clustered in a cone of radius R = 0.2 around each lepton. The plot on the right hand side stems from an analysis in the electron channel [38] which uses Monte-Carlo models to correct the leptons for all acceptances including the pseudorapidity range and minimal transverse momentum. Here we display the peak region of the transverse momentum of forward Z bosons with |y Z | > 2. The agreement between all three approaches and the measurement is outstanding. In the bins at p ⊥ < 10 GeV non-perturbative effects like the intrinsic transverse momentum of partons in a proton might play a role. Related Monte-Carlo models in SHERPA could be tuned to reach an even better agreement. Still, the Monte-Carlo prediction lies within the experimental error band over the full range. Two more measurements from the DØ experiment are displayed in Fig. 11. The pseudorapidity of the Z boson [39] was measured in the electron channel requiring electrons with p ⊥ > 15 GeV in the mass window 71 < m ee /GeV < 111. Again, all three Monte-Carlo approaches agree very well with the experimental data. The right hand plot shows the azimuthal correlation between the Z boson and the leading jet [40]. This is a measurement in the muon channel with the same selection cuts as described above. The distribution has been normalised using the inclusive Z cross section and the comparison shows that the three approaches underestimate the total rate for Z+jet production with respect to inclusive Z production by approximately 10%. This might hint at the need for NLO accuracy also in the Z+jet process. It is remarkable though that the inclusion of higher-order tree-level matrix elements significantly improves the shape of the distribution with respect to the POWHEG simulation. The observables presented so far are mainly sensitive to the correct description of the leading jet. For that reason even the POWHEG approach is well capable of providing sufficient accuracy in their prediction. We now proceed to observables sensitive to higher-order corrections. Figure 12 (left) shows the inclusive jet multiplicity [41] for jets constructed using the DØ improved legacy cone algorithm [42] with a cone radius of R = 0.5 and p ⊥ > 20 GeV. Jets were required to lie in |η| < 2.5 and to be separated from the leptons by ∆R( , jet) > 0.4. While POWHEG agrees with the data for the N jet = 1 bin it fails to predict the rate of events with more than one jet. The MENLOPS and ME+PS predictions impressively demonstrate the effect of higher-order corrections provided by tree-level matrix elements up to the third jet. They agree with the measurement within the error bands for N jet = 2, 3 and as expected fail to predict the correct four-jet rate because no matrix-element corrections have been applied at that multiplicity. Transverse momentum spectra of the three leading jets accompanying the Z boson were measured by DØ in [43]. The distributions in Fig. 12 (right) and 13 are normalised to the inclusive cross section for Z production and the jets have been constructed using the same settings as in the multiplicity measurement. Both MENLOPS and ME+PS deliver a very good description of these spectra while POWHEG fails to describe the rate and shape for the second and third jet. MC/data Figure 10: The transverse momentum of the reconstructed Z boson in Drell-Yan lepton-pair production at the Tevatron at √ s = 1.96 TeV. Experimental data stem from the DØ experiment [37,38] and are described in the text. MC/data Figure 13: Transverse momentum of the second and third jet [43] in Z+jets events at the Tevatron at √ s = 1.96 TeV.

W +jets Production
In this section we focus on the production of W -bosons and their subsequent decay into an electron-neutrino pair at the Tevatron at √ s = 1.8 TeV. The core process of the Monte-Carlo simulation is therefore qq → ν.
The separation criterion is set to Q cut = 20 GeV and up to three extra jets are taken into account. The electron-neutrino pair is required to have an invariant mass of m eν > 10 GeV. The W → eν decay is corrected for QED next-to-leading order and soft-resummation effects using the YFS approach [36]. Virtual matrix elements are supplied by BlackHat [31]. The left panel of Fig. 14 displays the transverse momentum of the W -boson as compared to data taken by the DØ collaboration [44], while the right panel shows the exclusive jet multiplicity of k ⊥ -clustered jets (D = 0.7) with at least 20 GeV. Although the event sample generated using the POWHEG technique only provides the best match to the central value of the data, all three event samples are well within the experimental uncertainties. On the other hand, already in the rate of single-jet events deviations between the POWHEG sample and both the MENLOPS and ME+PS samples are visible, with the latter two agreeing very well. Similarly, the POWHEG sample underestimates the amount of radiation into the central detector region, as exemplified in Fig. 15. The right panel of this figure shows that, since the POWHEG approach is capable of modelling the second hardest emission using the soft-collinear approximation of the parton shower only, its description of the angular separation of the the first two hardest jets is missing prominent features originating in the wide angle region. These features are of course present in the approaches having fixed-order matrix elements at their disposal. Figure 16 shows the differential jet rates d 01 , d 12 , d 23 and d 34 using the above k ⊥ -algorithm. While the first three of them, for the matrix-element merged samples, are described by matrix element to matrix element transitions, only the softer part of d 01 is described by such a transition for the POWHEG sample. The harder part of the d 01 receives corrections by matrix elements of higher jet multiplicity which are clustered into a single hard jet first. Of course, these corrections are missing in the POWHEG sample. Furthermore,   Figure 16: Differential jet rates d n n+1 in W production at the Tevatron at √ s = 1.8 TeV.

Higgs boson production
This section presents predictions for Higgs boson production via gluon fusion at nominal LHC energies of √ s = 14 TeV. As NLO corrections to the core process gg → h → τ + τ − are rather large, tremendous efforts have been made to perform fully differential calculations at NNLO [45] and several predictions have been presented which merged such fixed-order results with resummation at next-to-next-to-leading logarithmic accuracy [46]. In this publication, we have no means for an improvement of the resummed calculation, instead we are restricted by the limitations of the parton-shower model. However, the systematic inclusion of higher-order tree-level matrix elements through the MENLOPS method can yield a significant improvement of existing NLO predictions, thus partially closing the gap between full NNLO predictions and Monte-Carlo results. It was shown, for example, in [47] that the predictions from ME+PS algorithms are often competitive to NNLO results if only the shape, not the normalisation, of observable distributions is concerned. In our simulations we set m H = 120 GeV and we include the decay h → τ + τ − , however, the analysis focuses on the properties of QCD radiation associated with production of the Higgs boson. The invariant τ -pair mass is restricted to 115 < m τ τ /GeV < 125 at the matrix-element level. Virtual matrix elements are implemented according to [48]. The decay h → τ + τ − is corrected for QED soft-resummation and approximate next-to-leading order effects using the YFS approach [36]. Figure 17 shows the transverse momentum spectrum of the reconstructed Higgs boson. We observe that the POWHEG and MENLOPS samples are very consistent in the prediction of this rather inclusive observable. On the other hand, differences are observed in the results for individual jet transverse momentum spectra, cf. Fig. 18. They increase with jet multiplicity and with increasing transverse momentum, as can be expected, since the higher multiplicity jets are described by the uncorrected parton shower in the POWHEG method. Deviations are also found in the prediction of the dijet separation in η − φ space, which is shown in Fig. 19. However, it was previously found that the ME+PS result yields a prediction which is very similar to the NNLO result [47]. This feature is naturally retained in the MENLOPS simulation.

W + W − +jets Production
In this section we present predictions for the production of the W + [→ e + ν e ] W − [→ µ −ν µ ] final state at nominal LHC energies of √ s = 14 TeV. The lepton-neutrino pairs are required to have an invariant mass of m ν > 10 GeV each. The W → ν decays are corrected for QED next-to-leading order and soft-resummation effects using the YFS approach [36]. Virtual matrix elements are supplied by MCFM [49,50]. Again, this study focuses mainly on the properties of QCD radiation accompanying the diboson production process. Up to three additional jets at Q cut = 20 GeV are simulated in both, the MENLOPS and the ME+PS sample. It is known that high-multiplicity matrix elements in the ME+PS approach yield sizable effects on total event rates and shapes in this reaction [51], a feature which is inherited by the MENLOPS method. Setting the phase-space separation criterion to a rather low value compared to the average partonic centre-of-mass energy will thus always lead to sizable emission-rate differences, which might be an indication of potentially large higher-order corrections. A similar effect was observed in a recent analysis of Z-boson pair production in association with a hard jet [52]. While the NLO corrections to this process are comparably small at Tevatron, they can be rather large at nominal LHC energies. Restricting the available final-state phase space by a jet veto, the corrections were again limited to smaller values, which makes the importance of the ZZ+2 jets final state explicit. As we include up to three additional jets in our simulation of W + W − production, we observe similar effects. Figure 20 displays the invariant mass of the lepton pair and the scalar sum of transverse momenta of the jets, leptons and the missing transverse energy, H T . While the former is described very well by the next-to-leading order calculation used in the POWHEG sample and receives only mild corrections from higher-order matrix elements, H T receives sizable corrections at rather low values already. The reason for this is easily found in the sensitivity of H T to any jet activity and thus to higher-order matrix element corrections of the parton shower. This can be seen in comparison to Fig. 23, where a veto on additional jet activity was applied. We exemplify in Fig. 21 that the ME+PS part of the MENLOPS simulation predicts significantly harder radiation than the POWHEG subsample. The corresponding corrections naturally amplify the deviations between the respective predictions of H T . We show the impact of a jet veto on this distribution in the right panel of Fig. 24. Figure 22 presents predictions for the azimuthal separation of the leptons and the two hardest jets. Again, the former receives only comparably small corrections, while higher-order matrix-element corrections have  Transverse momentum of leading jet dσ/dp ⊥ (jet

1)
[pb/GeV] Transverse momentum of leading jet dσ/dp ⊥ (jet 1) [pb/GeV]  Figure 24: Azimuthal separation of the electron and the muon (left) and transverse momentum of the hardest jet in W + W − production at nominal LHC energies (14 TeV) after vetoing events with more than one jet with p T > 20 GeV.

Conclusions
In this publication, a parallel development and independent implementation of the MENLOPS algorithm, first discussed in [3], has been presented. This new algorithm combines the so far most advanced methods to include higher-order corrections to a given core process: The POWHEG technique, which allows to produce inclusive samples for that process with next-to-leading order accuracy, and the ME+PS technique, which allows to generate inclusive samples with a leading-order cross section, but with the production of additional hard radiation corrected by higher-order tree-level matrix elements. Until the work of Hamilton and Nason [3] and the work presented here, these two approaches were considered orthogonal and thus used independent from each other, in the regime of their respective strengths and validity.
With the recent efforts on combining them, the shortcomings of each method, i.e. the description of higher jet multiplicities in POWHEG and lack of the correct NLO cross section in ME+PS, have been expunged. We fully confirm the findings of Hamilton and Nason concerning both the formalism and the relative improvement in the simulation obtained through it. This is even more emphasised here, since the implementation of the MENLOPS method as presented in [3] seems to suffer from the choice of tools. As already indicated in the introduction, the omission of truncated showering in the program used to simulate the ME region may have caused a few of the uncertainties. We are convinced that in total, the superior quality of the ME+PS part of the simulation in SHERPA, including the truncated showering, are the only reason behind the improved simulation here -the formalism is identical in both publications. The drastically reduced uncertainties stress the great improvement by the MENLOPS method.
Our results and the ones presented by Hamilton and Nason in fact show a significant improvement of many aspects of previous simulations in a variety of processes, including here e + e − annihilation to hadrons, hadronic final states in DIS, jets in association with single vector bosons and with vector boson pairs, and the production of Higgs bosons through gluon fusion.
In the future, the description of many more processes with this combined NLO matching and multijet merging will become feasible. This is possible, because both the POWHEG and the ME+PS part of the implementation are fully automated in SHERPA.
We would like to also point out that the methods developed so far will naturally serve as a starting point to promote the ME+PS idea to full NLO, in the sense that merging sequences of multijet matrix elements at NLO into one inclusive sample becomes feasible. A first attempt to achieve this from a somewhat different angle has been presented in [53].