Fractal based observables to probe jet substructure of quarks and gluons

New jet observables are defined which characterize both fractal and scale-dependent contributions to the distribution of hadrons in a jet. These infrared safe observables, named Extended Fractal Observables (EFOs), have been applied to quark–gluon discrimination to demonstrate their potential utility. The EFOs are found to be individually discriminating and only weakly correlated to variables used in existing discriminators. Consequently, their inclusion improves discriminator performance, as here demonstrated with particle level simulation from the parton shower.


Introduction
A hadronic jet is produced from an initial parton via a sequence of perturbative QCD branching interactions (the parton shower), followed by the non-perturbative conversion of partons to the hadrons we observe in experiments (hadronization). A Markov chain description of the parton shower suggests the spatial distribution of partons will exhibit some fractal character [1][2][3][4][5][6], and this will be inherited by the final hadron distribution (invoking local parton-hadron duality [7]). However, true scale invariance of the hadron distribution within a jet is broken by the running of the branching probability, termination of the shower due to hadronization, and finite detector resolution. Here we define new observables to characterize jet branching structure, named Extended Fractal Observables (EFOs), which accommodate deviations from fractal structure through simple parametrizations. The idea is to apply box-counting techniques, used widely in the study of dynamical systems and scale invariant objects, to the substructure of QCD jets. Box counting has previously been employed in particle physics to calculate the fractal dimension of electromagnetic showers [8] for highly granular a e-mail: philip.coleman.harris@cern.ch calorimetric reconstruction. Here, we extend the generality and information content of this technique in our characterization of QCD jets.
The motivation for this study is twofold. Firstly, we would like to characterize the spatial substructure of jets into a set of new observables. Secondly, we would like to demonstrate the use of such observables in the discrimination of quark and gluon jets. Quark and gluon discrimination has long been used as a tool to enhance the sensitivity of signatures with additional quarks [9][10][11][12]. In particular, weak boson fusion induced Higgs-production is enhanced due to the distinct signature of two additional hard quark jets in the gluon-dominated forward region of the detector [9,[13][14][15][16][17][18][19][20][21]. Quark and gluon tagging are also expected to be useful for physics searches beyond the Standard Model, including the detection of supersymmetric particles [22,23]. Additionally, if well designed, these taggers can be further extended to the subjets of boosted boson signatures [24]. We demonstrate that modest improvements can be made to existing quark-gluon taggers by incorporating the new jet observables defined in this paper.
Finally, our construction of pixel-based jet observables resonates with the recent development of the jet image paradigm [25,26], in which the energy measured in each detector cell is interpreted as the intensity of a pixel in a 2D image. Within this approach, powerful machine-learning algorithms for classifying images have been brought to bear on a range of jet classification problems. This has included tagging boosted weak bosons [26,27], boosted top quarks [28], and heavy-flavors [29,30].
We define EFOs in the following section. In Sect. 3 we analyze the performance of these observables in quark-gluon discrimination, before concluding. Fig. 1 An illustration of the iterated box-counting procedure used to calculate fractal-based quantities on a set of points. The filled blue circles are the (η, φ) angular coordinates of the hadrons within a particular sample jet (in particular, this jet has total p T = 157 GeV, and 30 con-stituent hadrons). The box-counting is illustrated for four sample scales, corresponding to successively finer values of 0.2, 0.1, 0.067 and 0.05. The cells registering particle hits are highlighted with red shading

Extended fractal observables
The computation of the EFOs is performed on a jet by jet basis using a variation of the Minkowski-Bouligand (boxcounting) dimension, as follows.

Variable definitions
To define our variables we implement a two-stage recipe: firstly, the jet cone is divided in the familiar (η, φ) angular coordinates into a square grid of cells, each cell having sidelength . For a given scale , we count the number of cells N hits ( ) which register particle hits with a total transverse momentum greater than some pixel-level soft cutoff, in this study chosen to be p T > 1.0 GeV. This low energy cut represents a limiting threshold due to detector resolution. This counting is iterated over a range of scales, as is illustrated in Fig. 1. The second stage is to fit smooth functions to the variation of y = log N hits ( ) with x = log (1/ ), and to extract the parameters of the fit as a set of (correlated) jet observables, which we call Extended Fractal Observables (EFOs). This is a generalization of the traditional box-counting method, in which only linear functions y = mx + c are fitted, with the gradient m identified as the fractal dimension [8].
Indeed, in Fig. 2 there is no distinct region of linear scaling, as would be needed to extract a fractal dimension. Rather, log N hits ( ) levels off smoothly from large to small scales as saturation is approached, motivating a non-linear fit to extract whatever information this curve might encode about the jet. In particular, the hadronization region (i.e. at small ) obviously carries non-perturbative information sensitive to the flavor of the jet. The observed curves are distinct between quarks, gluons and b-quarks, as summarized in Figs. 2 and 3. This scaling is a fundamental property of QCD resulting from the differences in the splitting of quarks and gluons. Further measurements of this scaling allows for an alternative approach to extract QCD properties such as the strong coupling constant [32,33].
The generic plateauing curves in Fig. 2 can be fitted by almost any non-linear function (given a suitably restricted range in x), so we studied fit functions with at most three parameters, for speed and robustness of fitting. Fits were carried out simply by a binned χ 2 minimization of the chosen function. Example fit functions included the following: between quarks, b-quarks and gluons are predicted to be energy independent [31]. The small but non-zero slopes in this plot reflect the fact that box-counting at a given angular scale probes spatial information in addition to the rate of splitting at the corresponding energy scale 1. logarithmic fits of the form y = p 0 + p 1 x + p 2 log x. 2. quadratic fits: y = p 0 + p 1 x + p 2 x 2 . 3. hyperbolic tangent fits: The values of the best fit parameters { p i } for each fitting function constitute three possible sets of EFOs. For a polynomial in x = log(1/ ), like the quadratic fit function, the fit reduces to a matrix inversion and thus has a well-defined convergence. The other two parametrizations are not polynomials, hence we perform a χ 2 minimization.
Functions which actually saturate, such as the hyperbolic tangent parametrization above, are more physically motivated because they can model the saturation itself (asymptoting to the jet multiplicity). However, for the range of box scales used in our study (of width ≥ 0.05, -see 2.2 below), and for all but the lowest p T jets, the non-saturating fit functions also provide adequate models for the observed scaling. For the purpose of quark-gluon discrimination (see Sect. 3), the logarithmic fitting function was found to give the best discrimination performance of the three functions above (see Fig. 6 to compare the performance between the logarithmic and hyperbolic tangent fitting functions).

The range of box-counting scales
The range of angular scales has been chosen by paving the jet cone with a square grid of N × N cells, where the splitting scale N ranges in integer steps from 3 to 16. For each N , the angular scale is = 2R/N , where R is the jet radius, in this study R = 0.4. The coarsest scale chosen, corresponding to N = 3, is essentially the coarsest scale carrying potentially discriminating information (for N = 2 the jet cone would be divided into four quarters, all of which will register a hit for realistic jet shapes). The finest scale chosen is min = 0.8/16 = 0.05, because this is approximately the angular detector resolution in both LHC experiments, CMS and ATLAS [34,35]. For the p T ≥ 100 GeV jets studied here, the number of hits is just beginning to saturate at this scale (see Fig. 2), so we are probing into the hadronization region prior to the flat plateau.
Finally, we would like to highlight that these fractalbased observables are similar in spirit to calculating subjet rates of jets [15,36], given subjets clustered using the p T -independent Cambridge-Aachen algorithm [37]. Both observables compute p T -independent branching information on a succession of angular scales down to some threshold. And both observables perform what is essentially a further clustering on the substructure of the jet to extract this information pertaining to the branching history of the jet. In light of this, the EFO approach could be extended to utilize subjet counts (instead of hit grid cell counts) to assign scaledependent multiplicities N ( ).

Infrared and collinear safety
Preserving infrared and collinear safety ensure calculability in perturbative QCD. An observable is infrared (collinear) safe if its value is unchanged by the emission of soft (comoving) particles. The EFOs, as defined in 2.1 with a pixellevel soft cutoff, are fully IRC safe.
Firstly, the box counting procedure is intrinsically collinear safe: if one particle splits into two particles with the same (η, φ) coordinates, we still count just one cell hit by both daughter particles, at any finite scale of probing. Hence collinear splittings will not affect the number of cells N hits ( ) to register particle hits at any choice of scale. On the other hand, infrared safety of the EFOs can only be engineered by imposing some low momentum cutoff to cleanse the jet of its soft constituents. However, this soft cutoff must be implemented consistently with collinear safety. If we simply discarded all soft hadrons with, say p T < 1 GeV, this would spoil collinear safety. To see this, consider the following pathological example: if a particle with p T = 1.5 GeV splits into two comoving particles with p T = 0.8 GeV and p T = 0.7 GeV, then both would be discarded by a particle-level soft cut, and so N hits ( ) would not be invariant under this collinear splitting. This is remedied by defining a pixel-level (rather than particle-level) sort cutoff. That is, we only consider a cell to register a hit if it measures a total p T greater than our soft cutoff of 1 GeV. This way, if the troublesome 1.5 GeV particle in the example above splits collinearly into any number of daughters, the pixel still measures a total p T of 1.5 GeV, and so registers a hit regardless of these splittings. Thus, boxcounting with a pixel-level soft cutoff is fully IRC safe. In addition, a pixel-level rather than particle-level cut is more naturally realized experimentally since a pixel hit is consistent with an LHC detector cell.
Numerically, the performance of a quark-gluon discriminant built using the EFOs was found to be essentially insensitive to varying the value of this p T cut (over values between 0.1 GeV and 1.0 GeV), suggesting the variables are not strongly shaped by the IR emission, at least in simulations. In the following section, a p T cut of 1 GeV is used throughout. Finally, we acknowledge that pixel-level cutoffs have been used previously in the context of jet images analyses (for example in [25]) to ensure IRC safety in the same context.

Performance in Quark-Gluon discrimination
We now investigate whether these observables might be a useful new tool in the important and challenging problem of distinguishing light quarks from gluon jets.

Event generation and setup
In this study, we use QCD dijet samples at a center-of-mass energy of 13 TeV. Because previous quark-gluon studies have revealed that discrimination performance varies a lot between the different generators [9-11,14,38], 1 we here produce and shower events (at leading order) using both Her-wig++ (version 2.7.0 with tune UE-EE-5C ) [39,40] and Pythia 8 (version 8.185 with tune CUETP8M1) [41], with order 150k events in each. Jets are clustered with the anti-k T algorithm using the final state particles following showering and hadronization; a cone size of R = 0.4 and the FastJet code package [42] are used for the jet clustering. The EFOs (here computed using the logarithmic fitting function), along with a set of other established jet observables, have been computed for the highest p T jet in each event. We define the flavor of that jet by matching to the highestp T parton within R < 0.3 of the jet axis, and classify the event as signal (background) if matched to a light quark (gluon). 2 As a baseline for comparison, we shall consider the variables currently used by the Compact Muon Solenoid (CMS) quark-gluon tagger, which are [10]: (i) the total number of reconstructed particles in the jet (the multiplicity) [43]; (ii) the p T D variable (C β=0 1 ) [44],  where i sums over the constituents of the jet, which describes the distribution of transverse momentum between the particles in the jet; and (iii) σ 2 , the (p T -weighted) semi-minor axis of the jet in the (η, φ) plane [10], defined by where λ 2 is the smaller eigenvalue of the 2 × 2 symmetric matrix with components T,i Δφ 2 i , and M 12 = −Σ i p 2 T,i Δη i Δφ i . Throughout this study, we build multi-variable quark-gluon discriminants using a boosted decision tree (BDT), implemented using the Toolkit for Multivariate Analysis (TMVA) via adaptive boosting. The p T of the quark and gluon samples are reweighted to match the exact same kinematics in both cases, so as to avoid selection biases induced by kinematic differences in the simulation.

Results
We first compare the discriminator performance of single variables and the correlations between them, before going on to compare multi-variable taggers built with and without inclusion of the new EFO observables.
We can measure discriminator performance by receiver operator characteristic (ROC) curves, which plot background rejection against signal efficiency. Roughly speaking, the more convex the curve, the better the performance. The left plot of Fig. 4, made using the Herwig samples, shows that ROC curves for BDT discriminators constructed from various combinations of observables, as indicated by the legend, for events showered using both Herwig and Pythia with jet p T ≥ 100 GeV. The discrimination is superior in Pythia. We see in both event generators that including the EFOs rather than multiplicity (which is used in the CMS tagger) yields a marginally better performance the EFOs 3 are individually well-discriminating, particularly if we seek high signal efficiency. Their performance is significantly better than that of the jet multiplicity variable.  Fig. 6 Left: the relative gain for the three-variable taggers with respect to a baseline tagger using just p T D and σ 2 , for the Herwig events (which yield more conservative discrimination estimates). The gain is also plotted for EFOs computed with the hyperbolic tangent fitting function specified in Sect. 2.1, for which the performance is worse. Right: for Pythia events. Note the wider range of the y-axis, to accommodate the larger gains found in Pythia The right plot of Fig. 4 presents the linear correlation coefficients (calculated using the TMVA toolkit) between the EFOs and the existing CMS quark-gluon tagger variables: multiplicity, p T D and σ 2 . We also include a computation of the fractal dimension, which has been calculated from a linear fit over a small range of box scales. Strong correlations are present amongst the EFOs, as is natural given they are parameters derived from the same fit. However, their correlations with the other variables are no greater than 43% (for either quarks or gluons). 4 Interestingly, the EFOs are most highly correlated with σ 2 , not multiplicity as might have been expected. This evidence suggests the discrimination power of the EFOs is not simply a result of higher multiplicities in gluon jets, and therefore that the addition of these parameters to a quark-gluon discriminator might improve performance.
We find that replacing the multiplicity variable in the existing CMS quark-gluon tagger with the EFO variable yields a gain in discriminator performance, albeit only a modest one. This gain is seen using both Herwig and Pythia event generators (with the setup described above) in the ROCs presented in Fig. 5, which are for jets with p T ≥ 100 GeV. We see the performance in Pythia is significantly better than Herwig for each combination of variables, consistent with previous studies [9][10][11]14].
Moreover, the incremental gain upon replacing multiplicity with the EFOs is larger in Pythia than Herwig, so Herwig gives the more conservative estimate of the impact of including the EFOs. We see the gain in performance (relative to a baseline tagger using just p T D and σ 2 ) more clearly in Fig. 6, with the left panel for Herwig and the right for Pythia. The gain is at the level of 1-2% in the more conservative Fig. 7 Performance of a possible new quark-gluon tagger (using p T D, σ 2 , and the EFOs), in three p T bins, for Herwig-generated dijet events. Quarks and gluons are found to be easier to distinguish using this tagger at higher p T Herwig setup, and slightly larger in Pythia (note the different scaling of the y-axis). To emphasize a previous point, these gains were found to be stable across different values of the soft p T cut. Finally, we investigated how the performance varies with energy scale, by performing the analysis in p T bins of 50-100, 100-200, and 200-500 GeV. Discrimination was found to increase with p T in both Herwig and Pythia (see Fig. 7 for the Herwig results).
Combining all four variables (multiplicity, p T D, σ 2 and the EFOs) was seen to give no further improvement. This suggests all the information from multiplicity is captured by the EFOs, 5 while the converse is not true. In summary, we have presented evidence in this study that the Extended Fractal Observables provide an additional handle that captures the salient features of jet multiplicity, incorporates new information from showering and hadronization, and which is also better behaved under IRC emission (see 2.3).

Conclusions
In this study we defined new jet observables, the Extended Fractal Observables, by a generalization of the box-counting method used in the study of fractal systems. Defined with a pixel-level low momentum cutoff, these observables are infrared and collinear safe. We have then sought to apply the EFOs to improve quark-gluon discrimination. At the generator level, we find some modest improvement in discrimination by gluon rejection when we replace multiplicity with the EFOs in the existing CMS tagger, across both Herwig++ and Pythia 8. Extending the performance of these new variables to include detector effects can naturally be performed in the LHC environment with the CMS Particle Flow algorithm [45] in conjunction with the PUPPI algorithm [46] to reconstruct particle candidates in the presence of high pile-up.

Outlook
This method of studying jet substructure is a new approach. As such, there are many directions in which we would like to proceed, including: 1. Exploring particle hits in a 3-dimensional coordinate space spanned by η, φ and z −1 , where z is the fractional transverse momentum of the jet constituent. 2. Applying the EFOs beyond Quark-Gluon discrimination, for example to the identification of pile-up jets, or initial state radiation. 3. These box-counting methods extend very naturally from the substructure of a single jet to a whole-event analysis. Such a novel approach may provide new insight into searches for new physics topologies such as those in supersymmetry or top quark pair production [47]. 4. Furthermore, box-counting analyses could provide a useful characterization of event shapes in heavy ion collisions, where studies of jet properties beyond jet reconstruction are traditionally difficult, but well motivated [48][49][50]. 5. Finally, we would like to emphasize that the calculation of EFOs on quark and gluon jets probes parton shower scaling that results from the QCD color factor ratio. Calculating EFOs on cosmic ray air shower profiles [51] could therefore help discriminate QCD-induced air showers from more interesting signals; of particu-lar interest, showers induced by electroweak sphalerons. Experimentally, the calculation of EFOs in this air shower context is conceptually appealing: the 1660 individual Cerenkov detectors (spread over 3000 km 2 ) of the Pierre Auger Observatory in Argentina [52] would naturally function as the finest-scale cells in our box-counting algorithm. These techniques could therefore be useful in probing physics at energies far beyond that of the LHC.