Quark-gluon discrimination in the search for gluino pair production at the LHC

We study the impact of including quark- and gluon-initiated jet discrimination in the search for strongly interacting supersymmetric particles at the LHC. Taking the example of gluino pair production, considerable improvement is observed in the LHC search reach on including the jet substructure observables to the standard kinematic variables within a multivariate analysis. In particular, quark and gluon jet separation has higher impact in the region of intermediate mass-gap between the gluino and the lightest neutralino, as the difference between the signal and the standard model background kinematic distributions is reduced in this region. We also compare the predictions from different Monte Carlo event generators to estimate the uncertainty originating from the modelling of the parton shower and hadronization processes.


Introduction
The current LHC search for strongly interacting supersymmetric particles in a multi-jet final state primarily relies on kinematic discriminants to separate the signal from very large standard model (SM) backgrounds [1][2][3][4]. The signal from heavy squarks or gluinos decaying to a light neutralino lies in the high visible and missing momentum tail. The hadronic jets in the supersymmetry (SUSY) signal come from the decay of a gluino or a squark to one or more quarks and a neutralino. On the other hand, for the dominant SM background of a single weak boson and multiple jets, the jets originate from the initial state radiation of both quarks and gluons. A natural question therefore is whether the difference in the substructure of the decay jets in the signal process and the radiation jets in the SM background can be utilized to further improve the searches. This difference is related to the discrimination of quark-and gluon-initiated jets [5][6][7][8][9][10][11][12][13][14][15][16], a topic being actively explored by the ATLAS and CMS collaborations [17][18][19].
One goal of the studies in this direction by ATLAS and CMS is to derive template distributions from data for observables that can separate quark-and gluon-initiated jets [17][18][19]. Such a data-driven approach can avoid uncertainties coming from the Monte Carlo (MC) modelling of the low energy hadronization component, and to a lesser extent of the parton shower implementation. Although the data based templates are still in an early stage of development, especially when employing multiple observables, which requires large statistics not yet reached at the LHC, it is worthwhile to briefly discuss the method and its comparison to MC predictions. The only input from MCs in this approach is the quark JHEP01(2017)044 and gluon-initiated jet fractions in two different processes, 1 for example, in the dijet and γ+jet or Z+jet process, computed at the Born level including parton shower effects. With this definition, depending upon the jet transverse momenta, the dijet event sample consists of 50 − 60% gluon-initiated jets, while the γ+jet or Z(→ e + e − )+jet events contain 70 − 80% quark-initiated jets. The observed normalized distribution in the data for a given observable in the two samples can then be used to derive the normalized distribution for a "pure" quark and gluon jet by solving a pair of linear equations in two variables (for each bin of the normalized distribution, and repeated for different transverse momenta and rapidity intervals). While the uncertainties coming from the parton distribution functions and MC implementation of the Born process and initial and final state radiation are small, the largest systematic uncertainty for such studies arises from the dependence of the templates on the processes being used to derive it [19]. This key aspect of process dependence requires further studies before data-driven templates can be employed in the long run for physics searches without bringing in additional large systematics.
Comparison of the templates derived from data and from the MCs following the same procedure described above shows that while the MC predictions for quark-initiated jets agree reasonably well with the distributions extracted from data, the distributions for the gluon-initiated jets differ [17][18][19]. The data based templates for gluon jets fall in most cases in between the predictions of different MCs (Pythia [20,21] and Herwig [23] to be specific). Such difference in MC prediction for the distribution of quark-gluon tagging observables have also been observed in phenomenological studies [8][9][10][11][12][13][14]16]. With this in mind, the usefulness of quark-gluon discrimination in physics searches at the LHC can be studied using existing event generators, and future use of data-driven templates is expected to lead to a performance somewhere in between the Pythia-and Herwig-based predictions. If promising improvements are found irrespective of the MC used, and after folding in additional systematic uncertainties in the background prediction from substructure variables, the physics case to pursue quark-gluon discrimination as a tool to find new physics at the LHC would be well motivated.
The goal of this study is to evaluate the expected improvement in the search for gluino pair production at the LHC by including the quark-gluon tagging observables to the standard supersymmetry search strategy in the multijet and missing transverse momentum channel. After including initial and final state parton shower effects to leading order matrix elements, it is estimated that while the third and fourth highest transverse momentum jets in gluino-pair events are expected to be quark-initiated, in the dominant V +jets (V = Z, W ) backgrounds, they are more likely to be gluon-initiated. This leads to a considerable improvement in the signal to background ratio, when jet substructure based observables are utilized. Moreover, including both the kinematic and the jet substructure observables within a multivariate analysis is found to enhance the search prospects further, especially when the mass difference between the gluino and the neutralino lies in an intermediate region. The projected improvement over standard kinematics based searches is observed independent of the MC generator used, though to a different degree.

JHEP01(2017)044
In section 2, we describe the quark-gluon separation variables used to define a multivariate discriminant, our Monte Carlo simulation of the signal and background processes as well as the kinematic selection of the signal region. We begin section 3 by first describing the expected quark-gluon fraction of jets in the signal and background processes based on truth level MC information. This is followed by a discussion on the distribution of relevant kinematic variables. The multivariate analysis procedure is described next, followed by the results on the boosted decision tree based separation of the signal and background jet substructure. Combining the information from both kinematics and jet substructure we obtain the signal and background likelihood distributions, which are then used to estimate the expected LHC search reach using different methods in the gluino-neutralino mass plane. We summarize our findings in section 4.

Overview of quark-gluon tagging variables
Based on the difference in splitting probabilities in a parton-shower picture, different possible variables have been proposed for quark-gluon discrimination, which essentially rely on the fact that a gluon produced with similar kinematics leads to a larger multiplicity of soft emissions compared to a quark, and a gluon-initiated jet is wider than a quark-initiated one. These differences follow from the higher colour charge-squared of the gluon, C A = 3, versus C F = 4/3 for a quark. As demonstrated in previous studies, based on both perturbative methods as well as MC simulations, it is found that the following variables lead to a better quark-gluon separation: 1. The number of charged tracks inside the jet cone (n ch ), with each charged track having p T > 1 GeV, where p T denotes its transverse momentum. Even though it is difficult to model this observable accurately by MC generators, the recent ATLAS studies on the charged track multiplicity distribution using 8 TeV LHC data shows reasonable agreement for a set of MC tunes upto very high jet transverse momenta [24]. We shall utilize such tunes in our study for both Pythia and Herwig MCs, as discussed in section 2.3.
2. Energy-energy-correlation (EEC) angularity [9] variables, for example, the observable denoted by C can be defined in terms of the charged track momenta as Here, the sums over i and j run over all the tracks associated to the jet with j > i, while β is a tunable parameter. As determined in previous studies [9], from perturbative calculations and MC simulations, β = 0.2 is found to be an optimal choice that maximizes the quark-gluon separation. The distance in the rapidity-azimuthal angle plane between the tracks i and j is denoted by ∆R(i, j).
3. Jet mass (m J ) scaled by its transverse momentum m J /p T,J .

JHEP01(2017)044
4. In addition to the above set of variables, as discussed in our previous study [10], the input for the number of softer reconstructed jets (associated jets) around a primary hard jet can also improve quark-gluon separation, since it captures additional information from radiation outside the jet radius not included in the above variables.
In this study we shall use n ch , C (β) 1 and m J /p T,J as the inputs to a multivariate discriminant for quark-and gluon-like jets. While the inclusion of associated jets can be helpful, it is challenging to do so in a multijet environment, as one needs to remove overlap with ISR jets. We leave the investigation of such overlap removal methods to a future study.
Since the variables n ch and C (β) 1 are defined in terms of charged tracks with p T > 1 GeV, it is expected that they would be less sensitive to pile-up effects, as long as all such charged tracks can be traced back to the primary interaction vertex. However, the jet mass variable is sensitive to pile-up contamination, and using groomed jets can be useful in this regard. A detailed study of pile-up subtraction or the impact of different grooming algorithms is, however, beyond the scope of our work.

Kinematic selection of signal region
The ATLAS and CMS searches define multiple signal regions determined in terms of kinematic selection criteria that can separate a SUSY squark or gluino production process from the SM backgrounds in the multijets+E T / channel. Even though this is a challenging analysis in an hadronic environment, for high squark-gluino masses the hard scale of the signal process is higher than the hard scale of most SM processes. This latter fact is reflected in the high values of sum of jet transverse momenta (H T ) or effective mass (M eff = H T + E T / ) demanded in the signal regions. Following the ATLAS search strategies for 14 TeV LHC [25], we first make a pre-selection of events based on the following cuts: Cut-1.
1. The number of jets, n j ≥ 4, with p j 1 T ≥ 160 GeV and p j 2 ,j 3 ,j 4 T ≥ 60 GeV. For all other jets we demand p j T ≥ 20 GeV. The rapidity coverage of the jets is determined by ATLAS calorimeter design, where the forward calorimeter covers the pseudo-rapidity range of |η| < 4.9. However, the tracker covers only upto |η| < 2.5, and therefore it is not possible to obtain the information on the number of charged tracks inside jets in the forward region. Since the quark-gluon discrimination variables can be more accurately determined in terms of charged track momenta, we therefore count n j only within |η| < 2.5.
3. Missing transverse momentum in the event E T / > 160 GeV.
The jet p T cuts and the E T / cut are applied at the matrix element (ME) level while generating the background events, which is modelled by Z(→ νν+)jets. Furthermore, in order to obtain a large statistics of events with a high M eff cut, we have generated several JHEP01(2017)044 different samples of the Z+jets events, one with each value of the M eff cut. As discussed in detail in section 2.3, we normalize our total Z+jets event rate by comparison with the number of events reported in the ATLAS simulation after the cuts in 4jm category [25] (defined as Cut-1 followed by E T / /M eff > 0.25 and M eff > 3200 GeV). With this, we are able to reproduce with a reasonable accuracy the ATLAS projected sensitivity in the Mg − Mχ0 1 plane for 14 TeV LHC with 300 fb −1 data.
In addition to the above basic set of cuts, in order to compare with the search reach of ATLAS 14 TeV projections [25], we have computed the signal and background event yields in seven different signal regions (4jl, 4jm, 4jt, 5j, 6jl, 6jm, 6jt) as defined in the ATLAS study [25], essentially differing in the values of the M eff , E T / /M eff and E T / / √ H T cuts.

Monte Carlo simulation of signal and background processes
For both the signal and background processes, the parton level matrix elements are computed, and the events generated using MG5aMC@NLO [26]. The parton level events are passed onto both Pythia 6.4.28 (with the P2012-RadLo tune) [20], and Herwig++ 2.7.1 (with the default tune) [23], for simulating parton shower, hadronization and underlying events. The above choice for the Pythia tune is based on better data-model agreement in a recent ATLAS study comparing the charged track multiplicity distribution in the data with MC predictions [24]. The parton shower and hadronization effects are simulated using two different MCs to estimate the uncertainty in quark-gluon tagging coming from MC modelling. The signal cross-section is normalized to predictions including the resummation of soft-gluon emission at next-to-leading logarithmic accuracy, matched to next-to-leading order supersymmetric QCD corrections [28].
We use the CTEQ6L1 [29] parton distribution functions from the LHAPDF [30] library, and the factorization and renormalization scales are kept at the default event-by-event choice of MG5aMC@NLO. Detector effects have been simulated using Delphes3 [31], where the jet clustering is performed with FastJet3 [33]. Jets are reconstructed using the anti-k T clustering algorithm [33,35] with radius parameter R = 0.4. We have implemented the variables used for studying quark and gluon jet tagging in the Delphes3 framework.
As the signal process, we consider gluino pair production, followed by its three-body decay with 100% branching ratio to a pair of quarks and the lightest neutralino, via intermediate off-shell squarks. In general, depending on the squark mass, on-shell squark production will also contribute to the same final state. However, for studying the usefulness of quark-gluon tagging in SUSY searches, a simplified model with only the gluino and the lightest (bino-like) neutralino is adequate, and the rest of the MSSM particles are assumed to be decoupled. The final state of interest will then be ≥ 4-jets and missing transverse momentum.
It is well-understood that the primary background to such a multi-jet and missing momentum search comes from Z+jets production (with Z decaying to neutrinos), followed by a similar contribution from W +jets (where the charged lepton from the W boson decay falls outside the tracker acceptance, and therefore is not reconstructed as a lepton). The fractional contribution of tt+jets and single top production is reduced at higher M eff regions, but it can also become comparable to the individual weak-boson contributions JHEP01(2017)044 depending upon the signal region of interest. A strong cut on the E T / variable reduces the QCD multijet background, especially by ensuring that the jet direction and the E T / vector direction are not correlated. For a comparison of different SM background contributions, see, for example, the recent ATLAS note on squark-gluino search at the 13 TeV LHC with 13.3 fb −1 of data [2]. Both the recent 13 TeV ATLAS analysis and the ATLAS projection results for 14 TeV LHC with 300 fb −1 data show that the total SM background in our signal region of interest (i.e., after Cut-1 and with M eff > 1.8 TeV) is always less than twice the Z+jets contribution. The kinematic and quark-gluon fraction properties in Z+jets and the subdominant W +jets processes are nearly identical. Therefore, we perform the MC simulations using only the Z+jets process, and take the total SM background as twice the Z+jets prediction, which is a conservative estimate.
Since we shall focus on a multivariate analysis (MVA) strategy especially for the quarkgluon separation, the statistics of MC events required to perform the boosted decision tree (BDT) training is very high, especially if the number of input variables to the BDT training is large (eventually we shall use a ten variable BDT). Furthermore, these event samples are all required to pass a pre-selection of Cut-1 and different values of high M eff cuts. Therefore, generating such a large statistics of events with matrix element (ME) -parton shower (PS) matching is beyond the scope of our computational resources. On the other hand, as is well-known, to obtain accurate predictions for the jet p T s in processes such as Z+jets, ME-PS matching is important. However, since we are primarily interested in four relatively hard and central jets, the expectation is that events based on Z + 3−jets or Z + 4−jets matrix elements followed by PS can cover the relevant phase space region, and therefore the normalized differential distributions should be well-predicted by these event samples. In order to check this fact, we generated three different samples of Z+jets events and compared all the kinematic and jet-substructure distributions between them. The three samples are: (1) Z+jets, ME-PS merged upto 4−jets, (2) Z + 3−jet ME followed by PS and (3) Z + 4−jet ME followed by PS. We find that all the distributions have very similar shape in the three samples (as shown in the appendix). Thus it is possible to obtain accurate normalized distributions by just using the Z + 3−jet ME (followed by PS) event sample, for which generating a large enough statistics is least resource intensive among the three. For the overall normalization, as discussed earlier, we normalize our Z+jets event yield to the number reported in ATLAS simulation [25], and take the total SM background as two times the Z+jets contribution.

MC truth level quark-gluon fraction
As discussed in section 2.3, in the signal process of gluino pair production, with gluino dominantly decaying via (onshell or offshell) squarks, the decay jets are all quark-initiated. In addition, there are additional jets in the signal events from initial state radiation (ISR), which may reduce the difference between the signal and background likelihoods if a gluoninitiated ISR jet is harder than the decay jets and also lies in the central region of the detector. At Born level, the dominant background of Z+jets has a higher gluon fraction JHEP01(2017)044  Table 1. Quark fraction (f q ) at the MC truth level for the first four highest-p T jets ingg+jets and Z+jets processes. All events are selected after passing the jet-p T , E T / (Cut-1) and M eff > 1.8 TeV cuts, at the 14 TeV LHC. See text for details on the determination of f q .
in the third and fourth highest p T jets (denoted by j 3 and j 4 respectively). It is thus expected that the maximum discriminating power in the likelihood would come from j 3 and j 4 , rather than the first and second highest p T jets (denoted by j 1 and j 2 ).
To define the MC truth level quark and gluon jet fraction, we adopt the following method. Assume that we are looking for quark jets in an event. In the first step we find quarks in the matrix element, and a quark of flavour f is denoted by f i . Next, in the parton history related to the mother parton i, we find the parton P i with the same flavour as f i (we choose the parton with the highest transverse momentum if there are multiple quark partons of flavour f ). Finally, if the distance between the jet J and the parton P i is less than the jet cone size, ∆R(J, P i ) < R = 0.4, we define the jet J as a quark jet. If not, then J is defined as a gluon jet. As pointed out in ref. [36], such a definition is not infrared-safe to all orders and therefore has some higher-order dependence on the infrared cutoff in the MC. However, we emphasize that in the actual study of signal-background discrimination, this definition does not play any role, since in that case, we compare the likelihood of an event being signal-like or background-like, based on an MVA with the discriminating variables as inputs.
For illustration, we show in table 1 the parton level quark fraction of the first four jets, as defined above. A representative signal point with Mg = 2000 GeV and Mχ0 1 = 1000 GeV has been chosen for table 1, and the quark fractions are shown after the preselection of Cut-1 and with M eff > 1.8 TeV. The parton shower MC used for this figure is Pythia 6.4.28. In general, we see from this table that among the first four hardest jets, most signal events contain 3 − 4 quark jets, while most Z+jets events contain 1 − 2 quark jets.

Inclusive and exclusive kinematic variables
In the dominant SM background processes of Z/W +jets, the jets come from initial state QCD radiation, which exhibits a strong ordering of the jet p T s for a given H T value, primarily because of the enhancement in the soft gluon emission probability given by the QCD splitting functions. On the other hand, for the decay jets coming from gluino decay, the jet transverse momenta are not in general so strongly ordered, as in this case the p T s of the jets are determined by the mass-gap between the gluino and the lightest neutralino and the mass of the lightest neutralino itself. Admittedly, this is then a SUSY parameter dependent statement as to how the ordered jet p T distributions would differ between the signal and the background. Nevertheless, for certain ranges of the gluino and neutralino masses, the transverse momentum of the first four jets, ordered according to their p T s, can carry additional information not entirely captured in the M eff or H T distributions. We use JHEP01 (2017)  the nomenclature of exclusive kinematic variables to refer to the ordered jet p T s, while we shall refer to M eff , H T and E T / as inclusive kinematic variables.
We show in figure 1, the normalized (to unit area) distributions for the kinematic variables used as inputs in defining the combined signal and background likelihood functions, after the event pre-selection of Cut-1 and an M eff cut of 1.8 TeV. For the signal, we show JHEP01(2017)044 the distributions at a benchmark point with Mg = 2000 GeV and Mχ0 1 = 1000 GeV, and for illustration results from only the Pythia MC are presented. Since only events passing Cut-1 and M eff > 1.8 TeV are included, the H T and p T j distributions have a non-standard shape (first rise to a peak value and then fall). As we can see from this figure, for this signal benchmark point, the exclusive kinematic variables also provide discriminating power over the Z+jets background. For the gluino pair production events, we have also checked that including additional jets in the matrix element and using ME-PS matching, the kinematic distributions do not show any significant difference.

Multivariate analysis
Using the quark-gluon separation variables described in section 2.1, namely, n ch , C (β) 1 and m J /p T,J as inputs, we first develop an optimized discriminant using a multivariate analysis. This has been carried out by employing a Boosted Decision Tree (BDT) algorithm with the help of the TMVA-Toolkit [38] in the ROOT framework [40,41]. The training of the BDT classifier has been performed using the Z + q and Z + g processes at the Born level. The MC samples for these processes are generated such that we obtain a uniform statistical coverage across the entire jet p T range of interest, and the BDT training is performed for different p T ranges taken as different categories.
Following the above method, for the signal and background processes, we compute the BDT score B i for each of the first four jets ordered according to their p T . This procedure has been carried out using both the Pythia6 and Herwig++ MCs to simulate the parton shower and hadronization aspects. In figure 2, we show the distribution of the BDT scores for the first four highest p T jets in the gluino pair signal and the Z+jets background processes (for illustration, the distributions are shown using Pythia6). As expected from the truth level quark-gluon fractions discussed in section 3.1, significant separation in the BDT scores for the third and fourth highest p T jets (B 3 and B 4 ) are observed, for which the signal jets are mostly quark-initiated, and the background ones are mostly gluon-initiated.
As a final ingredient to our analysis, we perform a further MVA study with ten input variables containing: {M eff , H T , p T,j1 , p T,j2 , p T,j3 , p T,j4 , B 1 , B 2 , B 3 , B 4 }. This defines a signal and background likelihood with all the kinematic and jet substructure information of the event. The BDT score cut is chosen to maximize the exclusion (or discovery) significance for a given model point. For illustrating the separation power from each subset of variables, we show in figure 3 the BDT score distributions obtained with the inclusive kinematic variables (M eff and H T ), the exclusive kinematic variables (p T,j1 , p T,j2 , p T,j3 and p T,j4 ), and the jet substructure based BDT variables (B 1 , B 2 , B 3 and B 4 ). We also show in the bottom right panel of figure 3 the signal-background separation with all ten variables included together in the MVA. The results are shown for the signal point (Mg = 2000 GeV and Mχ0 1 = 1000 GeV) and with Pythia MC. Based on the final BDT score distribution with ten observables, we can now obtain the ROC curve which shows the signal acceptance ( S ) versus background rejection (1 − B ) efficiencies as a function of the BDT cut. In figure 4, the red, green and cyan curves show the ROC curves for the MVA analyses based on the inclusive, exclusive and inclu- sive and exclusive variable sets combined respectively. These two sets carry independent information, and therefore the background rejection using the combined set increases compared to the ones using the individual sub-sets, and the S and 1/ B values on the cyan curve is roughly given by the product of their corresponding values on the red and green curves. For example, the efficiencies ( S , −1 B ) for the red and the green curves pass through (0.4, 10) and (0.4, 1.4 × 10) respectively, and the cyan curve passes through the product (0.4 2 , 1.4 × 10 2 ). The black solid and dashed curves show the performance of the MVA analysis with all the variables taken together, using Pythia6 and Herwig++ respectively. We find that, by adding the jet substructure variables to the MVA, the background rejection factor increases by about a factor of 4 in the Pythia results and by a factor of 2 − 2.5 in the Herwig results, for S ∼ 0.1. As we shall see in the next subsection, this latter improvement has a considerable impact while considering the exclusion (discovery) reach in the Mg − Mχ0

Projected reach in Mg − Mχ0
1 plane By varying the BDT score cut with all or a subset of observables as input, we can choose the cut that maximizes the search significance for a given SUSY parameter point. Here, we use the profile-likelihood method [42] to determine the 95% C.L. exclusion region in the Mg − Mχ0 1 plane. The likelihood function is defined as follows: where S and B denote the expected number of signal (for a given point in the parameter space) and background event yields with a particular integrated luminosity, and N obs represents the number of events observed in the corresponding search with the same luminosity. For determining the exclusion contours, we set N obs to be equal to the mean value of the expected number of background events B. The systematic uncertainty in the background prediction is taken into account by convoluting the Poission likelihood function with a Gaussian with mean B and variance σ B . Combining the different components of the systematic uncertainty, we set σ B = (δ Incl + δ Excl + δ Jsub ) × B, where δ Incl , δ Excl and δ Jsub are the fractional systematic uncertainties in JHEP01(2017)044 the background prediction coming from inclusive, exclusive and jet substructure observables respectively. For our significance computation we have set δ Incl = δ Excl = δ Jsub = 10%, and to obtain a conservative estimate of the reach we have added these uncertainties linearly, making the total systematic uncertainty in the background yield prediction to be 30%, when all the variables are included together. We introduced a nuisance parameter B to deal with the systematic uncertainty, which is profiled out by maximizing the likelihood function by varying B in the interval 0 ≤ B ≤ ∞.
In figure 5, we show the projected 95% C.L. exclusion contours in the Mg − Mχ0 1 plane at the 14 TeV LHC with an integrated luminosity of 300 fb −1 . The orange curve is the ATLAS projected sensitivity with standard kinematic cuts (as reproduced by us), while the red, green and blue solid lines show the reach with each subset of variables described in the previous subsection. As expected from the ROC curve in figure 4, each of the subsets individually can lead to similar reach in this parameter space. We recall that all the curves include the effect of the pre-selection cuts on the jet p T 's and E T / (Cut-1) as well as a high M eff cut. Thus these improvements are within a high mass signal region. It is further observed that on including the information of the ordered jet p T s of the first four jets the reach improves to a good extent (cyan solid curve). Finally, if we now include the jet substructure information as well, the reach in the Mg − Mχ0 1 plane (black solid line) shows considerable improvement over the standard analysis. It should be noted in particular that especially in the region where the mass difference between the gluino and the neutralino falls in an intermediate range, the jet substructure observables provide stronger separation power. We also note that the signal benchmark point used to show JHEP01(2017)044 the various distributions in this study, namely, (Mg = 2000 GeV, Mχ0 1 = 1000 GeV) can be excluded at 2σ level only when the jet substructure variables are included in the MVA. Since we have also included additional systematic uncertainties in the background rate coming from the modelling of both the exclusive and jet substructure observables (upto 30% in total systematic uncertainty), our estimates for the improvement in the LHC reach should be conservative. It is thus promising that utilizing quark-gluon discrimination within an MVA including kinematic observables can considerably improve the LHC search prospects of strongly interacting SUSY particles.
In order to understand the uncertainty in the predictions from the MC modelling of jet substructure, we have performed the full analysis using both the Pythia6 and Herwig++ MCs. In figure 6 we show the 95% C.L. exclusion contours predicted by the two MCs using either only the jet substructure subset (blue curves) or the full variable set (black curves). For reference, the exclusion contours based on ATLAS cuts [25] are also shown (orange curves), and they are almost identical for Pythia6 and Herwig++. The Pythia6 exclusion contours (solid lines) show a better reach than the Herwig++ ones (dashed lines), and the difference between the two essentially comes from the jet substructure modelling, which, as remarked earlier, differs significantly for gluon jets. It is however encouraging that both MCs predict significant improvement over the standard analysis. Thus to the extent these two MCs provide an estimate of the uncertainty in prediction, our results show that irrespective of such differences, an improvement is expected in the LHC reach of gluino pair production, especially in the intermediate mass gap region, when we include the quark-gluon separation information within the MVA analysis. Future availability of databased templates and improved MC tunes are expected to lead to more reliable predictions and a reduction of the systematics in the application of quark-gluon discrimination. Figure 6. The 95% C.L. exclusion contours predicted by Pythia6 (solid lines) and Herwig++ (dashed lines) using either only the jet substructure subset (blue curves) or the full variable set (black curves). For reference, the exclusion contours based on ATLAS cuts [25] are also shown (orange curves), and they are almost identical for Pythia6 and Herwig++.

Summary and outlook
Quark-gluon discrimination is becoming a topic of growing interest, both in the theoretical and Monte Carlo front with improved jet substructure based observables being designed to capture the detailed pattern of QCD radiation, and on the experimental front with the development of data-based templates for tagging observables as well as validation of existing MC tunes. It is thus an ideal juncture when the importance of quark-gluon jet separation methods in the search for physics beyond the standard model should be thoroughly explored. With this goal in mind, in this paper, we studied the impact of including quarkand gluon-initiated jet discrimination in the search for gluino pair production events at the LHC. As seen in table 1, when ordered according to their transverse momenta, the third and fourth jets are more likely to be quark-initiated for the signal process, while for the dominant background of Z/W +jets, they are more likely to be gluon-initiated. With the quark and gluon separation variables of the number of charged tracks, energy correlation functions C β 1 , and jet mass (m J /p T,J ) as inputs to a multivariate analysis, we first develop a BDT-based quark-gluon discriminant across a large range of jet p T using the Z + q and Z + g processes as the training samples. In addition to the standard "inclusive" kinematic variables of E T / , H T and M eff , we also observe that for a given H T value, there is a strong ordering of the jet p T s for the Z+jets background process, while for the signal process the jets are not so strongly ordered. This is of course a parameter point dependent statement, as the gluino-neutralino mass splitting and the mass of the lightest neutralino determines the ordering of the p T s of the decay jets. However, in certain regions in the Mg − Mχ0 JHEP01(2017)044 of the inclusive, exclusive and jet substructure observables as MVA input variables to understand the importance of each category, and find that all three sub-categories, when added individually to a set of pre-selection cuts and a minimum effective mass cut (chosen according to the working point in the (Mg, Mχ0 1 ) plane), lead to a similar improvement in S/B. Consequently, compared to an optimized kinematic-category based search (as currently carried out by the ATLAS and CMS collaborations), inclusion of the quark-gluon discrimination variables improves the reach in the Mg − Mχ0 1 plane, especially in a region where the difference between Mg and Mχ0 1 falls in an intermediate range. This is because for such intermediate mass gaps, the H T and M eff distributions in the signal can become similar to the SM background ones. Given the fact that the jet substructure based variables, as well as the inclusive and exclusive kinematic distributions can bring in additional systematic uncertainties in the background rate determination, we have included a total systematic uncertainty of 30% on our background event yield, which should be a reasonable estimate.
As discussed in the introduction, there exist differences in the Monte Carlo prediction of the quark-gluon separation observables, and the data-based templates for gluon-initiated jets tend to lie in between the predictions of the Pythia and Herwig MCs, while for quarkinitiated jets the data-based templates largely agree with the MCs. With this observation in view, we carry out our complete analysis using both the MC event generators, in order to get an understanding of the variation in signal and background rates from MC modelling of parton shower and hadronization processes. This translates into a variation in the expected reach in the Mg −Mχ0 1 plane as well. While the expected improvement in reach does depend upon the MC generator used, the generic patterns remain the same. The reach based on different sets of kinematic variables are similarly predicted by both the event generators, as largely expected, since the low energy hadronization component does not enter in the jet transverse momentum distributions, while the effect of parton shower variation is weaker if we focus on high-p T jets only. Therefore, the MC variation almost entirely originates from the modelling of the jet substructure. It is however encouraging that independent of the MC generator used, the inclusion of quark-gluon discrimination leads to an improvement in probing the gluino pair production process, especially in the intermediate mass-gap region. This fact, combined with the future prospect of obtaining data-driven multivariate templates that do not rely on the MC modelling of the hadronization component (and possible improvements in the MC tunes as well), makes the utilization of quark-gluon discrimination in new physics searches sufficiently promising. We therefore expect that it would be explored in further detail by the LHC experimental collaborations in the future search for strongly interacting supersymmetric particles.

A MC simulation of the Z+jets background process
In this appendix, we discuss the details of our simulation of the Z+jets background process. As discussed in section 2.3, due to the necessity to generate several large statistics event samples as an input to the MVA after Cut-1 and different values of high M eff cuts, we use the Z + 3−jets ME (followed by PS) event sample, since it accurately reproduces JHEP01 (2017) Figure 7. Normalized distribution of the inclusive and exclusive kinematic variables using three different event samples: (1) Z+jets, ME-PS merged upto 4−jets, (2) Z + 3−jet ME followed by PS and (3) Z + 4−jet ME followed by PS . The distributions are presented after Cut-1 and an M eff cut of 1.8 TeV, and for this reason the H T and p T j distributions have a non-standard shape.
normalized differential distributions for all the input variables, and is less resource intensive. The overall normalization of the Z+jets background is fixed by comparison with the ATLAS simulation results in the 4jm category [25]. This method is also cross-checked by reproducing to a good accuracy the ATLAS projected exclusion contour [25]. In figures 7 and 8, we show the normalized distributions of the inclusive, exclusive and jet substructure variables utilized in this study for the three event samples of (1) Z+jets, ME-PS merged upto 4−jets, (2) Z + 3−jet ME followed by PS and (3) Z + 4−jet ME followed by PS. As we can see from this figure, the difference in shape between these three event samples is negligibly small. Normalized distribution of the jet substructure based BDT variables, using three different event samples: (1) Z+jets, ME-PS merged upto 4−jets, (2) Z + 3−jet ME followed by PS and (3) Z + 4−jet ME followed by PS. The distributions are presented after Cut-1 and an M eff cut of 1.8 TeV.