Investigation of Monte Carlo uncertainties on Higgs boson searches using jet substructure

We present an investigation of the dependence of searches for boosted Higgs bosons using jet substructure on the perturbative and non-perturbative parameters of the Herwig++ Monte Carlo event generator. Values are presented for a new tune of the parameters of the event generator, together with the an estimate of the uncertainties based on varying the parameters around the best-fit values.


Introduction
Monte Carlo simulations are an essential tool in the analysis of modern collider experiments. These event generators contain a large number of both perturbative and nonperturbative parameters which are tuned to a wide range of experimental data. While significant effort has been devoted to the tuning of the parameters to produce a best fit there has been much less effort understanding the uncertainties in these results. Historically a best fit result, or at best a small number of tunes, are produced and used to predict observables making it difficult to assess the uncertainty on any prediction. The "Perugia" tunes [1,2] have addressed this by producing a range of tunes by varying specific parameters in the PYTHIA [3] event generator to produce an uncertainty.
Here we make use of the Professor Monte Carlo tuning system [4] to give an assessment of the uncertainty by varying all the parameters simultaneously about the best-fit values by diagonalizing the error matrix. This then allows us to systematically estimate the uncertainty on any Monte Carlo prediction from the tuning of the event generator. We will illustrate this by considering the uncertainty on jet substructure searches for the Higgs boson at the LHC.
As the LHC takes increasing amounts of data the discovery of the Higgs boson is likely in the near future. Once we a e-mail: peter.richardson@durham.ac.uk b e-mail: d.e.winn@durham.ac.uk have discovered the Higgs boson, most likely in the diphoton channel, it will be vital to explore other channels and determine if the properties of the observed Higgs boson are consistent with the Standard Model. For many years it was believed that it would be difficult, if not impossible, to observe the dominant h 0 → bb decay mode of a light Higgs boson. However, in recent years the use of jet substructure [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] offers the possibility of observing this mode. Jet substructure for h 0 → bb as a Higgs boson search channel, was first studied in Ref. [5] building on previous work of a heavy Higgs boson decaying to W ± bosons [16], highenergy W W scattering [21] and SUSY decay chains [22], and subsequently reexamined in Refs. [8,15]. Recent studies at the LHC [23][24][25] have also shown this approach to be promising.
The study in Ref. [5] was carried out using the (FOR-TRAN) HERWIG 6.510 event generator [26,27] together with the simulation of the underlying event using JIMMY 4.31 [28]. In order to allow the inclusion of new theoretical developments and improvements in non-perturbative modelling a new simulation based on the same physics philosophy Herwig++, currently version 2.6 [29,30], is now preferred for the simulation of hadron-hadron collisions.
Herwig++ includes both an improved theoretical description of perturbative QCD radiation, in particular for radiation from heavy quarks, such as bottom, together with improved non-perturbative modeling, especially of multiple parton-parton scattering and the underlying event. In FOR-TRAN HERWIG a crude implementation of the dead-cone effect [31] meant that there was no radiation from heavy quarks for evolution scales below the quark mass, rather than a smooth suppression of soft collinear radiation. In Herwig++ an improved choice of evolution variable [32] allows evolution down to zero transverse momentum for radiation from heavy particles and reproduces the correct soft limit. There have also been significant developments of the multiple-parton scattering model of the underlying event [33,34], including colour reconnections [35] and tuning to LHC data [36].
The background to jet substructure searches for the Higgs boson comes from QCD jets which mimic the decay of a boosted heavy particle. Although Herwig++ has performed well in some early studies of jet substructure [25,37,38], it is important that we understand the uncertainties in our modelling of the background jets which lie at the tail of the jet mass distribution.
In addition we improve the simulation of Higgs boson decay by implementing the next-to-leading-order (NLO) corrections to Higgs boson decay to heavy quarks in the POWHEG [39,40] formalism.
In the next section we present our approach for the tuning of the parameters, which effect QCD radiation and hadronization, in Herwig++ together with the results of our new tune. We then recap the key features of the Butterworth, Davison, Rubin and Salam (BDRS) jet substructure technique of Ref. [5]. This is followed by our results using both the leading and next-to-leading-order matrix elements in Herwig++ with implementation of the next-to-leading-order Higgs boson decays and our estimate on the uncertainties.

Tuning Herwig++
Any jet substructure analysis is sensitive to changes in the simulation of initial-and final-state radiation, and hadronization. In particular the non-perturbative nature of the phenomenological hadronization model means there are a number of parameters which are tuned to experimental results. Herwig++ uses an improved angular-ordered parton shower algorithm [29,32] to describe perturbative QCD radiation together with a cluster hadronization model [29,41]. The Herwig++ cluster model is based on the concept of preconfinement [42]. At the end of the parton-shower evolution all gluons are non-perturbatively split into quarkantiquark pairs. All the partons can then be formed into colour-singlet clusters which are assumed to be hadron precursors and decay according to phase space into the observed hadrons. There is a small fraction of heavy clusters for which this is not a reasonable approximation which are therefore first fissioned into lighter clusters. The main advantage of this model, when coupled with the angularordered parton shower is that it has fewer parameters than the string model as implemented in the PYTHIA [3] event generator yet still gives a reasonable description of collider observables [43].
To tune Herwig++, and investigate the dependency of observables on the shower and hadronization parameters, the Professor Monte Carlo tuning system [4] was used. Professor uses the Rivet analysis framework [58] and a number of simulated event samples, with different Monte Carlo parameters, to parameterise the dependence of each observable 1 used in the tuning on the parameters of the Monte Carlo event generator. A heuristic chi-squared function is constructed where p is the set of parameters being tuned, O are the observables used each with weight w O , b are the different bins in each observable distribution with associated experimental measurement R b , error b and Monte Carlo prediction f b (p). Weighting of those observables for which a good description of the experimental result is important is used in most cases. The parameterisation of the event generator response, f (p), is then used to minimize the χ 2 and find the optimum parameter values. There are ten main free parameters which affect the shower and hadronization in Herwig++. These are shown in Table 1 along with their default values and allowed ranges.
The gluon mass, GluonMass, is required to allow the non-perturbative decay of gluons into qq pairs and controls the energy release in this process. PSplitLight, ClPowLight and ClMaxLight control the mass distributions of the clusters produced during the fission of heavy clusters. ClSmrLight controls the smearing of the direction of hadrons containing a (anti)quark from the perturbative evolution about the direction of the (anti)quark. Al-phaMZ is strong coupling at the Z 0 boson mass and controls the amount of QCD radiation in the parton shower, while Qmin controls the infrared behaviour of the strong coupling. pTmin is the minimum allowed transverse momentum in the parton shower and controls the amount of radiation and the scale at which the perturbative evolution terminates. PwtDIquark and PwtSquark are the probabilities of selecting a diquark-antidiquark or ss quark pair from the vacuum during cluster splitting, and affect the production of baryons and strange hadrons respectively.
Previous experience of tuning Herwig++ has found that Qmin, GluonMass, ClSmrLight and ClPowLight to be flat, and so it was chosen to fix these at their default values [29].
To determine the allowed variation of these parameters Professor was used to tune the variables in Table 1 to the observables and weights found in Appendix A in Tables 5,  6, 7 and 8. The dependence of χ 2 on the various parameters, about the minimum χ 2 value, is then diagonalized.
The variation of the parameters along the eigenvectors in parameter space obtained corresponding to a certain change, χ 2 , in χ 2 can then be used to predict the uncertainty in the Monte Carlo predictions for specific observables. In theory, if the χ 2 measure for the parameterised generator response is actually distributed as a true χ 2 , then a change in the goodness of fit of one will correspond to a one sigma deviation from the minima, i.e. the best tune. In practice, even the best tune does not fit the data ideally and nor is the χ 2 measure actually distributed according to a true χ 2 distribution. This means that one cannot just use Professor to vary the parameters about the minima to a given deviation in the χ 2 measure without using some subjective opinion on the quality of the results.
We simulated one thousand event samples with different randomly selected values of the parameters we were tuning. Six hundred of these were used to interpolate the generator response. All the event samples were used to select two hundred samples randomly two hundred times in order to assess the systematics of the interpolation and tuning procedure. A cubic interpolation of the generator response was used as this has been shown to give a good description of the Monte Carlo behaviour in the region of best generator response [4]. The parameters were varied between values shown in Table 1. The quality of the interpolation was checked by comparing the χ 2 /N df , where N df is the number of observable bins used in the tune, in the allowed parameter range on a parameter by parameter basis for the observables by comparing the interpolation response with actual generator response at the simulated parameter values. Bad regions were removed and the interpolation repeated leaving a volume in the 5-dimensional parameter space where the interpolation worked well. Figure 1 shows the χ 2 /N df distributions for two hundred tunes based on two hundred randomly selected event samples points for the cubic interpolation. The spread of these values gives an idea of the systematics of the tuning process showing that we have obtained a good fit for our parameterisation of the generator response.
The line indicates the tune which is based on a cubic interpolation from six hundred event samples. It is this interpolation which was used to vary χ 2 about the minimum to assess the uncertainty on the measured distributions. During the tune it was discovered that PSplitLight was relatively insensitive to the observables used in the tune. As such, PSplitLight was fixed at the default value of 1.20 during the tune and subsequent χ 2 variation.
Professor was used to vary χ 2 about the minimum value, as described above, determining the allowed range for the parameters. As five parameters were eventually varied, there are 10 new sample points-one for each of the parameters and one "+" and one "−" along each eigenvector direction in parameter space.
We follow the example set by the parton distribution function (PDF) fitting groups in determining how much to allow χ 2 to vary. Our situation is different to the PDF fitters in that we are using leading-order calculations with leadinglog accuracy in the parton shower, where they fit to next-toleading order calculations which gives better overall agreement with the observables used. Generally, PDF groups fit to fully inclusive variables, where as we have fitted to more exclusive processes and by nature, these are more model dependent, in particular hadronization.
In Refs. [45,46] these issues are explored in terms of PDFs and the allowed variation is related to a tolerance parameter T , where A tolerance parameter of T ≈ 10 to 15 is generally chosen for the PDF groups, where they are fitting to around 1300 data points. As our fit is likley to have a higher χ 2 than their fit due to the aforementioned reasons, and that we fit to a greater number of parameters, we will have a higher tolerance parameter. In our fit, we have 1665 degrees-of-freedom and we examined various changes in χ 2 , whilst considering the effects of the precision data from LEP. A variation of χ 2 /N df = 5, equivalent to T ≈ 90, seems, subjectively to keep the LEP data within reasonable limits while a variation of χ 2 /N df = 10, i.e. T ≈ 130 is too large. Anything less T ≈ 40 had very little variation and was therefore   Tables 2  and 3 respectively. The Professor tune was then compared with the internal Herwig++ tuning procedure [29] as not all analyses that are in the internal Herwig++ tuning system are available in Rivet and subsequently accessible to Professor. Looking at Fig. 4 it is found that PSplitLight at a value of 0.90 is favoured and gives a significant reduction in the χ 2 /N df . It was therefore decided to use the values obtained from minimisation procedure, but using the value of 0.90 for PSplitLight to maintain a good overall description of the data. The new  These error tune values can now be used to predict the uncertainty from the tuning of the shower parameters on any observable. In the next section we will present an example of using these tunes to estimate the uncertainty on the pre-

Jet substructure boosted Higgs
The analysis of Ref. [5] uses a number of different channels for the production of the Higgs boson decaying to bb in association with an electroweak gauge boson, i.e. the production of h 0 Z 0 and h 0 W ± . Reference [5] uses the fact that the Higgs boson predominantly decays to bb in a jet substructure analysis to extract the signal of a boosted Higgs boson above the various backgrounds. Their study found that the Cambridge-Aachen algorithm [47,48] with radius parameter R = 1.2 gave the best results when combined with their jet substructure technique. For our study, we used the Cambridge-Aachen algorithm as implemented in the Fast-Jet package [49]. Three different event selection criteria are used: (a) a lepton pair with 80 GeV < m l + l − < 100 GeV and p T > p min T to select events for Z 0 → + − ; (b) missing transverse momentum / p T > p min T to select events with Z 0 → νν; (c) missing transverse momentum / p T > 30 GeV and a lepton with p T > 30 GeV consistent with the presence of a W boson with p T > p min T to select events with W → ν; where p min T = 200 GeV. In addition the presence of a hard jet with p T j > p min T with substructure is required. The substructure analysis of Ref. [5] proceeds with the hard jet j with some radius R j , a mass m j and in a mass-drop algorithm: 1. the two subjets which were merged to form the jet, ordered such that the mass of the first jet m j 1 is greater than that of the second jet m j 2 , are obtained; 2. if m j 1 < μ m j and y = min(p 2 where R 2 j 1 ,j 2 = (y j 1 − y j 2 ) 2 + (φ j 1 − φ j 2 ) 2 , and p Tj 1,2 , η j 1,2 , φ j 1,2 are the transverse momenta, rapidities and azimuthal angles of jets 1 and 2, respectively, then j is in the heavy particle region. If the jet is not in the heavy particle region the procedure is repeated using the first jet.
This algorithm requires that j 1,2 are b-tagged and takes μ = 0.67 and y cut = 0.09. A uniform b-tagging efficiency of 60 % was used with a uniform mistagging probability of 2 %. The heavy jet selected by this procedure is considered to be the Higgs boson candidate jet. Finally, there is a filtering procedure on the Higgs boson candidate jet, j . The jet, j , is resolved on a finer scale by setting a new radius R filt = min(0.3, R bb /2), where from the previous massdrop condition, R bb = R 2 j 1 ,j 2 . The three hardest subjects of this filtering process are taken to be the Higgs boson decay products, where the two hardest are required to be b-tagged.
All three analyses require that: • after the reconstruction of the vector boson, there are no additional leptons with pseudorapidity |η| < 2.5 and p T > 30 GeV; • other than the Higgs boson candidate, there are no additional b-tagged jets with pseudorapidity |η| < 2.5 and p T > 50 GeV.
In addition, due to top contamination, criterion (c) requires that other than the Higgs boson candidate, there are no additional jets with |η| < 3 and p T > 30 GeV. For all events, the candidate Higgs boson jet should have p T > p min T . The analyses were implemented using the Rivet system [58].
The plots shown in Fig. 5 use the leading-order matrix elements for the production and decay of Higgs boson but the W , Z and top [50] have matrix element corrections for their decays. The plots shown in Fig. 6 have leading-order tt production, leading-order vector boson plus jet production (with the same matrix element corrections as the LO matrix elements) but the NLO vector boson pair production [51] and NLO vector and Higgs boson associated production [52]. In addition we have implemented the corrections to the decay h 0 → bb in the POWHEG scheme, as described in Appendix B. The signal significances are outlined in Table 4.  5 Results for the reconstructed Higgs boson mass distribution using leading-order matrix elements. A SM Higgs boson was assumed with a mass of 115 GeV. In addition to the full result the contribution from top quark pair production (tt ), the production of W ± (W + Jet) and Z 0 (Z + Jet) bosons in association with a hard jet, vector boson pair production (VV) and the production of a vector boson in association with the Higgs boson (V + Higgs), are shown Fig. 6 Results for the reconstructed Higgs boson mass distribution using leading-order matrix elements for top quark pair production (tt ), and the production of W ± (W + Jet) and Z 0 (Z + Jet) bosons in association with a hard jet. The next-to-leading-order corrections are included for vector boson pair production (VV) and the production of a vector boson in association with the Higgs boson (V + Higgs) as well as in the decay of the Higgs boson, h 0 → bb. A SM Higgs boson was assumed with a mass of 115 GeV The uncertainties due to the Monte Carlo simulation are shown as bands in Figs. 7 and 8. As there are correlations between the different processes the uncertainty is determined for the sum of all processes. Whilst it would be possible to show the envelope for the individual processes, this would not offer any information on the envelope for the sum of the Fig. 7 Results for the reconstructed Higgs boson mass distribution using leading-order matrix elements. A SM Higgs boson was assumed with a mass of 115 GeV. The envelope shows the uncertainty from the Monte Carlo simulation Fig. 8 Results for the reconstructed Higgs boson mass distribution using leading-order matrix elements for top quark pair production, and the production of W ± and Z 0 bosons in association with a hard jet. The next-to-leading-order corrections are included for vector boson pair production and the production of a vector boson in association with the Higgs boson as well as in the decay of the Higgs boson, h 0 → bb. A SM Higgs boson was assumed with a mass of 115 GeV. The envelope shows the uncertainty from the Monte Carlo simulation Table 4 The significance of the different processes for the leadingand next-to-leading-order matrix elements. The significance is calculated using all masses in the range 112-120 GeV processes which is the result of interest. In addition the uncertainty on the significance is shown in Table 4.

Conclusions
Monte Carlo simulations are an essential tool in the analysis of modern collider experiments. While significant effort has been devoted to the tuning of the parameters to produce a best fit there has been much less effort understanding the uncertainties in these results. In this paper we have produced a set of tunes which can be used to assess this uncertainty using the Herwig++ Monte Carlo event generator. We then used these tunes to assess the uncertainties on the mass-drop analysis of Ref. [5] using Herwig++ with both leading-and next-to-leading-order matrix elements including a POWHEG simulation of the decay h 0 → bb.
We find that while the jet substructure technique has significant potential as a Higgs boson discovery channel, we need to be confident of our tunes to investigate this with Monte Carlo simulations.
The error tunes and procedure here can now be used in other analyses where the uncertainty due to the Monte Carlo simulation is important.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribu-tion, and reproduction in any medium, provided the original author(s) and the source are credited.

Herwig++
The weights and observables used in the Professor tuning system are outlined in Tables 5, 6, 7 and 8. Table 5 Observables used in the tuning and associated weights for observables taken from [53] Observable Weight  Table 6 Observables used in the tuning and associated weights for observables taken from [44] Observable Weight   The NLO differential decay rate in the POWHEG [39] approach is wherē Here B(Φ m ) is the leading-order Born differential decay rate, V (Φ m ) the regularized virtual contribution, D i (Φ m , Φ 1 ) the counter terms regularizing the real emission and R(Φ m , Φ 1 ) the real emission contribution. The leadingorder process has m outgoing partons, with associated phase space Φ m . The virtual and Born contributions depend only on this m-body phase space. The real emission phase space, Φ m+1 , is factorised into the m-body phase space and the phase space, Φ 1 , describing the radiation of an extra parton. The Sudakov form factor in the POWHEG method is Fig. 9 The two real-emission processes contributing to the NLO decay rate where k T (Φ m , Φ 1 ) is the transverse momentum of the emitted parton.
In order to implement the decay of the Higgs boson in the POWHEG scheme in Herwig++ we need to generate the Born configuration according to Eq. (5) and the subsequent hardest emission according to Eq. (6). The generation of the truncated and vetoed parton showers from these configurations then proceeds as described in Refs. [29,52,54,55].
The virtual contribution for h 0 → bb was calculated in Ref. [56]. The corresponding real emission contribution, see Fig. 9, is gives Together with the definition, x 3 = x ⊥ cosh y, we obtain the Jacobian , (13) for the transformation of the radiation variables. We can then generate the additional radiation according to Eq. (10) using the veto algorithm [3]. To achieve this we use an overestimate of the integrand in the Sudakov form factor, f (p T ) = c p T , where c is a suitable constant. We first generate an emission according to Once the trial p T has been generated, y and φ are also generated uniformly between [y min , y max ] and [0, 2π], respectively. The energy fractions of the partons are obtained using the definition x 3 = x ⊥ cosh y, and x 2 using energy conservation. As there are two solutions for x 1 both solutions must be kept and used to calculate the weight for a particular trial p T . The signs of the z-components of the momenta are fixed by the sign of the rapidity and momentum conservation. Any momentum configurations outside of the physically allowed phase space are rejected and a new set of variables generated. The momentum configuration is accepted with a probability given by the ratio of the true integrand to the overestimated value. If the configuration is rejected, the procedure continues with p max T set to the rejected p T until the trial value of p T is accepted or falls below the minimum allowed value, p min T . This procedure generates the radiation variables correctly as shown in Ref. [3].
This procedure is used to generate a trial emission from both the bottom and antibottom. The hardest potential emission is then selected which correctly generates events according to Eq. (10) using this competition algorithm.