Signal mixture estimation for degenerate heavy Higgses using a deep neural network

If a new signal is established in future LHC data, a next question will be to determine the signal composition, in particular whether the signal is due to multiple near-degenerate states. We investigate the performance of a deep learning approach to signal mixture estimation for the challenging scenario of a ditau signal coming from a pair of degenerate Higgs bosons of opposite CP charge. This constitutes a parameter estimation problem for a mixture model with highly overlapping features. We use an unbinned maximum likelihood fit to a neural network output, and compare the results to mixture estimation via a fit to a single kinematic variable. For our benchmark scenarios we find a ~20% improvement in the estimate uncertainty.


Introduction
Machine learning techniques have already proven useful in particle physics, especially for separating signal from background events in analyses of LHC data. More recently, deep learning methods, such as multi-layer neural networks, have been shown to perform very well, due to their ability to learn complex non-linear correlations in high-dimensional data [1][2][3]. In this paper we study the performance of a deep neural network classifier, but rather than classifying signal vs. background we focus on estimating the mixture of different signal classes in a dataset. This is motivated by the not-unlikely scenario where a new (and possibly broad) resonance is discovered in future LHC data, but limited statistics makes the interpretation difficult, in particular the question of whether the signal is due to multiple degenerate states. a e-mail: steffen.maeland@uib.no b e-mail: inga.strumke@uib.no In such a scenario it will clearly be important to squeeze as much information as possible from the available data.
While the approach studied here is general, we take a Two-Higgs-Doublet Model (THDM) as our example scenario. In these models the Higgs sector of the Standard Model (SM) is extended with an additional SU (2) doublet, predicting the existence of a pair of charged scalars (H ± ) and three neutral scalars (h, H, A), one of which should be the observed 125 GeV Higgs. Several more extensive frameworks for New Physics predict a Higgs sector with the structure of a THDM, the prime example being the Minimal Supersymmetric Standard Model (MSSM). A further motivation for THDMs comes from the fact that the extended scalar sector can allow for additional sources of CP violation and a strongly first-order electroweak phase transition, as required for electroweak baryogenesis [4][5][6][7]. For a recent study of this, see [8].
We associate the light scalar h with the observed 125 GeV Higgs and take the heavier scalars H, A and H ± to be mass degenerate. The focus of our study is on the ditau LHC signal from decays of the neutral states H and A, which in this case are indistinguishable save for their opposite CP charges. Searches for heavy neutral Higgses in ditau final states are carried out by both the ATLAS and CMS collaborations, see [9,10] for recent results.
The remainder of this paper is structured as follows. In section 2 we motivate why it is reasonable to expect a certain level of mass-degeneracy among the new scalars in THDMs and present our example THDM scenario. The technical setup for our analysis is given in section 3. Here we define our signal models, describe the procedure for Monte Carlo event generation and detail the neural network layout and training. In section 4 we demonstrate H/A signal mixture estimation using the arXiv:1804.07737v3 [hep-ph] 13 Dec 2018 method of fitting a single kinematic variable. The result serves as our baseline for judging the performance of the deep learning approach. Our main results are presented in section 5. Here we estimate the signal mixture via a maximum likelihood fit to the output distribution from a network trained to separate H and A ditau events. The results are compared to those from section 4. We state our conclusions in section 6.

Theory and motivation
The starting point for our study is a THDM scenario where m H ≈ m A . Our main motivation for this choice is to obtain a challenging test case for signal mixture estimation. However, there are also physical reasons to expect the H and A states to have similar masses. After requiring that the scalar potential has a minimum in accordance with electroweak symmetry breaking, we are left with a model with only two mass scales, v ≈ 246 GeV and a free mass parameter µ, to control the four masses m h , m H , m A and m H ± . From the point of view of the general THDM parameter space, the least fine-tuned way to align the light state h with SM predictions, as favoured by LHC Higgs data, is to move towards simultaneous decoupling of the three heavier states by increasing µ, leaving v to set the scale for m h = 125 GeV [11]. This points to a scenario where |m H −m A | 100 GeV, and quite possibly much smaller, depending on the quartic couplings of the scalar potential. 1 Further motivation for a small H-A mass difference can be found in less general realisations of THDMs. For the type-II THDM in the MSSM the quartic couplings are fixed by the squares of the SM gauge couplings, resulting in the tree-level prediction that m H − m A 10 GeV for m A ∼ 400 GeV and tan β ∼ 1, and decreasing further with increasing tan β or m A [13]. Another well-motivated scenario predicting closely degenerate H and A states is the SO(5)-based Maximally Symmetric THDM [14].
When mass degenerate, the H and A appear identical except for their CP charge. If the properties of the light h deviates from SM predictions, this difference in CP charge can manifest as non-zero ZZ and W W couplings for H, while for the CP -odd A the Zh coupling is available. However, these couplings all vanish in the perfect SM-alignment limit we assume here. Yet the CP nature of H and A is still expressed as spin correlations in fermionic decay modes, impacting the kinematics of subsequent decays. Here we study the channels H → τ τ and A → τ τ . Methods for reconstructing spin correlations in ditau decays of the 125 GeV Higgs have been investigated in detail [15][16][17][18], providing a good baseline for comparison. The use of neural networks to optimize CP measurements for the 125 GeV state is studied in [19,20].

Benchmark scenario
Two-Higgs-Doublet Models are classified in different types based on the structure of the Yukawa sector. We choose a benchmark scenario within the CP -conserving lepton-specific THDM, with m H = m A = m H ± = 450 GeV. In this model, the quarks couple to one of the Higgs doublets and the leptons to the other. This enables large branching ratios for H/A → τ τ , even for masses above the 350 GeV threshold for H/A → tt.
For comparison, in Appendix A we also show the result of a similar scan of the type-I THDM. In this model all fermions couple to only one of the two Higgs doublets. Compared to the lepton-specific THDM, the ditau signal in type-I THDM suffers a much stronger suppression from the H/A → tt channel.
As further described in sections 3 and 4, the mixture estimation techniques we study require each tau to decay through the τ ± → π ± π 0 ν channel, which has a branching ratio of 25%. However, the neural network method we employ can be extended to include other tau decay modes as well, by implementing the "impact parameter method" in [18] in addition to the "ρ decayplane method" used here.
If we only assume the τ ± → π ± π 0 ν decay channel and an acceptance times efficiency of 5%-10% for the signal selection, our example scenarios predict no more than ∼ 100 signal events for the anticipated 300 fb −1 dataset at the end of Run 3. However, as the model scan in Appendix A shows, considering slightly lower benchmark masses can provide an order of magnitude increase in the predicted cross-section. Also, extending the method to include more tau decay channels can greatly increase the statistics available to the analysis discussed here. Still, the large backgrounds in the ditau channel, e.g. from "fake QCD taus", implies that a signal mixture estimation study for the THDM benchmark scenario we present here likely will require the improved statistics of the full High-Luminosity LHC dataset.
We do not include a third mixture component representing ditau backgrounds for our benchmark study. Clearly, the inclusion of backgrounds will increase the uncertainty in the estimated H/A → τ τ signal mixture. However, as we discuss in more detail in section 5.1, the mixture estimate obtained from the neural network approach we study here is likely to be less affected by backgrounds than traditional mixture estimation from fitting a single kinematic variable.
For our further discussions we define the parameter α as the ratio of the A → τ τ signal strength to the total ditau signal strength, This is the parameter we seek to determine in our signal mixture estimation. 2 The parameter region of our benchmark scenario predicts values of α between 0.5 and 0.7. To allow for some further variation in the assumptions, we will in our tests use α values of 0.5, 0.7 and 0.9.
3 Analysis setup

Event generation
We generate 13 TeV pp Monte Carlo events for this study using Pythia 8.219 [29,30]. Only gluon-gluon fusion and bottom-quark annihilation are considered, as these are the dominant H/A production modes at the LHC. 3 For our analysis we select opposite-sign taus decaying to π ± π 0 ν, which is the decay mode with the highest branching ratio. In order to roughly match recent LHC searches for H/A → τ τ , taus are required 2 In linking this theory quantity directly with the H/A event mixture in the datasets we simulate, we make the approximation that the acceptance times efficiency is equal for H → τ τ and A → τ τ events. 3 The magnitudes of the up-type and down-type Yukawa couplings have the same tan β dependence in both the leptonspecific and the type-I THDM. Gluon-gluon fusion through a top loop is therefore by far the most important production channel for the scenarios considered here.
to have visible transverse momentum p T larger than 40 GeV and pseudorapidity less than 2.1. Further, we require the taus to be separated by ∆R = (∆φ) 2 + (∆η) 2 > 0.5, and that there are no more than two taus in the event which pass the p T selection. Events with muons or electrons with p T > 20 GeV are rejected. Detector effects are taken into account by randomly smearing the directions and energies of the outgoing pions, following the procedure described in [18]: Each track is deflected by a random polar angle θ, which is drawn from a Gaussian distribution with width σ θ , so that the smeared track lies within a cone around the true track direction. For charged pions a value of σ θ = 1 mrad is used, while the energy resolution is ∆E/E = 5 %. For neutral pions, we use σ θ = 0.025/ √ 12 rad and ∆E/E = 10 %. To gauge the impact of such detector effects on our results, we repeat the main analyses in sections 4 and 5 for simulated data with and without detector smearing.

Network input features
For the neural signal mixture estimation in section 5, we train a network to separate H → τ τ events from A → τ τ events. The four-momenta of the visible tau decay products (π ± and π 0 ) constitute the most basic kinematic input features to our network. The momenta are boosted back to the visible ditau rest frame (the zero-momentum frame for the four pions) and rotated so that the visible taus are back-to-back along the zaxis. The system is then rotated a second time, now around the z-axis, so that the x-component of the π + is zero. This is done in order to align all events to a common orientation, as the azimuthal angle around the z-axis carries no physics information.
In addition to the pion momenta, the network is trained on missing transverse energy (E miss T ); the invariant mass of the four-pion system (m vis ); the transverse mass (m tot T ); the impact parameter vectors of the charged pions, which help constrain the neutrino directions; the pion energy ratios Υ ± , defined as and the angle ϕ * between the tau decay planes. For ϕ * we follow the definition in [18], 4 which uses the direction p The quantities in (a) and (b) are momentum components of the π − and π 0 from the τ − decay, after each event has been boosted back to the visible ditau restframe and rotated such that the taus are back-to-back in the z direction and the x-component of the π + momentum is zero. Panel (c) shows the transverse mass m tot T , defined in [31]. The observables Υ (d) and O * (e), defined in eq. (2) and eq. (4), respectively, are required for the computation of ϕ * , along with the momentum vectors of the tau decay products. The distribution of ϕ * is shown in (f). The green graph below each plot shows the ratio of the A-event and H-event distributions.
corresponding π ± , to form an intermediate observable ϕ ∈ [0, π) and a CP -odd triple correlation product O * , From these, we can define an angle continuous on the interval [0, 2π): The distribution of ϕ depends on the sign of the product Υ + Υ − ; in the case of Υ + Υ − ≥ 0, the distribution is phase-shifted by π relative to the case of Υ + Υ − < 0.
To incorporate this into a single consistent CP -sensitive observable, we define ϕ * as Before being input to the network, all feature distributions are standardised to have zero mean and unit variance. A selection of the feature distributions in the training data is shown in fig. 1. The univariate feature distributions are severely overlapping for H and A events, indicating that the classification task is very challenging. The one feature which stands out here is ϕ * , which is the basis for the single-variable mixture estimation described in section 4.
For the results presented in section 5 we use a network trained on all features discussed above. However, features such as ϕ * and m T are derived from the basic pion momenta that the network also has access to. These "high-level" features can in principle be inferred by the network itself from the "low-level" pion momenta. To briefly investigate this we repeat the network training with varying subsets of the input features, starting with only the pion four-momenta and sequentially adding ϕ * , Υ ± and the remaining features. For all networks we obtain ROC AUC scores of ∼ 0.630. While a full statistical comparison of the resulting networks is beyond the scope of our study, this indicates that the network is itself able to extract the relevant information from high-dimensional correlations between the pion momenta, making the explicit inclusion of the high-level inputs mostly redundant. We note that this observation is in agreement with the results of [1,2].
It is still interesting to investigate how much of the discriminatory power can be captured by the high-level features alone. For this we train several classifiers on high-level features only, adding a new set of features for each classifier. The first classifier is trained only on ϕ * and achieves a ROC AUC score of ∼ 0.605. When Υ + and Υ − are included as input features the performance improves to a score of ∼ 0.618. This improvement can be understood qualitatively from the fact that the difference between the Υ ± -conditional ϕ * distributions for H and A events increases with |Υ ± |. Adding E miss T , m T , π ± impact parameter vectors and O * raises the ROC AUC score to ∼ 0.620, and finally including m vis further increases the score to ∼ 0.623, which seems to be the limit for our network when trained on high-level features only. This indicates that ϕ * and Υ ± together capture most of the sensitivity, but that the neural network is able to extract from the pion four-momenta some additional information which is not contained in the high-level quantities. Similar behaviour was seen in [19] in a study focusing on the CP -nature of the 125 GeV Higgs.

Network layout
In this study we employ a fully-connected feed-forward network. The input layer has 26 nodes, followed by 500 nodes in the first hidden layer, 1000 nodes in the second hidden layer, and 100 nodes in the final hidden layer. These have leaky ReLU [32] activation functions, and dropout [33] is applied with a dropping probability of 0.375. No further regularisation is imposed. All network weights are initialised from a normal distribution, following the He procedure [34]. The output layer has a softmax activation function, and we apply batch normalisation [35] between all layers. The weights are optimised using Adam [36] with cross-entropy loss and an initial learning rate of 0.03. 20% of the training data are set aside to validate the model performance during training. If there is no improvement of the loss on the validation data for ten consecutive epochs, the learning rate is reduced by a factor ten. The network is trained for 100 epochs or until no improvement is observed during 15 epochs, whichever occurs first. The neural network implementation is done using the Keras [37] and TensorFlow [38] frameworks.

The ϕ * method
Traditional approaches for separating CP -even and -odd decays are based on the angle ϕ * between the tau decay planes, as defined in eq. (6). The ϕ * distribution for H and A events can be seen in fig. 2a. The CP -sensitive parameter in this distribution is the phase of the sinusoidal curve, which is shifted by π radians between the H and A hypotheses. We note that the distributions overlap across the full ϕ * range, hence no absolute event separation is possible based on this variable. Using the simplified notation p(ϕ * |A) ≡ p A (ϕ * ) and p(ϕ * |H) ≡ p H (ϕ * ), the ϕ * distribution for H/A signal data can be expressed as a simple mixture model, = α a cos ϕ * + c + 1 − α a cos(ϕ * + π) + c = α a cos ϕ * + (1 − α) a cos(ϕ * + π) + c, where we fix the amplitude a and offset c to a = 0.041 and c = 0.159, obtained from a separate fit to H and A training data. This leaves us with a model for the ϕ * distribution where α is the only free parameter. Given a dataset {ϕ * i } with N events, we can now obtain an estimateα for α by maximising the likelihood function We demonstrate this method in fig. 2 for a dataset with 100 H/A Pythia events, generated using a model with a true α of 0.7. The pdfs p H (ϕ * ) and p A (ϕ * ) are shown in fig. 2a, while the fit result is shown in fig. 2b. For this example the best-fit α estimate comes out atα = 0.74. To demonstrate the statistical performance of this estimator we repeat the fit using 10 000 independent test sets with 100 Pythia events each, generated with true α values of 0.5, 0.7 and 0.9. The resulting distributions of α estimates are shown in fig. 4a, where the purple (green) distributions depict results without (with) detector effects. By fitting a Gaussian to each distribution we find the spread in the estimates to be σ α = 0.27 (σ det α 0.45) when detector smearing is omitted (included). Further, the estimator is mean-unbiased for all three cases. Note that to demonstrate the unbiasedness we have allowed the fit to vary α beyond the physically valid range of [0, 1].

The neural network method
When estimating some parameter θ using collider data we ideally want to make use of the multivariate density p(x|θ) for the complete set of event features x. 5 However, it is typically infeasible to evaluate this density directly for a given x. A common approach is then to construct a new variable y(x) and base the parameter estimation on the simpler, univariate distribution p(y(x)|θ), as exemplified by the ϕ * fit in section 4.
The performance of such a univariate approach depends on how well the distribution p(y(x)|θ) retains the sensitivity to θ found in the underlying distribution p(x|θ). In the special case where the map y(x) is the output from a trained classifier, it can be shown that using p(y(x)|θ) to estimate θ in the ideal limit is equivalent to using the full data distribution p(x|θ). Here we briefly review this argument before applying the classifier approach to our mixture estimation problem.
After training on θ-labeled data, a classifier that minimizes a suitably chosen error function will approximate a decision function s(x) that is a strictly monotonic function of the density ratio p(x|θ)/p(x|θ ) [39]. 6 As shown in [40], the monotonicity of s(x) ensures that density ratios based on the multivariate distribution p(x|θ) and the univariate distribution p(s(x)|θ) are equivalent, If we now take θ to be a fixed value such that the support of p(x|θ ) covers the support of p(x|θ), 7 the maximum likelihood estimator for θ based on p(x|θ) can be rewritten as follows [40]: Hence, if the classifier output y(x) provides a reasonable approximation of s(x) we can expect the maximum likelihood estimator based on p(y(x)|θ) to exhibit similar performance to an estimator based on p(x|θ).
The main drawbacks of this approach are the complications associated with training the classifier, and that the physics underlying the parameter sensitivity may remain hidden from view. We now apply this classifier approach to our H/A mixture estimation problem. The maximum likelihood estimator for the mixture parameter α is then given by where we have expressed the overall network output distribution p(y|α) as a mixture of the pure-class distributions p(y|A) ≡ p A (y) and p(y|H) ≡ p H (y). We use a network trained on a balanced set of H and A events. The network is trained to associate outputs y = 0 and y = 1 with H and A events, respectively. By applying this network to another labeled dataset of equal size to the training set, we construct templates for the probability densities p H (y) and p A (y) in eq. (11) using a nonparametric kernel density estimation method (KDE) [41]. The resulting templates are shown in fig. 3a. We note that the pdfs do not span the entire allowed range y ∈ [0, 1]. This is expected, since the CP nature of a single event cannot be determined with complete certainty. Proper determination of the pdf shapes in the extremities -where the sensitivity is highestrequires a sufficient amount of data, which is why we devote a similarly sized data set to the template creation as to the network training.
Given a set of unlabeled data we can now estimate α by carrying out the maximization in eq. (11) as an unbinned maximum-likelihood fit. The resulting fit to the same example dataset as used for the ϕ * fit in fig. 2b is shown in fig. 3b. The best-fit α estimate in this case isα = 0.67.

Results
We can now compare the performance of the neural network method with that of the ϕ * method of section 4. To this end, we apply the network method to the same test sets as used in fig. 4a, i.e. 10 000 datasets of 100 Pythia events each, for each of the three scenarios α = 0.5, 0.7, 0.9. The analysis is repeated with network training and test sets with and without detector smearing. The results are given in fig. 4b, for easy comparison with the corresponding results of the ϕ * method in fig. 4a. We fit each distribution of α estimates with a Gaussian and summarize the fit parameters in table 1.
As for the ϕ * method, we find that detector smearing significantly impacts the width of the α distribution, which increases from σ α = 0.21 to σ det α = 0.37 upon inclusion of detector effects. Yet, the network approach consistently outperforms the ϕ * method, as σ α and σ det are reduced by ∼ 22% and ∼ 18%, respectively, compared to the ϕ * results. So while the absolute widths of the α distributions in fig. 4 illustrate that a dataset of 100 events is probably too small to obtain an accurate α estimate, the comparison with the ϕ * results indicates that the relative performance gain offered by the network approach is relatively robust against detector smearing.
Similar to the ϕ * method, the network method provides a mean-unbiased estimator. In order to demonstrate this we allow α to vary outside the physical range [0, 1] in our fits. However, for α > 1, the combined mixture model p(y|α) = αp A (y) + (1 − α)p H (y) will become negative for y values that satisfy p A (y)/p H (y) < (α − 1)/α. This we do not allow in our fits, and in such cases we lower the α estimate until p(y|α) is nonnegative everywhere. This choice explains the slight de-viation from Gaussianity in the region aroundα = 1.2 in the bottom right plot. 8 Figure 5 shows the distributions of α estimates for the cases of 20 events per test dataset (top row) and 500 events per test dataset (bottom row), where all sets have been generated with α = 0.7 and no detector smearing has been included. Compared to the results with 100 events per set, σ α for both fit methods increase (decrease) by approximately a factor √ 5 for the case with 20 (500) events per set, as expected from the factor 5 decrease (increase) in statistics. Thus, the relative accuracy improvement of the neural network approach over the ϕ * method remains approximately the same: 30% for the 20-events case, and 25% for the 500-events case. However, the absolute spread of estimates in the 20-events case shows that this is clearly not enough statistics to obtain a useful estimate of α.
As a cross-check of the behaviour of the network fit method, we plot in fig. 6a the distribution of the log-likelihood ratio −2 ln(L(α = 0.7)/L(α)) for all test datasets of our benchmark point with α = 0.7. According to Wilks' theorem [42], the distribution of this statistic should tend towards a χ 2 distribution with one degree of freedom. By overlaying a χ 2 distribution in fig. 6a we see that this is indeed the case. Thus, confidence intervals constructed from the log-likelihood ratio for a neural network fit should have the expected coverage. In fig. 6b we show the log-likelihood ratio curves for the example dataset used in figs. 2b and 3b. The narrowing of the log-likelihood parabola for the network method again illustrates the increase in precision over the ϕ * method.
For this study we focus only on the separation of two signal classes, not the separation of signal from background. Of course, a realistic dataset is likely to contain a significant fraction of background events. For the signal scenario studied here, the most important backgrounds are due to "fake taus" from QCD production, single Z production (pp → Z → τ τ ), double Z and W production (pp → ZZ/W Z/W W → τ τ + X) and top pair production (tt → W bW b → τ τ + X). While such backgrounds will degrade the absolute accuracy in the signal mixture estimate, it is likely to impact the ϕ * method more severely than the neural network method. With one or several background components in the mixture model, the network's ability to extract information from the many-dimensional kinematic space should allow it to differentiate the background components from the signal components better than what is possible with the ϕ * variable alone. We therefore expect a similar or better relative performance of the network method in the presence of background, compared to the results we have presented here. There are two ways to extend the network method to take into account additional components in the mixture model: either by implementing a multi-class classifier, or by training multiple binary classifiers on pairwise combinations of the model components. Based on [43] we expect the latter approach would give the best performance.

Conclusions
Estimating the component weights in mixture models with largely overlapping kinematics is a generic problem in high-energy physics. In this paper we have investigated how a deep neural network approach can improve signal mixture estimates in the challenging scenario of a ditau LHC signal coming from a pair of heavy, degenerate Higgs bosons of opposite CP charge. This is a theoretically well-motivated scenario within both general and more constrained Two-Higgs-Doublet Models.
We have studied a benchmark scenario with degenerate H and A states at m H = m A = 450 GeV. For this case we find that the neural network approach provides a ∼ 20% reduction in the uncertainty of signal mixture estimates, compared to estimates based on fitting the single most discriminating kinematic variable (ϕ * ). However, the improved accuracy of the neural network approach comes with a greater computational complexity.
The network method we have studied here can be extended to include additional mixture components, such as one or several background processes, either by training a multi-class classifier or by training multiple binary classifiers. To increase the available statistics, the method can also be extended to work with a wider range of tau decay modes, for instance by using the "impact parameter method" described in [18].
The code used to generate events, train the network and run the maximum likelihood estimates will be made available on gitlab.com/BSML after publication.  Fig. 4 Comparison of the distributions of α estimates using (a) the ϕ * method and (b) the neural network method, for test sets generated with α = 0.5 (top), α = 0.7 (middle) and α = 0.9 (bottom). The slight deviation from Gaussianity seen around α = 1.2 in the bottom right plot is due to the fact that we let α vary beyond [0, 1] in our fits, but still demand that that the mixture model p(y|α) is always non-negative. See the text for further details. Comparison of the distributions of α estimates using (a) the ϕ * method and (b) the neural network method, for test sets generated with α = 0.7. The top row shows results for test sets containing 20 events each, while the bottom row corresponds to test sets with 500 events each. The deviation from the Gaussian distribution seen at high α in the upper right plot is due to the same effect as discussed for fig. 4. Table 1 Summary of α estimation on 10 000 independent test sets with 100 events in each set, using the ϕ * fit and the neural network (NN) fit methods.