Reweighting Parton Showers

We report on the possibility of reweighting parton-shower Monte Carlo predictions for scale variations in the parton-shower algorithm. The method is based on a generalization of the Sudakov veto algorithm. We demonstrate the feasibility of this approach using example physical distributions. Implementations are available for both the parton-shower modules in the Herwig 7 event generator.


Introduction
Monte Carlo simulations [1][2][3][4][5][6] have become essential tools in both the analysis of data from energy frontier particle physics experiments and the design of future experiments. The last ten years have seen a dramatic improvement in the accuracy of these simulations with the development of techniques to improve the description of high-multiplicity jet production at leading order 1 [8][9][10][11][12][13][14][15][16][17][18], simulations accurate at next-to-leading order, including the description of the hardest emission at leading order, [19,20], and more recently multiple jet production at nextto-leading order (NLO) [21][22][23][24][25][26][27]. This unprecedented increase in accuracy means that often the results of modern Monte Carlo event generators are the main, or even only, way in which theoretical predictions are compared to the latest results of the LHC experiments. This means that it is vital that wherever possible we must be able to assess the uncertainty on the predictions of event generators, as well as the central value of the prediction. These uncertainties come from a number of different sources: missing higher order corrections in the calculation of the hard matrix elements and shower evolution, normally estimated by varying the factorization and renormalization scales [28,29]; -uncertainties from the perturbative and non-perturbative modelling in the event generator, usually estimated by using different tunes of the event generator parameters [30]; -uncertainties from the fitting of the parton distribution functions (PDFs), now normally estimated using the recommendations of the PDF4LHC working group [31,32].
As event generators have become more sophisticated, in the calculation of the perturbative physics the time taken 1 See Ref. [7] for a recent review of older techniques.
for the calculation of the hard partonic configurations has increased which, together with the time taken for any simulation of the detector, means that it is often unfeasible to rerun the event generator for each scale choice, set of parton distribution functions, and non-perturbative parameters needed to fully assess all the sources of uncertainty in the Monte Carlo simulation.
In the calculation of the hard process the calculation of the uncertainty from the variation of the factorization and renormalization scales, together with the PDF uncertainty, can be more efficiently calculated. This is achieved by calculating the effect of changing the scale, or PDF, as a weight with respect to the central values at the same time as computing the central value. These weights can then be used to reweight the result of the simulation to obtain the uncertainties without requiring additional runs of the event generator. While this does increase the run time of the event generator, over that of simply performing the calculation for one PDF and scale choice, it is expected to be much more efficient than fully simulating events for all the required choices of scales and PDFs.
In contrast, currently the effect of varying the scale in the parton shower, or any of the perturbative parameters controlling the parton shower, or non-perturbative parameters for the simulation of the underlying event and hadronization, can only be calculated by running the full event simulation for each scale or parameter choice. Given the number of binary choices which are made in both the parton shower and the various non-perturbative models it is far from obvious that the variation of the parameters of these models can be achieved using a reweighting procedure.
In this paper we will present a generalization of the Sudakov veto algorithm used in most modern Monte Carlo event generators to generate the parton shower. This modification will allow us to calculate the effect of changing parameters in the parton shower via a reweighting of the central result, rather than a full resimulation of the events. The first approach [33] to calculating these weights on an event-by-event basis required calculating a weight for each variation which was complicated as this weight was not calculated via the veto algorithm making it difficult to implement in a full event generator, particularly for weights involving the PDFs. Performing the changes in the Sudakov veto algorithm was introduced for final-state radiation in [34]. Related work on modifying the Sudakov veto algorithm to address a number of applications has been presented e.g. in [35,36], while detailed studies regarding negative splitting kernels and effects of the infrared cutoff have been addressed in [37]. In the next section we will present the full details of the algorithm. This is followed by various checks that the reweighting procedure correctly reproduces the results of resimulating the events with different parameters and the calculation of the uncertainty for a number of physical distributions. In this first proof-of-concept paper we will concentrate on the effect of varying the scale used in the strong coupling and PDFs, although other changes can be simulated using the same approach. Finally we present our conclusions and directions for future work.

Standard Veto Algorithm
The standard veto algorithm proceeds, given the starting scale Q, to generate the scale of the next emission q and the d additional splitting variables 2 x according to the distribution where x µ is a parameter point associated to the cutoff µ, the splitting kernel is P (q, x) and the Sudakov form factor is The distribution S P is normalized to unity. The standard veto algorithm proceeds by taking an overestimate of the kernel R(q, x) such that Normally, the overestimate is chosen to be integrable and invertible so q, x can easily be generated according to the overestimated distribution with a Sudakov form factor The generation of the splitting scale and variables starting at scale k = Q proceeds as follows: 1. A trial splitting scale and variables, q, x, are generated according to S R (µ, x µ |q, x|k). 2. If the scale q = µ then there is no emission and the cut-off scale, µ, and associated parameter point x µ are returned. 3. The trial scale and splitting variables are accepted with probability otherwise the process is repeated with k = q.
For a proof that this correctly reproduces the distribution in Eqn. 1, see for example [7,37].

Weighted Algorithm
We can generalise the veto algorithm to include weights while simultaneously relaxing the requirements so that P is not required to be positive and removing the restriction on R, Eqn. 3. In this case S P is still normalized to unity. In order to achieve this we need to introduce an acceptance probability (q, x|k, y) such that 0 ≤ (q, x|k, y) < 1.
In this case we start with a weight w = 1. The generation of the splitting scale and variables together with the calculation of the weight proceeds as follows: 1. A trial splitting scale and variables, q, x, are generated according to S R (µ, x µ |q, x|k). 2. If the scale q = µ then there is no emission and the cut-off scale, µ, and associated parameter point x µ are returned with weight w. 3. The trial splitting variables q, x are accepted with probability (q, x|k, y) and the returned weight is 4. Otherwise the weight becomes and the algorithm continues with k = q.
We stress that, in general, the acceptance probability can depend both on the point under consideration for a veto and the previously vetoed point, allowing the algorithm to be biased to traverse certain sequences more often than others. In general the algorithm is not guaranteed to terminate, however this is not an issue for the applications we are considering.

Proof of the Algorithm
In order to prove that this algorithm gives the correct result we note that the probability density for the algorithm to traverse a sequence (q, x|q n , x n |...|q 1 , x 1 ) of n − 1 veto steps to give a result q, x from an initial condition where we have introduced an arbitrary parameter point x Q at the start of the algorithm, which may be chosen to improve on the efficiency of the algorithm. The weight accumulated through such a sequence is The density produced by the algorithm is therefore the difference R(q, x) − P (q, x) exponentiates when performing the sum as for the standard veto algorithm hence i.e. the correct distribution is produced.

Competing Channels
Often we have to deal with the case of competing processes, for example the branching of a gluon, g → gg and g → qq. This can be correctly handled using the competition algorithm, 3 i.e. generating a trial emission for all the possible processes and then selecting the one with the highest emission scale. In this case the proof follows the standard one for the competition. We can continue to use competition to generate the different branchings and still take the emission with the highest scale but using the weighted Sudakov veto algorithm the weight is the product of the weights for all the trial emissions, including those which are rejected.

Applications
There are many potential uses of this weighted Sudakov Veto algorithm. One important use which we will not consider here is handling more complicated kernels, which can be negative, while still generating emissions, albeit weighted, using standard techniques.
A second use, which is the main aim of this paper, is that it gives us a method of performing the parton shower for a default splitting kernel P (q, x) while at the same time calculating the weights for different choices of the kernel. Here we will consider the simplest possible choice, i.e. changing the scale used in the strong coupling and PDFs, but the method allows for any variation which can be expressed as a change of the kernel.
In this case we will choose the acceptance probability for the default choice of kernel P (q, x). Using this choice the unweighted result reproduces the result of the standard veto algorithm for the default kernel. While the weighted results will produce the result for different choices of the kernel. This choice ensures that in our case the weighted Sudakov veto algorithm will terminate. Variations of the splitting kernels can now be introduced by changing P → P in Eqns. 8a and 8b, while keeping the acceptance probability given in Eqn. 14.

Results
In this section we will show that using the weighted Sudakov veto algorithm allows us to correctly reproduce the effect of varying the scale of the strong coupling and PDFs in the parton shower kernels without the need for multiple runs of an event generator using different scales. We will use this new approach to study the uncertainty from scale variations for a few example physical distributions. We will use Herwig 7 [2] for these studies. Herwig 7 provides the option of using two different parton-shower algorithms: the default angular-ordered shower [38] using 1 → 2 branchings together with global momentum reshuffling to ensure momentum conservation; and a dipole based approach using local recoils [39]. The implementations of both of these showers use the veto algorithm to generate the emissions, in the case of the angular-ordered shower using simple overestimate functions and vetos while for the dipole shower the ExSample library [40] is used to adaptively sample the Sudakov distribution using the veto algorithm.
Herwig 7 allows us to compare the results of two physically different parton-shower algorithms. At the same time we also can check that the weighted results are correct using two algorithms which differ both in the physical approach used and in the technical implementation of the veto algorithm in the program.
Figs. 1 and 2 show the differential distribution for 1 − T for e + e − → qq with √ s = 91.2 GeV for the angularordered and dipole showers, respectively. These results were obtained at the parton level after the parton shower without any corrections to describe hard radiation.
Figs. 3 and 4 show the differential distribution for transverse momentum of the Higgs boson for gg → h 0 with √ s = 13 TeV for the angular-ordered and dipole showers, respectively. These results were obtained at the parton level after the parton shower without any corrections to describe hard radiation.
As can be seen the results of calculating the scale variation using our reweighting approach are in excellent agreement with those obtained from directly simulating the events with the modified scale choice in both cases.

Discussion
There are two main issues which effect the practicality of using our reweighting approach to calculate the scale uncertainty in the parton shower: 1. The time taken to calculate the result of the scale variations using reweighting should be less than running the event generator for the different scale choices considered. In general this will always be the case if the other stages of the event generation, for example the hard process evaluation, take significantly longer than the generation of the parton shower, or if detector simulation is included. However, for simple processes without detector simulation the time taken for the two approaches can be comparable, at least for the angularordered parton shower, Table 1. 2. If the weight variation is large then a large number of events will have to be simulated in order for the reweighted result to converge on that generated by directly simulating the events with an acceptable error. This can be particularly problematic if there are regions of phase space which would be populated with varied scales which are not filled for the central value, and hence have infinite weight.
The difference in the time taken with the two different shower algorithms is due to the different technical implementations of the veto algorithm. The dipole shower which uses an adaptive-sampling approach in which only one acceptance probability is calculated shows a significant reduction in the time taken for the simulation using  Table 1. Time taken (s) to simulate 10000 gg → h 0 events with √ s = 13 TeV for the angular-ordered and dipole showers. The time taken to simulate three scale variations ((µR, µF )/ √ 2, (µR, µF ), √ 2(µR, µF )) by direct simulation and reweighting is shown together with the fractional difference in times, i.e. (T (direct)−T (reweight))/T (direct). Events with only the hard process were simulated as well as events with multiple-parton interactions (MPI) both with and without varying the factorization and renormalization scales for the secondary scattering processes. reweighting because only one additional weight needs to be calculated for each variation.
The situation is very different for the angular-ordered parton shower where Eqn.6 is split into a number of different components. For example for space-like evolution the splitting kernel is where q is the angular-ordered evolution variable, x is the momentum fraction of the branching parton and P AP is the Altarelli-Parisi splitting function. A simple overestimate can be written as where α over S , P over AP (z) and PDF over (z) are the overestimates of α S (z(1 − z)q), P AP (z, q) and xf (x,q) , respectively. The veto is then separately applied for the weights This calculation is organised so that the most time consuming piece, i.e. the evaluation of w 3 , is only performed if the emission is accepted after the tests on w 1,2 . However, using the reweighting approach all the weights have to be evaluated for each trial emission, both for the default scale and any variations, therefore the shower evolution can be slower if the time taken for the evaluation of these weights is significant compared to that for the rest of the shower evolution. This can be seen in Table 1 where the angular-ordered parton shower is faster when the scale is varied directly. However, due to the additional weights which need to be calculated it is slower than the dipole shower, where the more sophisticated sampling of the Sudakov form factor means fewer additional weights have to be calculated, when reweighting is used.
However, when all the necessary parts of the simulation are included even for the relatively simple hard scattering processes considered here the reweighting approach is significantly faster and this performance improvement will only increase when more complicated, and hence timeconsuming, processes are simulated.
Owing to the need to divide out the veto probability in the reweighting procedure, weight distributions of the reweighted results may broaden significantly for very efficient sequences of the veto algorithm, (q, x|k, y) ∼ 1 (see Eqn. 8). While, for the central value, a very efficient algorithm is desirable it may at the same time force us to use larger statistics to obtain convergent results for the reweighted distributions.
The situation can be improved by explicitly making the veto algorithm for the central prediction more inefficient than originally designed by introducing a 'detuning' parameter λ > 1 to increase the proposal kernel, R → λR. A faster convergence of the reweighted results can hence be obtained. Despite the increase in the run time, using reweighting is still expected to be faster than a full simulation of each variation. Detuning parameters are available for the reweighting mechanisms in both Herwig 7 parton shower algorithms.
We show an example of the improvements in weight distributions in Fig. 5. The negative weight events appearing for the 'down' variation are significantly reduced for a moderate increase in run time (with λ = 4 taking about twice as long as no detuning). This increase in run time needs to be compared with the longer time taken to obtain a similar statistical error for the reweighted distributions without detuning as a large increase in the number of simulated events (and hence run time) will be required. We note that for the case of the overestimate kernel being an overestimate for all variations, no negative weights appear. This is the case, for example, when varying the scale in strong coupling, but not the PDFs, in the angular-ordered parton shower.

Conclusions
We have presented a new algorithm to allow the inclusion of weights in the Sudakov veto algorithm. This new weighted Sudakov veto algorithm allows the computation of the weights for any variations in splitting kernel used in the veto algorithm at the same time as the central value is calculated, allowing efficient computation of the shower uncertainties.
This allows us to assess the uncertainty due to variations of the scales in the parton shower without resimu- lating the events for each scale choice of interest. This is significantly faster and makes the calculation of the scale uncertainty feasible as the time taken to simulate an individual event can be time consuming for complicated processes.
This new approach is available in the Herwig 7 (7.0.2) release and will be combined with the effect of varying the scales in the calculation of the hard process in a future release. This technique can be extended to include the effect of varying the parton distribution functions or changes to the splitting kernel.

Note Added in Proof
While this work has been finalized, a similar approach by the Pythia collaboration [41] came to our attention, and the Sherpa collaboration has also reported comparable functionality in [42].