Reweighting QCD matrix-element and parton-shower calculations

We present the implementation and validation of the techniques used to efficiently evaluate parametric and perturbative theoretical uncertainties in matrix-element plus parton-shower simulations within the Sherpa event-generator framework. By tracing the full $\alpha_s$ and PDF dependences, including the parton-shower component, as well as the fixed-order scale uncertainties, we compute variational event weights on-the-fly, thereby greatly reducing the computational costs to obtain theoretical-uncertainty estimates.


Introduction
The first operational run of the LHC collider during the years 2009-2013 was a tremendous success, clearly culminating in the announcement of the discovery of a Higgs-boson candidate by the ATLAS and CMS collaborations in July 2012 [1,2]. Through a large number of experimental analyses, focusing on a variety of final states and observables, the LHC experiments (re)established and underpinned to an unprecedented level of accuracy the validity of the Standard Model of particle physics (SM) [3].
When comparing theoretical predictions with actual collider data, Monte-Carlo event generators prove to be an indispensable tool. In particular parton-shower Monte-Carlo programmes like Herwig [4,5], Pythia [6] and Sherpa [7,8] provide simulations at the level of exclusive particle-level final states [9]. The cornerstones of these generators are their implementations of QCD parton-shower algorithms and their modelling of the non-perturbative parton-to-hadron fragmentation process. With the advent of sophisticated techniques to combine parton-shower simulations with exact higher-order QCD calculations at leading [10,11], next-to-leading [12,13] and even next-to-next-toleading order [14,15,16,17], Monte-Carlo simulations have developed into high-precision tools, encapsulating the best of our current knowledge of perturbative QCD.
With these simulations being widely used for making SM predictions, e.g. of the background expectation in searches for New Physics or the detailed properties of Higgsboson production final states, a comprehensive and efficient evaluation of associated theoretical uncertainties is of utmost importance. A comprehensive list of sources for generator uncertainties has been quoted in [18]. Following the categories identified there, when focusing on systematics related to the perturbative phases of event evolution, the following uncertainties might be distinguished: • Parametric uncertainties reflecting the dependence of the prediction on input parameters such as couplings, particle masses or the parton-density functions (PDFs).
• Perturbative uncertainties originating from the fact that perturbation theory is used in making predictions, to fixed-order in the matrix elements and resummed to all-orders with a certain logarithmic accuracy in the showers, thereby, however, neglecting higher-order contributions. Similarly, the use of the large-N c approximation in the showers belongs in this category.
• Algorithmic uncertainties corresponding to the actual choices made in the implementation of the shower algorithm, i.e. for the evolution variable, the inclusion of non-singular terms in the splitting functions, or the employed matching/merging prescription. Per construction, for sensible choices, these systematics also correspond to higher-order perturbative corrections, but might be addressed separately.
In addition to the listed categories, generically non-perturbative effects such as hadronisation or the underlying event are described through phenomenological models that feature various generator-specific choices and parameters, typically subject to tuning against experimental data, see for instance [19,20]. This publication focuses on the efficient evaluation of parametric and (some) perturbative uncertainties in matrix-element plus parton-shower simulations within the Sherpa event-generator framework. We present a comprehensive approach to fully trace the α s and PDF dependences in the matrix-element and parton-shower components of particle-level Sherpa simulations in leading- [21] and next-to-leading [22] order merged calculations based on the Sherpa dipole-shower implementation [23]. Furthermore, we provide the means to quickly evaluate the renormalisation-and factorisation-scale dependence of the fixed-order matrix-element contributions. Our approach is based on event-wise reweighting and allows us to provide with a single generator run a set of variational event weights corresponding to the predefined parameter and scale variations, that would otherwise have to be determined through dedicated re-evaluations. The alternative event weights can either be accessed through the output of a HepMC event record [24], or directly passed via the internal interface of Sherpa to the Rivet analysis framework [25].
The systematics of leading-order parton-shower simulations with Herwig 7 have recently been discussed in [18], a corresponding reweighting procedure has been presented in [26]. A similar reweighting implementation for the Pythia 8 parton shower has also appeared recently [27]. A discussion of uncertainty estimates for the Vincia shower model can be found in [28,29,30]. A comprehensive comparison of various generators is presented in [31]. The impact of PDFs in parton-shower simulations has been discussed in [32,33].
Our paper is organised as follows. In Sec. 2 we review the dependence structure of leading-order (LO) and next-to-leading-order QCD calculations on α s , the PDFs and the renormalisation and factorisation scales, and introduce the reweighting approach. In Sec. 3 we extend this to parton-shower simulations and in particular the algorithm employed in the Sherpa framework. In Sec. 4 we present the generalisation of the reweighting approach to multijet-merged calculations, based on leading and next-to-leading-order matrix elements matched to the parton shower. Our conclusions are summarised in Sec. 6. In App. A we present CPU time measurements that assess the reduction in computational time when the reweighting is used. The technical details on enabling and accessing the variations considered in Sherpa runs are listed in App. B.
Note, while fixed-order reweighting is already available with Sherpa-2.2, the general reweighting implementation described here, including parton showers and multijet merging, will be part of the next release, i.e. Sherpa-2.3.

Reweighting fixed-order calculations
In order to re-evaluate a QCD cross-section calculation for a new choice of input parameters, i.e. α s , PDFs or renormalisation and factorisation scales, it is necessary to understand and trace-out its respective dependences. This is a rather easy task at leading-order (LO) but is already more involved when considering next-to-leading order (NLO) calculations in a given subtraction scheme. However, these decompositions have been presented for Catani-Seymour dipole subtraction and the FKS subtraction formalism in [34,35].
In this section, we briefly review the dependence structure and discuss the corresponding reweighting equations for LO and Catani-Seymour subtracted NLO calculations within the Sherpa framework. With this paragraph we also introduce the notation used in the later sections, which explore the reweighting of more intricate QCD calculations, involving QCD parton showers and merging different final-state multiplicity processes.

The leading-order case
A LO parton-level calculation of some observable or measurement function of the finalstate momenta O is based on Born matrix elements B of O (α n s ). It exhibits explicit dependences on the PDFs f = f a (x, µ 2 F ), the running strong coupling α s = α s (µ 2 R ), the renormalisation scale µ R and the factorisation scale µ F : n trial,i , n trial denoting the number of attempts to generate an accepted event configuration, and Therein, the B contains all couplings, symmetry and flux factors, and PDFs, whereas B has the PDFs, here for assumed two incoming parton flavours a and b, and the strong coupling stripped off. Note that we have suppressed the event index i here. It is understood that B depends on the event kinematics and that µ R and µ F can be chosen dynamically, i.e. in a momentum (and flavour) dependent way. Changing the input parameters µ R →μ R , µ F →μ F , and the input functions f →f , α s →α s results in From eq. (2.3) we conclude that for PDF reweighting it is necessary to know the x a,b values of the event.
For an unweighted event generation, the event weights are uniform initially, i.e. B(Φ B ; α s , f ; µ R , µ F ) = w norm , eq. (2.1) thus simplifies to (2.4) Scale and parameter variations then work the very same way as for weighted events. Applying eq. (2.3) then, however, leads to a broader weight distribution and eq. (2.1) has to be used again. Partially unweighted events can be treated on the same footing. These conclusions hold irrespective of the type of event generation whenever (partially) unweighted event generation is possible, i.e. when the weight distribution is bounded from above and below. We therefore will not comment further on it.

The next-to-leading-order case
A full NLO parton-level calculation including real-emission and one-loop corrections of O (α n+1 s ) based in Catani-Seymour dipole subtraction [36,37] has the following structure where the new parts have the following dependences (2.6) Therein, VI combines the renormalised one-loop matrix element with the I-operator of the Catani-Seymour subtraction scheme. This operator gives the flavour-diagonal endpoint contribution of the integrated subtraction terms. VI is thus separately infrared finite and exhibits a common transformation behaviour. Thus, for α →α s , f →f , µ R →μ R with α s -and PDF-independent coefficients c (i) R and l R = log(μ 2 R /µ 2 R ). Again, VI is stripped of all coupling and PDF factors.
The KP-terms are defined as the remainders of the integrated dipole subtraction terms, containing all flavour changing and x a/b -dependent pieces, combined with the collinear counterterms. Here, x a/b are the ratios of the partonic momentum fractions in the respective dipole before and after radiation. Again, this combination is separately infrared finite and transforms as one unit. When evaluated for the modified set of input parameters, they read for a, b = {q, g}, respectively. Thereby, the sum over q includes all light-quark flavours, corresponding to all potential quarks emitting a gluon. We note that in order to obtain the reweighted expressions for the VI and KP contributions the additional book-keeping of the c F,a/b (altogether 18) 1 coefficients is required [34]. Due to its composite structure, the KP-terms do not possess a coupling-and PDF-stripped version 1 The two parameters c KP . Nonetheless, we formally introduce a still PDF-dependent version KP in eq. (2.8) for reference in later sections.
The remaining pieces of eq. (2.5) are the Born matrix element B, the real emission contribution R and the differential dipole subtraction terms D S,j . The latter defines an underlying Born configuration Φ B,j through its dipole-dependent phase-space map, employing the phase-space factorisation Φ R = Φ B,j · Φ j 1 . While the transformation of B under the exchange of input parameters was detailed in eq. (2.3), the transformation of R and the D S,j contributions works identically, merely having to adjust the power of the strong-coupling factor.

Validation
The reweighting approach outlined above has been implemented in the Sherpa framework for the two matrix-element generators Amegic [38] and Comix [39,40] in conjunction with the corresponding Catani-Seymour dipole-subtraction implementation [41]. The required decomposition of virtual amplitudes is generic and can be used for matrix elements from BlackHat [42,34], OpenLoops [43], GoSam [44], Njet [45], the internal library of simple 2 → 2 processes, or, via the BLHA interface [46].
Here we shall present the validation of the reweighting approach in particular of NLO QCD event samples. For that purpose we consider W-boson production in 13 TeV protonproton collisions at NLO QCD, and focus on the transverse-momentum distribution for the W and the lepton it decays to. In Fig. 1, the scale, α s and PDF uncertainty bands for the W p ⊥ and the lepton p ⊥ distributions are presented. All three bands have been produced for both observables using the internal reweighting of Sherpa from a single event generation run using a scale choice that has been motivated in [47]. For the PDFs the NNPDF 3.0 NLO set [48] has been used with α s (m 2 Z ) = 0.118. The running of α s (µ 2 R ) is calculated within Sherpa using its renormalisation group equation at NLO with parton thresholds as given by the PDF.
The treatment of partonic thresholds deserves a short discussion. While any flavour thresholds in the running of α s do not present any challenges to the reweighting algorithm as α s (µ 2 ) > 0 for all µ 2 > 0 and any loop order, this is different for the PDFs, where crossing a parton threshold results in a vanishing PDF for that flavour. Hence, the cross section component of the given partonic channel may be zero if no other non-zero contribution exists. Such an event will be discarded and, thus, cannot be reweighted. If now the respective parton threshold of the target PDF is smaller than the target factorisation scale while the one of the nominal PDF is larger than the nominal factorisation scale we are in a region of phase space where the reweighting must fail to reproduce a dedicated calculation. This could be remedied by storing events as well which vanish solely due to crossing PDF thresholds. However, as only observables sensitive to on-threshold production of light quarks (typically bottom quarks) are susceptible to these effects, they are of little relevance to the vast majority of LHC observables. 2 For the scale uncertainty band we employ a 7-point scale variation for µ R and µ F : Both scales are varied independently by factors of 1 ⁄2 and 2, omitting the variations with ratios of 4 between the two scales. The uncertainty is then taken as the envelope of all variations. The α s uncertainty band is generated by varying the numerical value of the starting point of the running coupling, α s (m 2 Z ), to the following five values: 0.115, 0.117, 0.118, 0.119 and 0.121. Note that this variation of α s should also enter the PDF fit, and hence the PDFs are varied consistently. This is expected to extenuate the effect of the α s variation in most cases, as the PDF of the varied α s is still fitted to describe the same data as the PDF of the nominal α s . This consistent α s +PDF variation is also part of the PDF4LHC recommendations for LHC Run II [49]. The envelope of these α s +PDF variations is taken as the respective uncertainty. The pure PDF uncertainty estimate is generated using the average and the standard deviation over the 100 PDF replicas provided by the NNPDF3.0 set (at a fixed value of α s = 0.118). This corresponds to the 68 % confidence level. This set-up is repeated for later reference in Table 1, along with a CT14 PDF variant, which is used in later studies.
Comparing the uncertainties for the W p ⊥ , we observe that the scale uncertainties are the largest, with relative deviations of O (10 %). The relative deviations related to the PDF and the strong coupling do not exceed ∼ 3 %. The scale uncertainty exhibits a minimum for 100 GeV < p W ⊥ < 200 GeV. The reason is that the variations of µ F alone cross the central value prediction in this range, such that only the µ R variation contributes to the overall scale uncertainty here.
Note that p W ⊥ = 0 at O (α 0 s ), and therefore only real-emission events contribute to the distribution. Hence, the observable is only described to leading-order. We introduce it here as a reference for our later validations including the parton-shower, which use this      at the LHC with independent variations of µ F,R (green), α s (red) and the PDF (blue). In the right-hand panels, the individual uncertainty bands, calculated via an on-the-fly reweighting, are compared to uncertainty bands from dedicated calculations (yellow). They are found to be equal.
observable. For the current validation, we complement the discussion of the W transverse momentum with the one of the lepton it decays to, as the region below m W /2 is already filled at O (α 0 s ), and therefore we have in part a true next-to-leading description for this observable. In fact, the scale uncertainties are much larger in that region, especially towards the m W /2 threshold, and at the lepton p ⊥ cut at 25 GeV. This gives a more realistic picture of the perturbative uncertainties than in the leading-order region above the threshold.
The small panels on the right of Fig. 1 compare the uncertainty bands calculated using the reweighting approach to uncertainty bands where dedicated calculations have been done for each variation. We observe that all bands overlap perfectly for both observables. This is because the reweighting as presented above is exact and for all runs the same phase-space points could be used: The reweighted and the dedicated predictions for each variation are therefore equal, and so are the uncertainty bands. 3

Reweighting parton-shower calculations
If parton-showering is added to a LO calculation, the value of the observable is not evaluated at Φ B any longer, but at PS(Φ B ), which denotes the phase-space point after showering. Applying this modification to eq. (2.1) yields Therefore the reweighting for B does not need to be altered, but the parton-shower emissions depend on the PDF, the strong coupling, their respective scale prefactors k αs and k f (detailed below) and the starting scale µ Q , i.e.
In order to reweight the parton-shower emissions, we first need to identify its exact dependence structure. Schematically, it acts on the phase-space element in the following way where the Sudakov form factor of the n-parton state, ∆ n , and its splitting kernel K n have been introduced. While the first term describes the no-emission probability between the starting scale t and the infrared cut-off t IR and therefore does not change the phase-space element, the second term describes the emission of a parton at scale t in the configuration dΦ 1 = dt dz dφ J(t, z) (the integration boundaries are to be understood Table 2: Definition of the evolution and splitting variables for each dipole type. The fifth column lists the splitting process as seen from the Born process, c and c refer to the flavour of the initial state before and after the splitting process, respectively. The variables y ij,k ,z i , x ij,a , x jk,a , x j,ab , u j and v j are defined in [36,37,23]. in this decomposition), leading to a configuration dΦ n+1 = dΦ n · dΦ 1 . The Jacobian J is not relevant to the discussion here and is subsequently absorbed in the splitting kernel K n . As the emissions are ordered in t, the Sudakov form factor in the second term ensures that the current emission is the hardest after starting the evolution at t . Additional emissions may occur at smaller t and are not resolved at this stage -they are described by the parton shower acting on the newly produced state Φ n+1 with the new starting scale t. In eq. (3.3) the dependences on α s , the PDFs, and their respective scale prefactors k αs and k f have been omitted for brevity. They directly carry over to the splitting kernel and the Sudakov form factor, according to When considering parton-shower emissions off NLO QCD matrix elements special emphasis has to be given to the first emission as described in Sec. 3.3 below.

Parton-shower dependence structure
The default parton shower of Sherpa, dubbed CSShower [23], is based on Catani-Seymour dipole factorisation [36,37]. Each branching of an emitter parton into two daughters is witnessed by a spectator parton, which takes the recoil, and ensures that on-shell states are transferred into on-shell states and energy-momentum conservation is respected simultaneously. The emitter and spectator partons reside either in the initial-state (I) or final-state (F), such that four dipole types need to be distinguished: II, IF, FI and FF. In this notation, the first letter refers to emitter, and the second to the spectator parton. The no-branching probabilities are given by the four corresponding Sudakov form factors They share the common form wherein the kinematics of the splitting are given by the default choice for t = Q 2 y z(1 − z) in the massless case while the K ij,k (t, z) denote the coupling and PDF stripped splitting kernels incorporating the remaining pieces of the K ij,k and the Jacobian J of the phasespace parametrisation. The precise definitions of the variables for each dipole type are given in Table 2. It directly follows that for FF-type dipole splittings the ratio of PDFs is simply unity. Eq. (3.6) further details the dependence on the α s and PDF scale factors k αs and k f . These multiplicative factors as well as their variations are assumed to be of order one, such that they do not induce spurious large logarithms. The generalisation to the massive case is straightforward and only involves generalised definitions of t, x, y and z, cf. [23].

Reweighting trial emissions
To numerically integrate Sudakov form factors typically the Sudakov Veto Algorithm is used [50,51,52,53,54,55]. Therein the integrands K found in the Sudakov form factors are replaced with integrable overestimatesK. This is balanced by only accepting a proposed emission with probability P acc = K/K. A multiplicative factor in K is therefore equivalent to a multiplicative factor in P acc [52]. This observation is for example used to apply matrix-element corrections [54], where the splitting kernels are replaced with a real-emission-like kernel R/B. This is done a-posteriori, i.e. the event weight is multiplied by (R/B)/K, the emission itself is unchanged. The same method is also used in the Vincia parton shower to calculate uncertainty variations for different scales, finite terms of the antenna functions, ordering parameters and sub-leading colour corrections [28].
Here we employ this technique to account for variations of the strong-coupling parameter and the PDFs in the shower evolution of LO and NLO QCD matrix elements.
As has been laid out in the previous section, the emission kernels K depend linearly on α s and on a ratio of parton densities f c (η c /x, k f t)/f c (η c , k f t). A change of PDFs f →f , the strong coupling α s →α s and the scale prefactors entering both, i.e. k αs →k αs and k f →k f , is equivalent to modifying the emission probability accordingly 4 : where the scale dependence and the definition of η c and x can be read off the Sudakov form factors given in eq. (3.6) and Table 2. In case of FF dipoles eq. (3.7) simplifies significantly as the ratios of PDF factors reduces to unity. It further follows, that the event weight for each accepted emission needs to be multiplied by the corresponding factor q acc in order to incorporate the new choice of α s , PDFs and the scales they are evaluated at. Accordingly, the probability to reject an emission is changed to Consequently, for each rejected emission the event weight receives a corrective weight of q rej . Proofs that this treatment indeed results in the correct Sudakov form factors can be found in [52,27,26].

Next-to-leading-order matching
To match NLO QCD parton-level calculations with subsequent parton-shower evolution Sherpa employs a variant of the original Mc@Nlo algorithm presented in [12], referred to as S-Mc@Nlo [54]. Schematically, such a S-Mc@Nlo calculation has the following structure: Here the real-emission contribution R of the NLO calculation has effectively been split into an infrared-singular (soft) and an infrared-regular (hard) part, the resummation kernel D A and the finite hard remainder H A , respectively, such that R = D A +H A [57,54].
The B-function has the following explicit parameter dependences (3.10) From the perspective of parameter reweighting, the resummation kernel D A behaves the same way as the subtraction term D S . In fact, in our reweighting implementation the (D A − D S ) contribution is treated as a single term, as indicated. It is only to note that their PDFs are evaluated at the partonic momentum fraction x a/b,j and external flavours a j and b j of their Φ B · Φ j 1 phase-space configuration rather than those of Φ B . The other parts of the B-function can then be reweighted as described in Sec. 2.2, leading to The H A -function then transforms as wherein each subtraction term D A,j has its own scales µ R,j , µ F,j defined on its underlying Born configuration Φ B,j . Writing eq. (3.9) as a Monte-Carlo sum over events with B-like and R-like structure, which are conventionally called S and H events in Mc@Nlo calculations, and with N = N S + N H , we obtain Thus, under µ R →μ R , µ F →μ F , α s →α s and f →f both B of the S-events and H A of the H-events transform as composite objects in terms of their constituents, as defined above. This leaves the S-Mc@Nlo parton shower, PS NloPs , defined through (3.14) It differs from the usual PS of eq. (3.3) with respect to the splitting kernel for the first emission and the associated definition of the Sudakov form factor, cf. [54,58]. However, for the purpose of reweighting, all trial emissions can be treated in the same way as in the standard parton shower as the parameter dependences are identical.

Validation
To validate the reweighting of scale and parameter dependences in CSShower and S-Mc@Nlo calculations within the Sherpa framework we perform closure tests between reweighting results and dedicated simulations.
Our implementation allows to constrain the maximum number of reweighted shower emissions per event. For a pure leading-order parton-shower run or H-like events in S-Mc@Nlo calculations this amounts to setting n PS ∈ {0, 1, 2, . . . , ∞}. When considering S-Mc@Nlo simulations in addition the parameter n NloPs ∈ {0, 1} can be used to disable the reweighting of the O(α s ) emission for S-events.
Of course the reweighting result will only coincide with a dedicated calculation if all emissions are reweighted, i.e. n NloPs = 1 and n PS = ∞. However, by subsequently enabling the reweighting of more and more emissions the relevance of their dependences for the determination of the full uncertainty can be studied. A finite value of n PS can also be useful in production, if the effect of reweighting higher-order emissions becomes negligible. The reduced amount of reweighting per event then allows for a faster event generation. An additional benefit would be that rare high-multiplicity shower histories do not spoil the statistical convergence of the reweighted result, even if their exact kinematics might be irrelevant for the studied observable.  The final-state only case: Thrust in e + e − → qq events To validate LoPs reweighting, we consider two observables, which are complementary in their sensitivity to parton-shower emissions. At first, we consider the event-shape variable thrust T [59] in hadronic events in e + e − -collisions at √ s = 91.2 GeV. In this case QCD emissions are restricted to the final-state. Accordingly, there appear no PDF factors in the shower reweighting, cf. eq. (3.7), and thus no factorisation scale dependence. Moreover, as we consider the leading-order matrix element for e + e − → qq only, the renormalisation scale is also absent in the hard-process component. Therefore, we can concentrate on the pure α s uncertainty in the parton shower here. Leaving the perturbative order of its running invariant it is defined by its value at the input scale m Z .
In Fig. 2, we compare α s uncertainty bands generated by reweighting the nominal prediction with the one generated by dedicated predictions for each variation. As in Fig. 1, the uncertainty band is defined as the envelope over the distributions with different α s (m 2 Z ) input values. The nominal value is taken as α s (m 2 Z ) = 0.120, and its up/down variations are 0.128 and 0.108, respectively. Reweighting bands are presented for n PS = 1, 4, 8, ∞. The n PS = 1 band underestimates the uncertainty, especially for T ≤ 2/3, where multiple hard emissions are required, and for T ≈ 1, the region sensitive to multiple soft emissions. For n PS = 4, the uncertainty is underestimated only for bins with T ≤ 2/3, and less so than for the n PS = 1 case. The difference between the two choices of n PS = 8 and ∞ is merely statistical and both reproduce the dedicated result very accurately. However, for low values of T , the statistical fluctuations of the reweighting results with higher n PS grow larger, corresponding to a widening of the distribution of reweighting factors. Low values of the thrust observable correspond to the emission of several hard partons, which is less probable in the parton-shower approximation, and more appropriately modelled in multijet-merged calculations, cf. Sec. 4. In this phase-space region it is difficult for the reweighting to compensate the multitude of accepted soft emissions off these hard legs, that turn unstable for P acc → 1, with rejected ones, cf. eq. (3.8). This issue can be addressed by introducing a prefactor for the over-estimator functionK in the reweighting runs, to ensure that P acc does not approach 1, cf. [52,27,26]. This renders the Sudakov Veto Algorithm somewhat less efficient, but is shown to reduce statistical fluctuations in the reweighting.  Table 1. We now use the CT14nlo PDF set, which uses a Hessian error representation at a 90 % confidence level [60]. 5 Therefore the PDF error band will be larger than before, as it now corresponds to nearly two standard deviations instead of only one.
Considering a hadronic environment, initial-state emissions are present, which means that our reweighting factors now include PDF double ratios. In Fig. 3, we compare LoPs uncertainty bands for scale, α s and PDF variations, including comparisons between reweighted and dedicated predictions for a varying maximum number of reweighted shower emissions n PS . Before discussing the bands, we observe that the tail of the p W ⊥ spectrum is not populated, in particular in comparison to Fig. 1. This is expected, as the LO configuration is restricted to p W ⊥ = 0, such that all other bins are filled through recoils against parton-shower emissions only. However, the phase space of parton-shower emissions is restricted to the soft region, and therefore the W boson can not build up a large recoil.   other maximum numbers of reweighted emissions n PS Figure 3: The same as in Fig. 1, but for LO + parton-shower (PS) generation. The uncertainty bands are calculated by reweighting the ME and up to n PS shower emissions. In the upper four plots, n PS = 3. In the lower plots, n PS is varied for comparison. The scale uncertainties do not change with n PS and are therefore not repeated.
We now turn to the scale uncertainty band-which is entirely due to factorisation scale variations, because the LO matrix element is independent of α s , and therefore the band underestimates the perturbative uncertainty. We also observe that the band is nearly flat. As we vary only the scales of the matrix-element calculation, the constant spread corresponds to the factorisation-scale uncertainty of the Born configuration at p W ⊥ = 0, merely propagating to higher p W ⊥ bins through the parton shower, which is unaware of the scale variations. In the matrix-element reweighting, we can guarantee the same phase-space points as in the dedicated run, such that we see perfect agreement between dedicated and reweighted predictions. We therefore omit comparisons for different n PS for the scale-uncertainty band.
Looking at the α s uncertainty band, we can see that the envelope constricts at the position of the peak of the distributions. This reflects that the variation of α s shifts the position of the peak, such that variations that are below the nominal distribution on the left side of the peak, are exceeding the nominal distribution on the right side, and vice versa. Comparing the reweighted prediction to the dedicated one, we find a flat band for n PS = 0, corresponding to restricting the reweighting to the fixed-order matrix element. As the LO calculation is independent of α s , this only reflects the change of the PDFs, which are fitted to α s (m Z ). The reproduction of the shape of the α s uncertainty improves a lot when reweighting up to one emission (n PS = 1), and slightly more when adding another emission on top (n PS = 2).
For the PDF uncertainty, we see that the reweighting with n PS = 0 underestimates it by at least 1-5 % for small transverse momenta, and overestimate it around the W mass. As for the α s uncertainty, this improves for n PS = 1, 2.
The last depicted step, i.e. n PS = 3, on the other hand, does not contribute further to the reproduction of the α s and PDF uncertainties. No significant differences with respect to the n PS = 2 case is observed. It can be concluded that it is sufficient to include up to two emissions to reproduce the uncertainty bands for this observable. This is to be expected, as the gauge boson recoils against the shower emissions and is therefore mostly affected by the few hardest branchings. These mainly originate from the incoming hard virtual partons, so the generally softer final-state emissions barely contribute. Although we do not reproduce this here, we confirmed this by entirely disabling final-state emissions, which showed no effect on the results.
In Fig. 4 we present the validation of NloPs predictions for the W-boson transversemomentum distribution. Overall, we get a similar picture as in the LoPs case. The main differences are the increased high-p W ⊥ reach and the significantly smaller scale uncertainties in the low p W ⊥ range, a consequence of including the complete set of O(α s ) corrections to the production process. For large p W ⊥ , the uncertainty increases again, because we fall back to a LO description again: Only the 2 → 3 matrix element still contributes.
To assess the quality of the reweighting, we consider again different settings for the parameters n NloPs and n PS . Assuming n NloPs = n PS = 0, only the scale variations of the hard process are considered and the parton-shower contribution to the O(α s ) correction is not reweighted. Furthermore, we present results for n NloPs = 1 and n PS = 0, 1, 2.  With these settings, the O(α s ) corrections get properly reweighted, but the number of subsequent shower emissions off the S-and H-like events treated correctly is varied. We observe a saturation for reproducing the dedicated calculations at n NloPs + n PS ≥ 2, with no further improvement when n PS is increased from 1 to 2. This confirms the findings made when considering the LoPs setup in Fig. 3: The gauge-boson transverse-momentum distribution is dominated by the few hardest emissions.

Reweighting multijet-merged calculations
In this section we address the reweighting of multijet-merged event generation runs. These approaches allow to combine LO or NLO QCD matrix elements of different multiplicity dressed with parton showers into inclusive samples. Accordingly, the production of jets associating a given core process can be modelled through exact matrix elements rather than relying on the logarithmic approximation of the parton shower only. In particular, when considering hard jet kinematics or angular correlations such techniques prove to be indispensable to properly describe experimental observations, see for instance [61,62,63,64].
To first approximation the reweighting as described in the previous sections can be used without change, only that the perturbative order p is no longer a constant across the sample, but varies for each event, corresponding to the considered matrix-element parton multiplicity. However, there are also new algorithm-specific intricacies which complicate the dependence on the input parameters and need to be dealt with to allow for a consistent reweighting. The LO and NLO merging techniques employed within the Sherpa framework are presented in [21] and [65,66], respectively. They rely on the reconstruction of parton-shower histories for multi-parton amplitudes that set the parton-shower initial conditions for their subsequent evolution. This is achieved by running a backward-clustering algorithm that identifies a corresponding core process and calculates hard-parton splitting scales that serve as predetermined shower branchings. In the Sherpa approach the actual parton shower then starts off the reconstructed core process and implements the predetermined hard splittings based on a truncated shower. Furthermore, it is the purpose of the truncated-shower evolution to implement possible Sudakov vetoes for shower emissions above the phase-space separation or merging scale Q cut . It should be emphasised here, parton-shower reweighting is vital when using modified input parameters in order to cancel the Q cut dependence to the accuracy of the parton shower. In case only the hard-process matrix element parameters get reweighted, the dependence on Q cut is cancelled to leading-logarithmic accuracy only, however, residual subleading contributions from the running coupling or the PDF evolution remain [65].
For the reweighting of the truncated-shower Sudakov veto probability the methods described in Sec. 3 can be applied. In what follows we will detail the specifics of the reweighting procedure for LO and NLO multijet-merging runs with Sherpa supplemented by an extensive validation of the implementation.

Preliminaries
Common to the LO and NLO merging techniques used in Sherpa, cf. [21,22,65,66,67], is the separation of the emission phase space into a soft and a hard region, defined through a suitable m-parton measure Q m and a separation criterion Q cut . For each parton configuration Φ m with Q m > Q cut a shower history that represents the event as a core process with subsequent 1 → 2 shower splittings is probabilistically build through backward clustering. The resulting sequence of cluster steps is characterised by tuples {a i , b i , x a,i , x b,i , t i }, recording the (possibly changing) initial-state flavours and momentum fractions as well as the evolution variable of each splitting. We allow for both QCD and EW splitting functions [21,68] to identify such splitting processes and veto recombinations that would lead to the reduction of configurations which are not present in the matrix elements 6 . Figure 5 details possible cluster histories for a given pp → Z + 4 jets configuration, depending on its kinematics, allowing for QCD splittings only (left) or both QCD and EW splittings (right). The sequence {t i } of reconstructed branching scales then may be either ordered or unordered, with an ordered history satisfying t j < t j−1 < . . . < t 1 < t 0 = µ 2 F,core . The recombination probabilities in each clustering step are determined by the forward-splitting probabilities and are therefore dependent on the parton shower and its parameters and choices. This is reflected, step-by-step, in the addition of one factor of α s (when appropriate) at the reconstructed splitting scale, a ratio of PDFs at the reconstructed initial flavours and their momentum fractions, and a Sudakov form factor describing the evolution of each step.
In the Sherpa implementations the α s and PDF factors are added explicitly onto the respective matrix elements and can therefore be reweighted directly. The Sudakov form factor, on the other hand, is implemented through a vetoed truncated parton shower [21,22]. The truncated shower itself, accounting for the possibility of soft parton-shower emissions between subsequent reconstructed hard emissions, i.e. with t m < t < t m−1 but Q < Q cut , can be reweighted with the methods described in Sec. 3. If, however, an emission with Q > Q cut occurs the event is vetoed. Practically, this is accounted for through increasing n trial of the next accepted event by n trial of the vetoed event. Thus, n trial becomes dependent on the parton-shower parameters.
A special remark concerning unordered histories, i.e. histories whose sequence of {t i } has at least one pair t k ≥ t k−1 , is in order. Such histories can be encountered in various configurations, e.g. when the last clustering step produces a splitting scale larger than the nominal starting scale of the core process 7 or the flavour structure only allows further clusterings at scales t k−1 lower than the last identified one t k 8 . As such configurations cannot be generated by a strictly ordered parton shower, for each unordered step neither the accompanying PDF ratio nor Sudakov factor is therefore present in the calculation. More than one unordering in a cluster history of a given event is possible and in fact likely at high multiplicities. PDF ratios and Sudakov factors then of course only occur in the ordered subhistories in between the unorderings. For the sake of clarity and brevity we will omit this case from the discussion of the following subsections. Its implications to the algorithm, and therefore to the reweighting, are straightforward. If ordered histories are enforced, core configurations beyond the standard 2 → 2 processes occur. Independent of the presence and number of unorderings, the renormalisation and factorisation scales µ R and µ F are always set in the way, as will be detailed below.

The leading-order case
We start the discussion with the simplest case, where all matrix elements used in the merging are given at leading order. A LO multijet-merged calculation, with Born matrix elements at O (α n+j s ), containing j additional partons relative to the core process, has the following structure Note that Φ j here denotes the entire final-state phase space of the process, including all particles of the core process. As before, Q j is a suitable infrared-safe distance measure of Φ j . The Θ-function thus realises a minimum separation of Q cut and acts as an infrared regulator. PS vt is the vetoed truncated parton shower derived from eq. (3.3). As the limit in the second line is well defined, it can be transposed with the summation over parton multiplicities. As the ingredient leading-order matrix elements need to incorporate the soft-collinear resummation properties of the parton shower, they have the following parameter dependences: The cluster steps {a i , b i , x a,i , x b,i , t i } denote the identified cluster history of the configuration Φ j , as discussed above. Therein, the a i , b i are the possibly changing initial-state flavours, the x a,i , x b,i their momentum fractions, and the t i are the reconstructed values of the parton-shower evolution variable at each splitting. Together with the α s and PDF scale prefactors k αs , k f of the parton shower, the cluster steps relate B merge j to the scale and PDF stripped Born matrix element B j encountered in Sec. 2.1, (4.3) In this notation, the core scale is t 0 = µ 2 F,core , it is therefore not multiplied by the prefactors of the parton shower. The partonic momentum fractions of the core process are x a,0 , x b,0 .
The scales of each single α s within the cluster history vary, but an effective global renormalisation scale can be defined through where i = 0 if the identified splitting process at branching i is of QCD-type, and 1 otherwise, e = j i=1 i . To consistently vary the µ R scale, we consider variations of the splitting scales t i and the core scale µ R,core on the right-hand side by a common factor, solving for the prefactor of the effective µ R to be used in the matrix-element calculation. Thus, while up to NLO accuracy µ R is varied by the same common factor, the full solution of this procedure results in slightly larger variations of the effective renormalisation scale.
Apart from the Sudakov form factors the soft-collinear structure of the B merge j is now identical to the emission of j partons off a B 0 configuration with the parton shower described in section 3. In case of final-state splittings the ratio of parton distribution functions is simply unity as neither the partonic x a/b,i and x a/b,i−1 nor the initial-state flavours a i , b i and a i−1 , b i−1 differ. In principle, with every ratio of PDFs there is also a ratio of flux factors. However, all such factors cancel except for the outermost ones, corresponding to Φ j and, hence, are regarded as part of B j .
Changing the scales µ R,core →μ R,core , µ F,core →μ F,core , k αs →k αs , k f →k f as well as α s →α s and f →f results in (4.5) The scaleμ 2 R is now calculated from eq. (4.4) usingμ 2 R,core andk αs as input. Eq. (4.5) describes what happens to the matrix-element part of a multijet-merged calculation. This leaves the vetoed truncated shower PS vt . While the truncated and standard shower part is described in section 3, the vetoed shower leads to vetoed events. As vetoed events correspond to events whose weights have been set to zero, their description is equivalent to increasing the number of trials, n trial , by one. Thus, when varying the parameters of the parton shower, also the probabilities of vetoing events are changed. Consequently, n trial acquires a dependence on the parameters of the variation. Thus, now explicitly stating the dependence on the shower's starting scale µ 2 Q and the merging scale Q cut , Thus, n trial corresponds to the survival probability between the unfolding of preexisting splittings when evolving from the n-parton core configuration to the (n + j)-parton configuration. Hence, when changing the parameters of the simulation the truncated showers probability to emit a parton with Q > Q cut must be re-evaluated following the substitutions k αs →k αs , k f →k f , α s →α s and f →f using the methods of Sec. 3. Note, all emissions produced by the truncated shower prior to the one that triggers the veto need to be reweighted as they impact the initial conditions for that emission.

The next-to-leading-order case
The merging of multijet matrix elements at next-to-leading order accuracy proceeds schematically similar as in the leading order case. The input quantity is now the NloPs matched (n + j)-parton configuration, thus are evaluated at the same phase-space point Φ j they share a common cluster history Hence, again suppressing any further Q cut -dependence which is not varied, transforms under the replacements k αs →k αs , k f →k f , α s →α s and f →f in the following way: (4.10) In addition to the transformation properties of the B-function of eq. (3.10), supplemented with the PDF ratios already encountered in the leading-order case, additional terms appear. They subtract the O (α s ) expansion of these ratios, in order to retain the NLO accuracy of the merged calculation. Again, please note t 0 = µ 2 F,core . The same does not hold, however, for the H-events. Their constituent real-emission matrix elements are defined on Φ j+1 , while each subtraction term D A,k has its own projection on a phase-space point Φ k j . Thus, wherein both R merge j and the D merge A,k,j separately transform as the leading-order counterpart B merge j . While the measure Q on Φ j+1 of R merge j is defined to act on the underlying Φ j after the first cluster step where the real emission configuration has been reduced to a Born configuration, it is defined directly on each Φ k j in each D A,k . Infrared safety is guaranteed through the infrared safety of their phase-space maps, the clustering algorithm and the measure Q.
Finally, we consider merging additional multiplicities up to j max described through leading-order matrix elements on top of a next-to-leading order merged calculation with up to j nlo max jets, where j max > j nlo max . The method of choice was outlined in [69,22,65,66,70] and is historically referred to as MeNloPs. Its methodology is defined as 12) with N = N S + N H + N LO . A differential K-factor k j nlo max is a applied to the highermultiplicity leading-order matrix elements in order to facilitate a smooth transition across Q cut . It has the form , (4.13) and therefore moulds the B j nlo max +1 into the same form as the H A,j nlo max it is replacing. The projection Φ j nlo max +1 (Φ j ) for j > j nlo max + 1 is defined through its cluster history, as is Φ m (Φ m+1 ) inside k m itself. When now changing the parameters of the calculation, α s →α s , f →f , µ R →μ R and µ F →μ F , k j nlo max transforms as a composite object in terms of its constituents, cf. Sec. 2 and 3. The scales are set directly by the B merge j process. Of course, in the interest of decreased computational costs one may decide to choose k m ≡ 1 throughout at the cost of larger merging systematics. Similarly, if an electroweak cluster history leads to a changed signature in Φ m+1 , k m ≡ 1 is chosen.

Validation
The reweighting for multijet-merged calculations as discussed in the previous sections has been implemented within Sherpa with the CSShower for leading-order matrix elements (MePs@Lo), next-to-leading order matrix elements (MePs@Nlo) and nextto-leading-order matrix elements with additional leading-order ones on top (MeNloPs). For the validation, we again perform closure tests between reweighted and dedicated predictions for the transverse momentum of the W boson in Figs. 6, 7 and 8. As before, the uncertainty bands are the ones defined in Table 1, with the CT14 PDF set. For all merged calculations, we employ a merging cut of Q cut = 20 GeV.
For the MePs@Lo validation in Fig. 6, we combine LO matrix elements for 0-, 1-and 2-jet multiplicities, obtained from Comix [39]. We can observe that we populate a much larger phase space than for a mere LoPs calculation in terms of p W ⊥ . Below the merging cut (i.e. p W ⊥ 20 GeV), the scale uncertainty band is equal to the one of the LoPs calculation. For higher p W ⊥ , the scale uncertainty increases corresponding to the larger uncertainty of the higher-multiplicity matrix elements, that contribute renormalisation scale uncertainties.
In Fig. 7, we consider the MeNloPs case. We combine an NLO matrix element for the 0-jet multiplicities with LO matrix elements for the 1-and 2-jet multiplicities. The scale uncertainty for low p W ⊥ values now features the reduced scale uncertainty, that we already have seen in the NloPs validation.
The same is true in the MePs@Nlo case depicted in Fig. 8. A direct comparison of the scale uncertainties to the MeNloPs case is not straightforward though, as we combine NLO matrix elements for the 0-and the 1-jet multiplicity, where the virtual amplitudes are obtained from BlackHat [42]. Hence, the 2-jet multiplicity is described at leading order through the 1-jet H-events. As such, the set-up is not a simple upgrade from our MeNloPs calculation.
In all multijet-merging validations, we find a similar behaviour with respect to the imprint of including emissions in the reweighting. For n NloPs + n PS = 2, the dedicated  In the lower plots, n PS is varied for comparison. Again, we find a saturation when reproducing dedicated calculations for n PS ≥ 2, with no further improvement when n PS is increased from 2 to 3.  In the upper four plots, n NloPs = 1 and n PS = 2, thus up to three emissions are reweighted. In the lower plots, both n are varied for comparison. Again, we find a saturation when reproducing dedicated calculations for n NloPs + n PS ≥ 2, with no further improvement when n PS is increased from 1 to 2.  In the lower plots, both n are varied for comparison. Again, we find a saturation when reproducing dedicated calculations for n NloPs +n PS ≥ 2, with no further improvement when n PS is increased from 1 to 2.
calculations are well reproduced, and no further improvement is found for n NloPs +n PS = 3. It is noteworthy, that for the MeNloPs case we find a worse reproduction for n NloPs = 1 and n PS = 0 compared to the NloPs and the MePs@Nlo cases. This originates in the fact that in the latter two cases, we enable the reweighting of emissions off S-events at all involved multiplicities, whereas in the MeNloPs case only the first of the three multiplicities is affected, because the other two are at LO and therefore do not have S-events. Thus, the overall importance of the S emission reweighting gets restricted to the region below Q cut of the 1-jet configuration in the MeNloPs case.

Consistent variations
In general, the renormalisation and factorisation scales, α s and the PDFs should be varied consistently throughout any of the presented calculations. While at fixed order the situation is clear, the matched and merged approaches allow for some degree of freedom regarding partial variations while still retaining their respective accuracies.
In the simplest case, LOPS, µ R and µ F of the short distance cross section and the parton shower may be varied independently as these variations can be expressed as higher-order terms in a perturbative expansion in the coupling parameter α s . This is not the case for α s itself and the PDFs as they are fixed through measured input values and parametrisations. Changes in these input values cannot be expressed as simple higher-order terms. Thus they need to be chosen consistently throughout.
Similarly, in NLOPS calculations, the renormalisation and factorisation scales may be varied in the matrix element (B and H A ) or the parton shower (PS NloPs and PS) separately without losing neither the fixed-order nor the resummation accuracy. As the pseudo-subtraction through the D A in any case employs different scales in PS NloPs and the B and H A functions, it always leaves remainders of O(α 2 s ). Hence, further scale variations in either one, the short-distance cross sections or the PS NloPs , do not worsen the nominal accuracy of the method. Retaining the logarithmic accuracy of the parton shower on the other hand requires identical renormalisation and factorisation scales throughout all resummation-relevant components, i.e. PS NloPs and PS. Again, variations in α s or the PDFs need to be consistent throughout the calculation.
The multijet-merged calculations impose further constraints since they treat multijet matrix elements and parton-shower emissions on the same footing. The notation of the scales already reflects this for µ R and µ F . In their definitions only the core scales remain as free parameters and may be varied independently. Again, the α s and PDF parametrisations need to be the same throughout.

Conclusions
In this publication we have presented the implementation and validation of reweighting techniques allowing for the fast and efficient evaluation of perturbative systematic uncertainties in the Sherpa event-generator framework. We have lifted the available techniques for the determination of PDF, α s and scale uncertainties in leading-and nextto-leading order QCD calculations to include the respective variations in parton-shower simulations. In turn we provide the means to perform consistent uncertainty evaluations for multijet-merged simulations based on leading-or next-to-leading-order accurate matrix elements of varying multiplicity matched with parton showers. The foundation for our reweighting method is the knowledge of the very dependence structure of the perturbative calculations on the parameters to be varied. For the fixed-order components this amounts to the corresponding decomposition of the Catani-Seymour dipole subtraction terms. This needed to be supplemented by the reweighting of the parametric dependences of the parton shower treated through the Sudakov Veto Algorithm.
With our extensive validation we have been able to prove on the one-hand-side the correctness of the implementation and have, furthermore, been able to illustrate the importance of parton-shower reweighting for reliable uncertainty estimates. With comparably little additional computational costs this allows for the on-the-fly determination of PDF, α S and scale uncertainties based on one single generator run, that, otherwise, would require explicit re-computations. The overall reduction in CPU time is by a factor of about 3 to 20, depending on the event-generation mode used, see App. A. The variational event weights provided are easily accessible through the HepMC event record and are furthermore consistently handed over to the Rivet analysis software by the corresponding Sherpa interface.
The methods presented in this publication are ideally suited for event-wise uncertainty estimates and can readily be used in arbitrary theoretical and experimental analyses. An extension to next-to-next-to-leading-order QCD calculations possibly dressed with parton showers, as presented in [16,17,71], is straightforward and planned for the near future.
The decomposition of the fixed-order part of QCD calculations employed here is also a necessary ingredient to produce cross-section grids as provided by the APPLgrid [72] and FastNLO [73,74] tools. These store the perturbative coefficients for a certain observable calculation discretised in Q 2 and x. Using interpolation methods, this allows for the a posteriori inclusion of PDFs, α s and variations of the renormalisation and factorisation scales. In turn, such techniques are well suited for (combined) fits of PDFs and α s that require a multitude of re-computations of the theoretical predictions. Over the last years, tools have been developed that automate the projection of arbitrary next-to-leading-order QCD calculations onto such grids, namely the aMCfast [75] and the MCgrid [76,77,78] packages. The first one produces APPLgrids with MadGraph5 aMC@NLO [79], the latter APPLgrids or FastNLO grids from Sherpa events projected on the observables through Rivet. The APFELgrid tool [80] provides an improved convolution method for use with APPLgrid files that furthermore speeds-up the re-evaluations.
However, none of these approaches includes generic parton-shower effects, i.e. the parametric dependence of the shower component on the PDFs and α s is ignored. With the methods presented in this publication we are confident that we can surmount this limitation and in the future provide interpolation grids that properly reflect showerresummation effects and allow for the inclusion of the affected phase-space regions in PDF determinations.
The benefit of reweighted calculations is given by the saving of CPU time. In order to evaluate the gain, we shall compare the event generation time of reweighted calculations with the sum of generation times for all corresponding dedicated computations. Here we consider both parton-level calculations, as well as runs including multiple interactions and hadronisation, the typical default in physics analyses applications. For the latter it can be expected that the gain in CPU time by using the reweighting approach is most considerable, as the CPU intense non-perturbative event generation phases do not need to be re-evaluated. In what follows we compare actual event-generation times, neglecting the set-up times of the individual runs. 9 In Fig. 9, we consider event generations using LoPs, NloPs, MePs@Lo and MePs@Nlo calculations for pp → W[e − ν ] at 13 TeV. The ratio of CPU time between the reweighting and the dedicated generations is shown for different maximum numbers of reweighted shower emissions n PS + n NloPs . Whether non-perturbative effects are included or not, the time needed for the reweighting calculation is below 10 % of the time needed for dedicated calculations if only the matrix element is reweighted (n PS = n NloPs = 0). The ratio then increases for larger numbers of reweighted emissions, as their reweighting needs additional time, asymptotically approaching the value when all parton-shower emissions are reweighted. For parton-level-only calculations, this ratio is around 0.35 for LoPs events, and around 0.3 for NloPs events. This reduction can be explained due to relatively smaller computational cost of the parton shower as a whole when the rest of the calculation is more complex. Also note that n PS for LoPs is only equivalent to n PS +n NloPs for S events. H events do not feature the S-Mc@Nlo-emission, and hence for them n NloPs does not contribute to their reweighting.
For the same reason, when non-perturbative effects are included, that ratio improves to about 0.1: The parton shower (and its reweighting) component plays a relatively smaller rôle in terms of CPU cycles, when multiple interactions and hadronisation are enabled.
If on top of the non-perturbative effects the events are also unweighted, the ratio does not change in the LoPs case, but in the NloPs case (by about 20 %). A reason might be, that only for NloPs a sizeable number of events gets rejected. For these, the jet evolution and non-perturbative phases are not performed at all, whereas the matrix-element calculation (and its reweighting) is always done, for accepted and rejected events alike. The same is true for the S-Mc@Nlo emission from S events. As a consequence, the relative cost of the reweighting grows slightly. A future improvement of the implementation would postpone these futile reweightings in unweighted calculations to a time point after the possible rejection. This of course requires that the dependence of the rejection probability is negligible. For the observables studied so far this was found to be true, at least to O (10 −4 ). Note that the effective gains will be lower than the results presented in this section, when we take into account the reduced statistical accuracy which comes with the partonshower reweighting. This requires more events to be generated in a reweighting calculation to reach the same statistical accuracy as in a dedicated calculation.

B. Configuring and accessing event-weight variations
Sherpa provides a list of pre-calculated alternative event weights, which are automatically output to the HepMC event record [24] or directly to an interfaced Rivet analysis [25]. For versions of Sherpa later than v.2.2.0, the variations to calculate can be specified with the following line in the (run) section of the Sherpa run card: VARIATIONS muR2fac1,muF2fac1,PDF1 muR2fac2,muF2fac2,PDF2 ...; Each variation is characterised by up to three arguments muR2fac a prefactor multiplying the nominal (squared) renormalisation scale muF2fac a prefactor multiplying the nominal (squared) factorisation scale PDF a parton density and its accompanying α s parametrisation.
This syntax works for all employed scale setters of Sherpa and both Sherpa's internal PDFs and PDFs interfaced through Lhapdf5/6 [81,82]. If trailing arguments are omitted from a variation, their default values are used, which is 1.0 for scale factors and the PDF set used by Sherpa for the nominal calculation.
In HepMC event records (v. 2.06 or later), the alternate weights can be accessed as named weights within the HepMC::WeightContainer of each event. The keys are given in one of the following formats:

MUR<muR2fac>_MUF<muF2fac>_PDF<ID> MUR<muR2fac>_MUF<muF2fac>_PDF<ID>_PSMUR<muR2fac>_PSMUF<muF2fac>
The parts in angle brackets are replaced with the respective scale factors and Lhapdf IDs. The second form is used, if a factor is applied to the renormalisation/factorisation scale of parton-shower splittings. This includes splittings within cluster histories determined by the multijet merging procedure, as discussed in section 4. If the scale reweighting with parton-shower splittings has been enabled (we discuss below how to do so), the scale factors for MUR, MUF and PSMUR, PSMUF are always equal, respectively, in the current implementation.
If the internal Rivet interface of Sherpa is used to analyse events during the generation, one histogram file per variation is written to disk, along with the nominal one. The file names follow a pattern resembling the HepMC weight-container keys as specified above.

Scale variations
The scale argument can also be specified by enclosing it in square brackets: [mu2fac]. This syntactic sugar implies both the given factor, its inverse and the default value. PDF and α s variations PDF and α s variations both work by specifying a PDF set through the PDF argument of a variation. This is because Sherpa per default uses the value for α s (m 2 Z ) given by the PDF set in use. Therefore an α s variation can be achieved by using PDF fits for different values of α s (m 2 Z ). To specify a specific member of a PDF set, its number is given as an additional argument separated by a slash. Thus, 1.0,1.0,CT14nlo/38 asks for the 38th member of the CT14nlo PDF set, without modifying the renormalisation and factorisation scales. If the slash and the number are not given, the central PDF member is used, i.e. CT14nlo is equivalent to CT14nlo/0. Sherpa can also be asked to do variations for all members of a PDF set by enclosing it in square brackets. Hence, 1.0,1.0,[CT14nlo] is equivalent to

Configuring how variations are calculated
The following options always affect all variations that are specified by arguments to the VARIATIONS keyword.
REWEIGHT SPLITTING ALPHAS SCALES (default: 0) If this is set to 1, the renormalisation scale factor is applied to the α s argument of individual splittings, instead of applying it only to the overall renormalisation scale, see section 4.2. This means that partonshower splittings are only included in the rescaling, if this option is enabled. In the notation of sections 3 and 4, this setsk αs =μ R /µ R .
REWEIGHT SPLITTING PDF SCALES (default: 0) If this is set to 1, the factorisation scale factor is also applied to PDF scale arguments within shower splittings (and intermediate cluster history PDF ratios), and not only to the core-process PDFs. In the notation of sections 3 and 4, this setsk f =μ F /µ F .
REWEIGHT MAXEM (default: -1) This option specifies the number of ordinary parton-shower emissions included in the reweighting per event. If this is set to 0, no emission is reweighted. The default value -1 means that all emissions should be reweighted.
REWEIGHT MCATNLO EM (default: 1) If this is set to 0, the single parton-shower emission within the Mc@Nlo contribution is not reweighted.
VARIATIONS INCLUDE CV (default: 1) If this is set to 0, the behaviour of the square bracket syntax is changed, such that the central-value variation is not included when expanding a parameter in square brackets. It is recommended not to disable it, such that one can do a closure test between the dedicated calculation and the reweighting. However, in CPU intensive applications, this setting can be used to omit this one obsolete variation while still making use of the convenient square-bracket syntax.