Matrix element corrections in the Pythia8 parton shower in the context of matched simulations at next-to-leading order

We discuss the role of matrix element corrections (MEC) to parton showers in the context of MC@NLO-type matchings for processes that feature unstable resonances, where MEC are liable to result in double-counting issues, and are thus generally not employed. By working with Pythia8, we show that disabling all MEC is actually unnecessary in computations based on the narrow-width approximation, and we propose alternative MEC settings which, while still avoiding double counting, allow one to include hard-recoil effects in the simulations of resonance decays. We illustrate our findings by considering tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t\bar{t}$$\end{document} production at the LHC, and by comparing MadGraph5-aMC@NLO predictions with those of POWHEG-BOX and standalone Pythia8.


Introduction
The exclusive simulation of processes that feature heavy unstable particles, such as the top, W , Z, and Higgs in the Standard Model, is complicated for a variety of reasons.Conceptually, the most straightforward approach is to choose a (set of) decay channel(s) for the unstable particle(s) and compute the process where the initial and final states only include light partons and the products of such decay channels.For example, in the case of t t production, one may focus on a dilepton channel, and thus simulate (in hadronic collisions) 1pp −→ b li ν i bℓ j νj , rather than pp −→ t t . ( From the viewpoint of computing resources, the process of eq. ( 1) is much more expensive than that of eq. ( 2).This does not change when one takes into account the fact that eq. ( 2) must be supplemented by the simulation of the leptonic decays of the tops, namely: since the complexity of such a simulation is generally smaller than that of either eq.( 1) or (2) 2 .
Having said that, one remarks that the combination of eqs.( 2) and ( 3) is equivalent to eq. ( 1) only in the limit of a vanishing top-quark width (the so-called narrow-width approximation); away from that limit, the agreement between the results of the two approaches might be degraded by the presence of non-resonant (i.e., non-t t mediated) contributions.Even when only resonant contributions are present, it is generally the case that the incoherent simulation of the production process (eq.( 2)) and of the decay processes (eq.( 3)) constitutes a much poorer physics description than that emerging from eq. ( 1), and this because of two mechanisms.Firstly, the kinematics of a decay product of a given unstable particle may be affected by that of partons/leptons in the event which are themselves not decay products of that unstable particle.These effects are called (production) spin correlations, and cannot be accounted for by the separate simulations of the individual decays, while they are correctly included by the matrix elements that underpin eq.(1).Secondly, the kinematics of the decay products is affected by the emissions of the extra partons that appear in the hard process when perturbative corrections are considered, as in the case of NLO+PS simulations.By construction, a parton shower mimics the effects of such extra emissions, however with an accuracy which decreases with the hardness of the recoils they induce.Therefore, while perturbative corrections to eq. ( 1) will automatically include hard-recoil effects, such effects are poorly described when the simulation of the decays of eq. ( 3) is followed by an ordinary parton shower.
As far as the first mechanism is concerned, spin correlations are nowadays routinely taken into account in the context of production+decay simulations (eqs.( 2) and ( 3)) by means of the procedures of ref. [2] (especially in NLO+PS approaches) or ref. [3]; we shall not discuss them any further in this work.Conversely, hard recoils, absent computations of the complete processes such as that of eq. ( 1), are approximated by Monte Carlo event generators (MC henceforth) by means of the so-called Matrix Element Corrections (MEC henceforth).The manner in which MEC are simulated depends on the specific MC one employs, but they all rely on (exact or approximated) tree-level matrix elements that describe the emission(s) of interest.Before the widespread adoption of matching and merging techniques, MCs (and specifically Pythia) used to add MEC to both the production (eq.( 2)) and the decay (eq.( 3)) processes.These mechanisms will henceforth be referred to as production and decay MEC, respectively; sample Feynman graphs involved in their computations, in the specific case of top-quark production and decay, are depicted in fig. 1.It must be stressed that MEC graphs for production and decay are not allowed to interfere; in other words, they are employed incoherently.

H H
In the context of an NLO+PS simulation 3 , MEC can potentially spoil the accuracy of the computation by double counting (essentially, by generating again an emission already accounted for at the hard matrix element level -this is the case for the contribution shown on the left-hand panel of fig.1).In fact, this is not an issue for decay MEC, even in the presence of spin-correlation corrections, since the latter do not feature (at present) any extra emissions.There is therefore a good motivation to always include decay MEC in one's simulation, since they may induce, and accurately account for, large shifts in observables sensitive to emissions off decay products 4 .In contrast, the case of production MEC must be considered more carefully.In MC@NLO-type [5] simulations, the MC counterterms that guarantee the absence of double counting are determined by considering the showers without any additional MEC; this implies that production MEC are certain to spoil the accuracy of an MC@NLO simulation.This is not a simple formal problem; since emissions in MC@NLO are not ordered in hardness, it is not unlikely that the effects of such double counting can be visible in phenomenologicallyrelevant quantities.On the other hand, in a POWHEG-type simulation the hardest emission is always generated before the MC shower is added; potential NLO effects driven by production MEC are thus restricted to the region where the POWHEG-generated emission did not occur, i.e. a small-p T region where the cross section is Sudakov suppressed.
The bottom line is that both MC@NLO-and POWHEG-based simulations are compatible with, and phenomenologically benefit from, decay MEC.Conversely, while production-MEC effects are essentially negligible in POWHEG and can be either included or discarded, they are responsible for double counting in MC@NLO, and must therefore not be employed there.Thus, the procedure to adopt in MC@NLO and POWHEG is in principle clear.In practice, unfortunately, the MC@NLO-type MC interfaces to date have disabled all MEC instead of just those necessary to avoid double counting.We have already mentioned that the phenomenological implications of such choices should be restricted to a certain narrow class of observables.However, observed differences between MC@NLO-and POWHEG-based t t predictions at the LHC are commonly attributed to the different underlying matching mechanisms, whereas it might be the case, and certainly for the observables mentioned above, that they stem from the absence as opposed to the presence, respectively, of decay-MEC effects.
The aim of this note is to document how the separation of production and decay MEC can be achieved 5 in Pythia8 [7], thus paving the way for MC@NLO-based simulations which are phenomenologically more complete than those carried out thus far.We shall present examples for t t production, where Pythia8 standalone runs are compared with POWHEG-based simulations performed with POWHEG-BOX [8], and MC@NLO-based simulations performed with MadGraph5_aMC@NLO [1] (shortened as MG5_aMC henceforth).

Matrix Element Corrections in the Pythia8 dipole shower
An introduction to the basics of MEC has been given in sect.1; here, we focus on their implementations within Pythia8, and explicitly consider the various settings that must be chosen in order to preserve the accuracy and precision of the underlying perturbative predictions.The Pythia8 standalone parton shower algorithm is based on a leading-logarithm DGLAP evolution, which is by default corrected by means of MEC to reproduce the appropriate tree-level matrix element behaviour, including in the cases where Pythia8 is interfaced with external calculations.Technically, this is done by populating all of the phase space by candidate emissions whose number is overestimated, and which are then possibly vetoed with a rejection rate given by the ratio of the matrix elements over their shower approximations.While the MEC provide continuity between the exact matrix-element based calculation and the parton shower in the hard emission limit, it also provides the soft emission limit necessary for particles with mass on the order of the soft emission scale.Various options are available that allow for all, some, or none of the MEC to be included in a given simulation.It is therefore important to understand which of these options need to be employed in the context of NLO+PS simulations, so as to avoid the double counting issue discussed in sect. 1.
We start by pointing out that the separation between MEC in production and decay, exemplified by the two graphs in fig. 1, is conceptually well defined in the zero-width limit (the width being that of the particle that decays), since in such a limit the interference between the two types of tree-level graphs vanishes; this is important, because MEC are effected at the level of amplitude squared.However, this is not the way in which MCs, and specifically Pythia8, operate, since the primary distinction there is that between emissions off initial-state and off final-state legs (conventionally referred to as ISR and FSR, respectively).Thus, production MEC are relevant to both ISR and FSR, whereas decay MEC are solely relevant to FSR.This implies that production MEC constitute a difficult problem, since in general at the level of amplitude squared one cannot unambiguously separate the effects due the graphs where the extra parton is emitted by an initial-state leg from those due to emissions by a final-state leg 6 .
The Pythia8 solution to the problem posed by production MEC is rather draconian.Namely, MEC are applied to ISR only for colour-singlet resonant production [9] (i.e. to qq → V , with V a heavy EW boson, and to g g/γγ → H), since in such a case the ambiguity mentioned before is absent.As far as FSR is concerned, one applies MEC by considering, for each possible emitting colour dipole, a matrix element deemed equivalent to it, according to table 1 of ref. [10] (for example, for t t production, the relevant matrix element is γ * → t t g).
According to what was said previously, the latter FSR mechanism applies to decay MEC as well.However, before resorting to that, priority is given to matrix elements which exactly match the decay one is considering 7 ; if those exist, they are used for the MEC.In the Standard Model, these include the t → W b g and the V → qq matrix elements; in particular, this implies that a top that decays hadronically is never seen as such, but rather as a decay chain of a top decaying to a W b pair, followed by a hadronic W decay.
Regardless of whether equivalent or exactly-matching matrix elements are employed for MEC, their effects can be iterated by Pythia8 for subsequent emissions as well (as opposed to limiting them to the first emission).This is done by applying MEC on the system obtained by discarding all radiated partons except that whose kinematical configuration one seeks to correct presently; this may require a momentum reshuffling (for example: if a t → W b g g configuration is obtained by means of two subsequent emissions off the b, such emissions can be both corrected, one after the other, by using the t → W b g matrix element; g → g g branchings are not corrected).For a decay chain (such as t → W (→ qq)b), the iteration of MEC is applied individually to each resonance that decays (here, the top and the W ).
The various strategies described above are controlled in Pythia8 by means of the following parameters, separately for ISR (X=SpaceShower) and FSR (X=TimeShower), which are to be set equal to either on (which is the default for all of them) or off: X::MEcorrections : If =on, use MEC whenever available, regardless of whether they entail the use of equivalent or exactly-matching matrix elements.
X::MEafterFirst : If =on, apply MEC also to all emissions after the first; ignored if X::MEcorrections=off.
Among other things, the above implies that if X::MEcorrections=on, then the MEC stemming from exactly-matching matrix elements are always applied, irrespective of the value of X::MEextended.Thus, presently one can generally switch production MEC off while keeping decay MEC on, but not the other way around.This is just as well, since as it was discussed in sect. 1 MEC in production have to be employed only when neither NLO matching nor multi-jet merged simulations are available for the process of interest (which nowadays is a virtual impossibility), whereas MEC in decay are important for better phenomenological predictions in the context of relatively cheap (CPU-wise) generations.Having said that, this discussion clarifies that a more flexible approach to MEC in Pythia8 with respect to that given by the three (per shower type) parameters above would be desirable, namely one that allows the MEC to be applied or disabled until a certain condition is met 8 .
In summary, the necessity of avoiding double counting in MC@NLO-type simulations led to the recommendation that all MEC be switched off, by means of: SpaceShower:MEcorrections = off TimeShower:MEcorrections = off An undesirable by-product of these settings is the absence of MEC for decays.However, we have seen that even with the public versions of Pythia8 (starting with Pythia8.219), this can be avoided by setting: SpaceShower:MEcorrections = off TimeShower:MEcorrections = on TimeShower:MEextended = off in the majority of the cases of interest including, but not limited to, the processes that feature a t t final state; both values of TimeShower::MEafterFirst are acceptable.
We stress again that this is not a clean separation between production and decay MEC in the context of a generic process, and therefore the above setting recommendations must not be seen as universal.As a cautionary tale, consider Z-mediated dijet production, qq → Z ⋆ → qq.This process is seen by Pythia8 as the production of a colour singlet followed by its decay, rather than as a unique 2 → 2 production process.As such, exactly-matching matrix elements exist for both production and decay (Z ↔ qq g), which would not be the case had the (2 → 2)process picture been adopted.It follows that the value of MEextended is irrelevant in this case, and MEC are applied or not depending solely on the value of MEcorrections.Thus, if the latter is equal to on, one has (does not have) double counting if the process is generated by MG5_aMC as a 2 → 2 (2 → 1 followed by tree-level Z decay).The bottom line is that, since avoiding any double counting is of paramount importance, in case of doubt the user is encouraged to contact the Pythia8 authors.

Monte Carlo event samples
We now proceed to evaluate the impact of MEC in resonance decays in realistic NLO+PSaccurate MC simulations of relevance for LHC phenomenology.For this, we consider different MC samples for the production of stable top-quark pairs in proton-proton collisions at s = 13 TeV.
The first set of samples is generated using the MG5_aMC [1] program at the leading-and the next-to-leading-order (the latter with MC@NLO matching) in the strong coupling constant.The renormalization and factorization scales are set equal to the sum of the transverse energies of the final-state particles at the hard-process level, namely the top, the antitop, and when present a light quark or a gluon.The top quarks are subsequently decayed either leptonically or semileptonically (by including only electrons and muons), and preserving the tree-level spin correlations, with MadSpin [11].A second sample is generated using standalone Pythia8.309(thus, at the LO), with scales set equal to the geometric mean of the squared transverse masses of the two outgoing particles.The last sample is generated at the NLO with the HVQ [12] program in POWHEG-BOX-V2 [8,13].In POWHEG-BOX the scales are set equal to m 2 t + p 2 T , with m t and p T the top quark mass and transverse momentum, respectively, evaluated by employing the underlying Born configuration; the h damp parameter is set equal to m t .The top decay is performed at the tree-level, and including the effect of spin correlations, using the approach of ref. [2] as is implemented in the POWHEG-BOX internal routines.In all samples the NNPDF31_nnlo_as_0118 [14] parton distributions are used, and the value of the top quark mass is set equal to m t = 172.5 GeV.Both the MG5_aMC LO sample and the Pythia8 standalone one are normalized, prior to any final-state cuts, to the NLO total cross section as is predicted by MG5_aMC.
All of the predictions are then interfaced to Pythia8.309 (employed with the Monash tune [15]) to include the effects of parton showering, multiple parton interactions and underlying event.For the standalone Pythia8 and for the POWHEG-BOX samples, MEC are included wherever available, for both production and decay.For each of the MG5_aMC LO and NLO event sets, we shower the events twice.Namely, we use the current recommended MEC settings, in which matrix-element corrections are completely switched off, as well as the newer settings presented in sect.2, in which MEC are included only in the decays.Finally, the showered particle-level events are passed through Rivet [16]-based analyses that implement the observables of interest.

Phenomenological results
The inclusion of MEC to resonance decays is in general expected to affect the kinematics of the reconstructed decay products.We shall show in this section a few illustrative observables that document the extent of these effects.In particular, the Pythia8 LO+PS and POWHEG-BOX +Pythia8 NLO+PS predictions that include MEC to production and decays are compared with MG5_aMC NLO+PS predictions where MEC are either switched off (the current default), or included only for the decay (the settings proposed here).Since the focus of this work is the impact of MEC, no effort is made to simulate hard radiation in production in a manner which is as mutually consistent as possible across the various programs we employ.We note that the corresponding phase-space region is in any case better described by merged approaches, which also feature smaller systematics w.r.t.those of matched simulations.In order to facilitate the disentangling of the effects due to decay MEC from those due to hard radiation in the context of simulations stemming from exactly the same assumptions, in appendix B we present comparisons between the LO-and NLO-accurate predictions of MG5_aMC, for the same observables as those considered in this section.
All of the figures of this section have the same layout, namely consist of a main frame and a lower inset.In the main frame the four predictions (and possibly the experimental data) are displayed in absolute value, as blue (MG5_aMC with decay MEC), green (MG5_aMC without MEC), orange (POWHEG-BOX +Pythia8), and red (Pythia8) histograms.The lower inset shows the ratios of the latter three predictions (and possibly of the data) over the one for MG5_aMC with decay MEC, using the same colour patterns as in the main frame.For ease of reading the plots, the labels indicate, where relevant, which kind of MEC are applied: MEC in decay (MEC dec ) for MG5_aMC, and MEC in both production and decay (MEC dec prod ) for POWHEG-BOX and Pythia8 standalone.
For the case of the semileptonic channel, we can expect differences in the reconstructed hadronic W -boson and top quark masses.These two distributions are shown in fig. 2. The MEC in decay lead to a difference of about 10% (20%) in the reconstructed W -boson (top-quark) mass shape in the vicinity of the resonance peak.After inclusion of MEC, the three generators considered are in a better than 5% agreement with each other around the resonances peak, where the validity of the narrow-width approximation is good, and the impact of production hard radiation is negligible.As one can see from fig. 6, hard radiation does have a visible (but still moderate) impact at large invariant masses, more pronounced in the case of the reconstructed top mass than for the W mass, which explains the small residual differences among our three benchmark results in that region.t − m 2 W ∼ 150 GeV, sensitive to the value of the top mass.We observe that decay MEC shift the position of the peak, and after their inclusion in MG5_aMC this simulation and the POWHEG-BOX one (and to a good extent that of Pythia8 as well) agree fairly well with each other around and below this peak.At larger m(l, b jet ) values the detailed description of the production mechanism, which is treated differently in the three codes, becomes important, and we observe relative differences of up to 30% in this region.Inspection of the left panel of fig.7 confirms that the vast majority of these discrepancies are indeed due to hard radiation, with a residual O(5%) effect stemming from decay MEC.The impact of decay MEC is also evident at small values of the b-jet transverse momentum (for p T smaller than about 50 GeV); thus, in this region decay MEC improve the agreement between the three generators.However, at variance with m(l, b jet ), for the present observable the separation between hard-radiation and decay-MEC effects is less clear-cut; this can be best understood by looking at the right panel of fig. 7. Having said that, it is at large p T that the impact of hard radiation is larger than that of decay MEC; hence, the inclusion of decay MEC in MG5_aMC does not help reduce significantly the O(10%) discrepancies between MG5_aMC, POWHEG-BOX, and Pythia8.
Since decay MEC modify the kinematic properties of the b-quark and of the B-hadron resulting from the hadronization of the former, they also have an impact on the radiation pattern inside its corresponding b-tagged jet.In order to illustrate this, in fig. 4 we consider the distribution of the b-jet profile, and of the scaled B-hadron energy spectrum.The b-jet profile r, a.k.a.b-jet shape, is defined as the average fraction of the jet transverse energy that lies inside an inner cone of radius r < R, with R the jet-radius parameter; jets are defined according to the anti-k T algorithm [17], with parameters that depend on the specific analysis one considers.The scaled B-hadron energy, defined as the ratio of the B-hadron energy over the energy of the b-jet containing it (with energies defined in the lab frame), is a proxy for the longitudinal b-quark fragmentation function.The decay MEC are found to narrow the distribution of energy around the b-jet axis (i.e.there are comparatively more events at small r values) , and to shift the peak of the scaled B-hadron energy spectrum to higher values.The effect of decay MEC is significant in the whole r range considered, while that of hard radiation is negligible (see the left panel of fig.8): thus, after their inclusion, the agreement among the MG5_aMC, POWHEG-BOX, and Pythia8 simulations is of O(2%).The radiation pattern is more complicated in the case of the scaled B-hadron energy (see the right panel of fig.8), with residual non-negligible hard-radiation effects at small values; still, in the bulk of the distribution the agreement among our three benchmark results is at the same level as for the jet energy profile.
We finally compare our predictions to experimental data of jet substructure in t t semileptonic events that use charged-particle tracks from the CMS collaboration [18].In fig. 5 we consider two observables that are sensitive to the radiation pattern inside b-tagged jets: the groomed subjet distance, ∆R g , and the jet width, λ 1  1 .As is clearly demonstrated by fig.9, these observables are essentially insensitive to hard radiation, and decay MEC dominate their behaviour.Given what has been observed so far, it is therefore not particularly surprising that the inclusion of decay MEC in MG5_aMC vastly improve its agreement with both POWHEG-BOX and Pythia8 results, whereas the previous MEC settings in MG5_aMC led to discrepancies of O(±15%) with these two codes.The description of the data is also significantly improved.

Conclusions
We have illustrated how matrix element corrections to resonance decays as implemented in the Pythia8 parton shower can be consistently included in MC@NLO-type NLO+PS simulations in order to improve their phenomenological accuracy, without resulting in any double counting.
We have discussed the impact of these corrections in a process of particular relevance for LHC physics, namely the production of top quark pairs.We have also verified that the same conclusions apply to the associated production of top-quark pairs and a Higgs boson -this must be expected, since this kind of effects are thought to largely factorize w.r.t. the hard process; we thus regard our findings as to be universally valid.We have found that decay MEC can have a relative impact on the shape of distributions of up to 20%.By comparing NLO+PS predictions stemming from MC@NLO-and POWHEG-type matching with standalone Pythia8 ones we find that a significant reduction in the spread among the results occurs if the MEC are included whenever possible in the various simulations.We thus encourage the usage of the MEC settings proposed in this paper for any practical applications, in order to improve the phenomenological accuracy of MC@NLO-type simulations, and to reduce systematic uncertainties.

C Control of Matrix Element Corrections
The Pythia8 setting TimeShower:MEextended=off disables the equivalent-matrix-element FSR MEC for all parton emissions in the production process.As discussed in the text, consistency with MC@NLO-type matching requires only that the MEC be disabled for the first emission.However, in the Pythia8 shower, the MEC also provide the soft emission limit nec-  CMS Data, L=35.9 fb 1 mg5_amc@nlo+py8 MEC dec mg5_amc@nlo+py8 mg5_amc@lo+py8 MEC dec mg5_amc@lo+py8 0.0 0.2 0.4 0.6 0.8 1.0  distribution as reconstructed from charged particle tracks in semileptonic t t events, to measured data by the CMS Collaboration [18].
essary for particles with mass larger than or of the order of the soft emission scale.These corrections lead to the dead-cone effect.
For the case of a heavy particle system radiating a gluon, such as t t → t t g, it can be argued that the MEC to heavy particle radiators beyond the first emission, e.g.t t g → t t g g, are unimportant -dipoles with the g as the radiator will dominate over the mass suppressed t dipoles.We have tested this explicitly by modifying the Pythia8 shower to disable the MEC for only a user-determined number of (QCD) emissions.We observe little or no impact of this choice on any of our numerical results.To simplify the discussion and allow users to experiment with the impact of the MEextended setting in a public version of Pythia8, we have not used this new capability in our results.This option will be publicly available in an upcoming Pythia8 release.

Figure 1 :
Figure 1: Sample graphs relevant to MEC in top production (left panel) and top decay (right panel).The top quark is depicted by means of a thicker fermion line.The hard production process H is indicated by the blue circle.

Figure 3 :
Figure 3: The distribution in the invariant mass of the (lepton,b-jet) pair (left panel) and the leading b-jet transverse momentum (right panel) in dileptonic t t events.

Figure 4 :
Figure 4: The differential b-jet shape distribution (left panel) and the scaled energy fraction of the B-hadron (right panel) in semileptonic t t events.

Figure 5 :
Figure5: The groomed subjet distance (left panel) and the jet width (right panel) distribution as reconstructed from charged particle tracks in semileptonic t t events, compared to measured data by the CMS Collaboration[18].

Figure 7 :
Figure 7: The distribution in the invariant mass of the (lepton,b-jet) pair (left panel) and the leading b-jet transverse momentum (right panel) in dileptonic t t events.

Figure 8 :
Figure 8: The differential b-jet shape distribution (left panel) and the scaled energy fraction of the B-hadron (right panel) in semileptonic t t events.

Figure 9 :
Figure9: The groomed subjet distance (left panel) and the jet width (right panel) distribution as reconstructed from charged particle tracks in semileptonic t t events, to measured data by the CMS Collaboration[18].