Machine Learning Templates for QCD Factorization in the Search for Physics Beyond the Standard Model

TL;DR

This paper introduces a novel GAN-based method for modeling jet mass distributions in high-multiplicity hadronic final states, improving background estimation for new physics searches at the LHC.

Contribution

It proposes a GAN-based approach, including a modified WGAN, for unbinned, conditioned background modeling, outperforming traditional histogram methods.

Findings

GAN methods outperform histograms in modeling jet mass distributions

Modified WGAN with circle constraint improves stability and performance

Approach enhances sensitivity of LHC searches for new physics

Abstract

High-multiplicity all-hadronic final states are an important, but difficult final state for searching for physics beyond the Standard Model. A powerful search method is to look for large jets with accidental substructure due to multiple hard partons falling within a single jet. One way for estimating the background in this search is to exploit an approximate factorization in quantum chromodynamics whereby the jet mass distribution is determined only by its kinematic properties. Traditionally, this approach has been executed using histograms constructed in a background-rich region. We propose a new approach based on Generative Adversarial Networks (GANs). These neural network approaches are naturally unbinned and can be readily conditioned on multiple jet properties. In addition to using vanilla GANs for this purpose, a modification to the traditional WGAN approach has been investigated where weight clipping is replaced with a naturally compact set (in this case, the circle). Both the vanilla and modified WGAN approaches significantly outperform the histogram method, especially when modeling the dependence on features not used in the histogram construction. These results can be useful for enhancing the sensitivity of LHC searches to high-multiplicity final states involving many quarks and gluons and serve as a useful benchmark where GANs may have immediate benefit to the HEP community.

Tables

Table 1. Table 1: Phasespace requirements for the different regions considered..

Region	$N_{\text{jet}}\left(p_T>200\text{GeV}\right)$	$p_{T,1}$	$\left|\Delta {\eta }_{1,2}\right|$	$M_J^{\Sigma }$
Control	$=3$	-	-	-
Validation	$=4$	$>400\text{GeV}$	$>1.4$	-
	$\ge 5$	-	$>1.4$	-
Signal	$=4$	$>400\text{GeV}$		$>1.0\text{TeV}$
	$\ge 5$	-		$>0.8\text{TeV}$

Equations

$$W({\mathbb{P}}_{\text{real}},{\mathbb{P}}_{\text{generated}})=\underset{\gamma \in \Pi ({\mathbb{P}}_{\text{real}},{\mathbb{P}}_{\text{generated}})}{\inf }{\mathbb{E}}_{(x,y)\sim \gamma }[\parallel x-y\parallel ],$$

$$\displaystyle x^{(k+1)}_{n}=\psi\bigg{(}\bigg{|}\sum_{m}x^{(k)}_{m}e^{i\theta^{(k)}_{nm}}+e^{i\phi^{(k+1)}_{n}}\bigg{|}\bigg{)},$$

$$S({\mathbb{P}}_1,{\mathbb{P}}_2)=\frac{1}{2}{\int }_X\frac{({\mathbb{P}}_1(x)-{\mathbb{P}}_2(x){)}^2}{{\mathbb{P}}_1(x)+{\mathbb{P}}_2(x)}\text{d}x,$$

$${\mathbb{P}}_{\text{real}}$$

$${\mathbb{P}}_{\text{generated}}$$

Full text

1 Introduction

Even though collimated sprays of particles (jets) produced from high energy quarks and gluons are ubiquitous at the Large Hadron Collider (LHC), analyzing their substructure has proven to be a powerful tool in the search for physics beyond the Standard Model (SM) Larkoski:2017jix ; Asquith:2018igt . Many theories of physics beyond the SM predict new particles with cascade decays that can result in large multiplicity final states. When many quarks and gluons are produced in these cascades, multiple large-radius jets with non-trivial substructure can be created Cohen:2012yc ; Hedri:2013pvl ; Hook:2012fd . As a result, one powerful method for searching for new particles in the all-hadronic channel is to look for events with a large $\sum_{j\in J}m_{j}$, where $m_{j}$ is the jet mass and $J$ is a set of jets in an event. The key challenge for such an analysis is to estimate the SM background, as high multiplicity multi-jet final states are difficult to accurately predict with current simulation tools.

Based on the approximate factorization of quantum chromodynamic (QCD) jet production at the LHC Collins:1989gx , the authors of Ref. Cohen:2014epa proposed an innovative background estimation technique. The idea of the procedure is to estimate the conditional probability $p(m_{j}|\text{jet kinematics})$ with an event selection suppressed in signal and then to convolve it with the jet kinematic spectrum in the signal region (from data). A comparison between the predicted $m_{j}$ and observed $m_{j}$ is then sensitive to the presence of new particles. The ATLAS collaboration has successfully applied this method in both Run 1 and Run 2 to set strong limits on potential gluino and squark production Aad:2015lea ; Aaboud:2018lpl .

The background estimation procedure described above has two major limitations111The extensive smoothing studies in Ref. Cohen:2014epa help to mitigate binning effects, but do not have an impact on the feature conditioning challenge.. First, $p(m_{j}|\text{jet kinematics})$ is represented as a histogram and each bin is uncorrelated so many events are needed for a precise determination. Second, there are physics and detector effects which change the distribution of the jet mass between the region it is constructed (‘trained’) and the region where it is applied (‘tested’). For example, the quark/gluon composition of the background may be different between the two regions. One way to mitigate this source of method bias is to condition on more features of the jet when constructing the conditional probability. To reduce the impact of changes in quark/gluon composition, one could add the number of charged-particle tracks inside the jet. Gluon jets tend to have more particles than quark jets due to their larger color factor. Since the templates $p(m_{j}|\text{jet kinematics})$ are binned, it is not simple to condition on more features as one needs many more bins and thus larger samples of events for training.

This paper proposes a solution to both of these limitations using modern machine learning. Deep neural networks are becoming popular tools for classification and regression tasks in high energy physics (HEP) data analysis, but there is a growing machine learning literature on neural network-based generative models as well. Training a generative model can be viewed as a regression task that maps noise to structure, mimicking the Jacobian from a pre-defined probability distribution to a target probability distribution. One of the most well-studied paradigms for such models is the Generative Adversarial Network (GAN)Goodfellow:2014upx (details in Sec. 2). GANs have also been studied in HEP and show great promise for accelerating simulations Paganini:2017hrr ; Paganini:2017dwg ; deOliveira:2017rwa ; Chekalina:2018hxi ; Carminati:2018khv ; Vallecorsa:2018zco ; Erdmann:2018jxd ; Musella:2018rdi ; Erdmann:2018kuh ; ATL-SOFT-PUB-2018-001 ; deOliveira:2017pjk and may also be useful for other tasks such as sampling from the space of effective field theory models Erbin:2018csv . In the context of QCD factorization studied in this paper, the GAN will learn the probability distribution of the jet mass given the jet kinematics and any other useful information.

This paper is organized as follows. Section 2 introduces GANs and how they can be used to exploit QCD factorization. The application of GANs to the phase space relevant to the Supersymmetry (SUSY) search from Ref. Cohen:2014epa ; Aad:2015lea ; Aaboud:2018lpl is presented in Sec. 3. The paper ends with conclusions and outlook in Sec. 4.

2 Machine Learning Architectures

2.1 Overview of GAN setup

The goal of this section is to introduce an approach to learn the conditional distribution of the jet mass given various kinematic features. Neural networks are chosen due to their flexibility and the resulting algorithms are naturally unbinned. There are multiple neural network-based approaches to generative modeling such as Variational Autoencoders (VAE) DBLP:journals/corr/KingmaW13 ; Rezende:2014:SBA:3044805.3045035 , Mixture Density Networks (MDN) mdn , and Generative Adversarial Networks (GAN) Goodfellow:2014upx . GANs are selected because they can readily model asymmetric distributions and accommodate conditional features.

Generative Adversarial Network training uses a pair of neural networks: one to map noise into structure (generator) and one to classify (discriminator) physics-based examples from the generator examples. The generator is structured as a densely connected feed-forward neural network that takes as input both jet kinematic features and noise and outputs a jet mass. The generated masses are then input to the discriminator network where they are compared against masses from a physics-based generator that match the kinematic quantities. The two networks ‘compete’ until the discriminator network is as bad as possible, which means that the generator is proficient at modeling the conditional jet mass distribution. This minimax structure uses the loss function described in Ref. Goodfellow:2014upx , constructed to minimize the Jensen-Shannon divergence between the distribution of the real data (in this case, a physics-based simulator) and the distribution of the generated data. A schematic of this GAN is shown in Fig. 1.

The GAN networks for the jet mass are relatively small compared to others in the literature, which are mostly used for modeling image data. For the vanilla GAN implementation, both the Generator and Discriminator networks are composed of three hidden layers between input (kinetic variables and noise/mass for generator/discriminator respectively) and output (generated mass/likelihood for generator/discriminator respectively). The first layer has 160 neurons, the second has 80 neurons, and the last hidden layer has 40 neurons. Additionally, the generator network is made more robust by adding 50% dropout JMLR:v15:srivastava14a . Network weights were chosen using the Adam optimizer DBLP:journals/corr/KingmaB14 with early stopping. These settings were chosen after a modest hyper-parameter scan.

All of the neural networks are built using Tensorflow tensorflow on Nvidia GeForce 1080 Ti GPUs. Since the jet mass distribution is skewed toward the high mass tail, the input noise distribution was varied according to a skew-normal distribution skewnormal : $f(x|\alpha)=2\phi(x)\Phi(\alpha x)$, where $\alpha$ is a hyper-parameter and $\phi$ and $\Phi$ are the probability density and cumulative distribution of the standard normal distribution, respectively. The best value of the skew parameter was identified to be $\alpha=5$ (labeled skew in the results plots). This will be compared to the standard no-skew case ($\alpha=0$; labeled noskew).

2.2 Modified WGAN

As an alternative to the vanilla GAN described in Sec. 2.1, a popular variant called the Wasserstein-GAN (WGAN) 1701.07875 was also studied (for another HEP application, see Ref. Erdmann:2018jxd ). The WGAN differs from the vanilla GAN in that it minimizes the Earth-Mover distance (also known as the Wasserstein Distance):

[TABLE]

where $\gamma\in\Pi(\mathbb{P}_{\text{real}},\mathbb{P}_{\text{generated}})$ is a joint distribution with marginal distributions $\mathbb{P}_{\text{real}},\mathbb{P}_{\text{generated}}$ respectively, and $(x,y)\sim\gamma$ means that the random variable $(x,y)$ is drawn from the distribution $\gamma$. This ’softer’ metric was introduced as a way to combat the vanishing gradients that often occured when training regular GANs 1701.07875 ; Arjovsky2017TowardsPM . In the algorithm suggested by the WGAN paper to minimise the Earth-Mover distance, the discriminator is a function $f$ that is taken from a space of trial functions that are all $K$-Lipshitz for some $K$. To enforce such conditions, the functions $f$ are constructed as feed-forward neural networks with weights $w$ that are clipped (after every update) to a compact space $[-\alpha,\alpha]$ for some fixed $\alpha$ which enforces $K$-lipshitz for some $K$.

An exact implementation of the WGAN approach resulted in weights that would often aggregate around the specific limit values $\alpha$ chosen, leading to vanishing gradients and generated mass distributions that did not match the physical mass distributions. The challenges surrounding the clipping operation to enforce the Lipshitz condition are discussed in the original WGAN paper, and explored in variations of WGAN where the clipping is replaced by gradient penalty DBLP:journals/corr/GulrajaniAADC17 or asymmetric clipping 2017arXiv170507164G .

In this paper, a modification to weight clipping is introduced that changes how the weights act on the outputs of neurons to enforce the Lipshitz condition. A schematic for this modification to the WGAN is shown in Fig. 2. The original weight-clipping operation enforces a Lipshitz function because the composition of Lipshitz functions is still a Lipshitz function; in particular, there is the weak assumption that all the activation functions used in the neural networks are themselves Lipshitz (true for the most popular activation functions including sigmoid, tanh, and ReLU). Then, the fact that there is a single constant $K$ for which all the functions $f$ that we are considering are $K$-Lipshitz is due to the weight-clipping – this restricts the function space to be compact.

Another natural way to enforce the Lipshitz constraint that eliminates boundary effects is to draw weights from a compact space with no boundaries. One way to achieve this is to draw the weights and biases from the unit circle222In principle, one could generalize this idea to modify the WGAN where all weights and biases are drawn from a generic compact Lie Group. The outputs of the neurons exist in a vector space with the Lie Group acting by a chosen representation. In this framework, a normal linear neural network uses the Lie Group $\mathbb{R}$ with canonical action on itself, and the modified WGAN shown in Eq. 2 uses $U(1)\approx S^{1}$. The point of this construction is that the compactness of the weight-space assures the convergence of the WGAN and is well-formed because Lie Group actions can be differentiated. Such applications of Lie Groups to Machine Learning have been explored elsewhere in the literature, for example in applications to 3D-classification problems DBLP:journals/corr/HuangWPG16 .:

[TABLE]

where $x^{(k)}_{m}$ represents the output of the $m$-th neuron in the $k$-th layer of the neural network and $\psi$ is the activation function. The bias is replaced by $b_{i}=e^{i\phi_{n}}$ and the weights are replaced by $w_{nm}=e^{i\theta_{nm}^{(k)}}$.

3 Templates for RPV-SUSY

3.1 Overview of RPV-SUSY and the Template Method

Supersymmetry (SUSY) Golfand:1971iw ; Volkov:1973ix ; Wess:1974tw ; Wess:1974jb ; Ferrara:1974pu ; Salam:1974ig is a well-studied extension of the Standard Model in which there is a new symmetry relating fermions and bosons. In models of SUSY that conserve $R$-parity FARRAR1978575 ; DIMOPOULOS1981150 ; PhysRevD.26.287 ; SAKAI1982533 ; DIMOPOULOS1982133 – an additional symmetry that requires SUSY particles to be produced in pairs – collider signatures often feature large missing transverse momentum (MET) carried by the lightest supersymmetric particle (LSP) which must be neutral under the electromagnetic and strong forces. However, $R$-parity is not present in all SUSY models and collider-based limits relying on MET are typically insensitive to such models. One $R$-parity violating (RPV) coupling (often denoted $\lambda^{\prime\prime}$) results in gluino/neutralino decay into three quarks, as illustrated in Fig. 3. When this coupling is present, traditional searches for SUSY are largely ineffective because there can be little MET and no charged leptons. Despite lacking standard handles for separating potential SUSY events from SM background processes, RPV signatures like those in Fig. 3, all-hadronic SUSY searches have been able to set strong limits on gluino and squark production in models with a large $\lambda^{\prime\prime}$ (see references within Ref Aad:2015lea ; Aaboud:2018lpl for other constraints on such models).

One approach to search for all-hadronic SUSY decays like those in Fig. 3 is to exploit the high-multiplicity of hard, well-separated partons in the final state. In the recent ATLAS analysis Aaboud:2018lpl , the main kinematic observable used in the search for $\lambda^{\prime\prime}$ RPV decays is the total jet mass $M_{\Sigma}^{J}$ defined as the sum of the masses of the four leading large $R=1.0$ jets. Multiple jets with a large mass are generated from well-separated high-energy partons that happen to be clustered within a single jet. The SM multijet background is estimated by a data-driven method Cohen:2014epa , whereby mass templates are constructed from control regions containing a low purity of potential signal. These templates are histograms that model the dependence of QCD mass on jet kinematic properties ($p_{T}$ and $\eta$). If the jets in the signal region are only due to QCD, then the mass distribution can be estimated by convolving the jet kinematics with the mass templates. This estimate is compared with the actual mass distribution and deviations would be an indication of BSM physics.

The goal of this study is to demonstrate that GANs may be a useful alternative method to simple histogram templates for learning the dependence of the jet mass on jet properties.

3.2 Simulation Setup

To demonstrate the potential of GAN-based templates, $pp\rightarrow\text{jets}$ at $\sqrt{s}=13\text{ TeV}$ are generated with Pythia 8.223 Sjostrand:2006za using the fast detector simulation Delphes 3.4.0 deFavereau:2013fsa with the detector card delphes_card_ATLAS.tcl. Following the ATLAS analysis strategy Aaboud:2018lpl , events are clustered Cacciari:2011ma into jets using the anti-$k_{t}$ algorithm Cacciari:2008gp with radius parameter $R=1$. These jets are groomed according to the trimming procedure Krohn:2009th where subjets with radius $R=0.2$ are created from the large-radius jet constituents and removed if their transverse momentum is below 5% of the parent jet’s $p_{T}$. The remaining large-radius jets are only considered if $p_{T}>200$ GeV and $|\eta|<2.0$ and are divided into control, validation and signal regions according to Table 1. The control region is used to construct the templates (train the GAN) and the signal region is where they are applied333We note that the original template proposal from Ref. Cohen:2014epa called for a smoothing procedure that was used in the earlier ATLAS result Aad:2015lea but not the later one Aaboud:2018lpl . Our analysis here more closely resembles the later approach. The smoothing should in principle help mitigate binning effects, though would not change the qualitative conclusions about feature dependence.. The validation region is expected to be sufficiently devoid of potential signals that it can be used to study the fidelity of the templates. Even though the Pythia is only leading order in the strong coupling constant, Ref. Aaboud:2018lpl used the same setup and found a good agreement with data so this setup is sensible for testing new methods.

As a baseline, mass templates (histograms) are constructed following the ATLAS study: jets in the control region divided into $4$ $|\eta|$ bins defined uniformly between [math] and $2$, and $15$ $p_{T}$ bins defined uniformly in $\log_{10}(p_{T})$. This results in $60$ mass histograms in total. No $b$-tagging is applied, though future extensions to include flavor tagging information are possible. Using the mass histograms, jets in the validation and signal regions can be dressed with random masses given their $p_{T}$ and $\eta$.

3.3 Machine Learning Results

After the selections described in the previous section, there were 1.1 million jets in the control region and 30k in the validation region. The validation region is used to test the efficacy of the neural network training, filling the role of the ‘test set’ and by construction is independent from the training set. The jets in the control region are split 50%-50% for the purpose of training the neural network and ‘validating’ the network to enforce the early stopping condition.

The accuracy of the generated mass templates was quantified444There is no unique way to monitor the GAN performance during training. For a typical GAN trained with images, this is qualitatively different than classifier training because the entire multi-dimensional probability is being modeled, not just the likelihood ratio. The 1D case here is not as extreme, but still different than classification or regression monitoring. One can use the full GAN loss, the discriminator loss, or any divergence that gives a ‘distance’ between probability distributions. This particular divergence is popular in HEP and is therefore used as a diagnostic here. using the separation power metric Harrison:1998yr ; Hocker:2007ht :

[TABLE]

where $\mathbb{P}_{1},\mathbb{P}_{2}$ are probability distributions over a space $X$ and Eq. 3 is normalized to be between [math] and $1$. For the RPV SUSY case, $X\simeq\mathbb{R}_{p_{T}}\times\mathbb{R}_{\eta}\times\mathbb{R}_{N}\times\mathbb{R}_{m}$, one real line for each of the jet properties $p_{T},\eta$, constituent track multiplicity555Due to their robustness to pileup and excellent angular resolution, charged-particle tracks are associated to jets and used as proxy for the number of particles inside the jet. ($N$), and $m$. The $\mathbb{P}_{i}$ are the real and generated distributions:

[TABLE]

where $G(m|p_{T},\eta,N)$ is the learned mass distribution by the GAN for a given $p_{T},\eta,N$ value, and $\mathbb{P}_{\text{real}}(p_{T},\eta,N)=\int\mathbb{P}_{\text{real}}(p_{T},\eta,N,m)\text{d}m$. Equation 5 explicitly encodes the QCD factorization of the jet kinematics and the mass given those kinematic properties. Neither the GAN or physics-based simulator provide $\mathbb{P}$ directly; instead, only examples are drawn from the distributions. Empirical distributions are constructed from samples and the separation power is approximated by first binning the jets into eight regions of equal statistics in the kinematic variables $p_{T},\eta,N$, by splitting the events into two collections based on the jet $p_{T}$, then splitting each of these collections evenly in $\eta$, and then in $N$. In each of these eight bins, the jet mass distribution is used to calculate the (binned) separation power for each dataset and then all eight sets666The estimated separation power, for our dataset size, is not sensitive to increasing bin number beyond $8$. are combined.

Figure 4 shows the separation power for various generative models as a function of the number of epochs used to train. GAN models that were initialized with a high separation error (above $0.6$) training much slower due to vanishing gradients, so only those GANs with an initialized value below 0.6 are considered for Fig. 4. Furthermore, to reduce the impact of fluctuations in the initialization, the average value over ten random initializations are used. By construction, the default template method does not involve neural networks and thus is constant. As desired, all of the GAN approaches converge to a separation power that is smaller than the template method, as they have access to more and unbinned information. For both the vanilla GAN and the WGAN, using a skew-normal distribution for the noise accelerates the training time. The Vanilla GAN also converges faster than the WGAN, though all GAN approaches have a similar final separation power.

The mass distributions in the validation region are presented in Fig. 5, in bins of jet $p_{T}$, $\eta$, and $N$. The average jet mass scales approximately as $\alpha_{s}\times p_{T}\times R$ and the width of the distribution also grows with $p_{T}$. Given the jet $p_{T}$, the jet mass should be approximately independent of $\eta$, aside from detector effects. Gluon jets have a large jet mass than quark jets and also a higher constituent multiplicity so there is a positive correlation between $N$ and mass.

Overall, the level of agreement between the GAN and the real mass distributions is better than for the template method and the real distributions. This is particularly true for $N$, where the real mass is shifted to lower values for low $N$ (more quark-like) and to higher values for high $N$ (more gluon-like). The GAN is typically well within 50% of the real distribution, while the template method can be much more than a factor of 2 off of the real distribution for low and high $N$. This is particularly important if the quark/gluon composition is different between the control, validation, and signal regions, either by chance or because some quark-tagging is applied to suppress the QCD background in such high multiplicity final states Sakaki:2018opq .

The modeling of the jet mass distribution in the validation region is used to determine systematic uncertainties on the templates. Figure 6 explicitly constructs the systematic uncertainty as the sum in quadrature of the non-closure from the validation region and the control region statistical uncertainty. These uncertainties are computed for the $\sum_{j\in J}m_{j}$ distribution (sum of the masses of the top four jets in the event), which is the main observable used in the RPV SUSY search Aad:2015lea ; Aaboud:2018lpl . In blue we have the deviation from the exact ratio $1$ for the template mass distributions and the GAN mass distributions. The GAN outperforms the template method in both the low mass and high mass limits; also note that when placing a cut on $N_{\text{track}}$, the performance of the GAN becomes even more pronounced. This is particularly encouraging because we expect such RPV-SUSY signals to be quark jet dominated, with a lower multiplicity on average than a gluon jet dominated background. For a 50% quark jet efficiency requirement on all four jets ($N_{\text{track}}<26$), the uncertainty for the GAN approach is about 20% in the high $\sum_{i\in J}m_{j}$ tail while it is well over 100% for the template appraoch.

An important part of any background estimation technique is the associated systematic uncertainty. One of the main sources of uncertainty here is the limited size of the training set in the control region. The authors of Ref. Cohen:2014epa suggested a bootstrapping technique to estimate the uncertainty by rerunning the template procedure on bootstrapped datasets. In principle, one could do the same procedure for the GAN training, with one GAN per bootstrap dataset. More sophisticated methods include modeling GAN weights and biases as nuisance parameters to be profiled by the data with prior distributions.

A challenge for assessing previous GAN applications in HEP is that they have been designed to model high-dimensional feature spaces that are difficult to visualize and study Paganini:2017hrr ; Paganini:2017dwg ; deOliveira:2017rwa ; Chekalina:2018hxi ; Carminati:2018khv ; Vallecorsa:2018zco ; Erdmann:2018jxd ; Musella:2018rdi ; Erdmann:2018kuh ; ATL-SOFT-PUB-2018-001 ; deOliveira:2017pjk . The mass distribution example presented here provides a concrete testing ground to study GAN approaches where quantitative agreement can be studied and achieved using existing techniques. While this study used only leading-order simulations of jet production, the methods are applicable more generally and can be applied to collision data from the LHC experiments.

4 Conclusions

Generative Adversarial Networks have been proposed as an alternative to histogram-based mass templates for the background estimation in LHC searches for RPV SUSY. These methods rely on the approximate QCD factorization whereby a jets type and kinematic properties are sufficient for determining the distribution of the jet mass. The neural network approaches are naturally unbinned and can be readily conditioned on multiple jet properties. In addition to using vanilla GANs for this purpose, a modification to the traditional WGAN approach has been investigated where weight clipping is replaced with drawing weights from a naturally compact set (in this case, the circle). Both the vanilla and modified WGAN approaches were able to outperform the histogram method, especially when modeling the dependence on features not used in the histogram construction. When training such generative models for physical applications, the usual limitations of the method apply such as the potential for overfitting, sensitivity to hyperparameters, and vanishing gradients slowing training down - though methods of circumnavigating these limitations have been studied in the last few years, such as using ’softer metrics’ such as the Wasserstein metric. These results can be useful for enhancing the sensitivity of LHC searches to high-multiplicity final states involving many quarks and gluons and serve as a useful benchmark where GANs may have immediate benefit to the HEP community.

Acknowledgements.

We would like to thank Tim Cohen, Luke de Oliveira, Mustafa Mustafa, Michela Paganini, Max Swiatloski, and Jesse Thaler for their helpful feedback on the manuscript. This work was supported by the U.S. Department of Energy, Office of Science under contract DE-AC02-05CH11231.