Measurement of $b$-quark fragmentation properties in jets using the decay $B^{\pm} \to J/\psi K^{\pm}$ in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

The fragmentation properties of jets containing $b$-hadrons are studied using charged $B$ mesons in 139 fb$^{-1}$ of $pp$ collisions at $\sqrt{s} = 13$ TeV, recorded with the ATLAS detector at the LHC during the period from 2015 to 2018. The $B$ mesons are reconstructed using the decay of $B^{\pm}$ into $J/\psi K^{\pm}$, with the $J/\psi$ decaying into a pair of muons. Jets are reconstructed using the anti-$k_t$ algorithm with radius parameter $R=0.4$. The measurement determines the longitudinal and transverse momentum profiles of the reconstructed $B$ hadrons with respect to the axes of the jets to which they are geometrically associated. These distributions are measured in intervals of the jet transverse momentum, ranging from 50 GeV to above 100 GeV. The results are corrected for detector effects and compared with several Monte Carlo predictions using different parton shower and hadronisation models. The results for the longitudinal and transverse profiles provide useful inputs to improve the description of heavy-flavour fragmentation in jets.


Introduction
1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of Δ ≡ √︁ (Δ ) 2 + (Δ ) 2 .

Data and Monte Carlo samples
The dataset used in this analysis comprises the data taken from 2015 to 2018 at a centre-of-mass energy of √ = 13 TeV. After applying quality criteria to ensure good ATLAS detector operation [39], the total integrated luminosity used in this analysis is 139 fb −1 [40,41]. The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [40], obtained using the LUCID-2 detector [42] for the primary luminosity measurements. The measurements of the fragmentation properties are unaffected by this uncertainty, given that they are normalised to the total cross section, thus cancelling out the contribution of the luminosity uncertainty. The average number of inelastic interactions produced per bunch crossing for the dataset considered, hereafter referred to as 'pile-up', is = 33.6.
Several different models of multĳet production are used in this analysis. The MC samples were generated using the P 8 [43,44], S [45], and H 7 [46][47][48] generators. These models differ in the matrix element (ME) calculation, the parton shower (PS) and the hadronisation model (HM). The main features of these samples are summarised in Table 1. The P 8 samples were generated using P 8.240. The matrix element was calculated at leading order for the 2 → 2 process. The PS algorithm includes initial-and final-state radiation based on the dipole-style T -ordered evolution, including →¯branchings and a detailed treatment of the colour connections between partons [43]. The renormalisation and factorisation scales were set to the geometric mean of the squared transverse masses of the two outgoing particles (labelled 3 and 4), i.e.
Two different sets of hadronisation and underlying-event parameter values (tunes) were used in the generation of the P 8 samples. Two of the samples make use of the ATLAS A14 tune [49], for which the CTEQ6L1 PDF set [50] was used for the ME generation, the PS, and the simulation of multi-parton interactions (MPI). Two additional samples make use of the Monash tune [51], interfaced to the NNPDF2.3 PDF set [52] for the ME generation, PS and MPI. The hadronisation is modelled using the Lund string model [53,54] in all samples. The two samples using the A14 tune make use of the Lund-Bowler parameterisation of the fragmentation function [55], differing in the value of the parameter , which controls the shape of the fragmentation function for -quarks: while the first of them uses the nominal value = 0.855, the second one uses the so-called A14-tune, for which = 1.05, as obtained from a combined fit to LEP and SLD data [20][21][22][23]. The two samples using the Monash tune make use of the Peterson [56] and Lund-Bowler [55] parameterisations for the fragmentation functions, where the parameter is set to = 0.855 for the Lund-Bowler sample. The P 8 sample making use of the A14tune is taken as the nominal MC sample, used across the analysis for unfolding and uncertainty estimation.
The S samples were generated using S 2.2.5. The matrix element was calculated at leading order (LO) for the 2 → 2 process, and the default S CSS dipole PS [57,58]. Matrix element renormalisation and factorisation scales for 2 → 2 processes were set to the harmonic mean of the Mandelstam variables , and [59]. The CT14 [60] PDF set was used for the matrix element calculation, as well as for the modelling of the PS and MPI. The different S samples use different hadronisation models. The first one makes use of the S AHADIC model for hadronisation [61], which is based on the cluster hadronisation algorithm [62], while the second one uses the Lund string model [43,53,54,63] for the modelling of the hadronisation.
Finally, two H 7 samples were generated at leading order using H 7.2.1. The matrix element for the 2 → 2 process was calculated at LO with the MMHT2014 PDF [64]. The renormalisation and factorisation scales were set as 2 r = 2 f = 2 2 + 2 + 2 , where , and are the Mandelstam variables [59]. The first sample uses the default H 7 angleordered PS, while the second sample uses a dipole-based PS [65]. For both H 7 samples, the MMHT2014 [64] PDF was used for the modelling of the MPI, and the hadronisation was modelled by means of the default H 7 cluster hadronisation algorithm.
The decays of the mesons, of particular importance for this analysis, were modelled using the E G 1.6.0 generator [66] for the totality of the samples described above. The samples simulated using P 8 A14-, P 8 A14, and S with string hadronisation were passed through the G 4-based [67] ATLAS detector-simulation program [68] since they were also used to correct the measurements for detector effects, as described in Section 6. They are reconstructed and analysed with the same processing chain as the data. These fully simulated samples include the effect of multiple interactions per bunch crossing, simulated using P 8.186 interfaced to the A3 tune [69], as well as the effect on the detector response of interactions from bunch crossings before or after the one containing the hard interaction. In addition, during the data-taking, some modules of the ATLAS hadronic calorimeter were disabled for some periods of time. The resulting non-functioning regions in the hadronic calorimeter were not necessarily included in the simulation for all the samples. However, the effect of removing jets pointing towards these regions has been found to be negligible.

Object and event selection
Events are selected using triggers optimised for / meson identification and selection in its decay into muon pairs [70]. These triggers select events with two muons with T > 6 GeV. During the 2017 and 2018 data-taking periods, the trigger also required the invariant mass of the dimuon pair to satisfy 2 GeV < < 9 GeV and the angular separation between the muons to satisfy Δ < 1.5. The trigger selection efficiency is about 60%.
The primary vertex of the event is reconstructed as the vertex maximising the value of 2 T for all tracks originating from it. The reconstruction of ± mesons is done by re-fitting a pair of oppositely charged muons with T > 6 GeV and | | < 2.5 to a common vertex, with a fit quality fulfilling 2 < 10. Both muons are reconstructed using information from both the inner detector and the muon spectrometer. They must pass the M quality requirements [71,72], including at least one hit in the precision chambers of the muon spectrometer, and are required to have an invariant mass consistent with the mass of the / meson, i.e. to lie within the range 2.6-3.6 GeV. Moreover, both muons must spatially match the muon objects used in the trigger selection within Δ = 0.01. Given the short / lifetime following the ± → / ± decay, a third track with T > 4 GeV and | | < 2.5 is fitted to a common vertex together with the two muon tracks, requiring 2 / dof < 2.0. In the fit, the invariant mass of the pair of muon tracks is constrained [73] to match the world average value for the / mass, / = 3096.900 MeV [74]. Moreover, each of the three tracks in the triplet must have at least one hit in the pixel detector and at least four hits in the silicon microstrip detector, while the invariant mass of the triplet is required to lie within the range 5.0-5.7 GeV. Finally, the pseudo-proper lifetime of the ± candidate, = / T , where is the PDG value of the ± mass, = 5279.320 MeV [74], and is the transverse flight distance from the primary vertex, is required to be above 0.20 ps. The efficiency of this cut has been shown to be well described by the MC simulations, and hence to have only a small impact on the results after the corrections for detector effects are applied. The mass of the charged kaon, = 493.677 MeV [74], is assumed for the third track for the sake of calculating the three-body invariant mass.
Jet reconstruction is performed using particle-flow objects [75,76] as inputs to the anti-algorithm as implemented in FastJet [77, 78] with a radius parameter = 0.4. Muon tracks are not considered in the particle-flow algorithm and, thus, the energy of muons from the ± decay chain is not accounted for in the jet energy. In order to correct for the presence of muons inside the jets, muons passing the same quality requirements as those used in the ± reconstruction are matched to the jet by means of the ghost-association method [79]. This procedure defines objects with infinitesimal energy and the same direction as the track momentum (ghosts), which are then used as inputs to the jet reconstruction algorithm. A track is considered matched to a given jet if the ghost associated with it is clustered inside the jet. The jet four-momentum is then corrected by vectorially adding to it the four momenta of the muons associated with the jet, after subtracting the muon energy loss in the calorimeter. After this correction is applied, the jet energy is calibrated and, in order to ensure full inner-detector acceptance, jets with T > 20 GeV and | | < 2.1 are preselected. The jet calibration procedure includes energy corrections for pile-up, as well as angular corrections. Effects due to energy losses in inactive material, shower leakage, the magnetic field and inefficiencies in energy clustering and jet reconstruction are taken into account. This is done using a simulation-based correction, in bins of and T , derived from the relation of the reconstructed jet energy to the energy of the corresponding particle-level jet. In a final step, an in situ calibration corrects for residual differences in the jet response between the MC simulation and the data using T -balance techniques for dĳet, +jet, +jet and multĳet final states. In order to reject pile-up jets, the so-called 'jet vertex tagger' (JVT) algorithm is used [80]. Moreover, in order to avoid overlaps between jets (and thus, incorrect matchings with mesons), jets at a distance Δ < 0.8 from any other jet with T > 20 GeV are discarded.
The reconstructed ± mesons are then matched to the corresponding jet by requiring both objects to be within Δ = 0.4 of each other. If there is more than one meson candidate for a single jet, the one with lowest 2 to the three-track vertex is selected. If any of the three tracks arising from the ± decay is found to lie at Δ > 0.4 from the jet, its four-momentum is added to the jet four-momentum. The basic object for this analysis is, therefore, a jet containing a ± meson. The selection criteria described above are met by 1 413 684 jets in 1 404 620 events.
To correct the simulation for differences in the muon identification and JVT acceptance efficiencies from those measured in the data, simulated events are weighted using dedicated corrections provided to that end [71,80]. The differences in the trigger efficiencies between data and the simulation are accounted for as a systematic uncertainty, as described in Section 7. Figure 1 shows the detector-level distributions of the jet T and the T of the ± meson within the jet, together with the MC predictions. The data distributions include the purity corrections described in Section 5. While the S sample with string-based fragmentation describes the jet T within 10%, it fails to describe the meson T , with discrepancies of over 50%. The P samples with the A14 and A14-tunes give similar descriptions of the jet T , with differences from data of up to 40% in the tail of the distribution. The differences between the two tunes are clearly observed in the distribution of the meson T , where A14 gives a more accurate description within 10% while A14-shows discrepancies of up to 30% in the tail of the distribution. For the fragmentation measurement, the data are binned in three intervals of the jet T , namely [50,70), [70,100) and T ≥ 100 GeV. This choice ensures that the bin width is larger than the resolution in the different T intervals [75]. For each of these intervals, the distributions are binned in intervals with a width of 0.07, while the rel T distributions use a variable bin width. This choice allows finer binning, while retaining sufficient statistical precision in each bin, as needed for the estimation of the signal purity described in Section 5. The lower limit of the distribution is extended towards lower values with increasing jet T , where a larger fiducial phase space becomes available. The and rel T distributions are nor-malised to unit area by dividing by the number of entries, allowing the systematic uncertainties to be reduced.
Since the momentum resolution for charged particles is much better than the jet energy resolution, the measured distribution extends beyond unity. The values larger than one are absorbed as an overflow into the last bin of the distribution. Similarly, values of below the lower limit of the first bin are absorbed into the first bin as an underflow, and values of rel T above the upper limit of the last bin are absorbed into it as an overflow.

Signal extraction
The selected sample of ± meson candidates does not only contain ± mesons, but also backgrounds arising from different sources. For each ( T , ) bin, with T being the jet T and = or rel T , only a fraction ( T , ) of the reconstructed entries, referred to as the purity, are real ± mesons. The number of ± mesons, , reconstructed in a given bin can thus be expressed as where is the total number of candidates in a given bin. In order to determine ( T , ), a binned maximumlikelihood fit to the invariant mass distribution of the ± candidates is performed. The probability density function (pdf) of the model describing the invariant mass distribution can be written as a combination of the signal and background pdfs as where F , F , F and F are the pdfs for each of the components, signal or background, and , , and are coefficients representing their relative fractions, and thus ( T , ) = ( T , ). The closure relation = 1 must be satisfied by these coefficients. Each of the fit components is described below.
• The signal F , arising from the real ± meson contribution, is modelled by a double-Gaussian where is the relative normalisation of the two Gaussian components.
• The misreconstructed background, F , arising from the decays ±/0 → / * ±/0 → / ( ) ±/0 and ±/0 → / ( ) ±/0 , creates a low-mass structure, displaced from the ± mass. It is modelled by the following function • The resonant background F , arising from the decays ± → / ± , creates a peaking structure displaced from the ± mass towards higher masses. It is modelled by the sum of a Gaussian and an asymmetric Gaussian function, where is a relative normalisation factor and • The combinatorial background F , arising from random combinations of real / mesons with additional tracks, is modelled by a first-order polynomial All the parameters in Eq.
(1) ( 1 , 2 ,ˆ1,ˆ2,ˆ3, ) are obtained from fits to MC simulations, performed using P 8 A14-, of ± → / ± in bins of the jet T and independently of or rel T . The normalisation of this background, , is fixed relative to the signal as the ratio R ( , ) of the branching fractions of signal and background [74], i.e.
The fits are performed individually for each ( T , ) and ( T , rel T ) bin in the mass range [5.0, 5.7] GeV, and each of them has ten free parameters: four for the signal ( , 1 , 2 , ), two for the misreconstructed background ( , ), two for the combinatorial background ( 0 , 1 ), and two for the free normalisations of each component, . Figure 2 shows the results of such fits in four representative bins of the longitudinal profile and transverse profile rel T for the lowest and highest jet-T selections.
The quality of the fits is quantified using 2 functions, defined as 2 = ( − ) 2 /(Δ ) 2 for each fit. Here, is the value of the data distribution in bin ; Δ is its statistical uncertainty, and is the corresponding value of the fit model. The values of the 2 divided by the number of degrees of freedom, dof = 30, range from 0.6 to 2.2. The fit results are studied as a function of the jet T and of the fragmentation variables and rel T . The results show that the signal component dominates the fit models with an average, global purity close to 69%, followed by the misreconstructed background, with a global value of about 15% and the combinatorial background, amounting to 12% of the sample. Finally, the ± → / ± background amounts to 3% of the sample. As a general trend, the purity of the sample increases as a function of , while it is relatively constant as a function of rel T , particularly at low T . Figure 3 shows the signal purity and the fractions of the different backgrounds resulting from the fit as a function of and rel T , for the lowest and highest jet-T selections. The hatched bands include the sum in quadrature of the statistical and systematic uncertainties of the fit, as discussed in Section 7.

Unfolding to particle level
The detector-level distributions, from which the background processes have been subtracted following the methodology described in Section 5, are corrected for detector effects such as inefficiencies and migrations caused by the finite resolution of the detector. The particle-level observables are defined using requirements equivalent to the detector-level selection described in Section 4. Candidate ± → / ± decays are selected if the muons from the / decay have T > 6 GeV and the kaon has T > 4 GeV. The kinematics of the ± hadrons are reconstructed before QED radiation. Jets are reconstructed using the anti-algorithm with = 0.4, using all particles with average lifetime > 10 ps, including muons and neutrinos. Jets with T > 20 GeV and | | < 2.1 are selected, and those overlapping with another jet with T > 20 GeV within Δ = 0.8 are discarded. While the measurement is performed in the fiducial phase space with T > 50 GeV, an underflow bin with 35 GeV < T < 50 GeV is used to take into account the migrations from phase space regions below the jet T threshold. 2 The ± mesons are matched to jets if the angular distance between them fulfils Δ < 0.4. Finally, the four-momenta of any decay products of the ± meson lying outside the jet (Δ > 0.4) are added to the jet four-momentum.
The unfolding procedure can be parameterised using a transfer matrix , which contains the information about the bin-by-bin migrations from the detector-level distribution to the particle-level distributions, and is written as The unfolding can thus be regarded as a system of linear equations, for which the solution is the particle-level distribution . The transfer matrix , the efficiency E and the purity 3 P are determined using the P 8 A14-MC sample, and are the values of the detector-level distribution. An iterative matrix inversion method based on Bayes' theorem [81] is used, as implemented in the R U program [82]. The number of iterations is chosen so that the sum in quadrature of the statistical uncertainty and the uncertainty associated with the mismodelling in the unfolding, described in Section 7, is minimised. This results in four iterations for the unfolding of both the longitudinal and transverse profiles. The unfolding is performed in such a way that the migrations in both dimensions, between the fragmentation variable ( or rel T ) and the jet T are taken into account. The correlations in the statistical uncertainties, as well as those between the particle-level and detector-level MC distributions, are treated by using pseudo-experiments [83] in the calculation of each component of the unfolding, including matrices, efficiencies and purities as well as in the detector-level data.
The efficiency E is estimated as the fraction of particle-level jets, associated with a meson, which are matched within Δ = 0.4 to a detector-level jet that contains a reconstructed meson. Its value ranges from 10% to 30% depending on and rel T . The efficiency decreases as a function of for all T bins, although this behaviour is more pronounced for the higher T bin. This is due to the fact that, in this boosted regime, the dimuon trigger cannot efficiently separate the trajectories of the two muons produced in the / decay, and the events containing ± mesons are therefore not recorded efficiently. The values of E as a function of rel T are significantly flatter.
The purity P is estimated as the fraction of detector-level jets which are matched to a particle-level 2 Data falling in this bin are also corrected for background sources as described in Section 5. 3 This purity is unrelated to the purity described in Section 5 and corrected for before the unfolding. jet within Δ = 0.4. Its value ranges from 60% to 100%, systematically increasing as a function of both and the jet T . As a function of rel T , the purity exhibits a decrease in the low T region. This is due to the fact that reconstructed mesons flying in the same direction as the jet, as well as hadrons with high T , are more likely to arise from a true meson. In the high T region, P is constant as a function of rel T .
The values of the matrix elements are estimated using particle-level -meson jets, geometrically matched to the corresponding detector-level objects within Δ = 0.4. Figure 4 shows the two-dimensional transfer matrices as a function of ( T , ) and ( T , rel T ).

Systematic uncertainties
The systematic uncertainties in the measurement are classified into four main categories. The first of them concerns the identification of ± mesons, and includes the systematic uncertainties in muon reconstruction as well as the uncertainties in the fit procedure for the determination of the signal purity described in Section 5; the second category concerns jet reconstruction and identification, including the jet energy scale and resolution, the jet angular resolution and the efficiencies of the JVT algorithm. The third category is related to the unfolding procedure, including both the uncertainty in the MC model used for the corrections and the uncertainty due to mismodelling of the data by the nominal MC simulation. The fourth category takes into account systematic effects in the MC description of the pile-up dependence of the measurement.
The systematic uncertainties from all experimental sources are estimated by using the detector-level MC events. The resulting shifted distributions are then corrected for detector effects following the unfolding procedure described in Section 6, normalised to unit area, and compared with the nominal MC distributions, also normalised to unit area. The results of this comparison constitute the systematic uncertainties of the unfolded results for the normalised differential cross sections, and exploit the correlations between the numerator and the denominator when applying the normalisation. The full set of systematic uncertainties are described below, together with a detailed description of how they are estimated.

meson reconstruction uncertainties
• The uncertainties in the purity corrections for the ± candidates are evaluated by considering other fit models as alternatives to those described in Section 5. The probability density function for each component of the fit is substituted by an alternative function. For the signal model, one of the Gaussian functions in the double-Gaussian pdf is replaced by an asymmetric Gaussian function. The hyperbolic tangent shape in the misreconstructed-background pdf, 1 − tanh( ), is replaced by . For the ± → / ± background, one of the Gaussian functions in the double-Gaussian shape is replaced by a Crystal Ball function. Finally, the first-order polynomial describing the shape of the combinatorial background is replaced by an exponential function. The total uncertainty in the purity corrections is defined as the sum in quadrature of the deviations of the varied fit results, evaluated separately, from the nominal fit model described in Section 5. The impact of this uncertainty ranges from 1% to a maximum of 17% for the measurement of , and from 1% to 8% for the measurement of rel T , systematically increasing with the jet T . • Muon momentum scale and resolution uncertainties [71] include variations in the smearing of the inner detector and muon spectrometer tracks, variations of the muon momentum scale, and additional charge-dependent variations related to the correction for the sagitta bias. These variations are applied to the muon selection and both when adding the muon four-momenta to the jets using the ghost-association method described in Section 4 and when reconstructing the mass and four-momenta of the ± mesons. These are applied in a fully correlated way. The impact of this uncertainty is below 2% in all regions of the phase space, for both the and rel T measurements. • Muon identification uncertainties [72] cover variations of the muon efficiency corrections described in Section 4. The variations include 1 shifts of the statistical and systematic uncertainties in the efficiency corrections. As in the muon reconstruction case, the variations are applied, in a fully correlated way, both when adding muons to the jets and when reconstructing the hadron momenta. The impact of this uncertainty is small, with values below 3% in all regions of the phase space.
• The uncertainty in the dimuon trigger efficiency [70] is evaluated by using correction factors to take into account the uncertainty in the mismodelling of the trigger efficiency. These correction factors are applied as per-jet weights depending on the T and rapidity of the muons, as well as the angular distance Δ between the two muons in the ± decay. The difference between the weighted distribution and the default one defines the systematic uncertainty. The values of this uncertainty range from a few per mille to a maximum of 4% in the lower tails of the distribution, and from a few per mille to a maximum of 2% as a function of rel T . • The uncertainty in the kaon reconstruction efficiency is estimated by randomly rejecting 2% of the kaon tracks from the reconstruction, following a uniform random distribution [84]. The impact on the fragmentation properties is at the per mille level for both the and rel T measurements.

Jet-related uncertainties
• The jet energy scale (JES) and jet energy resolution (JER) uncertainties are estimated as described in Ref.
[76]. The JES is calibrated on the basis of the simulation and in situ corrections obtained from data. The JES uncertainties are estimated using a correlation scheme comprising a set of 29 independent components, which depend on the jet T and . The total JES uncertainty in the T of individual jets ranges from 2% at T = 50 GeV to approximately 1% at T = 100 GeV, with a mild dependence on . The JER uncertainty is estimated using a correlation scheme involving 13 independent variations. Each of these variations ranges from about 0.8% at T = 50 GeV to about 0.5% at T = 100 GeV. In this measurement, the JES and JER uncertainties are propagated by varying the energy and T of each jet by one standard deviation of each of the independent components. The impact of the JES uncertainty varies from a few per mille for low rel T values at high jet-T to 10% for medium values of at low jet-T , and typically decreases with increasing jet T . The JER has a similar behaviour, varying from approximately 1% for low rel T values at high jet T to 30% at medium values of at low jet-T .
• The impact of the jet angular resolution (JAR) uncertainty is estimated by smearing the angular coordinates ( , ) of the jets by 10% of the angular resolution of jets reconstructed from topological clusters [85], estimated in MC simulation. The and smearing is applied with the T component of the jets held constant. The JAR uncertainty is negligible for the measurement of . On the other hand, due to a linear dependence of rel T on the sine of the angle between the meson and the jet, this uncertainty can have an impact of up to 3% for the measurement of rel T . • The uncertainty from the efficiency corrections for the JVT algorithm described in Section 4, used to mitigate the impact of pile-up jets, is evaluated by using alternative correction factors, shifted by the corresponding uncertainties. This variation has two components, each of which shifts the correction factors in opposite directions. The resulting uncertainty is very small, with values below 0.1%.

Unfolding-related uncertainties
• The effects of the mismodelling of the data by the MC simulation on the unfolding is accounted for as an additional source of uncertainty. This is assessed by reweighting the particle-level distributions so that the detector-level fragmentation variables predicted by the MC samples match those in the data. The modified detector-level distributions are then unfolded using the method described in Section 6. The difference between the modified unfolded distribution and the reweighted particle-level distribution is taken as the uncertainty. The resulting uncertainty ranges from 0.2% to 10%, and typically increases with increasing values of and decreasing values of rel T . • The uncertainty due to using a particular MC model in the unfolding is estimated from the impact of two independent variations. The first accounts for the different amounts of gluon splitting in samples generated with different tunes. It is evaluated by comparing the results of the unfolding using the A14 samples with those using the different gluon splitting fractions in the Monash samples. These two particular tunes, which differ in the value of s for final-state radiation, 4 are chosen because fits to the gluon splitting fraction in data are compatible with the uncertainty band spanned by them. The second is related to the description of the measured observables made by different MC models, after disentangling the effect of different gluon splitting fractions. It is evaluated by repeating the unfolding using two alternative samples. While the nominal result is obtained using the P 8 A14-sample, the alternative samples are simulated using P 8 A14 and S , in which the fragmentation is simulated by the string model. To avoid double counting of the differences arising from the description of the gluon splitting, the S and P 8 A14 samples are reweighted so that the fraction of the selected jets arising from gluon splittings →¯is the same as in the P 8 A14-sample. The envelope of the differences between the results using P 8 A14-and both alternative samples covers the second source of systematic uncertainty. The total modelling uncertainty is defined by the sum in quadrature of these two variations. This uncertainty ranges from 4% to 30% and it is, together with the JES uncertainty, dominant for these measurements.

Pile-up-related uncertainties
• The uncertainty due to pile-up is derived from the ratios of the unfolded results for low ( < 32) and high ( ≥ 32) pile-up subsamples to the results for the nominal measurement. The uncertainty is defined as the envelope of the two ratios, and is below 10% for both and rel T . The total uncertainty is then determined as the sum in quadrature of the all the systematic uncertainties described above. The total uncertainty is larger for low values of the jet T , where it can reach values up to 30% in some bins. For the higher T bins, the uncertainty is systematically smaller, mainly because of the smaller jet energy scale components. Figure 5 shows the values of the systematic uncertainties as a function of and rel T , for the lowest and highest bins of the jet T .
The dependence of the particle-level results on the choice of PDF was tested by using two alternative PDF sets for the S predictions with the string fragmentation model. Predictions using NNPDF3.0 and MMHT2014 were compared with the nominal prediction, which uses CT14 (see Table 1). A maximum deviation of 2% from the nominal results is observed. Moreover, the CT14 uncertainty, estimated using the full set of eigenvectors, produces a total variation of the order of a few per mille. The effect of varying the PDF is thus much smaller than the experimental uncertainties, and is therefore neglected.  Figure 5: Systematic uncertainties as a function of for all jet-T bins (left) and as a function of rel T for all jet-T bins (right).

Results
The particle-level results are presented and compared with the predictions described in Section 3. Figures 6  to 8 show the distributions of the longitudinal and transverse profiles for each T bin.
The results show important differences between the low and high T bins. In particular, the lower tails of the distributions contain a larger fraction of the data at high T , which translates into more high-rel T events at high T . This is understood to be due to the fact that gluon splittings, →¯, occur with a larger probability at high values of the jet T . The resulting -quarks are reconstructed in the same jet, while the ± meson originates from the fragmentation of one of them, thus leading to smaller values of and higher values of rel T .
The results for the longitudinal profile show reasonable agreement with the H 7 prediction with the angle-ordered parton shower, while large discrepancies are observed with the dipole parton shower. In particular, the prediction largely overestimates the data in the low tails at low T , with smaller discrepancies for higher values of the jet T . This is understood to be a result of the larger fraction of jets arising from gluon splittings in the H 7 sample with the dipole parton shower.
The S predictions give a reasonable description of the distributions in the low and medium T bins, although both predictions differ from data for very high values of , particularly for the cluster model. At high T , the description worsens for the cluster model, which tends to underestimate the data at low and significantly overestimate the data at high , while the Lund string model confirms the deviation from data at high .
The P 8 samples tend to give descriptions that are qualitatively similar, and these are compatible with data within the systematic uncertainties across the different jet T bins. However, the Monash +Peterson model overestimates the data at intermediate and low T , while it tends to conform with the rest at high T . The Monash+Lund-Bowler and A14+Lund-Bowler models both overestimate the data at very high values of in all T ranges.
The discrepancies observed for the distributions have their counterparts in the rel T distributions as follows. Due to the larger gluon splitting fractions, the H 7 sample with the dipole parton shower significantly overestimates the data for rel T values between 1.5 and 4.0 GeV at low T , while the differences are smaller with increasing T . The H angle-ordered parton shower gives a better description of the rel T distributions, although non-negligible discrepancies are also observed.
The S predictions, particularly the one from the sample interfaced with the cluster hadronisation model, show large discrepancies for low values of rel T , which increase when moving towards higher bins of jet T . The description of the rel T tails is within uncertainties at low T , but for higher T , both S samples tend to underestimate the data for high rel T . In general, the description from the sample interfaced to the Lund string model is observed to be better than the one provided by the cluster model.  In order to explicitly test the scale dependence of the longitudinal and transverse profiles, the average values of the longitudinal and transverse profiles, and rel T , are studied as a function of the jet T . The results, together with the MC predictions, are shown in Figure 9. All P samples describe the scale dependence reasonably well for both and rel T , although the samples using the A14 and A14-tunes predict slightly larger values of (and slightly lower values of rel T ) than measured in data. The S sample making use of the cluster hadronisation model fails to describe the rel T data, disagreeing by 10% to 25%, but describes the distribution reasonably well, except at high T . The S sample interfaced to the Lund string hadronisation model describes the data well, while showing small discrepancies for rel T , although much smaller than when using the cluster model. The H 7 sample making use of the angle-ordered parton shower describes the scale dependence very well, while showing discrepancies of up to 15% for the rel T data. Finally, the H 7 sample implementing the dipole-based parton shower fails to describe both the and rel T profiles, showing discrepancies of up to 10% for the former, and up to 20% for the latter. As expected, the incidence of gluon splittings in the MC simulation grows as the jet T increases. The P samples show similar fractions of →¯jets, although the samples making use of the Monash tune present a slightly higher fraction due to using a larger value of s than the others and, hence, having a larger amount of gluon radiation. Both S samples present larger fractions than the P samples for low and medium T . However, for the highest T bin, the S sample using the cluster model presents a lower fraction than P . The angle-ordered shower implemented in H 7 shows an amount of →¯splitting very similar to that of S , while the H 7 sample making use of the dipole-based parton shower presents a much larger fraction of jets arising from gluon splittings. As pointed out previously, these differences help to explain the discrepancies between this sample and the data, since the lower tails of the distributions, as well as the higher tails for rel T , are strongly influenced by the fraction of gluon splittings.

Summary and conclusions
The fragmentation properties of -jets have been measured by using the decay of charged mesons into / charmonia and charged kaons, using 139 fb −1 of collisions at √ = 13 TeV recorded by the ATLAS detector at the LHC. To this end, the decay channel of the mesons is fully reconstructed using a three-track fit to a common vertex. The ± candidates are matched to jets, reconstructed using the anti-algorithm with = 0.4, and the longitudinal and transverse profiles of the mesons over the jet momentum are measured. The measurement spans three bins of the jet transverse momentum, and the yield of mesons is extracted for each bin of the fragmentation variable ( or rel T ) and the jet T . The results are then corrected for detector resolution effects using an iterative Bayesian algorithm, and the systematic uncertainties due to the muon, jet and event reconstruction properties are evaluated.
The results are compared with different MC predictions, including P 8, H 7 and S samples using different models for the parton shower and hadronisation. Generally, the best description of the longitudinal profile is provided by the P 8 and S samples making use of the string hadronisation model, which provide similar descriptions for all values of the jet T . For medium and high values of the jet T , the S sample with cluster hadronisation shows large deviations from data at very high values of . The H 7 sample with the angle-ordered parton shower also gives a good description of the longitudinal profile for all values of the jet T . In contrast, the H 7 samples using the dipole-based parton shower show large deviations from the data in the lower tails of the distribution. For the transverse profile, similar comments apply. While the P 8 and S samples using the string hadronisation model give fair and similar descriptions of the data for all values of the jet T , H 7 fails to describe the data. The angle-ordered parton shower, however, provides a better description than the dipole-based parton shower. At medium and high T , the S sample with the cluster hadronisation model gives a poor description for low values of rel T . This analysis provides key measurements which help to better understand the fragmentation functions of heavy quarks. As has been shown, significant differences among different MC models are observed, and also between the models and the data. Some of the discrepancies are understood to arise from poor modelling of the →¯splittings, to which the present analysis has substantial sensitivity. Including the present measurements in a future tune of the MC predictions may help to improve the description and reduce the theoretical uncertainties of processes where heavy-flavour quarks are present in the final state, such as top quark pair production or Higgs boson decays into heavy quark pairs. [86] ATLAS Collaboration, ATLAS Computing Acknowledgements, ATL-SOFT-PUB-2021-003, : https://cds.cern.ch/record/2776662.