{'fraction_non_alphanumeric': 0.0279299113884896, 'fraction_numerical': 0.007787623501155407, 'mean_word_length': 4.653739327229709, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
QCD radiation in the production of heavy colored particles at the LHC

We present a study of effects of QCD radiation in the production of heavy colored states, employing inclusive multi-jet samples obtained by matching matrix elements and parton showers. We discuss several examples showing that matched samples are in general not only more accurate than a parton shower alone, but also often indispensable to make reliable predictions of beyond the Standard Model signals.


Introduction
Many models for physics beyond the Standard Model, notably Supersymmetry, Little Higgs models, various models for extra dimensions and several versions of technicolor models, contain new strongly interacting particles with masses below or near the TeV scale. These particles, if they exist, will be copiously produced at the LHC, and could provide the first and most important sign of new physics. By analyzing their production and decay modes, we might be able in only a few years to determine many of the new physics properties, including large parts of the particle spectrum, spin structure and the possible existence of new stable particles that can act as dark matter candidates.
The final state signatures characteristic of events where heavy colored particles are produced typically include a large energy measured in the detector, missing transverse energy and several hard jets coming from their decays. Such very spectacular events can be efficiently simulated by Monte Carlo generators that implement new physics models, including general-purpose generators like Pythia [1], Herwig [2] and Sherpa [3], and matrix element generators such as CalcHEP/CompHEP [4,5], MadGraph/MadEvent [6] and Whizard [7]. However, an additional difficulty in the simulation of particle production at hadron colliders is due to the presence of abundant QCD radiation, in particular initialstate radiation. This radiation, which is enhanced in the production of heavy and strongly interacting particles, can have important effects. It affects the event kinematics by giving a transverse boost to the heavy particle system, and can produce additional jets besides the jets originating from the decay of the heavy particles, thus complicating the reconstruction and indentification of the event.
This additional jet production has traditionally been simulated using Parton Shower (PS) Monte Carlo programs such as Pythia and Herwig, which describe parton radiation as successive parton emissions using Markov chain techniques based on Sudakov form factors. This description is formally correct only in the limit of soft and collinear emissions, but has been shown to give a good description of much data also relatively far away from this limit. However, for the production of hard and widely separated QCD radiation jets, this description breaks down due to the lack of subleading terms and interference. For that case, it is necessary to use the full tree-level amplitudes for the heavy particle production plus additional hard partons.
The Matrix Element (ME) description diverges as partons become soft or collinear, while the parton shower description breaks down when partons become hard and widely separated. In order to describe both these areas in phase space, the two approaches must be combined, without double counting or gaps between different parton multiplicities. An additional physical requirement is that such a procedure should give smooth distributions, and interpolate between the parton shower description in the soft and collinear limits and the matrix element description in the limit of hard and widely separated partons. Several algoritms have been proposed to achieve this, including the CKKW [8,9], Lönnblad [10] and Mangano [11,12] schemes. These different procedures are in substantial agreement and give consistent results at hadron colliders [13,14]. For Standard Model processes, in particular the production of jets and weak vector bosons, the matching schemes have been extensively used and compared to the available data from the Tevatron [15]. In addition, inclusive top quark pair production with matching has been compared to a next-to-leading order plus parton shower approach (MC@NLO) [16] in Ref. [12].
That jet matching is necessary for Standard Model backgrounds is obvious, since the only way many Standard Model processes (such as weak vector boson production) can emulate the production of new heavy states charged under QCD is by taking into account high-p ⊥ radiated jets. It is less clear, however, that jet matching should be important in new physics production, for several reasons: The new particles are expected to decay to high-p ⊥ jets making them fairly insensitive to the effects of QCD radiation; furthermore the parton shower formalism should work better the higher the mass is for the produced particles, since then the region in which emissions can be considered "collinear" is increased. As we will demonstrate in this paper, this is only partly true, and there are many scenarios where jet matching turns out to be surprisingly important also in the production of heavy new physics QCD states.
We will here use the production of squarks and gluinos within the framework of the Minimal Supersymmetric Standard Model, MSSM. The properties of initial state QCD radiation are however quite insensitive to the precise type of particle produced, and mainly depend on the mass scale and production mechanism of the new particles. All our results here are therefore readily applicable to any type of new physics producing colored heavy states decaying to jets.
The layout of the paper is as follows: In Sec. 2 we discuss the details of the matching schemes used, as well as how to solve the double counting problem that arises due to the decay of on-shell states in multiparton final states. In Sec. 3 we study the differences in QCD radiation between matched and unmatched simulations. The matching is found to reduce the uncertainties of the parton shower and hence to increase the predictivity of the simulations. This increased predictivity allows us, in Sec. 4, to study the effects of QCD radiation in the presence of jets from the decay of the new states. We analyse in detail a few examples where jet matching turns out to be crucial to get a qualitatively correct description of the new physics signatures. In Sec. 5 we take a look a the impact of using matched simulations in a more realistic experimental situation, where we include production of all relevant supersymmetric particles as well as Standard Model backgrounds. We draw our conclusion in Sec. 6.

Technical challenges
In this section we present the technical issues and the problems which can arise while merging matrix elements involving heavy colored states with a parton shower. Our approaches and their implementations are completely general and can therefore be applied equally well to the SM as to any Beyond the SM (BSM) construction. To be concrete we will focus on top pair production and on pair production of strongly interacting SUSY particles, such as squarks and gluinos. In this study we have employed the MadGraph/MadEvent matrix element generator [17,18,6] interfaced to Pythia 4 [1] for parton showering and hadronization. The MSSM implementation used was presented in Refs. [19,20,6].

Matching schemes
The goal of jet matching is to merge samples with different parton multiplicity obtained via matrix elements, correctly accounting for showering effects and avoiding double counting. The matching algorithms used in this study can be viewed as hybrids between the approaches currently employed by Sherpa [3] and Alpgen [21]. The phase space separation between the different multijet processes is achieved using the k ⊥ -measure [22]. No analytic Sudakov reweighting of the events is performed, but instead showered events are rejected if they are not matched to the matrix element-level partons. This method allows the use of the well-tuned showering and hadronization implementations of Pythia while retaining the advantages of matrix element production. The matching implementations in MadGraph/MadEvent can be used both for Standard Model and new physics processes.
The first matching scheme used in this study, the k ⊥ -jet MLM scheme, is the one used for MadGraph/MadEvent in Ref. [14]. The final-state partons in an event are clustered according to the k ⊥ -jet algorithm to find the "equivalent parton shower history" of the event. Here, the Feynman diagram information from MadGraph is used to allow only clusterings that correspond to diagrams existing in the generated matrix element. The smallest k ⊥ value is restricted to be above some cutoff scale Q ME cut . In order to closely mimic the behaviour of the parton shower, the k ⊥ value for each clustering vertex corresponding to a QCD emission is used as renormalization scale for α s in that vertex. As factorization scale, as well as renormalization scale for the central hard 2 → 1 or 2 → 2 process, the transverse mass m 2 ⊥ = p 2 ⊥ + m 2 of the particle(s) produced in the central process is used. This event is then passed to Pythia for parton showering. After showering, but before hadronization and decays, the final-state partons are clustered into jets, again using the k ⊥ jet algorithm, with a cutoff scale Q match > Q ME cut . These jets are then compared to the original partons from the matrix element event. A jet is considered to be matched to the closest parton if the jet measure k ⊥ (parton, jet) is smaller than the cutoff Q match . The event is rejected unless each jet is matched to a parton, except for the highest multiplicity sample, where extra jets are allowed below the k ⊥ scale of the softest matrix element parton in the event, Q ME softest . This matching scheme can be used with both the old (vituality-ordered) and the new (p ⊥ -ordered) shower implementations of Pythia.
In order to further study the systematics of the jet matching, a new matching scheme has been implemented, which we call the "shower k ⊥ scheme". In this scheme, which is here used for the first time, events are generated by MadGraph/MadEvent as described above, including the reweighting of α s . The event is then passed to Pythia and showered using the p ⊥ -ordered showers. For the p ⊥ -ordered showers, Pythia reports the scale of the first (hardest) emission in the shower, Q PS hardest . For events from lower-multiplicity samples, the event is rejected if Q PS hardest is above the matching scale Q match , while events from the highest multiplicity sample are rejected if Q PS hardest > Q ME softest , the scale of the softest matrix element parton in the event. This matching scheme is simpler and yet effectively mimics the workings of the k ⊥ -jet MLM scheme. However, it allows for the matching scale Q match to be set equal to the matrix element cutoff scale Q ME cut , and it more directly samples the Sudakov form factor used in the shower. Furthermore, the treatment of the highest multiplicity sample more closely mimics that used in the CKKW matching scheme.
Since we here study heavy QCD particle production at hadron colliders, we expect final-state radiation to be highly suppressed as compared to initial-state radiation. The treatment of final-state radiation from the (undecayed) produced particles is therefore of somewhat secondary importance, and does not strongly affect any results of this study. We have here chosen not to differentiate between shower particles originating from finalstate radiation and initial-state radiation in the jet clustering in the k ⊥ -jet MLM scheme, and correspondingly for the "shower k ⊥ " scheme, veto also events where the final-state radiation reports an emission scale above the matching scale.

Double counting from resonant diagrams
In the simulation of events where colored particles are produced that can decay into one another by emitting partons, a further problem arises, also related to double counting. Consider, for example, the contributions togq + 1 jet coming from various subprocesses. Among these, gg →gqq displays a peculiar behaviour. In Fig. 1 we show three representative diagrams out of a total of sixteen. Diagram (a) is a "genuine" correction to the 2 → 2 These processes are already taken into account in thegg andqq * channels, respectivetly, and have to be subtracted.
Born amplitude and it is correctly handled by the matching procedure. Diagrams (b) and (c), however, contain possibly resonant gluino and squark propagators. When integrated over the phase space these diagrams give rise either to gg →gg withg →qq or to gg →qq * withq * →gq, depending on the mass hierarchy. These contributions, however, are already taken into account by the Born level processes gg →gg or gg →qq * plus the corresponding decays and have therefore to be properly subtracted to avoid double counting. In fact, this kind of issue is not specific to SUSY, but appears every time a given final state can be reached through different decay cascades. 1 Several solutions have been proposed in the literature, with various degrees of approximations, only a few of which can be employed in a Monte Carlo framework. For example, subtractions at zero width, typically implemented in NLO calculations, are not suitable to a Monte Carlo approach. In Ref. [25] these diagrams were removed from the matrix element by hand. Even though this procedure in principle violates gauge invariance, it works very well in practice when the width of the resonances is sufficiently small, Γ/m ≪ 1. In order to deal with this kind of problem in full generality we have implemented and checked two independent types of solutions.
In the first, resonant diagrams are removed as in Ref. [25], but in an automatic way directly by MadGraph. This approach has the virtue of being very simple, though it is not gauge invariant and neglects the interference between the diagrams. The impact of these approximations therefore has to be carefully checked case by case.
The second solution is an hybrid between a proper gauge-invariant subtraction and the diagram removal. It is based on a general algorithm implemented in MadEvent to record the information of the presence of intermediate resonant propagators on an eventby-event basis. Such information is needed by the showering program to ensure that the QCD evolution does not shift or even cancel the Breit-Wigner peaks. When only resonant or only non-resonant diagrams are present, there is, in principle, no ambiguity and the information could be directly passed to the events. When both resonant and non-resonant diagrams are present such a unique assignment is not possible, due to interference effects. However, the interference terms are typically suppressed by Γ/m, and can therefore be neglected in the narrow width limit. In this case one can associate the presence of a propagator to an event on a statistical basis, using the relative size of the corresponding squared amplitudes. The above procedure exactly parallels that commonly employed to associate a color flow to an event [26]: The calculation of the matrix element is exact at all orders in the number of colors (as the amplitude-squared includes all diagrams and their interference), however the color flow assignement is approximated by ignoring interference terms which are proportional to 1/N 2 c . Such propagator/diagram mapping is possible in MadGraph/MadEvent due to the fact that automatic integration over space space is already based on the amplitudes corresponding to single diagrams. The phase space integration is split into independent channels, each corresponding to a single Feynman diagram [18]. Events generated in a resonant channel are provided with the resonant particle information and also one of the color flows consistent with the Feynman diagram itself (in so doing, we avoid possible mismatches between the propagator structure and the color flow, which would lead to an inconsistent shower evolution). Having at our disposal the above procedure, it easy to get rid of the resonant contributions, by simply dropping the events with resonant propagators before passing them to the Pythia interface. The results presented and discussed in the following are obtained using the latter method. However, we have checked that the two methods are consistent.

QCD radiation in ME/PS merged samples
QCD radiation in production of heavy states at hadron colliders is expected to be mainly due to initial state radiation. 2 The main reason for this is that the phase space for initial state radiation is typically large and that collinear radiation from the final state particles is suppressed by their large mass. The main impact of QCD radiation on the production of heavy states is therefore twofold: A transverse boost of the heavy particle pair, and production of additional possibly hard jets radiated from the incoming parton lines. The details of the initial state radiation are expected to depend mostly on the mass of the produced particles, which sets the scale for the radiation, and on the type of initial state partons (whether the production is from valence quarks, quark-antiquark or gluon fusion).
In this section we compare the QCD radiation pattern as obtained from a parton shower alone (Pythia) and from a matrix element plus parton shower (MadGraph/MadEvent + Pythia) approach, for some key processes with heavy colored final states. To simplify the discussion and reduce the model dependence, we study specific heavy final states which are left undecayed, and give all distributions after parton showering but without hadronization. We compare the systematic uncertainties involved in parton showering by using the two different showering implementations of Pythia, the "old" (virtualityordered) and "new" (p ⊥ -ordered) showers, with a range of shower parameters similar to that in [25]. Similarly, for the matching we use the same shower implementations and parameter variations, and furthermore use two different matching schemes, the k ⊥ -jet MLM scheme for virtuality-ordered showers and the "shower k ⊥ scheme" for p ⊥ -ordered showers, as discussed in Sec. 2.1.
For most of the studies below, the shower parameter that is varied is the starting scale of the shower. This is the most important parameter to determine the hardness of radiation allowed in parton shower emissions, and the default value has varied over the years. Originally, the default scale was set to the factorization scale µ 2 F of the process (typically m 2 ⊥ = p 2 ⊥ + m 2 of the produced particles), inspired by the notion of (initial state) parton showers being the deconvolution of the DGLAP equation for the parton density functions. However, with the so-called Tune A [27], based on W and Z boson production at the Tevatron, this was replaced by 4µ 2 F . Later (from v. 6.319), this upper bound was increased, for certain types of processes, to s pp , the center of mass energy of the proton-(anti)proton collisions, in order to further improve the description of extra jet radiation in hadron collisions. Following the nomenclature of Ref. [25], we will in the following call showers with a starting scale corresponding to µ 2 F "wimpy showers", and showers with the starting scale s pp "power showers".
To study the impact of QCD scale choices, we also vary the factorization and renormalization scales by a factor 1/2 and 2, in the matrix element as well as the parton showers.
We consistently simulate X+ jets for X =gg,q ( * )q( * ) ,gq ( * ) and tt, with up to two extra partons from the matrix element simulations.

Variation of matching parameters
Before starting the comparisons between unmatched (showered-only) and matched distributions, a series of "sanity checks" on the matched distributions were performed. These include the requirement of a smooth transition between the region of phase space described by the shower (below the matching scale Q match ) and the region described by matrix elements (above Q match ), and stability of distributions as the matching scale is varied. Since both the matching implementations employed here rely on the Durham k ⊥ measure to achieve the separation of the phase space, the most revealing distributions to study their features are the differential jet rates defined according to the same measure. In particular, in the k ⊥ -jet MLM scheme there is at parton level a sharp division in the jet rates between the shower and matrix element regions, making it very easy to see to which extent the transitions are smooth. For the shower k ⊥ scheme, as well as the cone jet MLM scheme implemented in Alpgen and the CKKW scheme implemented in Sherpa, the separation is less sharp, but the differential jet rates still tend to be the best variables to study the transition between parton showers and matrix elements.
Our guidelines for the choice of the scales Q ME cut and Q match are based on smoothness of distributions across the matching transition as well as efficiency. The higher the scale can be chosen, the higher will the proportion of lower-multiplicity events be in the combined sample. Since the computational effort is heavier for higher multiplicities, it is desirable to choose the matching scale as high as possible. Therefore for each choice of particle type, particle mass, and shower, the matching scale is chosen to be close to the highest scale that still gives a stable matched cross section and smooth differential distributions, in particular differential jet rate distributions. As it turns out, the "new", p ⊥ -ordered Pythia showers allow significantly higher choices of the matching scale than the "old", virtuality-ordered showers. The reason for this is that they give significantly harder emissions than the old showers, and therefore give distributions more similar to the matrix element distributions.
In the following, our default scale choices are: • For pair production of SUSY particles with heavy mass (600 GeV and above): For virtuality-ordered showers, Q ME cut = 40 GeV and Q match = 60 GeV, for p ⊥ -ordered showers Q ME cut and Q match = 100 GeV.
• For Standard Model tt production or light gluino par production: For virtualityordered showers, Q ME cut = 20 GeV and Q match = 30 GeV, for p ⊥ -ordered showers Q ME cut and Q match = 100 GeV.
In Fig. 2 we show the effects of varying the matching scale by a factor 1/2 and 2 from the default values, for 607 GeV gluino pair production and Standard Model top pair production, for the new and old Pythia showers. We have here used the default starting scales for the respective Pythia showers, corresponding to the factorization scale µ F for the p ⊥ -ordered showers and 2µ F for the virtuality-ordered showers.
In the right-hand column, the curve for unmatched Pythia showers with default parameters are shown, together with the pure matrix element prediction without any parton showering or matching and the matched curve. We see, as expected, that the matched curve smoothly interpolates between the unmatched Pythia curve below the matching scale, and the matrix element prediction for large scales.
We also see how the variation in the matched curves is small as the matching scale is varied, and the curves above and below the matching scale variation limits are quite stable.
The difference between the two matching schemes is visible in the plots, in the different behaviour of the parton multiplicity sample contributions. The lefthand column for each particle type shows the the k ⊥ -jet MLM matching scheme, with the contributions from the different parton multiplicity samples in grey. The matching scale cutoff is, in this scheme, done in the same variable that is plotted, the differential jet rate, and there is therefore a sharp cutoff between the 0-and 1-parton samples in DJR(0 → 1), and between the 1-and 2-parton samples in DJR(1 → 2), so that below the cutoff only the lowermultiplicity samples contribute and above the cutoff only the higher-multiplicity samples. In the middle column, the "shower k ⊥ " scheme is used (with the p ⊥ -ordered Pythia showers). This scheme cuts on the first emission of the parton shower rather than on the combined radiation of the whole shower, giving some smearing across the matching scale. This scheme therefore allows to use the same cut at matrix element level and matching level. The distributions for the p ⊥ -ordered showers have been double-checked using the k ⊥ -jet MLM matching method, with excellent agreement.
showers DJR(2    Interesting to notice is the differences in curve shapes depending on the choice of shower type. Below the matching scale, the shape of the curve is given completely by the shower, in particular for the 0 → 1 jet rate. Above the matching scale, however, the shape is mainly given by the matrix element. It is easy to see the reason for the different choices of matching scales for the different showers -the p ⊥ -ordered shower gives significantly harder Pt of the 2nd jet, unmatched t t Figure 3: p ⊥ spectrum for the first and second hardest radiated jet in (a)gg and (b) tt events. The right column in each group of plots shows the spread of Pythia predictions with different choices of the shower evolution variable (virtuality-and p ⊥ -ordered) and starting scale for the evolution (labeled as "wimpy" and "power" showers respectively). The left column presents the results obtained after matching in the same four shower scenarios. The grey band shows the systematic uncertainty associated with a variation of a factor of two of the renormalization and factorization scales with respect to their central values. The different curves have a normalization corresponding to their cross section. The gluino mass is here 607 GeV, while the top mass is set to 174 GeV. distributions than the virtuality-ordered shower, and is more similar to the matrix element curve, hence allowing a higher matching scale.

Parameter dependence in matched and unmatched generation
One of the advantages of the parton shower formalism, and also one of the arguments for using it, is that parton shower Monte Carlo generators have several parameters on which the behaviour of the shower depends, that can be tuned to the data. While this is certainly an advantage in general, it also means that the parton shower lacks predictability at least for some observables or areas of the phase space. It is not always clear that a tune done for one type of initial state will be applicable to other initial states, or that a tune done for a particular mass of a pair-produced particle will be applicable for other masses.
We first compare the matched and unmatched spectra of the first and second jet p ⊥ in Standard Model tt production events and SUSYgg production with mg = 607 GeV, for different shower parameters, Fig. 3. Here, we look at only jets coming from QCD radiation, leaving the heavy particles undecayed, in order to specifically study the differences in treatment of radiation only. The different curves correspond to the different shower parameters settings: the "old", virtuality-ordered shower and the "new", p ⊥ -ordered implementation, each with two different choices for the starting scale of the shower, µ F ("wimpy shower") and s pp ("power shower"). In order to study the intrinsic QCD uncertainty on the predictions, we also vary the factorization and renormalization scales for one of the parameter settings by a factor 2 up and down. Note that we here vary all scales together, both for the central process and for the parton shower/QCD radiation in the matching.
Several interesting features can be noted from Fig. 3: • The spread in predictions for the parton shower is very large and strongly affects the shapes of the distributions. This uncertainty due to shower parameters is almost completely removed when matching is applied.
• The region where the shower predictions start to diverge, and the rate of this divergence, is strongly correlated with the mass of the produced particles. This correlation is due to the choice of starting scale for the "wimpy" showers as the factorization scale, which is close to the mass of the produced particle.
• The "power" shower curves consistently overshoot the matched curves, and hence give too hard predictions, while the "wimpy" showers give too soft distributions.
• The uncertainty due to scale variations is considerable, but mainly affects the normalization and only to a small degree the shape of the curves.
The present Pythia default for the virtuality-ordered showers is close to the curve for the "wimpy" p ⊥ -ordered showers, which is also the default for the p ⊥ -ordered showers.

Impact of different inital states
The impact of different initial states is shown in Fig. 4, where we compare radiation for pair production of gluinos, up-type squarks and up-type squark-antisquark. For ease of comparison, the supersymmetric particle masses have all been set to 607 GeV. These different SUSY particles are produced by different initial states: gluino production is mainly through gluon fusion, squark-squark production from valence squarks exchanging a t-channel gluino, and squark-antisquark production from a combination of gluon fusion and quark-antiquark with t-channel gluio exchange. Besides the p ⊥ of the two hardest radiated jets, we also  Figure 4: H rad T and p jet ⊥ distributions for radiated jets ingg,qq andqq * production for mg ,q = 607 GeV. The grey band shows the spread of unmatched Pythia predictions with varying shower parameters, while the green band shows the corresponding matched predictions. The full black curve represents the default unmatched Pythia Q 2 -ordered shower, while the full red curve is the default matched curve forgg production, included in the other plots as a guide for the eye. Each curve has normalization relative to the Born cross section, which in turn is normalized to unity in order to allow comparison between the different samples.
show the total H T of radiated jets, defined as where the sum is taken over radiated jets with p ⊥ > 40 GeV, as a measure of the total QCD radiation activity. An interesting observation from Fig. 4 is that there is no choice of shower parameters that describe all three types of particle production. For squark pair production in particular, the shower always undershoots the matched prediction, except for the hardest shower choice (p ⊥ -ordered power shower), which is much too hard in the other productions. Indeed the matched curves are much more similar between the different productions than the shower curves, in particular in the tails (as illustrated by the inserted matched gg comparison curve). This indicates that the shower retains a too strong "memory" of and p jet ⊥ distributions for radiated jets ingg production for mg = 300, 600, 1200 GeV. The grey band shows the spread of unmatched Pythia predictions with varying shower parameters, while the green band shows the corresponding matched predictions. The full black curve represents the default unmatched Pythia Q 2 -ordered shower, while the full red curve is the default matched curve for mg = 600 GeV, included in the other plots as a guide for the eye. Each curve has normalization relative to the Born cross section, which in turn is normalized to unity in order to allow comparison between the different samples. the initial state, while the matrix element displays a larger independence from the initial state.
It is also interesting to note that the default virtuality-ordered Pythia shower undershoots all curves, but least so for gluino production and most notably for squark-squark production, where it undershoots the matched curve by at least 50% starting already below a p ⊥ of 200 GeV.

Produced particle mass dependence
Determination of the absolute masses of produced particles at the LHC is suprisingly nontrivial if the produced particles decay to an invisible stable massive particle such as a WIMP. This difficulty is due to the fact that only the kinematics of decays are due to mass differences only, and are relatively independent of the absolute masses involved. Many methods have been devised to extract this information, including the use of different trans-  Table 1: Benchmark scenarios employed in this work as modifications of SPS1a [34]. We always assume the squarks decaying 100% into quark+lightest neutralino.
verse mass quantities (see, e.g., [28,29]) and combined information on kinematics, cross sections and branching ratios [30,31], although especially the latter tend to be quite model dependent. Another possiblility would be to take advantage of the spectrum of QCD radiation, since this is directly sensitive to the absolute mass scale rather than mass differences. Studies have been done on the effect of QCD radiation in the decays of squarks and gluinos [32,33]. In a hadron collider however, the majority of visible radiation originates from the initial state. As can be seen from Fig. 5, the predictions for the shapes of the spectra of radiated jets are quite precise once matching is taken into account, while it is clear that this type of study cannot be done with Pythia radiation only.
For the lightest SUSY scenario (300 GeV gluinos), the distribution of the leading radiation jet is quite distinctive from the heavier scenarios, while the difference between the 600 GeV and 1200 GeV scenarios is small. Even barring the difficulties to distinguish radiation jets from jets from decays, it would be very difficult to differentiate between different high-mass scenarios using the distribution of radiated jets. One can however imagine scenarios where this could be a useful crosscheck, e.g. in the case where the mass splitting between the produced QCD particle and the stable WIMP is relatively small (such a scenario will be revisited in Sec. 4.4 below). More work on this would clearly be needed, but since this is beyond the scope for this paper, we choose to save it for a future publication.

Anatomy of E ⊥ + multi-jet final states
While the previous section dealt with jets from QCD radiation only, in models with new heavy QCD states produced at the LHC we expect jets also from their decays. It is therefore necessary to study the impact of jet matching in the context where we include also such jets. One could expect that jets from QCD radiation should be relatively unimportant compared to the hard decay jets from heavy particles, but as we will show, there are many situations when this is not true. Indeed, as it turns out, there are surprisingly many scenarios where radiation jets can significantly alter the analysis or interpretation of data.
In order to clarify the discussion, and keep our results conservative, in this section we will use a set of simplified supersymmetric benchmark scenarios, summarized in Table 1. In all the scenarios we assume all light-flavour squarks to have the same masses and that they all decay directly to the LSP, i.e. we ignore the existence of intermediate weak states. Introduction of cascade decays will have as main effect that jets from decays get softer, while the jets from QCD radiation are not affected, and will hence mainly further accentuate our results.
For ease of comparison between the scenarios, we have chosen to use the same masses, around 600 GeV, for the active heavy QCD states in all scenarios. Scenario I has a SUSY QCD spectrum similar to the SPS point 1a [34], with a gluino at 607 GeV which decays to squarks at 560 GeV, while the LSP is at 100 GeV. In scenario 2, the gluino has a mass of 607 GeV but all squarks are heavy, so that the gluino decays through offshell squarks to two quarks and the LSP. Scenario III has squarks at 560 GeV and the gluino too heavy to be produced at the LHC, and finally scenario 4 has a gluino at 607 GeV decaying through offshell heavy squarks, but the LSP mass is 500 GeV, only 100 GeV lighter than the gluino. These scenarios will act as "cartoons" to illustrate different effects of QCD radiation in the production of new heavy QCD states.
All plots and results in this section are generated using matched samples with the default Pythia parameter choices for virtuality-ordered showers unless otherwise stated.

H jet
T variables in gluino production It is something of common lore that, in order to study squark pair production one should select 2-jet observables, for associated gluino-squark production at 3-jet observables and for gluino pair production at 4-jet observables. While it is obviously true that gluinos must decay to two quarks and a color singlet (barring the exotic possibility that the dominant decay is the loop-mediated twobody decay to a gluon and a color singlet), it should be kept in mind that the visibility of these jets depends strongly on the mass hierarchy of QCD states. In Fig. 6 we show the H jet T (n), defined as for n = 2, 3, 4, for the scenarios 1 and 2. The sum in eq. (4.1) is taken over jets defined using the SISCone [35] algorithm with a radius of 0.5 and p ⊥ > 40 GeV. We also show the composition of H jet T (n) in terms of jets from the gluino decay and radiated jets ("ISR"), as well as show the average fraction of the H T coming from the decay. For scenario 2, where the gluinos decay through off-shell squarks to two quarks of similar energies and an LSP, the majority of events in the peak of H jet T (4) include only jets from the decay, while the tail of the distribution is dominated by 2-3 jets from the decay and one jet from radiation. For scenario 1 however, where the 600 GeV gluino decays into a fairly soft jet (with an energy around 50 GeV) and a squark which in turn decays to a hard jet and an LSP, the distributions are quite different. Here, H jet T (4) is dominated by events where at least one of the jets come from radiation, and we need to go down to H jet T (2) to be dominated by events with only decay jets across the whole H T range. The energy fraction of H jet T (4) coming from decay is still fairly high, even in the high-H T tail, since most of the transverse energy comes from the squark decay jets.
The immediate interpretation of this result is that in a scenario with a small mass splitting between gluinos and squarks, gluino production might be difficult to distinguish from squark production with additional QCD radiation. In Section 4.3 we will study this in greater detail.

Jet multiplicities for different scenarios
A question that the proper matching of jets is particularly apt to answer, is to specify the number of jets typically present in production of different particles. We here present a table with the jet multiplicities, for matched and unmatched (Pythia virtuality-ordered default) production, for the scenarios studied. In order to make the table as useful as possible, we have used jet cuts close to what is used in many preparatory analyses for squark and gluino searches: p jet 1 ⊥ > 180 GeV, p jet 2 ⊥ > 110 GeV, p jet j ⊥ > 50 GeV for j > 2 and |η jet | < 3 for all jets. We require all events to have at least two jets. The jet multiplicities are exclusive, and so add up to 100% of the events passing the 2-jet cut. The jet algorithm used is SISCone with a radius of 0.5 and p ⊥ > 40 GeV.
The difference between matched and unmatched generation (as seen in Fig. 4) is twofold: First matching tends to increase the transverse boost of the produced pair, and  Table 2: Contribution of events with N jets for matched and unmatched processes:gg,gq andqq in Scenario I (for squark production this is very similar to scenario III), andgg in Scenario II and IV. All numbers are in percent. "Signal efficiency" shows the percentage of events that pass the 2-jet cut. The jet cuts are described in the text.
hence the p ⊥ of the softest jets from the decays. Second, the ISR jets get harder, and more easily get above the threshold of 50 GeV. The matched production therefore in general populates higher jet number bins the unmatched ones. The first three lines in Table 2 represent Scenario I, which is very similar to the benchmark point SPS1a. In this scenario, the gluinos decay asg →qq →qqχ 0 (with m χ 0 ∼100 GeV), so there are typically two hard and two soft jets from the gluino decay.
Here the effect of the matching is large, due to the increase in p ⊥ for the soft jets as well as increased hardness of the ISR jets. For scenario II on the other hand, where gluinos decay in three body decay into two jets and a light neutralino, the sensitivity to the matching is much lower, since there, typically at least four reasonably hard jets are present from the decays. It is only in this and similar cases that the statement that gluino pair production generally corresponds to four hard jets in the event is true. In the last row in the table, the produced gluinos decay to a near-degenerate LSP, meaning that most hard jets are due to QCD radiation. Here, a large recoil against initial state jets is needed in order to even pass the 2-jet cut, hence the very low signal efficiency.
For squark-squark and gluino-squark production for Scenario I (row 2 and 3 in Table 2, there are only two hard partons from the decays (and one additional soft parton from the gluino decay) and additional QCD radiation. Again the main consequence of the matching is to increase the mean number of jets. It is interesting to note that in this particular case, the addition of matching to the generation of squark pair production gives very similar numbers as the unmatchedgq associated production (although the selection efficiency is different), indicating that the matching has a similar impact as the addition of one extra jet with p ⊥ ∼ 50 GeV.
We now look closer at two examples where an analysis based only on a parton shower approach might lead to the wrong conclusions.

Example 1: False gluino evidence
We here consider the case given in Scenario III, where the only observable SUSY particles are 560 GeV squarks, while the gluinos are too heavy to be produced. If such a scenario is simulated using only parton showers, the result will typically be a deficit in the number of events with multiple jets. This is illustrated in Fig. 7, where we show H jet T (2), H jet T (3) and H jet T (4) (defined in eq. (4.1)) for matched and unmatched simulation ofqq ( * ) production, using the Pythia virtuality-ordered shower with default parameters (the full black and full red curves, respectively). While the unmatched generation well reproduces the 2-jet H T , it increasingly falls below the matched curve for H jet T (3) and H jet T (4). In real data, this deficit could easily be interpreted as a sign for a missing production mode, such as associated production with gluinos and squarks, or gluino pair production.
Also shown in Fig. 7 are the curves for associated gluino production (dotted) and gluino pair production (broken red), with gluinos at 700 GeV, decaying to the squarks emitting an additional jet. While the cross section for the associatedgq production process is too large to be accommodated by the "data" (i.e. the matched squark production curve), there are many ways in which this cross section estimate could be wrong; if, for example, the decay modes are different for squarks and gluinos, or if the particles seen are not from the MSSM but from some other realization of new physics. We can therefore view the normalization of the different curves as free parameters, and look at how the curves could be combined to generate a fit to the "data" curve, i.e. the matchedqq ( * ) curve.
It turnes out that for these choices of masses, the data can be very well accommodated by adding the gluino-squark associated production, with the cross section reduced to 25% of the nominal cross section, to the unmatched squark production curve, as shown by the broken blue line in the figure.
Possible cross checks to avoid this type of false conclusions could be looking at detailed kinematic properties of the jets, or looking at different decay signatures. However, this example still illustrates the importance of using matching also for BSM signals, in particular when looking at observables where extra jets from QCD radiation may play an important role.

Example 2: Degenerate spectrum
Another case where jet matching is of considerable importance is when the mass splitting between the produced QCD particle and the LSP is small, as is the case in Scenario IV. There, 600 GeV gluinos decay through off-shell squarks to a neutral 550 GeV LSP and two quarks. The main problem with searches in this scenario is that no large missing E T is produced in the decay, since all decay products are soft. It doesn't help if the gluinos are produced at a large invariant mass, since they will then decay back-to-back, and the missing energy due to the boost of the individual gluinos is canceled between the two LSP's. While this type of scenario is not possible in typical unification scenarios (e.g. mSUGRA or mGMSB), a more model-independent approach necessitates taking them into account [36,37]. The only way to get a large missing energy in this scenario is through a transverse boost of the gluino pair center-of-mass system. This happens when the gluino pair system recoils against hard initial state radiation. In this case, a proper jet matching between parton showers and matrix element QCD emission is crucial to well describe the E ⊥ distribution as well as the resulting jet structure. This is particularly important since a E ⊥ of a certain magnitude as well as several hard jets is typically needed to pass the experiment triggers, so that a misrepresentation of these quantities will lead to very different trigger acceptences.
We show in Fig. 8 three representative quantities for this case, the E ⊥ and the p ⊥ of the two hardest jets. As can be seen from the figure, the virtuality-ordered Pythia default distribution severely underestimates the missing E T as well as the hardness of the leading jets, which mainly come from the initial state radiation. The effect, with typical cuts jet-E ⊥ cuts like E ⊥ > 200 GeV, p jet 1 ⊥ > 180 GeV and p jet 2 ⊥ > 110 GeV is a signal efficiency less then half that of the matched production. The p ⊥ -ordered shower rather overshoots the matrix element curves, up to the factorization scale (mg), where it falls off rapidly. Once again, the lesson here is that in order to get a description that is predictive and insensitive to the details of the shower parameterization, it is necessary to use jet matching.

Impact on BSM searches
We now consider the effects of an accurate simulation of QCD radiation in the typical observables employed in the BSM inclusive searches: high p ⊥ jets and high missing transverse energy. In order to see the effect of the matching on the sensitivity to shower after a smearing of the signal, and in order to be as complete as possible, we consider the production ofgg,gq,qq ( * ) andt 1,2t * 1,2 , withq defined asũ L,R ,d L,R ,s L,R ,c L,R ,b 1,2 . The signal is produced in both the matched (2→2, 3, 4) and unmatched (2→2) modes in the the SPS1a benchmark scenario [34]. For the background, we consider the most important processes leading to four hard jets and potentially large missing transverse energy: W ± → l ± ν + 4 jets, Z 0 → νν + 4 jets, W ± → τ jet ν + 3 jets and finally the inclusive tt+0,1,2,3 jets. We do not include QCD multijet production, since we have no means of realistically performing simulations of the missing energy distribution, which is due to decays of heavy quarks to neutrinos and jet mismeasurement in the detector. We instead base our analysis on cuts similar to those used in Refs. [38] and [39], and keep this contribution in mind. All background simulations are done using jet matching. Many comparisons have been done between matched and unmatched background simulations, which are well known to differ by up to several orders of magnitude for this type of multi-jet observables [15,40,12]. We therefore here look at the effects of including matching only in the signal simulation.
Detector simulation is performed using PGS 4 [41] with the MidPoint cone algorithm with a minimum p ⊥ of 40 GeV and a radius of 0.5. To be conservative, we use only kinematic variables associated to the jets and the missing transverse energy. The cuts used are The SPS1a scenario is affected by several of the difficulties described in Sec. 4. The gluino has a mass higher than, but close to, the squarks, and hence decays to a squark and a soft jet, making the jet counting complicated. Since gluinos and squarks are of similar mass, the QCD SUSY production includes associated squark-gluino production, gluino pair production and squark pair production (bothqq andqg * which are of similar  cross section, and to a negligible degreeq * q * ), in order of cross section. This means that the Pythia shower cannot simultaneously describe all production modes, as demonstrated in Fig. 4. Furthermore the separation of the different production modes based on jet counts or jet kinematics will be non-trivial. The production cross section is however in this scenario dominated by associated gluino-squark and gluino pair production, where the default Pythia description is reasonably good, and only undershoots the matched description by about 10-30% for the first couple of QCD radiation jets, so we expect the inclusive Pythia description to be reasonably close to the matched curve.
The result is illustrated in Fig. 9, which shows the H T = 4 i=2 p jet i ⊥ + E ⊥ for inclusive supersymmetric production of gluinos and squarks (includingt andb).
The effect of the matching is as expected -a significant reduction in the sensitivity to parton shower parameters, and a shift of the prediction as compared to the default virtuality-ordered Pythia shower by about 10-30%. Even with the smearing due to the detector simulation and more complex decays than the simplified scenarios used in Sec. 4, the "power" p ⊥ -ordered shower continues to overshoot the result obtained with the matching whereas the "wimpy" virtuality-ordered showers undershoot the matched curve.
Since we require four hard jets, the strongest impact is on the squark-(anti)squark pair production, with an efficiency increase close to 40% when passing from the unmatched default virtuality-ordered Pythia shower to the matched production. This happens for two reasons; as described in Sec. 4.2, only two hard jets are produced by the decay of the squarks, which means that two jets from QCD radiation are needed. Second, as shown in Fig. 4, the difference between the unmatched and matched radiation is particularly large for squark production. On the other hand, the effect is much smaller for gluino pair production (around 5%), which has to do mainly with the large fraction of events with at least one top quark in the decay, giving rise to multiple hard jets and hence a small sensitivity to the matching. The effect on the associated gluino-squark production lies between these two extremes, with an efficiency increase close to 10%.t andb production is dominated byt 1t * 1 pairs (due to a lowt 1 mass, around 390 GeV), which decay to top quarks and the LSP. For those, the effect of the matching is negligible.

Conclusions
Discovering new physics at the LHC will most probably be quite challenging. Apart from a few very clean but theoretically not very motivated signatures, such as narrow resonances decaying into lepton pairs, the vast majority of the theoretical constructions which are still viable nowadays gives a rich and complex phenomenology that will be hard to decipher. For instance, many models that address the hierarchy problem, predict new heavy colored states at the TeV scale that can be easily produced at hadron colliders. To what extent the decay products of such states will be identifiable on top of a large QCD background, including vector boson(s) + jets, tt jets and multi-jets, is currently subject of many theoretical and experimental investigations.
In this work we have for the first time addressed the issue of quantifying the effects of extra QCD radiation in the production of heavy colored states employing inclusive multi-jet samples obtained by matching multiparton matrix elements and parton showers. Previous work has either considered radiation in Standard Model tt events only [12] or in SUSY but at parton level [25].
Our results can be briefly summarized as follows. First the extension of matching techniques to beyond the Standard Model scenarios, such as SUSY or UED, while posing no problems of principle, requires dealing with new technical issues. For instance, the combination of multi-parton samples for production of different resonances decaying into jets leads to problems of double counting that need to be addressed. We have proposed two working solutions, both of which have been implemented in MadGraph/MadEvent.
The main result of this work is that a matched matrix element plus parton shower approach for heavy particle production is in general much more accurate and predictive than a parton shower alone. We find that, contrary to the "common lore" that showers alone provide a good description of extra radiation when the produced states are heavy, there are many cases where matrix element corrections are indispensable. Not only are the p ⊥ spectra of the extra jets in the parton shower approach extremely sensitive to the shower starting scale and the shower algorithm, while the matched simulations are not, but also that there is in general no tuning of the shower that can simultaneously reproduce the matched samples for all initial states. We have presented examples where simulations based only on the parton shower could be misleading. Another possibly relevant issue, that we leave to future studies, is matching of radiation in the decays of the heavy states. Though already possible in our current implementation, we have not included it in this study mainly because it is computationally expensive. Next-to-leading order effects in decays have previously been shown to be possibly important in precision studies, e.g. of the spins of new particles [33]. If found relevant, also matching of radiation in decays could be included in future simulations.
In conclusion, we recommend that matched samples should be used not only for backgrounds but also for beyond the Standard Model physics signals. Having the best available simulations could be important not only in designing better and more solid strategies to make discoveries but even more to identify the nature of the new physics. For precision measurements of new physics properties, as well as reliably distinguish between different scenarios, fully differential next-to-leading order simulations will most probably be needed. Given the large number of theoretical possibilities still open, it is clear that an automatic approach to NLO computations, and to their matching with a parton shower, will be certainly welcome.