NLO QCD predictions for $Z+\gamma$ + jets production with Sherpa

We present precise predictions for prompt photon production in association with a $Z$ boson and jets. They are obtained within the Sherpa framework as a consistently merged inclusive sample. Leptonic decays of the $Z$ boson are fully included in the calculation with all offshell effects. Virtual matrix elements are provided by OpenLoops and parton shower effects are simulated with a dipole parton shower. Thanks to the NLO QCD corrections included not only for inclusive $Z\gamma$ production but also for the $Z\gamma$ + 1-jet process we find significantly reduced systematic uncertainties and very good agreement with experimental measurements at $\sqrt{s}=8$ TeV. Predictions at $\sqrt{s}=13$ TeV are displayed including a study of theoretical uncertainties. In view of an application of these simulations within LHC experiments, we discuss in detail the necessary combination with a simulation of the $Z$ + jets final state. In addition to a corresponding prescription we introduce recommended cross checks to avoid common pitfalls during the overlap removal between the two samples.


Introduction
The production of a Z boson is one of the standard candle processes at hadron colliders like the LHC. The massive boson is often produced in association with photons, which are typically low-energetic or collinear with charged final state particles. Cases where one of the photons happens to be well-isolated and high-energetic can be regarded as an individual important final state, Zγ production.
Zγ production plays an important role both as a signal and as a background process at the LHC. The absence of couplings of the photon to the uncharged Z boson in the Standard Model can be probed by measurements in this channel, resulting in differential cross sections and limits on anomalous couplings from LEP experiments [1], Tevatron experiments [2], and by ATLAS [3] and CMS [4].
The Zγ process also constitutes an irreducible background in the search for the Higgs boson [5] or new massive gauge bosons decaying to Zγ [6] or in more inclusive searches in final states containing a photon and missing transverse momentum [7,8].
Theoretical predictions for Zγ production can be divided into onshell and offshell calculations. Results for onshell Zγ production leave out the decays of the Z boson or include them only in a narrow-width or pole approximation. Beyond the leading-order results [9], the first onshell higher-order calculations included NLO QCD [10] and, more recently, also NNLO QCD corrections [11]. In reality, the Z boson is unstable and is thus never produced as an onshell final-state particle. Recent calculations take this finite width into account and provide predictions for the offshell γ final state. The most accurate offshell predictions contain NNLO QCD [12] and NLO EW corrections [13].
Experimental analyses at the LHC rely heavily on theoretical predictions for their signal and background processes. Fixed-order predictions as listed above can provide such an input only to some extent. While they describe the dominant features of the given final state objects at the parton level, they are not made to simulate a realistic behaviour of the hadronic final state. Monte-Carlo event generator programs on the other hand combine fixed-order predictions and an approximate all-order resummation of QCD corrections to enable a full simulation at the hadron level. Different approaches and programs are available, but the simulation currently in use in LHC experiments for Zγ production is generated with leading-order multi-leg generators using approaches like CKKW(-L) or MLM merging [14]. For the related process of W γ production, an implementation of a NLO QCD calculation of the inclusive process matched to a parton shower exists within the Powheg framework [15].
This article applies a NLO accurate multi-leg merging formalism to the processes of Zγ and Zγ+jet production 1 for the first time. After a review of the relevant methods in Section 2 we present our computational setup and results comparing LO and NLO accurate merged predictions to each other and to experimental data in Section 3. A special aspect relevant for the application of multi-jet merged Zγ samples in experimental analyses is a combination of Zγ+jets and Z+jets. Recommended techniques and cross checks for such a combination are implemented and discussed in Section 4. 1 Even though the process is denoted with the shorthand Zγ, the calculations throughout this paper include the full offshell γ final state.

Matching and merging with Sherpa at NLO
To obtain NLO accurate multi-jet-merged predictions the "MEPS@NLO" formalism [16] is applied to Zγ production within the Sherpa framework. It combines two essential ingredients, NLO + parton-shower matching and multi-jet merging, which are briefly summarised in the following. For the combination of NLO-accurate matrixelement calculations with a parton shower (PS), a matching procedure to avoid double counting of the QCD emission effects at O(α s ) is needed. Here, the prescription proposed in [17] is applied, both to pp → Zγ and pp → Zγ+jet simulations. It is based on the original MC@NLO method [18] but extends it to a fully colour-correct formulation also in the limit of soft emissions.
Both simulations and further PS emissions are then consistently merged into one inclusive pp → Zγ+0,1j@NLO sample using the "MEPS@NLO" method [16]. Shower emissions above a pre-defined separation criterion, Q CUT , are vetoed in the lower multiplicity contribution, and Sudakov factors are applied where appropriate such as to make this contribution exclusive and allow the combination with a higher jet multiplicity. This applies not only to the emissions from the parton shower but to all contributions of the NLO-matched emission, including the hard remainder ("H events"). In the application of the Sudakov factors care has to be taken to remove the O(α s ) contribution which is already present in the NLO(-matched) emission.
Beyond the processes simulated at NLO accuracy, higher-multiplicity processes can be added to the simulation at LO accuracy to improve the modelling of high jet multiplicities beyond the parton shower approximation.
Similar to merging methods at LO, an appropriate scale choice for the evaluation of multi-jet configurations is obtained by statistically identifying a parton shower history in the matrix-element final state. To that end, the parton shower is run in reverse mode, i.e. the closest parton pair is identified according to the shower splitting probabilities and then recombined using the kinematical properties of the shower. When applied recursively, this clustering results in a core process, e.g. pp → Zγ, and in an ordered history of shower emissions. The factorisation scale is then determined by a typical momentum transfer within the core process (core scale), and the renormalisation scale is calculated from the n identified branchings to resemble where k ⊥,i is a suitably scaled relative transverse momentum encountered in the ith splitting. For events with very hard QCD emissions, the clustering may include electroweak combinations to preserve the ordering in the history ("inclusive clustering"). In such events, the core process can also contain jets and contribute to the renormalisation scale accordingly.

Soft photon resummation with YFS
The YFS-algorithm [19] describes a possibility to resum soft logarithmically enhanced photon radiation to all orders in a process-independent manner. This algorithm is implemented in Sherpa for both, leptonic decays of W / Z bosons and hadron decays. Details of this implementation are given in [20] and briefly summarised here. The decay width for a decay of a particle i with mass M into a final state f , corrected for the radiation of n R real and n V virtual photons, reads (1) Here, M n V + 1 2 n R n R describes a decay matrix element with additional n R real and n V virtual photons and Φ denotes the corresponding phase space. As shown in [19], the soft limits of all these (virtual and real) matrix elements can be resummed and factorized out. The corresponding infrared divergences cancel order by order and result in the finite YFS form factor Y (Ω). Y (Ω) includes the soft limits of all virtual and real contributions. However, only the divergent part of the real contributions is retained and separated by the symbolic cut off Ω from the non-divergent one. Eq. (1) can then be approximated with The eikonals S(k) include all contributions for the emission of a real photon with momentum k. Since the divergent, real part is already included in Y (Ω), these eikonals will only be integrated in the non divergent phase space. The number of additional, resolved photons is denoted with n γ , while β 0 0 is the undressed matrix element without any real or virtual photons and Φ the corresponding phase space. The YFS algorithm can in principle be improved order by order with exact, process dependent real and virtual matrix elements.
These possible corrections are incorporated in the factor C, which is equal to one in case of no corrections. In Sherpa, higher order corrections are available either in an approximative way using collinear splitting kernels or as exact matrix elements. The first option is independent of the process whereas the second one is limited to a few cases, including the decay of a vector boson into two fermions.
It is worth noting that in a Monte Carlo code photons radiated by YFS are -in contrast to e.g. QED showers -unordered.

Isolated photons
Matching theoretical predictions with experimental measurements including isolated photons has turned out to be a non-trivial problem when going beyond leading order in QCD. It is no longer possible to isolate the photon completely from all kind of hadronic activity since this would constrain the phase space of soft gluon radiation and destroy the cancellation of infrared singularities. As a consequence, many experiments allow a small fraction of hadronic energy within the isolation cone. However, this relaxed cone criterion requires special attention when calculating theoretical cross sections. An observable constructed with such an isolation criterion includes a divergence when the photon gets collinear to a massless quark. This divergence is of QED origin and does not cancel within the pertubative QCD calculation. In principle there are two common ways to solve this problem, either by absorbing the divergences into fragmentation functions or by using a smooth isolation criterion. Such a smooth isolation criterion has been proposed in [21] and is also used in this paper. This criterion suppresses the divergent collinear contribution by limiting the maximal transverse energy E max ⊥ close to the photon axis, ( Here, p γ ⊥ is the transverse momentum of the photon and , n and R are parameters which define the final shape of the smooth cone. E max ⊥ is defined as the sum of the transverse energies of all partons present at matrix element level within a cone of radius r around the photon axis. The condition in Eq. (3) has to be fulfilled for all cones with r = (∆η) 2 + (∆φ) 2 < R (4) and ensures that E max ⊥ converges smoothly to zero for r → 0.

Setup
All results in this publication are obtained with the Monte Carlo event generator Sherpa 2.2.2 [22] using merged calculations. Two setups will be compared in the following:
Therein, jet refers to an additional parton in the matrix element and LO/NLO denotes the accuracy of the corresponding multiplicity. The matrix elements are calculated by the internal matrix element generators Amegic++ [23] and Comix [24]. Virtual diagrams are calculated by OpenLoops 1.3.1 [25], using CutTools [26] and OneLoop [27]. In all MEPS@NLO setups Amegic++ is used only for Born-like processes. All real-subtracted contributions and the leading order diagrams of higher multiplicity are calculated by Comix. The merging cut Q CUT is set to 30 GeV.
Where not explicitly stated otherwise, scales are determined by the inclusive clustering algorithm described in Sec. 2.1 (STRICT METS). The core scale is calculated according to the core process as This corresponds to the default core scale implementation in Sherpa 2.2. The electroweak couplings are evaluated using a mixed scheme as recommended in [28]. First, all couplings are calculated using the G µ scheme. In this scheme, the coupling constant is evaluated as function of the Fermi constant G µ and the masses of W and Z, This behaviour effectively resums contributions which arise when evolving the electroweak coupling to the electroweak scale and is a common choice for processes involving heavy W or Z bosons. However, in the V γ processes an additional external -i.e. on-shell -photon is present. Taking this into account one electroweak coupling should be evaluated at α(0). This is achieved by a global reweighting with k = α(0)/α Gµ . Unstable particles are described using the complex mass scheme [29]. Following [13], on-shell masses and widths are converted to pole values and result in These values are used for all calculations.
As parton distribution functions NNPDF3.0 [30] sets are used, taking the NLO set for MEPS@LO calculations and the NNLO set for MEPS@NLO. The running of α s and its value at M Z are thereby set according to these PDF sets, resulting in α s (M Z ) = 0.118 and a running at two(three) loop order when using the NLO(NNLO) sets. YFS is set active and includes matrix element corrections for further photon emissions.
For the event generation, all parton level cuts are set to be more inclusive than the respective analysis cuts. As described in Section 2.3, it is not possible to use the experimental isolation criterion, instead the smooth cone criterion is used and validated for two different sets of parameters.
When comparing to experimental data the simulation is performed including the default multiple interactions [31] and hadronisation models [32]. The final state analyses are done within the Rivet framework [33].

Merging cut variation
Using the ME+PS merging method defined in Section 2.1 a new parameter is introduced, the merging cut Q CUT . It separates the different phase space regions for the parton shower and higher multiplicity matrix elements. Since this parameter is unphysical, physical observables should be independent of its exact value as long it is chosen in a reasonable range. An interesting observable for checking this behaviour are the splitting scales as defined by the k ⊥algorithm [34]. All final state partons 2 are clustered to jets according to this algorithm. The splitting scale d (n−1)(n) is than defined as the jet measure which describes the cluster step of a n-parton final state to a  (n − 1) final state. Following this, d 01 gives the p ⊥ of the hardest jet and d 12 either describes the production of a second jet or the second splitting of the first jet in case there is only one.
In the context of matching and merging, jets can emerge either from higher multiplicity matrix elements or from the parton shower. Since the k ⊥ algorithm uses a jet measure which is very similar to the jet criterion used by the merging algorithm, its splitting scales are very sensitive observables to study the interplay between LO/NLO matrix elements and partonic showers.
In Figure 1, these differential jet rates are evaluated for the first two splittings while varying the merging cut between 20 and 40 GeV. Figure 1a shows the hardest splitting scale. This jet rate is sensitive to the transition from the zero-jet NLO matrix element 3 with a shower emission at Q emission < Q CUT to a one jet NLO matrix element with Q emission > Q CUT .
By contrast, in Fig. 1b the subleading splitting scale is shown. This splitting probes the transition from the one-jet NLO matrix element with an unresolved emission to a LO matrix element with two resolved jets.
In both cases the uncertainties coming from the merging cut variation are less than 10%. As mentioned above these splitting scales are shown since they are expected to be very sensitive to Q CUT . Indeed, other observables have a much smaller merging cut uncertainty. The E γ ⊥ spectrum is given in Figure 1c, here the merging cut variation is found to be completely negligible in contrast to the scale uncertainties studied in Section 3.2.3.
3 A n-jet matrix element refers to a matrix element with n additional, well separated partons in the matrix element. Table 1: This table summarizes all cuts which define the differential cross section analysis for 13 TeV. The leading photon has to be isolated from all other particles by fulfilling the requirement ∆R<0.4 E < iso E γ where the sum includes all particles which have an angular distance of 0.4 or less to the photon axis. In addition, the photon is required to not come from a hadron decay. The leptons are dressed with all photons not coming from a hadron decay and within ∆R < 0.1.

QCD core scale choice
After validating the stability of the MEPS@NLO method, here and in the next section the scale choices and variations will be discussed for a centre of mass energy of 13 TeV. As before, all predictions are performed at shower level, meaning that hadronisation and multiple interactions are switched off explicitly. All phase space cuts for the 13 TeV analysis are summarized in Table 1.
Different scale choices have been employed for V γ production in the literature. Two of them are  and The former was used for the NNLO-QCD calculation of V + γ in [12], the latter is inspired by [35] and was used e.g. in the NLO QCD+EW calculation [36] in a slightly modified manner.
Since this is a merged calculation it is not possible to use these scale definitions directly as described in Sec. 2.1. However, the STRICT METS scale setter also allows to use custom scales for the core process. In order to use this possibility, the clustering has to be restricted such that it results in a Zγ core process. This can be enforced by using an exclusive cluster mode which exactly reconstructs a possible shower history using QCD splittings only. If an allowed history is found, the core process is always pp → e + e − γ. It might also happen, that no possible shower history can be found, e.g. if all possible histories are unordered in terms of the shower variable. In such a case the remaining pp → e + e − γ + n partons process is retained as core process. In both cases the remaining core process is finally evaluated using the core scales defined above.
Altogether, four different scale choices are compared. The first one is "inclusive", it uses the default settings of Sherpa as described in Section 3.1. By contrast, the three remaining setups make use of the exclusive cluster mode introduced above and use different core scales as defined in Fig. 2.
In Figure 2 the differential E γ ⊥ distribution and the p ⊥ spectrum of the leading jet are compared for the four different scale schemes.
In the E γ ⊥ spectrum all scale choices are in good agreement and differ by not more than 10%. Switching from the default to the exclusive clustering algorithm does not make a large difference for this observable. This also holds for all other observables which where measured by ATLAS and are discussed later on in Section 3.3, those measurements do not allow us to discriminate between the different scale choices.
By contrast, the difference is much higher when looking at the p ⊥ spectrum of the leading jet. Here, all scale and cluster choices are in good agreement for low p ⊥ but differ as soon as p ⊥ exceeds 100 GeV. At large values at order of 1000 GeV the difference reaches almost a factor of two. There, the highest cross section corresponds to the core scale defined in Eq. (8) in combination with the exclusive cluster model, whereas the lowest cross section is given by the default settings. This is not surprising since such a high p ⊥ region probes configurations where the jet is harder than the typical scale of this process, e.g. the Z mass. Such configurations are unordered in terms of the parton shower evolution variable and are thus very sensitive to the clustering definition. If here a pp → Zγ +n partons core process is determined but the core scale is evaluated solely based on Z and γ, this scale will underestimate the physical scale and thus overestimate the strong couplings, resulting in a larger cross section. However, based on the available information, there is no clear way to decide which cluster / scale settings are best suited for this process and the inclusive cluster settings are retained. Measuring the leading jet p ⊥ distribution in Zγ events could greatly help to improve this situation.

Scale and PDF variation uncertainties
In addition to the core scale variations studied above, independent variations of the µ R and µ F scales are performed. All replicas of the NNPDF set are used to estimate the PDF uncertainty. These variations include only contributions from the matrix elements but not from the shower. All figures are structured as follows. Each figure has three ratio plots. The main plot and the first ratio plot point out the difference between MEPS@NLO and MEPS@LO. Here, the MEPS@NLO prediction is choosen as reference. By contrast, the additional subplots show the size of all performed scale variations for each method as a ratio with respect to the corresponding nominal prediction. Figures 3 and 4 show predictions for E γ ⊥ and M llγ . The corrections between MEPS@LO and MEPS@NLO are almost flat at a level of 20%. As expected, the scale uncertainties are reduced when moving from MEPS@LO to MEPS@NLO. The factorisation scale dependency vanishes almost completely for p γ ⊥ 60 GeV and m llγ close to the Z peak, compared to 10% in the MEPS@LO case. The renormalisation scale uncertainty is reduced by roughly a factor of two for both observables, it shrinks from 20% to 10% for high E γ ⊥ . Figure 5 shows H ⊥ , which is defined as the sum of the transverse momenta of all jets fulfilling the conditions defined in Table 1, H ⊥ = jets p jet ⊥ . Here both the factorisation and renormalisation scale uncertainties are reduced significantly for small H ⊥ but have almost the same size if H ⊥ exceeds 200 GeV. The leading jet p ⊥ is depicted in Fig. 6 and shows a very similar behaviour. Again the renormalisation scale uncertainty is reduced from 15% to 5% for low p jet ⊥ but does not change for values from 200 GeV onwards. Finally, Fig. 7 shows the jet multiplicity. In the zero-jet bin the factorisation scale dominates when using MEPS@LO, this uncertainty vanishes almost completely when moving to MEPS@NLO. However, the one jet bin is dominated by the renormalisation scale unertainty. This uncertainty is reduced from 20% to 10% when moving to MEPS@NLO.
All of these observables only show minor improvements in the multi-jet regions. This is not surprising since observables which are sensitive to a high number of hard jets are hardly improved by the MEPS@NLO applied here. Only the zero-and one-jet matrix elements are calculated at next-to-leading order but the two-and three-jet calculations still have leading order accuracy. Both H ⊥ and p jet ⊥ are dominated at large values by multi-jet configurations and are thus described at leading order accuracy only. This effect is reflected by the ratio between the MEPS@NLO and the MEPS@LO

Comparison with √ s =8 TeV measurements
Both, the stability and the reduction of the perturbative uncertainities have been demonstrated in the 13 TeV results in the last section. Now, the focus is on the comparison with recent experimental data and a study of the interplay between the smooth isolation criterion with the experimental one. This section relies on a measurement of the ATLAS collaboration at 8 TeV [37], using 20.3 fb −1 of data. In this measurement, final states with llγ and up to three jets are studied. Here, the focus lies on the e + e − γ + jets final state. All cuts which define the extended, differential fiducial cross section are summarized in Table 2.
The first topic to be studied is the isolation criterion used by this analysis. As described in Section 2.3 our prediction uses the smooth cone isolation criterion. By contrast, the experimental isolation is based on antik ⊥ -jets with ∆R = 0.4. These jets include all particles except neutrinos and muons and are not required to fulfil the cuts described in Table 2. A photon is defined as isolated if either the nearest jet has an angular distance ∆R > 0.4 to the photon axis or this jet's transverse en- Even if one assumes that both the final jet E ⊥ and its direction are in perfect agreement with the closest parton E ⊥ at matrix element level, this criterion differs from the one used in our calculation and described in Section 2.3 when going to lower angular distances. As a consequence, two different parameter sets for the smooth isolation criterion are compared. The first set is based on the 2013 Les Houches report [28] which recommends the usage of the smooth isolation criterion for fixed order calculations if the parameters are matched to the experiment. Following this, the parameters are R = 0.4, n = 1 and = 0.5. However, since this is not a fixed order calculation, the smooth isolation criterion is used only at matrix element level and in the subsequent final state analysis the experimental cut has to be passed additionally.
The second parameter set is thus chosen more inclusively in R, here R = 0.1, n = 2 and = 0.1. Such a setup also reflects the requirement that experiments want to generate event samples to be as universal as possible, usable not only as signal process but also as background for many other measurements.
In Figure 8 both these parameter sets are compared with the E γ ⊥ spectrum measured by ATLAS. Both predictions are in good agreement with the data but the more inclusive set gives a slightly higher cross section. This is most obvious in a p γ ⊥ region of around 70 GeV, there the difference reaches almost 10%. Although it is expected that the first set with R = 0.4 may miss some contributions due to the smoothing of the cone, it is not guaranteed that the second set gives a more accurate prediction. A more inclusive parton-level isolation always allows configurations which come closer to the collinear, non-perturbative region. This region cannot be described without fragmentation functions or QED parton shower matching [38]. However, this is not expected to happen if the angular distance is large enough and thus the inclusive parameter set is used in the following.
All measured observables are shown in Figure 9-10. The plots are structured in the same way as in Section 3.2.3. Both methods, MEPS@LO and MEPS@NLO, are compared to the data and directly to each other. In the main plot and the first ratio plot both the MEPS@LO and the MEPS@NLO predictions are com-pared to the measured data. In addition, two further ratio plots show the impact of the described perturbative variations. In contrast to the first ratio plot the respective nominal prediction is chosen as reference here since this simplifies a direct comparison of both methods.
In almost all observables the MEPS@NLO prediction is in excellent agreement with the data. A small deviation is found in the invariant mass prediction with zero jets in a region around 250 GeV. The MEPS@NLO differential cross sections are about 20% larger with respect to the MEPS@LO results at small scales. At large scales the difference gets smaller since the contribution of the additional LO jets is increasing.
As already seen in the 13 TeV section, the uncertainties estimated by the scale variations are reduced when moving from MEPS@LO to MEPS@NLO. At MEPS@LO, for lower values of E γ ⊥ and M eeγ the factorisation scale is the dominant source of uncertainty and reaches up to 10% . This is reflected in the zero jet    [37]. Leptons are dressed with all photons having an angular distance of ∆R < 0.1.
bin of Figure 11, too. By contrast, in the MEPS@NLO case this uncertainty is removed almost completely.
At higher values of E γ ⊥ or m eeγ the renormalisation scale uncertainty takes over in all inclusive observables and reaches values of 10-20 % in the MEPS@LO case. This uncertainty is reduced for MEPS@NLO to 5-10% in both large E γ ⊥ and the lower jet multiplicity bins. Both, the size of the corrections and their uncertainties behave very similarly between 8 and 13 TeV.

Motivation
When predictions for V γ production are used in experimental searches to determine background contributions, they have to be combined with predictions for V +jets production in several cases. A jet from the V +jets sample can be misidentified as a photon at the detector level and thus contribute to the V γ event selection. Another example is a selection requiring multiple leptons, if the photon is misidentified as an electron. At the same time the two types of MC samples are not exactly complementary: the simulation of QED final state radiation (FSR) from the leptons in the V +jets sample generates a fragmentation contribution also contained in the FSR-like diagrams of the V γ process. It is obvious that this overlap has to be removed before the samples can be used for background estimation.
The overlap removal is a conceptually straightforward requirement, which is complicated by two facts. The QED FSR photons in the V +jets simulation are produced at the hadron level and can thus not simply be subjected to parton-level cuts matching the ones in the V γ simulation. Furthermore the photon cuts in a multi-jet merged sample of V γ+jets require an isolation of the photon with respect to partons from the multijet matrix elements. This constraint has to be respected when defining the complementary cuts for the V +jets sample.
An implementation of such an overlap removal at the event generation level is discussed in this section using the example of V = Z.

Implementation of overlap removal
In order to combine Z and Zγ events 4 the phase space is split into two regions. The Zγ process includes photons directly in the matrix elements. The phase space of this region should be as large as possible but is limited since the matrix elements diverge when the photon is soft or collinear either to a massless lepton or quark 5 . By contrast, photons generated by YFS in Z events do not have these limitations, but can not describe initial state radiation which usually gives most of the contribution to hard photons.
In principle, the phase space slicing is defined by three components. First, a p γ ⊥ cut, secondly a lepton photon isolation and finally a photon hadron isolation. These cuts exclude a region where collinear or soft divergences are present and no fixed order calculation in QCD is possible. Thus, an Zγ event has to pass all these cuts while an Z event has to fail at least one of them.
In case of Zγ events, there are already cuts at matrix element level present and it would be desirable to use them directly for the overlap removal. Unfortunately, this is not possible since the generation of additional final state photons via YFS happens technically after the parton shower. The shower can shift the kinematics of all particles, thus cutting once before and once after the shower would result in a mismatch. As a consequence, the slicing cuts are applied to both Z and Zγ events at hadron level.
Hadron level cuts are not supported by Sherpa outof-the-box. For this study, a support for custom modules was implemented which makes it possible to veto events at the hadron-level very flexibly. This feature will be available within the next Sherpa release.
Special care has to be taken when selecting the photon and leptons which take part in the slicing procedure. Non-prompt leptons and photons can easily be produced by the decay of hadrons and there is no requirement that these further particles are softer than the particles coming directly from the hard interaction or the YFS algorithm. However, it has to be guaranteed that all divergences which are present at matrix element level are covered by the slicing cuts since otherwise the result would still depend on the matrix element level cuts.
As a consequence, the hardest photon which does not come from a hadron decay is used for the definition of the slicing variables. The photon-lepton isolation is applied only to prompt leptons. For the hadronic isolation all particles excluding the prompt leptons and the photon are taken into account.
In the following, isolated photons are required to fulfil p ⊥ > 10 GeV, a photon-lepton isolation of ∆R > 0.4 and the hadronic isolation using a smooth cone isolation with R = 0.4, n = 1 and = 0.5.

Results
The validation of the overlap removal algorithm proceeds with analyses in three different phase space regions. Both the Zγ and Z prediction should not be altered by the overlap removal in their regions of validity. The first condition is checked using the Zγ analysis introduced in Section 3.2. It covers the explicit Zγ phase space and defines region I. By contrast, the Z phase space is probed by the default inclusive Z analysis provided by Rivet. Here, a reconstructed Z boson with an invariant mass between 65 and 115 GeV is required. The leptons are dressed with all photons having an angular distance of 0.2 or smaller. These cuts define the phase space region II.
In addition, a special region III is defined where final state radiation contributions via YFS are supposed to give similar contributions as the direct production from matrix elements. Having such a region it is directly possible to study the interplay between both components of the overlap removal and compare their sum with a pure YFS or direct sample.
A region dominated by final state radiation is defined by requiring the lepton pair to have an invariant mass below the Z peak, 30 GeV < m ll < 87 GeV. An event is accepted, if both leading leptons have electron flavour but opposite charge. As photon candidate the leading photon (p γ ⊥ > 5 GeV) is chosen, it has to be isolated from the selected leptons by requiring ∆R γ,e ± > 0.05. In addition, the total energy of all remaining par-ticles (excluding the selected leptons) in a cone with ∆R < 0.4 around the photon axis has to be less than 0.5 · E γ ⊥ . Details of all analyses are summarized in Table 3. All cuts and observables are implemented as a user module using the Rivet framework.
A comparison is performed using four sets of separately generated event samples; pure Zγ, pure Z, sliced Z and sliced Zγ. The latter two are summed up to give the total prediction after overlap removal. YFS is set active for all samples.
For event generation, the matrix element level cuts are selected to be more inclusive than the analysis cuts. For the generation of the sliced direct part, the matrix level cuts are additionally chosen to be more inclusive than the phase space slicing parameters. All events are generated with Q CUT = 30 GeV and up to 3 jets at leading-order accuracy 6 . The slicing parameters which are used for this test are summarized in Table 4.
In Figure 12, the inclusive jet multiplicity and the inclusive E γ ⊥ spectrum for the Zγ phase space (region I) are shown. In these and all further plots the inclusive Z, the direct Zγ and the summed overlap-removed predictions are shown, together with the corresponding components in the overlap removal. For better readability the statistical uncertainties of the latter have been omitted.
Here, the overlap removal result is in very good agreement with the pure Zγ prediction in both plots. The dominating contribution is the Zγ component, giving about 90% of the cross section for low p ⊥ and almost 100% if p ⊥ exceeds 80 GeV.
By contrast, the inclusive Z phase space (region II) is dominated by the Z component. The Z mass and Z p ⊥ distributions are shown in Figure 13. The overlap removal result is in excellent agreement with the pure Z prediction. The only region of phase space where the direct component of the overlap removal is sizeable is the low mass region of less than 85 GeV, there the direct component gives around 10% of the cross section.
Despite the good agreement between overlap removal and the respective reference, one might want to construct the overlap removal to require only the Zγ component when looking at Zγ analyses. This would require to take into account the analysis cuts for slicing the phase space and is thus not possible in a generic sample. An overlap-removed contribution with inclusive Z production allows more flexibility as needed in general purpose experiments.
Finally, in Figure 14 and 15 some observables of the FSR dominated phase space (region III) are shown. In     this region neither the pure Z nor the pure Zγ predictions are guaranteed to give an accurate result. The former one includes only final state radiation and will therefore miss contributions especially in the high E ⊥ region. By contrast, the latter one includes all contributions at a fixed order but leaves the perturbative region if the photon is soft or very close to the lepton. Thus, both these predictions can only be interpreted as lower and upper bounds for the overlap removal in the context of this validation. In Figure 14, different E γ ⊥ spectra are shown. The first subplot covers the whole phase space while the two remaining ones cover only the regions where the photon is either very close to (0.05 < ∆R < 0.5) or separated from (0.5 < ∆R < 3) the closest lepton. While in the inclusive plot both components of the overlap removal procedure give a very similar contribution if E γ ⊥ exceeds the slicing cut, the two remaining plots reveal the nature of the overlap removal procedure much better. Whereas the low ∆R region is entirely dominated by the Z component, the high ∆R region is dominated by the Zγ component as soon as the cut off is exceeded. Figure 15 shows the azimuthal angle between the photon and its closest lepton and the invariant mass of the photon and both leptons. Both observables show an interesting behaviour. As expected, the cross section of the combined overlap removal result always interpolates between the cross section of the pure Zγ and the pure Z sample. When looking at the azimuthal distance, the Z component dominates at lower and the Zγ component at higher values, while both components are equal in a large region of phase space. The invariant mass spectrum is dominated at lower values (< 105 GeV) by the Z component and at higher values by Zγ.

Conclusions
Precise Standard Model predictions for Zγ+jets production are crucial for the search for new particles or anomalous couplings in measurements of this final state at the LHC.
With the presented simulation within the Sherpa framework using the MEPS@NLO algorithm we provide a simulation which is at the same time precise and realistic: NLO QCD corrections for the Zγ and Zγ+jet processes are included and reduce the uncertainties in relevant observables significantly. At the same time, the matching and merging with the parton shower allows a realistic simulation of the full final state at the hadron level, and the inclusion of all off-shell effects allows to place realistic experimental cuts on the prompt leptons without approximations.
Comparing to data from experimental measurements at √ s = 8 TeV we find very good agreement. On that basis we make predictions at √ s = 13 TeV and identify the dominant theoretical uncertainties and the phase space regions affected by them.
To further the application of these precise Zγ+jets predictions in experiments we also demonstrate how they can be combined with event generator predictions for Z+jets including final-state photon radiation. As a validation we introduce a number of cross checks based on different phase space regions which can be repeated for any specific application of such samples in the experiments.