|Home | About | Journals | Submit | Contact Us | Français|
Accurate Monte Carlo simulations for high-energy events at CERN’s Large Hadron Collider, are very expensive, both from the computing and storage points of view. We describe a method that allows to consistently re-use parton-level samples accurate up to NLO in QCD under different theoretical hypotheses. We implement it in MadGraph5_aMC@NLO and show its validation by applying it to several cases of practical interest for the search of new physics at the LHC.
The search of new physics is one of the main priorities of the LHC. The recent observation of an anomaly in the di-photon spectra [1, 2] gives hope that we might have a first evidence of Beyond Standard Model (BSM) physics very soon. In that case, we would only be at the beginning of a long program of investigations of what the underlying physics is. In any case, searches of new particles or modifications of the interactions among the SM particles will continue as well as progress associated to our ability to provide precise predictions to be compared with data.
In the recent years, efforts have focussed on providing accurate theoretical predictions for a large number of BSM models at Leading Order (LO), in the form of event generators. First, various programs such as FeynRules , LanHep  or Sarah  have automated the extraction of the Feynman rules from a given Lagrangian. Secondly several matrix element based generators like MadGraph5_aMC@NLO  (referred to as MG5_aMC later on), Sherpa  or Whizard  have extended the class of BSM model they support with extensions in various directions: high spins, high color representations and any kind of Lorentz structure [9–11]. More recently, automated Next-to-Leading Order (NLO) prediction (in QCD) for BSM models are available thanks to the NLOCT  package of FeynRules which adds in the model the additional elements (R2 and UV counter-terms) required by loop computations.
It is now possible to generate Monte Carlo sample for a large class of BSM theories at LO and for an increasing number at NLO accuracy. Even though technically possible, producing samples for many models and benchmark points down to full detector level at the high luminosity expected at the LHC would require an unmanageable number of computing and storage resources. However, the stages of a simulation (parton-level generation, parton-shower and hadronisation, detector simulation, and reconstruction) are independent and factorise. Therefore changes in local probabilities happening at very short distance, i.e. from BSM physics, decouple from the rest of the simulation stages. This is particularly interesting since the slowest part of the simulation is the full simulation of the detector.
A logical possibility therefore arises: one can generate large samples under a SM or basic BSM hypothesis and then continuously and locally deform the probability functions associated to the distributions of parton-level events in the phase space by changing the “weight” of each event in a sample to account for an alternative theory or benchmark point. Under a not-too-restrictive set of hypotheses which are easy to list, such an event-by-event re-weighting can be shown to be exactly equivalent (at least in the infinite statistic limit) to a direct generation in the BSM. Note that such an event-by-event re-weighting is conceptually different from the very common yet very crude method where events are re-weighted using a pivotal one-dimensional distribution. Event-by-event re-weighting is a common practice in MC simulations, yet currently it has been only publicly available at LO [13, 14] or available at NLO for very specific cases (e.g. ) or in methods where NLO accuracy is far from ensured [16, 17]. It is the aim of this work to show that a consistent (and practical) re-weighting of events can also be done at NLO accuracy.
The plan of this paper is as follows. Before introducing the NLO re-weighting method, we will focus on the LO case in order to explain the intrinsic limitations of such types of methods (Sect. 2). In Sect. 3, we present three types of NLO re-weighting, two of them correspond to methods already introduced in the literature [14, 18]. The third one is the NLO accurate re-weighting method introduced here for the first time. In Sect. 4, we present some validation plots performed with MG5_aMC. We then present our conclusions in Sect. 5.
As stated in the introduction, the re-weighting method consists in attaching a new weight to every parton-level event as corresponding to a different scenario. The new weights allow to predict accurately (up to statistical precision) all the LO differential distributions at the parton level, leading also to the possibility of performing a single shower and detector simulation for all the models under consideration. At LO accuracy the new weight (Wnew) can be easily obtained from the original one (Worig) by simply multiplying it by the ratio of the matrix-elements estimated on that event for both models (noted respectively |Morig|2 and |Mnew|2) [13, 14]:
In practice, in a weighted Monte Carlo generation, the weights are simply given by1
where fi(xi, μF) is the parton-distribution function estimated on the Bjorken fraction xi at the factorization scale μF. ΩPS is the phase-space measure of the phase-space volume associated to the events.2 From this equation it is clear that Eq. 1 is the correct procedure since the weight is exactly multiplicative. This property is preserved by the unweighting procedure making Eq. 1 to hold for both weighted and un-weighted samples (an actual proof is presented in Appendix).
A few remarks are in order regarding the range of validity of this method. First, even if the method returns the correct weight, it requires that the event sampling related to Worig covers appropriately the phase-space for the new theory. In particular, Worig must be non-zero in all regions where Wnew is non-vanishing. Though obvious, this requirement is in fact the most important and critical one. In other words, the phase-space where the new theoretical hypothesis contributes should be a subset of the original one. For example, re-weighing can not be used for scanning over different mass values of the final state particles3, yet it is typically well-suited for probing different types of spin and/or coupling structures. More in detail, when the new theory has large contribution in a region of the phase-space where the original sample has only few events – since the original is sub-dominant in that part of the phase-space –, the statistical uncertainty of the re-weighted sample becomes very large and the resulting predictions unreliable. To appreciate quantitatively such an effect, we can use a naive estimator assuming a gaussian behavior. In that case one can write the estimated uncertainty as
where and Std(w) are respectively the mean and the standard deviation of the ratio of the weights and , are an observable and the associated statistical uncertainty. As a consequence, the relative uncertainty can be enhanced if the weights have a large variance. In Appendix, we introduce, as a proof of principle, a second method on how to estimate the statistical uncertainty from the distribution of the weights.
Second, the parton-level configuration feeder to parton-shower programs not only depends of the four-momenta but also of additional information, which is commonly encoded in the LesHouches Event File (LHEF) [19, 20]. Consequently, re-weighting by an hypothesis that does not preserve such additional information is not accurate. In general, such informations are related to:
Selected results obtained with this re-weighting are presented in Sect. 4.
In this section, we will present three re-weighting methods for NLO samples. First we will present a LO type of re-weighting that we dubbed “Naive LO-like” re-weighting introduced in VBFNLO (i.e. REPOLO ) and MadSpin [22, 23]. As it will become clear later, this method is not NLO accurate and should be used only if the difference between the two theories factorizes from the QCD production. The second method that we propose is original and consists in a fully accurate and general NLO re-weighting. Finally, we present the “loop-improved” re-weighting method  to perform approximate NLO computation for loop-induced processes when the associated two-loop computations are not available.
where R, S, C, SC, MC correspond respectively to the contributions of the fully-resolved configuration (the real), of its soft (including the Born matrix-element), collinear, soft-collinear limits (the counter-events) and the Monte Carlo (MC) counter-term. The (𝕊) (for standard) part corresponds to events generated with the Born configuration (N particles in the final state), while the (ℍ) (for hard) part corresponds to events generated with the real configuration (N + 1 particles in the final state). The MC counter-term (shower dependent) assures the coherent treatment with the parton-shower (no double counting) while preserving the NLO accuracy of the computation.
The Naive LO-like re-weighting computes the weights based on the multiplicity of the events before parton shower. i.e.,
, are respectively the weights for Born/real topology events for the hyppothesis (where is either the orig or new label). is the Born matrix element squared ( ) while is the real matrix element squared ().
As this method does not consider the dependence of the virtual contributions, it fails to be NLO accurate. To ensure NLO accuracy, it requires that the effect of the new theory factorises out, i.e., when
where is the finite piece of the virtual contribution (the interference term between the Born and the loop amplitude). Such relation should hold over the full phase-space with a universal constant since the MC counter terms connect the born and the real in a non local way. Nevertheless, as we will see later, the effect of the MC counter terms are quite mild, as expected since their contribution to the total cross-section are exactly zero by construction. This allows the Naive LO-like method to nicely approximate the NLO differential cross-section for many processes/theories where the last equation needs to be valid only phase-space point by phase-space point (i.e. when the ratio of the real matches the ratio of the Born and of the virtual in the soft and/or collinear limit).
In order to have an accurate NLO re-weighting method, one should explicitly factorise out the dependence in the (various) matrix elements (i.e. in the Born squared matrix element – ℬ – , the real squared matrix element – ℛ – and in the finite piece of the virtual – 𝒱 – ). We use the decomposition of the differential described in 4 introduced in the context of the evaluation of the systematics uncertainties:
where the α index is either R, S, C, SC, MC (see previous section). Q is the Ellis-Sexton scale and dχα is the phase-space measure.
The expression of the , , are given in the appendix of  and are not repeated here. All those expressions have linear dependencies in the Born, the virtual, the real, the color connected Born ℬCC (this term is defined Eq. (3.28) of ) and the reduced matrix element ℬRM (Eq. (D.1) of ). This allows us to decompose the corresponding expressions as:5
where the β index is either 0, R or F. The are expressions which do not depend of either the PDF/scale or the matrix-element. From this expression we define the following three terms:6
By keeping track of the at the generation time and writing it in the final event, one can perform an NLO re-weighting by:
The final weight associated to the event can then be calculated by combining those various pieces as it is done for the estimation of the systematics uncertainty (see Appendix of ). One can notice that both the color-connected Born and the reduced matrix-element are simply re-weighted by the ratio of the Born which can lead to a breaking of the NLO accuracy of the method. In the case of the color-connected Born, this does not consist in an additional limitation of the method since the re-weighting factors should differ only if the two theories present a difference in the relative importance of the various color-flows (a case already not handled at LO accuracy). The case of the reduced matrix-element is actually different since the contribution related to this matrix-element vanishes after integration over the azimutal angle . The infra-red observables are therefore not sensitive to such contribution and consequently neither on the re-weighting used for such contribution.
More generally, the possible drawbacks and limitations on the statistical precision of the method are the same as for the LO case. However, for NLO calculations in MG5_aMC we face one additional source of statistical uncertainty due to the method used to integrate the virtual contribution. This method reduces the number of computations of the virtual by using an approximate of the virtual contribution based on the Born amplitudes times a fitted parameter κ. It performs a separate phase-space integration to get the difference between the virtual and its approximation (full description of the method is presented in Section 2.4.3 of ). Schematically it can be written as:
If it exists a value of κ such that κℬ ≈ 𝒱, the second integral is approximately zero and does not need to be probed as often as the first integral (thanks to importance sampling ), reducing the amount of time used in the evaluation of the loop-diagrams. However the re-weighting proposed in Eq. 14 will highly enhance the contribution of the second integral since each term of the integral will be re-weighted by a different factor, having a direct impact on the statistical uncertainty.
To reduce this effect, we propose to use a slightly more advanced re-weighting technique. We split the contribution proportional to the Born () in two parts: and . is the part, proportional to the Born, related to the one of the counterterms, while includes all of the other contributions (the Born itself and the approximate virtual). We then apply the following re-weighting:
Both the virtual and the approximate virtual are re-weighted by the same pre-factor which should allow to limit the enhancement of the second integral. The demonstration that such re-weighting is NLO accurate is presented in Appendix. It can be intuitively understood considering (ℬ + 𝒱) as a single block which is re-weighted accordingly.
A third type of re-weighting was originally introduced in the context of multiple Higgs production [18, 28–30], which we now briefly describe. In this case the idea is to perform the NLO computation in the infinite top-mass limit and then re-introduce the finite top-mass effects via re-weighting. Equation 16 is directly applicable if the exact finite virtual part is known. If not, one can still use an approximate method:
Both this method and the Naive LO-like method are not NLO accurate. However one can expect that the loop improved method has a better accuracy than the other one due to the correct treatment of the various counter terms.
The various methods of re-weighting discussed in the previous section have been implemented in MG5_aMC and are publicly available starting from version 2.4.0. At the LO, the default re-weighting mode is based on the helicity information present in the event (Eq. 4), while for NLO samples, the default re-weighting mode is the NLO accurate one (Eq. 16). Fixed-order NLO generation can not be re-weighted since no event generation is performed in this mode. A manual of the code is available online at the following address: https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Reweight.
In this section, we will present four validation examples covering the various types of re-weighting introduced in the previous section. Since the purpose of this section is mainly to validate our method, the details of the simulation used (cuts, type of scale, ...) are kept to a minimum. Unless otherwise stated, the settings used correspond to the default value of MG5_aMC (version 2.4.0). In particular the minimal transverse momentum on jet is of 20 GeV at LO and of 10 GeV at NLO.
For the first validation, we will use the effective field theory (EFT) in the Electro-Weak sector . We will focus on the associated production of the W and Z boson for the following dimension six operator:
and gW is the weak gauge coupling, τI are the pauli matrices and is the gauge Field of SU(2).
In Fig. 1 we present the differential distributions for the transverse momenta of the Z boson at LO accuracy. Starting from a sample of Standard Model events (black solid curve), we have re-weighted our sample to get the SM plus the interference term with the dimension six operator for two values of the associated coupling: c = 50 TeV-2 (dashed blue) and c = 500 TeV-2 (dashed green). This second value is clearly outside the validity region for the EFT approach as the differential distributions turns to be negative at low transverse momentum. Nevertheless, having such large effects is interesting for the validation of the re-weighting method. The same differential distributions are generated with MG5_aMC (solid green and blue) and validates the re-weighting method.
The ratios between the differential curves obtained with each method are presented in the second inset. This inset contains also the statistical uncertainty (yellow band) for the ratio of two independent SM samples. The compatibility of those two ratio plots with the expected statistical fluctuation validates our approach/code implementation. The first inset presents the ratio between the EFT and SM predictions. It shows that the method works correctly for quite small and quite large modifications of the differential distributions.
One can note that in the context of EFTs, the weight is linear in the dim-6 coupling7 therefore it is trivial to predict the weight from any value of the coupling as soon as the weights for two different values of the coupling are known. This property can be used to further speed up the computation of the weight.
For our first NLO validation, we consider the associated production of a Z and H boson in the EFT as implemented in the Higgs Characterisation framework/model . We use two of the benchmarks introduced in : HD and HDder. In more details, the effective Lagrangian relevant for this example is
where Λ is the high energy scale (set to 1TeV), κHWW, κH∂Z, κH∂W are dimensionless couplings (set to one). H is the Higgs doublet field and Vμν = ∂μVν - ∂νVμ; V = Z, W-, W+.
In Fig. 2 we present the differential cross-section for the transverse momentum of the Higgs and for its rapidity. In both cases, we present the curve for the SM, HD and HDder benchmarks. For the transverse momentum, we start from an HDder sample of events and perform the re-weighting to the other scenarios. While for the rapidity we present the plot where the original sample is the HD theory. Each re-weighted curve is then compared with a dedicated generation and the associated ratio plot is displayed below with the statistical uncertainty expected for the generation of two independent samples. The agreement between the two is excellent for both the NLO accurate re-weighting and the Naive LO-like re-weighting. In this case the NLO QCD effects factorise from the BSM ones and therefore the NLO accuracy of the Naive LO-like approach can only be spoiled by MC counter terms – which are as expected quite mild – .
From the comparison of the two methods for the HD curve in the plot of the transverse momenta (top plot), we can observe that the statistical fluctuations are more pronounced for the curve obtained by re-weighting. This is an example of enhancement of statistical uncertainty due to the re-weighting as discussed around Eq. 3 since in the high pT region, the HDder is suppressed compare to the other theories under consideration (HD and SM).
In this second NLO example, we will use the EFT framework in the context of the top-quark  and focus on the chromomagnetic operator:
where Q is the third generation left-handed quark doublet, φ and t are respectively the Higgs and top quark fields, gs is the SM strong coupling constant, yt is the top-Yukawa coupling and TA is the SU(3) generator.
In Fig. 3, we present the transverse momentum of the Z boson in the associated production of this boson with a top/anti-top quark pair. We present the result for both the full matrix element squared (labelled σ(2)) and for the SM contribution plus the interference with the dimension 6 operator only (labelled σ(1)).
As in the previous section, we present our prediction both via the Naive LO-like re-weighting method (RWGT-LO) and via the NLO accurate one (RWGT-NLO), The ratio to the SM curves are presented in the first inset while the ratio between our prediction and the direct computation in MG5_aMC for σ(2) is presented in the second inset. The green band represents the expected statistical uncertainty for the ratio of two MG5_aMC samples. It is not possible to extract in automatic way the contribution of σ(1) from MG5_aMC and therefore we do not provide any comparison for this curve. As before, we observe a case where the statistical uncertainties are enhanced by the re-weighting approach and where both the Naive LO-like and the NLO re-weighting provides similar results. In this case the theory do not factorise and the ratio of the virtual and of the Born are not expected to match. The nice agreement is explained by the small contribution of the virtual and, once again, by the mild effect of the MC counter terms.
In this last example, we will present results for the associated production of a SM Higgs with one jet. In Fig. 4, we present the transverse momentum of the Higgs at both LO and NLO accuracy. For the LO case, we present three curves. The first one is the curved obtained within the heft model  featuring the dimension five operator obtained by integrating out the top quark (HEFT LO). The second line (SM LO/RWGT) is the one obtained by re-weighting the previous curve by the full one loop matrix element squared which contain the complete top-quark mass dependence. The last LO curve is the one obtained via direct integration of the one-loop amplitude squared by MG5_aMC  (SM LO). At NLO accuracy, we have the curve in the infinite top mass limit (HEFT NLO) using the Higgs characterization model . This sample is then re-weighted by the full-loop (Loop Improved) following the loop-improved method presented in the previous section. It is so far not possible to compute the NLO contribution directly in order to compare the accuracy of such method.
The first inset presents the ratio at LO and NLO of the infinite top mass limit over the full theory. For the NLO case, the full theory is approximated by the loop-improved method. The two ratios are very similar showing that the loop-improved method re-introduces the top-mass effects in a sensible way. The second inset presents the ratio between the re-weighting and the direct approach in the LO case, the statistical uncertainty of the ratio of two independent SM sample is presented by the yellow band. Its bumpy shape is due to the use of multiple samples with different cuts to decrease the statistical uncertainty. This ratio plot fully validates the re-weighting in the case of the LO curves.
We have presented the implementation of several methods that can be used for re-weighting LO and NLO samples and discuss the associated intrinsic limitations. We have released a new version of MG5_aMC that allows the users to employe the various re-weighting methods presented in this paper in a fully automatic and user-friendly way. In particular we have introduced for the first time an NLO accurate re-weighting method and compared it with the approximate methods available in the literature. Other re-weighting methods like the Naive LO-like and the loop-improved are for the first time available in a public code.
The comparison between the various methods shows that the approximate method (the Naive LO-like re-weighting) performs a satisfactory job. This indicates that the non locality of the MC counter terms is often more a theoretical problem than a contribution spoiling the NLO accuracy of the Naive LO-like re-weighting. Therefore the Naive LO-like re-weighting should be a good approximation in a quite large class of model/observable either when the virtual contribution is sub-dominant and/or when the effect of the BSM physics factorises. Consequently, we recommend phenomenologist to first test the Naive LO-like re-weighting and in case of loss of accuracy move forward to the slower NLO method. On the other hand for mass production at the LHC, where the samples are often used for more than one study, we recommend to always use the NLO accurate method.
The framework introduced here is flexible enough to accommodate different types of re-weighting approaches. In the near future we plan to capitalise on this to allow different type of functionalities. First we plan to implement a standard systematics uncertainty computation module as it is done in [16, 25, 38–42]. Compared to the existing module of MG5_aMC , this new module will allow to perform this determination independently of the event generation which will be extremely useful to evaluate the effect of a new PDF set/test a new scale scheme on existing samples. In a second stage, we plan to be able to compute the systematics uncertainty for the re-weighted BSM sample at the time of the re-weighting.
I would like to thank all the authors of MG5_aMC for their discussions, help and support at many stages of this Project. I would like also to thank C. Degrande, R. Frederix and F. Maltoni to have read and comment on this manuscript, E. Vryonidou, F. DeMartin, I. Tsinikos, V. Hirschi for their help during the validation of this implementation. O.M. is supported by a Durham International Junior Research Fellowship. This work is supported in part by the IISN “MadGraph” convention 4.4511.10, by the Belgian Federal Science Policy Office through the Interuniversity Attraction Pole P7/37, by the European Union as part of the FP7 Marie Curie Initial Training Network MCnetITN (PITN-GA-2012-315877), and by the ERC Grant 291377 LHCtheory: Theoretical predictions and analyses of LHC physics: advancing the precision frontier.
In order to have a formal proof that the unweighting procedure can commute with the re-weighting, we first have to formalize the procedure. Following the convention adopted in the previous sections, a standard Monte Carlo integration is:
To get an unweighted sample, we first need to multiply and divide this expression by :
Finally, the term can be re-interpretted as a probability to accept/reject the phase-space point.8
By randomly selecting a sub-sample of phase-space points with that probability, we reduce significantly the sample size. Additionally, all the remaining events have the same weight () and the associated distribution of events follows the physical distributions.
where Acci is either 0 or 1 depending on whether the event was kept or rejected following the probability distribution.
Let’s now proof that the re-weighting works on a unweighted sample, by doing the same for a second theory. But instead of multiplying and dividing by we will use the maximum weight of the original theory:
Since (See Eq. 23), this is equal to
We recover in that equation the same ratio which was used to unweight the original theory. We can therefore select the same sub-sample of events and just re-weight them by the ratio of the matrix element squared.
One can notice that the estimated uncertainty can not be obtained via re-weighting for an unweighted sample due to the non linear dependence in the matrix element squared.
We will show in this section what needs to be done in order to build an estimator of the variance from a re-weighted sample. Following the idea of the unweighting procedure, we can rewrite the standard estimator of the variance by:
As for the unweighting case, we can re-interpret the ratio as the probability to keep the event during the unweighting procedure. Therefore after the event unweighting the equation can be read as:
In this case, a dependence remains in the unweighting probability as well as in the number of generated and accepted events. If those informations were kept during the unweighting procedure it would be possible to construct the above estimator of the variance. The re-weighting of such information is then possible and one can construct such an estimator for any re-weighted sample:
Note that in presence of multi-channel integration such information need to be provided for each channel individually.
This method is currently not implemented in MG5_aMC but we plan to include it in a near future and study the accuracy of such an estimator.
In order to proof that the re-weighting proposed in Eq. 16 is correct we first need to formalise the loop integration method. We will use in this section a simplified notation such that
Where B, V, C represents respectively the Born, the virtual and the counter terms contribution. Since the counter terms do not play any role in this optimisation procedure (and have a natural re-weighting) we will focus on the pieces: In this simplified formalism the phase-space optimisation method can be written has (see Eq. 15):
In those equations, we first (Eq. 38) add and subtract the approximant of the virtual: , while in the second equation we integrate on different statistics the two pieces of the sum. We run k times less phase-space point in the second and therefore have to multiply it by the factor k.
If we want to use the re-weighting on the sample generated via this method, we have to apply the same method with the same value of κorig
Inspired by Eq. 16, we will multiply all those terms by the identity factor :
We can rewrite the expression as the expected re-weighting formula plus some rest-over
If the same phase-space sampling is used for both parts (k = 1) then the second and third lines cancel. The remaining lines correspond to the re-weighting of Eq. 16. If both integral are sampled in a different way (k ≠ 1), then the cancellation is not exact but should still occur for large enough samples. Therefore this optimization method introduces a new contribution to the statistical uncertainty.
1For the simplicity of the discussion, we will always consider that the sum of the weights is equal to the total cross-section of the sample.
2The normalisation choice implies that the phase-space factor ΩPS is proportional to N-1 where N is the number of phase-space points used to probe the phase space.
3For intermediate particle a small variation of the mass – order of the width – is reasonable.
4We also use the same (MC) counter terms as described in that paper.
5Due to the presence of multiple couter terms, the kinematic configuration on which the matrix-element is evaluated is not unique: an implicit sum over such kinematical configurations is assumed here and in the rest of the paper.
6One can notice that for β = R, F due to the use of the Ellis-Sexton scale .
7There would also be quadratic contribution if we include the squared matrix element associated to the dimension six operator.
8For non definite positive quantity the same idea holds by using .