Jet energy spectrum and substructure in $e^+e^-$ collisions at 91.2 GeV with ALEPH Archived Data

The first measurements of energy spectra and substructure of anti-$k_{T}$ jets in hadronic $Z^0$ decays in $e^+e^-$ collisions are presented. The archived $e^+e^-$ annihilation data at a center-of-mass energy of 91.2 GeV were collected with the ALEPH detector at LEP in 1994. In addition to inclusive jet and leading dijet energy spectra, various jet substructure observables are analyzed as a function of jet energy which includes groomed and ungroomed jet mass to jet energy ratios, groomed momentum sharing, and groomed jet radius. The results are compared with perturbative QCD calculations and predictions from the SHERPA, HERWIG v7.1.5, PYTHIA 6, PYTHIA 8, and PYQUEN event generators. The jet energy spectra agree with perturbative QCD calculations which include the treatment of logarithms of the jet radius and threshold logarithms. None of the event generators give a fully satisfactory description of the data.


I. INTRODUCTION
Jets, sprays of particles originating from the fragmentation and hadronization of a fastmoving quark or gluon, are some of the most useful tools for studying the Quantum Chromodynamics (QCD) and searching for new physics beyond the Standard Model (SM) in high energy colliders. Since the end of the Large Electron Positron Collider (LEP) operation in 2000, significant progress has been made on jet definitions and jet clustering algorithms.
Information about the parton shower is not only carried by the four-momentum of jets but also by their internal structure. Therefore, active developments on jet substructure observables have been carried out in order to fully exploit jets [1][2][3][4]. However, those jet reconstruction and jet substructure algorithms, which are widely used in the data analyses of proton-proton [5] and heavy-ion collisions [6][7][8], have not been used in the most elementary electron-positron e + e − collision system. Phenomenological models, such as Pythia [9], Sherpa [10], and Herwig [11], are tuned with hadron level or hadronic event shape observables from LEP experiments and are employed to predict or to describe the jet spectra and jet substructure in more complicated hadron-hadron collision environments.
Studies of jets, using identical clustering algorithms as those used in high energy hadron colliders such as the Large Hadron Collider (LHC) and Relativistic Heavy-Ion Collider (RHIC), are of great interest. Unlike hadron-hadron collisions, e + e − annihilations do not have beam remnants, gluonic initial state radiation, or the complications of parton distribution functions. The e + e − data provide the cleanest test for analytical QCD calculations and phenomenological models that are tuned with hadronic event shapes. Unlike the smooth and steeply falling jet transverse momentum spectra at RHIC and LHC, the inclusive jet and leading jet energy spectra in e + e − annihilation are peaked near the half of the collision energy. Therefore, leading and inclusive jet spectra are sensitive to the jet energy loss from wide-angle radiation that is not clustered into the jet. These spectra provide new tests of the calculation of leading-jet cross-sections, which is sensitive to the treatment of logarithms of the jet radius and threshold logarithms [12]. Moreover, fully reconstructed jets provide us with an opportunity to inspect quark and gluon fragmentation in great detail on a shower-by-shower basis. Finally, studies of jet substructure observables and their comparison to modern event generators are of great interest since jet substructure observables are novel tools for jet flavor identification, electroweak boson and top tagging, and studies of the Quark-Gluon Plasma properties at hadron colliders [13]. Those results can also serve as references for understanding proton-proton collision data and future electron-ion collisions (EIC) [14].
In this paper, the first measurement of inclusive jet momentum spectra, jet splitting functions, (groomed) jet mass, and groomed jet radius distribution in hadronic Z 0 boson decays are presented, with two types of jet selections. The inclusive observables include all jets above an energy threshold and inside a defined acceptance. They are sensitive to higher-order corrections of jet spectra in perturbative QCD. They include low momentum jets in events that are typically associated with soft gluon radiation or random combinatorial jets with hadrons from different partons. On the other hand, the leading dijet observables consider the leading and subleading jet in the event. This type of observable focuses more on the dominant energy flow and is less sensitive to soft radiation.

II. ALEPH DETECTOR
The ALEPH detector is described in ref. [15]. The central part of the detector is designed for the efficient reconstruction of charged particles. The trajectories of them are measured by a two-layer silicon strip vertex detector, a cylindrical drift chamber, and a large time projection chamber (TPC). Those tracking detectors are inside a 1.5 T axial magnetic field generated by a superconducting solenoidal coil. The charged particle transverse momenta are reconstructed with a resolution of δp T /p T = 6 × 10 −4 p T ⊕ 0.005 (p T in GeV/c).
Electrons and photons are identified in the electromagnetic calorimeter (ECAL) situated between the TPC and the superconducting coil. The ECAL is a sampling calorimeter, made by lead plates and proportional wire chambers segmented in 0.9 • × 0.9 • projective towers. They are read out in three sections in-depth and have a total thickness of around 22 radiation lengths. Isolated photons are reconstructed with a relative energy resolution of The iron return yoke constructed with 23 layers of streamer tubes is also used as the hadron calorimeter (HCAL) for the detection of charged and neutral hadrons. The relative energy resolution for hadrons is 0.85/ √ E (E in GeV). Muons are identified by their pattern in HCAL and by the muon chambers, made by two double-layers of streamer tubes outside the HCAL.
The information from trackers and calorimeters is combined in an energy flow algorithm [16]. This particle flow algorithm outputs a set of charged and neutral particles, called the energy flow objects. They are used in the jet reconstruction in this analysis.

III. DATASET AND EVENT SELECTION
The archived e + e − annihilation data at a center-of-mass energy of 91.2 GeV were collected with the ALEPH detector at LEP in 1994. To analyze these data, an MIT Open Data format [17] was created and validated in the two-particle correlation function analysis [18].
Hadronic events are selected by requiring the sphericity axis [19] to have a polar angle in the laboratory reference frame (θ lab ) between 7π/36 and 29π/36 to ensure that the event is well contained within the detector. At least five tracks having a minimum total energy of 15 GeV are also required to suppress electromagnetic interactions [20]. The residual contamination from processes such as e + e − → τ + τ − is expected to be less than 0.26% for these event selections [20]. Approximately 1.36 million e + e − collisions resulting in the decay of a Z 0 boson to quarks are analyzed.
High-quality tracks from particles are selected using requirements identical to those in previous ALEPH analyses [20] and are also required to have a transverse momentum with respect to the beam axis (p lab T ) above 0.2 GeV/c and | cos θ lab | < 0.94 in the lab frame. Secondary charged particles from neutral particle decays are suppressed by V 0 reconstruction in the energy flow algorithm [20]. In addition, the total visible energy from energy flow objects is required to be smaller than 200 GeV in order to remove laser calibration events during the run. Event thrust distributions [21] published by the ALEPH Collaboration using a similar dataset [19] were successfully reproduced within uncertainties, affirming that the archived data is analyzed properly.

IV. SIMULATED SAMPLES
Archived Pythia 6.1 [9] Monte Carlo (MC) simulation samples which were produced with the detector conditions during the 1994 run are used to correct for detector effects and to compare with the data. A set of new Pythia events are also generated with Pythia version 8.303 [22] at a center-of-mass energy of 91.2 GeV. The Monash 2013 tune [23] is used, with the weak boson exchange and weak single boson processes turned on. In order to select hadronic Z 0 decays, pure electroweak events are rejected by requiring at last one hadron in the final state particles. Sherpa samples are produced with version 2.2.5 [10], with electron positron events generating 2 to 5 outgoing partons, which are then showered into jets. The coupling constant α s (M Z ) is set to 0.1188 in the event generation. A set of Herwig [24] samples is generated using version 7.2.2. In this sample, 2 to 3 outgoing partons are specified in the hard process, and leptonic decays of Z 0 boson are turned off to increase the fraction of hadronic events. The order of α coupling is set to 2, whereas the colored α s order is set to 0 in the sample generation. In order to demonstrate the potential modification of the jet spectra and jet substructure observables from jet quenching [25][26][27], a sample is generated with the Pyquen generator [28] (version 1.5.3). The strength of jet quenching is set to be equivalent to that in a minimum bias sample of PbPb collisions at 5.02 TeV and wide-angle radiation is turned off during the sample generation.

V. JET RECONSTRUCTION AND PERFORMANCE
Experimental jets are clustered with the anti-k T algorithm [29] with a distance parameter R = 0.4, using the FastJet package [30] (version 3.3.2). This resolution parameter was chosen since it is widely used in the jet analyses in proton-proton and heavy-ion collisions [31][32][33][34][35] carried out at the hadron colliders. Moreover, the chosen value also gives us an opportunity to examine the shower from quarks in detail. The distance measure d ij between jets i and j and d iB between jet i and the beam is defined using the opening angle θ ij : where E i is the energy of the i-th jet and θ ij is the opening angle of the i-th and j-th jets.
The termination of the clustering process is defined by the distance between the jet and the beam direction, d iB . In the data analysis, energy flow objects which are reconstructed using the tracker and calorimeter information, are used for jet reconstruction. The mass of the objects are assumed to be pion for charged hadron candidates, and massless for photon candidates. Generator-level jets are clustered by considering all visible final state particles by the ALEPH detector (i.e., excluding neutrinos). Reconstructed-level jets are clustered with all energy flow candidates reconstructed by the tracker and calorimeters. In order to avoid overlapping with the beam pipe, only jets between 0.2π < θ jet < 0.8π are considered in this measurement, where θ jet is the angle between jet momentum and beam axis.
The jets are calibrated following a multi-stage strategy. In the first stage, simulated jet responses are corrected by using generated jets, clustered with visible final state particles with a distance parameter of R = 0.4, as the reference. The correction factor is derived as a function of jet 3-momentum p jet in 18 bins of jet direction θ jet . The closure in the simulated sample is seen to be better than 0.5%.
The jets in data are corrected with the first stage correction from simulation before deriving residual corrections characterizing the difference between data and simulation. In the second stage, there are two steps in the residual correction: relative residual correction, which equalizes the jet response between the two halves of the detector, and absolute residual correction, which calibrates the jet to an absolute scale. The jet relative correction is performed using the leading jet distribution in different bins of θ jet and comparing the two sides of the detector. This procedure aims to unify the jet energy scale in data between the two sides of the detector. The correction is typically of scale up to 0.5% in data, while the result on the simulated sample with the same procedure is consistent with 0.
A cross-check is carried out using the momentum balance of dijets in bins of third jet activity related to the leading dijet: . The residual scale is derived by extrapolating the third jet activity to zero. Due to the larger uncertainty induced by the extrapolation procedure, this procedure is not adopted as nominal. The corrections between the cross-check and the nominal one are consistent within the uncertainty.
The absolute correction of jets is performed using multi-jet events with an invariant mass around the Z 0 peak. Both the simulation correction and the residual relative correction are applied before the derivation of the absolute correction. The overall scale is fitted using up to N high energy jets (with N ranging from 3 to 9), and with different requirements on the energy of the (N + 1)-th leading jet from 1 to 5 GeV as a handle to control potential bias from the method. If there are less than N jets reconstructed in the event, the (N + 1)-th jet has 0 energy and the event is included. The correction is derived in bins of θ jet and is up to 0.3%.
The jet energy resolution in simulation is measured using an empirical function, in bins of θ jet : where a 0 to a 3 are parameters to be fitted. Depending on the direction of the jet, a 0 ranges from 0.005-0.010 GeV 2 , a 1 from 0.58-0.72 GeV 3 , a 2 up to 0.62 GeV 4 , while a 3 is negligible.
The difference between jet energy resolution between data and simulation is measured using the momentum balance of the leading dijet in an event. The magnitude of the thirdleading jet is used as a handle to understand potential systematic bias from the method, and it is constrained to be less than 3 GeV in the final result. The difference is characterized as a scale factor, in bins of θ jet . The resolution in data is found to be up to 5% worse than that in simulation, depending on the direction of the jet [36].

A. Jet Grooming
The soft drop algorithm [37,38] is proposed as a way to define a set of theoretically wellcontrolled observables that are sensitive to structures inside a jet and has been used in various QCD measurements in proton-proton and heavy-ion collisions [39][40][41][42][43][44]. It is characterized by two free parameters z cut and β. The constituents of a jet are first reclustered using the Cambridge-Aachen algorithm [35,45], and then the clustering history is traced, and the algorithm terminates when the soft drop condition is satisfied: where E 1 and E 2 are the energies of the two branches, and θ 12 is their opening angle. Note the definition is adapted to using the energy and opening angle (a true Great-Circle Distance), as opposed to using the transverse energy and ∆R ≡ (∆η) 2 + (∆φ) 2 as its counterpart in a hadronic collider environment. When the algorithm terminates, the z is defined as z G [46], and the opening angle is defined as R G .
In this measurement, in addition to the jet energy spectrum, the z G and R G spectra in bins of jet energy for z cut = 0.1, β = 0 are reported. For the z G and R G distributions, the distributions are self-normalized in order to decouple effects coming from jet energy migration and the rest.
In order to mitigate smearing effects from the detector, an unfolding is performed on the measured spectra. The unfolding is done with a Bayes' iterative method based on ref. [47], as implemented in the RooUnfold package. The method applies an iterative method with a corrected error calculation which decreases sensitivity to the prior knowledge used in the Bayes' formulation.
The ideal number of iterations for the unfolding is obtained with a two-step procedure.
The generator level distribution from the simulation is first combined with the smearing matrix to produce the idealized detector level distribution. Pseudodatasets are then generated using this distribution with the number of statistics compatible with that of the data. The ideal number of iterations can then be obtained by comparing the unfolded distribution with the generator level distribution. In the second step, data is first unfolded using the number of iterations found from the first step, and the whole procedure is repeated by using the unfolded distribution in place of the generator level distribution. The numbers range from 7 to 17 and the exact numbers can be found in ref. [36].
The response matrix for the jet energy measurement is shown in Figure 1. For all cases, the width of the pull distribution is consistent with 1.

C. Leading dijet selection
There is no detector acceptance close to the beamline, and therefore in order to measure only the event-wide leading and subleading jet, referred to as "leading dijet" in the later text, an additional selection is designed. The total energy visible either inside the acceptance 0.2π < θ < 0.8π or within radius R of any jet above 5 GeV is calculated, and events with total visible energy above 83 GeV are selected to ensure a purity of 99% of the events having leading dijet inside the acceptance. The efficiency of such a selection is about 60%.
A correction is derived using simulated events to correct for effects on the spectra due to this extra selection.

VII. SYSTEMATIC UNCERTAINTIES
The majority of systematic uncertainties are evaluated by unfolding the data with alternative versions of smearing matrices that encompass the systematic variations we would like to examine, with a few exceptions, as discussed individually below. Different sources are added in quadrature. These sources are discussed in the following paragraphs: Systematics related to hadronic event selections

Systematics related to jet reconstruction
The residual jet energy scale ("JEC") is varied by ±0.5% which is the maximum magnitude of the residual corrections. This uncertainty source ranges from a few percent up to 20% in jet energy spectra, and it is the dominating source of uncertainty at large E jet . For jet substructure results, due to the self-normalization, jet energy scale uncertainty, which manifests itself mostly as an overall normalization in each energy bin, is small.
The jet energy is smeared by ±2.5%, the difference between data and MC simulation extracted from dijet balance studies, to construct alternative response matrices in order to For the angular and energy resolution of the subjets, the associated uncertainties are determined in the following way. The directions of the subjet momenta are smeared randomly based on the observed simulated resolution. Although this choice is conservative, the overall effect is small. The energy of the subjets is also smeared based on twice the nominal energy resolution scale factor (+5%) obtained for the jets. The total uncertainty from these two sources amounts to 3% for z G and R G , and 5% for jet mass. The total uncertainties from substructure modeling and MC reweighting are shown in the "Model" category in the summary plots.
Contributions from fake or combinatorial jets are estimated by matching reconstructed jets in simulation to generator-level jets. The fraction of unmatched jets indicates the amount of potential fake jets. It ranges from 0-5% depending on the jet energy for inclusive jets.
A check on the results has been carried out by comparing the measured spectrum from each side of the detector (0.2π < θ jet < 0.5π and 0.5π < θ jet < 0.8π). The spectra are consistent with each other.

Systematics related to unfolding
We examine the sensitivity to the choice of prior which is part of the input to the unfolding. The difference between a flat prior and a prior with the generator level distribution from the simulation is taken as the uncertainty. It is grouped ("Unfold") with the uncertainty associated with the number of iterations done in the unfolding process. The number of iterations is varied from the determined ideal number of iterations to that optimized on a separate sample (ie., through Monte Carlo studies), and the difference in the unfolded distribution is taken as uncertainty.
In order to assess for potential bias in the unfolding method itself, an alternative unfolding method, the singular value decomposition (SVD) method [48], is used. The SVD unfolding method employs a regularization parameter, and the ideal parameter is determined using a method similar to that used to determine the ideal number of iterations in the nominal unfolding method. The difference between the SVD and the nominal is quoted as a systematic uncertainty.
A final source of systematic uncertainty related to the unfolding is the observed nonclosure on simulated events. This source of uncertainty is negligible for jet E and sizable, but still subdominant, for z G and R G when E is small.

Systematics related to leading dijet selection
The uncertainty on the extra leading dijet selection is evaluated by changing the thresholds on the total visible energy to those corresponding to a leading dijet purity of 98% and 99.5%. The difference to nominal is taken as the systematic uncertainty. The effect from potential imperfect modeling of the resolution of the total visible energy is investigated by smearing the quantity event by event and repeating the analysis. The difference is added in quadrature to the scale uncertainty.
Uncertainty that arises from the correction factor is evaluated through a reweighting procedure on the simulated sample. The events are weighted by the ratio of jet spectra for each jet between the nominal unfolded data and simulation, and the weighted simulated sample is used to derive the alternative correction factor. The difference is quoted as uncertainty.  The spectrum decreases rapidly as one moves to lower jet energy, reaching a minimum at around 20 GeV. Below that, the spectrum increases as we go to even lower energy due to the contributions from soft radiation and combinatorial jets. These features in the jet energy spectrum are captured by most of the event generators. Pythia6, which was tuned to describe ALEPH data, gives the best description of the energy spectrum. At low jet energy, Herwig gives the worse description of the results. It over-predicts the jet spectrum at low jet energy. Pyquen generator, which includes a large jet quenching effect equivalent to that in the inelastic PbPb collisions, predicts a large reduction of the population at around 45 GeV and a significant increase in the number of jets at low jet energy.
The same data could also be compared to a next-to-leading-order perturbative QCD calculation and a calculation based on a next-to-leading logarithmic (NLL') threshold and jet radius resummation (NLL'+R) [12]. The NLO calculation predicts a sharper peak at large jet energy. The NLL'+R calculation gives a better description of the inclusive jet energy spectrum, which shows the importance to include the jet radius resummation.
In order to focus on the dominant energy flow and suppress the contributions from soft radiation and combinatorial jets, the leading dijet energy is shown in Figure 13. Due to the leading dijet selection, the rise at low jet energy is suppressed. The results are consistent with predictions from Pythia6 within the quoted systematical uncertainties. Pythia8, Herwig, Sherpa tend to overpredict the spectrum at low jet energy. The dijet sum energy ( Figure 14), which is equivalent to the collision energy minus the sum of all the lower energy jets, is more sensitive to the modeling of subleading jets. The levels of (dis)agreement between the simulation and data for leading dijet energy and the leading dijet sum energy are similar.
Now we change our focus to the characterization of the substructure inside the anti-k T jets. The groomed momentum sharing z G spectra and the opening angle of subjets R G are presented in bins of jet energy. The distributions are self-normalized for each energy bin to factor out effects from an imperfect jet energy spectrum in the simulated samples. As shown in Figure 15, the z G spectrum is falling as a function of z G value, reaching a minimum at z G ∼ 0.5, which is similar to the data from proton-proton and heavy-ion collisions [39].
The jets which failed the grooming procedure are populated in the first bin z G = 0 in the plot. At high jet energy E > 40 GeV, Herwig predicts a flatter z G spectra than data compared with pQCD calculations at NLL'+R resummation [12] and over-predicts the jets with z G close to 0.5 while Pythia6 under-predicts the population at z G ∼ 0.5. This is similar to the results from pp collisions at a higher jet transverse momentum [39]. At low jet energy E < 40 GeV, the data are consistent with predictions from HERWIG 7. In most of the jet energy intervals, Pythia6, Pythia8 and Sherpa under-predict the number of jets at large z G by up to roughly 10% while they tend to over-predict the z G spectra at low z G .
The R G is the distance between the groomed subjets. The agreement between R G spectra and event generators, as shown in Figure 16, is worse than that observed in the z G spectra.
At high jet energy, event generators predict a slightly narrower core of the R G spectra compared to data. At low jet energy, most event generators predict on average a larger separation between subjects (larger R G ) than the unfolded data. The PYQUEN generator dijet energy spectra agree better with the predictions from the archived Pythia6 MC, which was tuned to describe the ALEPH data, none of the event generators gives a satisfactory description of the data. While the inclusive jet spectrum can not be described by partonlevel perturbative QCD calculation at next-to-leading order, calculations based on next-toleading logarithmic threshold and jet radius resummation could give a good description of the leading and inclusive jet energy spectra. These results provide the cleanest test of the perturbative QCD calculation of jets and serve as an important reference to the studies of jets in proton-proton, heavy-ion, and future electron-ion collisions.