Azimuthal correlations of prompt D mesons with charged particles in pp and p–Pb collisions at sNN=5.02TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\sqrt{s_\mathrm{NN}}} = 5.02\ \hbox {TeV}$$\end{document}

The measurement of the azimuthal-correlation function of prompt D mesons with charged particles in pp collisions at s=5.02TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s} =5.02\ \hbox {TeV}$$\end{document} and p–Pb collisions at sNN=5.02TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s_{\mathrm{NN}}} = 5.02\ \hbox {TeV}$$\end{document} with the ALICE detector at the LHC is reported. The D0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{D}^{0}$$\end{document}, D+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{D}^{+} $$\end{document}, and D∗+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{D}^{*+} $$\end{document} mesons, together with their charge conjugates, were reconstructed at midrapidity in the transverse momentum interval 30.3GeV/c\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_\mathrm{T} > 0.3\ \hbox {GeV}/c$$\end{document} and pseudorapidity |η|<0.8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\eta | < 0.8$$\end{document}. The properties of the correlation peaks appearing in the near- and away-side regions (for Δφ≈0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \varphi \approx 0$$\end{document} and Δφ≈π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \varphi \approx \pi $$\end{document}, respectively) were extracted via a fit to the azimuthal correlation functions. The shape of the correlation functions and the near- and away-side peak features are found to be consistent in pp and p–Pb collisions, showing no modifications due to nuclear effects within uncertainties. The results are compared with predictions from Monte Carlo simulations performed with the PYTHIA, POWHEG+PYTHIA, HERWIG, and EPOS 3 event generators.


Introduction
Two-particle angular correlations allow the mechanisms of particle production to be investigated and the event properties of ultra-relativistic hadronic collisions to be studied. In particular, the azimuthal and pseudorapidity distribution of "associated" charged particles with respect to a "trigger" D meson is sensitive to the charm-quark production, fragmentation, and hadronisation processes in proton-proton (pp) collisions and to their possible modifications in larger collision systems, like proton-nucleus (pA) or nucleus-nucleus (AA) [1]. The typical structure of the correlation function, featuring a "near-side" (NS) peak at ( ϕ, η) = (0, 0) (where ϕ is the difference between charged-particle and D-meson azimuthal angles ϕ ch − ϕ D , and η the difference between their pseudorapidities η ch − η D ) and an "away-side" See Appendix A for the list of collaboration members. e-mail: alice-publications@cern.ch (AS) peak at ϕ = π extending over a wide η range, as well as its sensitivity to the different charm-quark production mechanisms, are described in details in [2].
In this paper, results of azimuthal correlations of prompt D mesons with charged particles at midrapidity in pp and p-Pb collisions at √ s NN = 5.02 TeV are presented, where "prompt" refers to D mesons produced from charm-quark fragmentation, including the decay of excited charmed resonances and excluding D mesons produced from beautyhadron weak decays. The study of the near-side correlation peak is strongly connected to the characterisation of charm jets and of their internal structure, in terms of their particle multiplicity and angular profile. Probing the near-side peak features as a function of the charged-particle transverse momentum ( p T ), possibly up to values of a few GeV/c, gives not only access to the transverse-momentum distribution of the jet constituents, but can also provide insight into how the jet-momentum fraction not carried by the D meson is shared among the other particles produced by the parton fragmentation, as well as on the correlation between the p T of these particles and their radial displacement from the jet axis, which is closely related to the width of the near-side correlation peak. This study provides further and complementary information with respect to the analysis of charm jets reconstructed as a single object through a track-clustering algorithm and tagged by their charm content [3][4][5].
The azimuthal-correlation function of D mesons with charged particles is largely sensitive to the various stages of the D-meson and particle evolution, as hard-parton scattering, parton showering, fragmentation and hadronisation [6]. Its description by the available Monte Carlo event generators like PYTHIA [7,8], HERWIG [9][10][11], and EPOS 3 [12,13] or pQCD calculations like POWHEG [14,15] coupled to event generators handling the parton shower, depends on several features, including the order of the hard-scattering matrixelement calculations (leading order or next-to-leading order), the modelling of the parton shower, the algorithm used for the fragmentation and hadronisation, and the description of the underlying event. The azimuthal-correlation func-tion of D mesons with charged particles in pp collisions at √ s = 7 TeV measured by ALICE is described within uncertainties by simulations produced using PYTHIA6, PYTHIA8 and POWHEG+PYTHIA6 event generators [2]. However, more precise and differential measurements are needed to set constraints to models and be sensitive to the differences among their expectations.
The validation of Monte Carlo simulations for angular correlations of heavy-flavour particles in pp collisions is also useful for interpreting the results in nucleus-nucleus collisions, for which the measurements in pp collisions are used as reference. The temperature and energy density reached in nucleus-nucleus collisions at LHC energies are large enough to produce a quark-gluon plasma (QGP), a deconfined state of strongly-interacting matter [16,17]. The interaction of heavy quarks (charm and beauty) with the QGP should affect the angular-correlation function [1, 18,19]. First measurements performed at RHIC and the LHC showed modifications of the correlation function in nucleus-nucleus collisions when the trigger was a heavy-flavour particle, where a suppression of the away-side correlation peak and an enhancement of the near-side correlation peak for associated particles with p T < 2 GeV/c was observed [20,21]. A comparison of the results in nucleus-nucleus collisions to those in pp collisions, along with a successful description by models, would allow the modifications of the correlation function to be related to the in-medium heavy-quark dynamics [18,22,23].
In proton-nucleus collisions, several cold-nuclear-matter effects can influence the production, fragmentation and hadronisation of heavy-flavour quarks. They are induced by the presence of a nucleus in the initial state of the collision and, possibly, by the high density of particles in its final state. The most relevant effect is a modification of the parton distribution functions due to nuclear shadowing [24], which can consequently affect the heavy-flavour production cross section. Measurements of the nuclear modification factor of D mesons and of electrons from heavy-flavour hadron decays in p-Pb collisions at √ s NN = 5.02 TeV [25,26] point towards a small influence of cold-nuclear-matter effects on the heavy-flavour quark production at midrapidity. Nevertheless, nuclear effects could still affect the fragmentation and hadronisation of heavy quarks. These can be investigated by measuring potential modifications of the shape of the angular correlation between heavy-flavour particles [27] or, more indirectly, between heavy-flavour particles and charged particles.
Additionally, the search and characterisation of collectivelike effects in high-multiplicity proton-proton and protonnucleus collisions are a crucial topic, due to the observation of long-range, ridge-like structures in two-particle angularcorrelation functions at RHIC [28,29] and the LHC [30][31][32][33][34][35], resembling those observed in Pb-Pb collisions. The mechanism leading to these structures in small collision systems is not straightforward to identify. Possible explanations include final-state effects due to a hydrodynamic behaviour of the produced particles [36,37], colour-charge exchanges [38,39], initial-state effects, such as gluon saturation as described within the Color-Glass Condensate effective field theory [40,41], or gluon bremsstrahlung by a quark-antiquark string [42]. In addition, a positive ellipticflow coefficient was observed also for heavy-flavour particles, from the analysis of their azimuthal correlations with charged particles, by the ALICE [43][44][45], ATLAS [46][47][48], and CMS [49,50] Collaborations. This approach generally assumes that the jet-induced correlation peaks do not differ in low-and high-multiplicity collisions, i.e. nuclear effects have the same impact on the heavy-quark fragmentation and hadronisation at different event multiplicities. This assumption can be tested by looking for modifications of the azimuthal-correlation function.
The results presented in this paper significantly improve the precision and extend the kinematic reach, with respect to our previous measurements [2] in both pp (at a different energy) and minimum bias p-Pb collisions. Correlations with associated particles at higher p T probe the angular and p T distribution of the hardest jet fragments, which retain more closely the imprint of the hard-scattering topology. The properties of the away-side peak are also studied for the first time. The paper is structured as follows. In Sect. 2, the ALICE apparatus, its main detectors and the data samples used for the analysis are presented. In Sect. 3 the procedure adopted for building the azimuthal-correlation functions, correcting them for experimental effects, and extracting physical quantities is described. Section 4 describes the systematic uncertainties associated to the measurement. The results of the analysis are presented and discussed in Sect. 5. The paper is briefly summarised in Sect. 6.

Experimental apparatus and data sample
The ALICE apparatus consists of a central barrel, covering the pseudorapidity region |η| < 0.9, a muon spectrometer with −4 < η < −2.5 coverage, and forward-and backward-pseudorapidity detectors employed for triggering, background rejection, and event characterisation. A complete description of the detector and an overview of its performance are presented in [51,52]. The central-barrel detectors used in the analysis presented in this paper, employed for chargedparticle reconstruction and identification at midrapidity, are the Inner Tracking System (ITS), the Time Projection Chamber (TPC), and the Time-Of-Flight detector (TOF). They are embedded in a large solenoidal magnet that provides a magnetic field of 0.5 T, parallel to the beams. The ITS consists of six layers of silicon detectors, with the innermost two composed of Silicon Pixel Detectors (SPD). It is used to track charged particles and to reconstruct primary and secondary vertices. The TPC is the main tracking detector of the central barrel. In addition, it performs particle identification via the measurement of the particle specific energy loss (dE/dx) in the detector gas. Additional information for particle identification is provided by the TOF, via the measurement of the charged-particle flight time from the interaction point to the detector. The TOF information is also employed to evaluate the starting time of the event [53], together with the time information provided by the T0 detector, an array of Cherenkov counters located along the beam line, at +370 cm and −70 cm from the nominal interaction point.
The results reported in this paper were obtained on the data samples collected during the 2016 LHC p-Pb run at √ s NN = 5.02 TeV and the 2017 LHC pp run at √ s = 5.02 TeV, corresponding, after the event selection, to integrated luminosities of L int = (295 ± 11) μb −1 and L int = (19.3 ± 0.4) nb −1 , respectively. The events were selected using a minimum bias (MB) trigger provided by the V0 detector [54], a system of two arrays of 32 scintillators each, covering the full azimuthal angle in a pseudorapidity range of 2.8 < η < 5.1 (V0A) and −3.7 < η < −1.7 (V0C). The trigger condition required at least one hit in both the V0A and the V0C scintillator arrays. This trigger is fully efficient for recording collisions in which a D meson is produced at midrapidity [2]. The V0 time information and the correlation between number of hits and track segments in the SPD were used to reject background events from the interaction of one of the beams with the residual gas in the vacuum tube. Pile-up events, whose probability was below 1% (0.5%) in pp collisions (p-Pb collisions), were rejected with almost 100% efficiency by using an algorithm based on track segments, reconstructed with the SPD, to detect multiple primary vertices. The remaining undetected pile-up events are a negligible fraction of the analysed sample. In order to obtain a uniform acceptance of the detectors, only events with a reconstructed primary vertex within ±10 cm from the centre of the detector along the beam line were considered for both pp and p-Pb collisions. In p-Pb collisions, the √ s NN = 5.02 TeV energy was obtained by delivering proton and lead beams with energies of 4 TeV and 1.58 TeV per nucleon, respectively. Therefore, the proton-nucleus center-of-mass frame was shifted in rapidity by y NN = 0.465 in the proton direction with respect to the laboratory frame. The azimuthal correlations between D mesons and charged particles in p-Pb collisions were studied as a function of the collision centrality. The centrality estimator is based on the energy deposited in the zerodegree neutron calorimeter in the Pb-going direction (ZNA). The procedure used to define the centrality classes and to determine the average number of binary nucleon-nucleon collisions for each class is described in [55].
Some of the corrections for the azimuthal-correlation functions described in Sect. 3 were evaluated by exploiting Monte Carlo simulations, which included a detailed description of the apparatus geometry and of the detector response, using the GEANT3 package [56], as well as the luminous region distribution during the pp and p-Pb collision runs. For the evaluation of the charged-particle reconstruction efficiency, pp collisions were simulated with the PYTHIA8 event generator [8] with Monash-2013 tune [57], while p-Pb collisions were simulated using the HIJING 1.36 event generator [58] in order to describe the charged-particle multiplicity and detector occupancy observed in data [59]. For the corrections requiring the presence of a D meson in the event, enriched Monte Carlo samples were used, obtained by generating pp collisions containing a cc or bb pair in the rapidity range [−1.5, 1.5], employing PYTHIA 6.4.21 with Perugia-2011 tune. For p-Pb collisions, an underlying event generated with HIJING 1.36, was superimposed to each heavy-quark enhanced PYTHIA event.

Data analysis
The analysis largely follows the procedure described in detail in [2]. It consists of three main parts: (1) reconstruction and selection of D mesons and primary charged particles (see [60] for the definition of primary particle); (2) construction of the azimuthal-correlation function and corrections for detectorrelated effects, secondary particle contamination, and beauty feed-down contribution; (3) extraction of correlation properties via fits to the average D-meson azimuthal-correlation functions with charged particles.

Selection of D mesons and primary charged particles
The analysis procedure begins with the reconstruction of D mesons (D 0 , D * (2010) + , and D + and their charge conjugates), defined as "trigger" particles, and primary charged particles, considered as "associated" particles. The D mesons are reconstructed from the following hadronic decay channels: D 0 → K − π + (BR = 3.89 ± 0.04%), D + → K − π + π + (BR = 8.98 ± 0.28%), and D * + → D 0 π + → K − π + π + (BR = 2.63 ± 0.03%) [61] in the transverse-momentum interval 3 < p T < 24 GeV/c. The D-meson selection strategy, described in detail in [25,62], exploits the displaced topology of the decay and utilises the particle identification capabilities of the TPC and TOF to select on the D-meson decay particles. A dedicated optimisation on the selection variables was done, where the selections were tightened to increase the signalto-background ratio of the D-meson invariant mass peaks. A gain up to a factor 5 at low p D T was obtained with respect to the selection defined in [25,62], at the expenses of a reduction of the raw yield. This allowed reducing the impact of the D-meson combinatorial background, whose subtraction induces the largest source of statistical uncertainty on the The D-meson raw yields were extracted from fits applied to the invariant mass (M) distributions of D 0 and D + candidates, and to the distribution of the mass difference M = M(Kππ) − M(Kπ) for D * + candidates, for several subranges in the interval 3 < p T < 24 GeV/c. The fit function was composed of two terms, one for the signal and one for the background. The signal was described by a Gaussian, while the background was modelled by an exponential term for D 0 and D + mesons, and by a threshold function multiplied by an exponential for the D * + meson, as detailed in [2]. Examples of the invariant mass distributions in pp and in p-Pb collision systems are shown in Fig. 1 for D 0 , D + , and D * + mesons in different p T intervals.
Associated particles are defined as charged primary particles with p assoc T > 0.3 GeV/c and with pseudorapidity |η| < 0.8. As additional requirement, for this study only pions, kaons, protons, electrons and muons are considered as associated particles. The associated-particle sample does not include the decay products of the trigger D meson. Reconstructed charged-particle tracks with at least 70 space points out of 159 in the TPC, 2 out of 6 in the ITS, and a χ 2 /ndf of the momentum fit in the TPC smaller than 2 were considered. The contamination of non-primary particles was largely suppressed by requiring the distance of closest approach (DCA) of the track to the primary vertex to be less than 1 cm in the transverse (x y) plane and along the beam line (z-direction). This selection identifies primary particles with a purity varying from 95% to 99% (increasing with p assoc T ) and rejects a negligible amount of primary particles. In particular, less than 1% of the primary particles originating from decays of heavyflavour hadrons are discarded. For the D 0 mesons produced in D * + → D 0 π + decays, the lowp T pion accompanying the D 0 was removed from the sample of associated particles by rejecting tracks that, combined with the D 0 , yielded a M consistent within 3σ with the D * + mass peak. It was verified with Monte Carlo simulations that this selection rejects more than 99% of the pions from D * + decays in all D-meson p T intervals considered and has an efficiency larger than 99% for primary particles with p assoc T > 0.3 GeV/c. The selection criteria described above provided an average track reconstruction efficiency for charged particles with p assoc T > 0.3 GeV/c of about 83% (82%) in pp (p-Pb) collisions in the pseudorapidity interval |η| < 0.8, with an increasing trend as a function of p assoc T up to ≈ 1 GeV/c, followed by saturation at about 90%. As the track reconstruction efficiency has a sudden drop below ≈ 0.3 GeV/c, caused by the TPC requirements in the track selection, this transverse momentum value was chosen as the minimum p assoc T for the analysis.

Evaluation and correction of the azimuthal-correlation functions
Selected D-meson candidates with an invariant mass in the range |M −μ| < 2σ (peak region), where μ and σ denote the mean and width of the Gaussian term of the invariant mass fit function, were correlated to the primary charged particles selected in the same event. A two-dimensional angularcorrelation function C( ϕ, η) peak was evaluated by computing the difference of the azimuthal angle and the pseudorapidity of each pair. The azimuthal-correlation functions were studied in four D-meson p T intervals: 3 < p D T < 5 GeV/c, 5 < p D T < 8 GeV/c, 8 < p D T < 16 GeV/c, and 16 < p D T < 24 GeV/c and in the following p T ranges of the associated tracks: p assoc T > 0.3 GeV/c, 0.3 < p assoc T < 1 GeV/c, 1 < p assoc T < 2 GeV/c, and 2 < p assoc T < 3 GeV/c, significantly extending both transverse momentum coverages with respect to the previous measurements reported in [2].
The two-dimensional correlation functions are affected by the limited detector acceptance and reconstruction efficiency of the associated tracks (A assoc × assoc ), as well as the variation of those values for prompt D mesons (A trig × trig ) inside a given p D T interval. In order to correct for these effects, a weight equal to 1/(A assoc × assoc ) × 1/(A trig × trig ) was assigned to each correlation pair, as described in detail in [2]. A weight of 1/(A trig × trig ) was applied also to the entries in the D-meson invariant mass distributions, used for the evaluation of the amount of signal S peak and background B peak triggers in the peak region.
The two-dimensional correlation function C( ϕ, η) peak also includes correlation pairs obtained by considering Dmeson candidates from combinatorial background as trigger particles. This contribution was subtracted by evaluating the per-trigger correlation function obtained selecting D mesons with an invariant mass in the sidebands, 1/B sidebands × C( ϕ, η) sidebands , and multiplying it by B peak . The term B sidebands is the amount of background candidates in the sideband region, i.e. 4σ < |M − μ| < 8σ (5σ < M − μ < 10σ , for D * + mesons) of the invariant mass distributions weighted by the inverse of the prompt Dmeson reconstruction efficiency.
The event-mixing technique was used to correct the correlation functions C( ϕ, η) peak and C( ϕ, η) sidebands for the limited detector acceptance and its spatial inhomogeneities. The peak and sideband region event-mixing functions ME( ϕ, η) peak and ME( ϕ, η) sidebands were evaluated as explained in [2]. The inverse of these functions was used to weight the functions C( ϕ, η) peak and C( ϕ, η) sidebands , respectively.
The per-trigger angular-correlation function was obtained by subtracting the sideband-region correlation function from the peak-region one, as follows: (1) The division by S peak provides the normalisation to the number of D mesons. In our notation per-trigger quantities are specified by theC symbol. In Eq. 1, p prim ( ϕ) is a correction for the residual contamination of non-primary associated particles not rejected by the track selection (purity correction). This was evaluated with Monte Carlo simulations based on PYTHIA6 (Perugia-2011 tune) by quantifying the fraction of primary particles, among all the tracks satisfying the selection criteria. The correction was applied differentially in ϕ, since from Monte Carlo studies it was verified that this contamination shows a ϕ modulation, typically of about 1-2%. The largest value of the contamination was found in the nearside region, for the lowest p T range of the associated tracks, where p prim ( ϕ) approaches 95%. Statistical fluctuations prevented a ( ϕ, η)-doubledifferential study of the correlation peak properties. Therefore, the per-trigger azimuthal-correlation functionC inclusive ( ϕ) was obtained by integratingC inclusive ( ϕ, η) in the range | η| < 1.
A fraction of reconstructed D mesons originates from the decay of beauty hadrons (feed-down D mesons). It was verified with Monte Carlo simulations that azimuthal correlations of prompt and feed-down D mesons with charged particles show different functions. This is a result of the different fragmentation of beauty and charm quarks, as well as of the additional presence of beauty-hadron decay particles in the correlation function of feed-down D-meson triggers. The contribution of feed-down D-meson triggers to the measured angular-correlation function was subtracted using tem-plates of the azimuthal-correlation function of feed-down D mesons with charged particles, obtained with Monte Carlo simulations at generator level (i.e. without detector effects and particle selection), as detailed in [2].
Before performing this subtraction,C inclusive ( ϕ) has to be corrected for a bias which distorts the shape of the nearside region of the feed-down contribution, induced by the D-meson topological selection. For feed-down D-meson triggers, indeed, the selection criteria are more likely to be satisfied by decay topologies with small angular opening between the trigger D meson and the other products of the beautyhadron decay. This induces an enhancement of correlation pairs from feed-down D-meson triggers at ϕ ≈ 0 and a depletion at larger ϕ values. This bias was accounted for as a systematic uncertainty in [2]. In this paper, instead, a ϕ dependent correction factor (c FD−bias ( ϕ)) was determined by comparing Monte-Carlo templates of feed-down D mesons and associated particles at generator level and after performing the event reconstruction and particle selection as on data. This correction factor ranges between 0.6 at ϕ ≈ 0 and 1.3 at ϕ ≈ π/4, decreasing then to 1, and was applied to the feed-down contribution toC inclusive ( ϕ) as follows, to restore this contribution to an unbiased value: (2) is the value of the per-trigger correlation function of prompt (feed-down) Dmesons with associated particles, and the term A total NS ( ϕ) is the value of the per-trigger correlation function considering both prompt and feed-down components. The terms A prompt NS ( ϕ) and A feed−down NS ( ϕ) were evaluated from an analysis on reconstructed Monte Carlo events, where the reconstruction was performed as on data. The fraction of prompt D mesons in the raw yields, f prompt , was evaluated as detailed in [63]. It typically decreases from 95% to 90% with increasing p D T in the studied transverse-momentum intervals, independently of the collision centrality. The maximum effect of the correction when applied onC inclusive ( ϕ) is about 5%, at ϕ ≈ 0, for the lowest D-meson p T range and the highest p assoc T interval. The correction becomes negligible for p D T > 8 GeV/c. After performing this correction, the feed-down contamination was subtracted as described above.
As a result, the fully-corrected, per-trigger azimuthalcorrelation function of prompt D mesons with associated particles was obtained, denoted as 1/N D × dN assoc /d ϕ from Fig. 2 onwards.

Average and fit to the correlation functions
The correlation functions obtained from D 0 , D + , and D * + mesons were averaged using, as weights, the inverse of the quadratic sum of the statistical and D-meson uncorrelated systematic uncertainties, discussed in Sect. 4, since the three functions were found to be consistent within uncertainties. Since the correlation functions are symmetric around ϕ = 0 and ϕ = π , they were reflected in the range 0 < ϕ < π to reduce statistical fluctuations. In order to quantify the properties of the average D-meson azimuthalcorrelation function, it was fitted with the following function: (3) The fit function is composed of a constant term b describing the flat contribution below the correlation peaks, a generalised Gaussian term describing the NS peak, and a Gaussian reproducing the AS peak. In the generalised Gaussian, the term α is related to the variance of the function, hence to its width, while the term β drives the shape of the peak (the Gaussian function is obtained for β = 2). The function in Eq. 3 is a generalisation of that adopted in [2], where a Gaussian function was used for the near-side, corresponding to the case β = 2. The new parametrisation allowed to improve the χ 2 /ndf value in all the kinematic ranges studied, especially in the high p T ranges of both the D mesons and associated particles, where the standard Gaussian fit systematically underestimates the near-side peak yields (widths) up to 10% (20%) with respect to the generalised Gaussian. In both collision systems, the β parameter decreases monotonically from ≈ 2.2 at low p D T and p assoc T to ≈ 1 in the highest p D T and p assoc T intervals. By symmetry considerations, the means of the Gaussian functions were fixed to ϕ = 0 and ϕ = π . Figure 2 shows examples of fits to the azimuthalcorrelation functions of D mesons with associated particles, The integrals of the functions describing the near-and away-side peaks, Y NS and Y AS , correspond to the associatedparticle yields (i.e. the average number of associated particles contained in the peak), while the widths of the correlation peaks are described by the square root of the variance of their fitting terms, α √ (3/β)/ (1/β) and σ AS , for the nearand away-side, respectively. The baseline b represents the physical minimum of the ϕ function, and depends on the average charged-particle multiplicity. To reduce the effect of statistical fluctuations on the estimate of the associated yields, b was fixed to the weighted average of the points in the transverse region, defined as π/4 < | ϕ| < π/2, using the inverse of the point squared statistical uncertainties as weights.

Systematic uncertainties
The systematic uncertainty induced on the correlation function from the evaluation of S peak and B peak , obtained by fitting the D-meson invariant-mass distribution, was evaluated by varying the fit procedure. In particular, the fit was repeated modelling the background distribution with a linear function and a second-order polynomial function instead of an exponential function (for D 0 and D + mesons only), varying the fit range, fixing the mean of the Gaussian term describing the mass peak to the world-average D-meson mass [61], or fixing the Gaussian width to the value obtained from Monte Carlo studies. A systematic uncertainty ranging from 1 to 3% (1 to 2%), depending on the p D T , was estimated from the corresponding variation of the azimuthal-correlation function for pp (p-Pb) collisions. No dependence on ϕ was observed and the same uncertainty was estimated for all Dmeson species.
An uncertainty ranging from 1 to 3%, depending on p D T and on the D-meson species, was assigned in both pp and p-Pb collisions for the possible dependence of the shape of background correlation function on the invariant-mass value of the trigger D meson. This source of uncertainty was determined by evaluatingC( ϕ, η) sidebands , defining a different invariant-mass sideband range, and also considering, for D 0 and D + mesons, only the left or only the right sideband for the evaluation ofC( ϕ, η) sidebands . No significant dependence on ϕ was obtained for this uncertainty. A systematic effect originating from the correction of the D-meson reconstruction efficiency, due to possible differences of the topological variable distributions between Monte Carlo and data, was evaluated by repeating the analysis applying tighter and looser topological selections on the D-meson candidates, with a corresponding variation of the D-meson reconstruction efficiencies larger than ±25%. An uncertainty up to 2.5% (2%), increasing for smaller p D T values, was assigned in pp (p-Pb) collisions. No significant dependence on ϕ was observed. The same uncertainty was estimated for the three D-meson species.
The systematic uncertainty originating from the evaluation of the associated track reconstruction efficiency was estimated by varying the quality selection criteria applied to the reconstructed tracks, removing the request of at least two associated clusters in the ITS, or requiring a hit on at least one of the two SPD layers, or varying the request on the number of space points reconstructed in the TPC. An uncertainty Table 1 List of systematic uncertainties for the azimuthal-correlation functions in pp and in p-Pb collisions. If not specified, the uncertainty does not depend on ϕ System pp p-Pb D-meson species D 0 , D * + , D + D 0 , D * + , D + Signal, background normalisation ± 1-3% ± 1-2% Background ϕ function ± 1-3% ± 1-3% Associated-track reconstruction efficiency ± 2.5-4.5% ± 3% Primary-particle purity ± 1-2% ± 1.5-3% D-meson efficiency ± 1-2.5% ± 1-2% Feed-down subtraction up to 5%, ϕ-dependent up to 3%, ϕ-dependent Bias on topological selection up to 2%, ϕ-dependent up to 2%, ϕ-dependent up to 4.5% (3%), was assessed for pp (p-Pb) collisions. No significant trend in ϕ was observed.
The uncertainty on the evaluation of the residual contamination from secondary tracks was determined by repeating the analysis varying the selection on the DCA in the x y plane from 0.1 to 1 cm, and re-evaluating the purity of associated primary particles for each variation. This resulted in a 2% (3%) maximum systematic uncertainty on the azimuthalcorrelation functions in pp (p-Pb) collisions, decreasing with increasing p assoc T and with negligible ϕ dependence. The uncertainty on the subtraction of the beauty feeddown contribution was quantified by generating the templates of feed-down azimuthal-correlation functions with different event generators (PYTHIA6 with the Perugia-2010 tune, PYTHIA8 with the 4C tune) and by varying the value of f prompt within its uncertainty band, as described in details in [62]. The resulting uncertainty was found to be dependent on ϕ, with a maximum value of 5% (3%) in pp (p-Pb) collisions, and was applied point-by-point on the correlation functions.
As discussed in Sect. 3, Monte Carlo studies revealed the presence of a bias on the near-side region of the correlation function for feed-down D-mesons triggers, induced by the topological selections applied to the D mesons. The correction applied to remove this bias relies on a proper description of the azimuthal-correlation functions of prompt and feed-down D-meson triggers by the Monte Carlo simulations. A ϕ-dependent, symmetric systematic uncertainty of ±δC( ϕ)/ √ 12 was introduced to account for underor overestimation of the correction, where δC( ϕ) is the point-by-point shift of the correlation function induced by the correction. The largest value of the uncertainty was 2%, at ϕ ≈ 0, for both pp and p-Pb collisions.
In Tab. 1, the minimum and maximum values of the systematic uncertainties affecting the azimuthal-correlation functions, depending on the kinematic range, are listed for both collision systems. Only the uncertainties deriving from the feed-down subtraction and from the correction on the bias of feed-down D-meson correlations are ϕ dependent. All the other contributions define a ϕ-independent systematic uncertainty, which acts as a scale uncertainty for the correlation function. In both pp and p-Pb collisions the total scale uncertainty ranges from ±4% to ±5%.
The systematic uncertainties on the near-and away-side peak yields and widths, and on the baseline height, obtained from the fits to the azimuthal-correlation functions, were evaluated as follows. The main source of uncertainty arises from the definition of the ϕ transverse region used to determine the baseline height (term b of Eq. 3). The impact on the physical observables induced by the baseline value was estimated by considering different ϕ ranges for determining the baseline position and performing the fits again using Eq. 3. Moreover, the fits were repeated by moving the points of the correlation functions upwards and downwards using the corresponding value of the ϕ-dependent systematic uncertainty. The total systematic uncertainty was calculated by summing in quadrature the aforementioned contributions. For the associated yields and for the baseline, whose values depend on the normalisation of the correlation function, also the ϕ-independent systematic uncertainties affecting the correlation function (i.e. the first five contributions listed in Table 1), which act as a scale factor, were summed in quadrature.
In p-Pb collisions, the presence of long-range correlations among the particles produced in the collision can have an impact on the values of the quantities extracted from the fits, in particular for the analysis as a function of centrality. This effect was studied by fitting the functions with a v 2 -like modulation [43], in place of a flat baseline. The v 2 values adopted for D mesons, ranging up to 8% for the lowest p D T range in 0-20% central events, were estimated employing the available results for heavy-flavour particle v 2 in p-Pb collisions from CMS [49], ALICE [43], and ATLAS [46,47], while those for associated particles were estimated based on di-hadron correlation measurements by ALICE [30]. For the centrality-integrated analysis and for the case when p assoc T > 0.3 GeV/c, considering a v 2 -like modulation reduced the near-side peak yields by about 16% (5%) The uncertainties from the subtraction of the baseline are displayed as boxes at ϕ > π for 3 < p D T < 5 GeV/c (5 < p D T < 8 GeV/c) and the awayside peak yields by about 20% (3%) for 3 < p D T < 5 GeV/c (5 < p D T < 8 GeV/c). A smaller variation was observed for the peak widths and for the baseline value. For the analysis as a function of the event centrality, the largest effect was obtained for the 0-20% centrality class, where for p assoc T > 0.3 GeV/c a decrease of 27%, 17%, and 5% was found for 3 < p D T < 5 GeV/c, 5 < p D T < 8 GeV/c, and 8 < p D T < 16 GeV/c, respectively. Smaller variations were found for the near-side peak width and the baseline. This systematic uncertainty was summed in quadrature with the others to obtain the total uncertainty.

Comparison of results in pp and p-Pb collisions
The averaged azimuthal-correlation functions of the D 0 , D + , and D * + mesons with associated particles in pp and p-Pb collision systems are compared, after baseline subtraction, in Fig. 3, for four D-meson transverse momentum ranges, 3 < p D T < 5 GeV/c, 5 < p D T < 8 GeV/c, 8 < p D T < 16 GeV/c, and 16 < p D T < 24 GeV/c. The functions are presented for p assoc T > 0.3 GeV/c as well as for three sub-ranges, 0.3 < p assoc T < 1 GeV/c, 1 < p assoc T < 2 GeV/c, and < 3 GeV/c. The qualitative shape of the correlation function and the evolution of the near-and away-side peaks with trigger and associated particle p T are consistent within uncertainties in the two collision systems. In particular, an increase of the height of the near-side correlation peak is observed for increasing values of the D-meson p T . This reflects the production of a higher number of particles in the jet accompanying the fragmenting charm quark, when the energy of the latter increases. A similar, though milder, effect can be observed also for the away-side peak.
A more quantitative comparison of the near-and awayside peak features and p T evolution in the two collision systems can be obtained by fitting the azimuthal-correlation functions and evaluating the peak yields and widths, as it was explained in Sect. 3. Figure 4 compares these observables for the near-side correlation peaks in pp and p-Pb collisions, as a function of the D-meson p T , for p assoc T > 0.3 GeV/c and in three p assoc T sub-ranges. For both yields and widths, the values measured in the two collision systems are in agreement. The increase of associated particle production inside the near-side peak with p D T , qualitatively observed in Fig. 3, is present for all the associated particle p T intervals, and is similar in the two collision systems. A tendency for a narrowing of the near-side peak with increasing p D T is also observed in most of the p T ranges, though a flat behaviour cannot be excluded with the current uncertainties.
The away-side peak yields and widths measured in pp and p-Pb collisions are compared in Fig. 5 as a function of the D-meson p T , with the common associated-particle p T ranges analysed. For pp collisions, specific kinematic regions where the χ 2 /ndf of the fit was much larger than unity, or where the uncertainties on the peak observables were larger than 100%, were excluded from the results. As in the near-side analysis, the away-side yields show an increasing trend with p D T , and overall have similar values in the two collision systems. In the intermediate D-meson transverse momentum range, there is a hint for larger yields in p-Pb than in pp, but not a statistically significant one (about 2.2σ for the combined range 5 < p D T < 16 GeV/c for all p assoc T ranges). The away-side peak widths show consistent values in pp and p-Pb collisions in all kinematic ranges. No significant impact from cold-nuclear-matter effects on the fragmentation and hadronisation of charm quarks appears from the comparison of the results in the two collision systems, within the current precision of the measurements. This result complements the observation, emerged from the measurements reported in Refs. [25,26], that cold-nuclear-matter effects have a small impact on the production of charm quarks at midrapidity in p-Pb collisions.

Results in p-Pb collisions as a function of the event centrality
The correlation functions of D mesons with associated particles for p-Pb collisions in the 0-20%, 20-60%, and 60-100% centrality classes are compared in Fig. 6, for nine kinematic ranges with 3 < p D T < 16 GeV/c and p assoc T > 0.3 GeV/c. No results are shown for the 60-100% centrality class, for The baseline-subtracted correlation functions do not show significant differences among the three centrality intervals studied. Figure 7 shows the near-side yields and widths extracted by a fit to the correlation functions, for the three centrality intervals. A similar increase of the near-side peak yields, as a function of the D-meson p T , is observed for the three centrality ranges, with the absolute values of the yields also being generally in agreement. The only exception is for the 3 < p D T < 5 GeV/c, p assoc T > 0.3 GeV/c interval, where the yield for the 60-100% centrality class is lower than for the 0-20% and 20-60% centrality classes, with a statistical significance of 1.4σ and 2.1σ , respectively. This effect could be due to statistical fluctuations of the correlation function data points (see Fig. 6). The near-side peaks also have consistent widths among the three centrality ranges, for all the kinematic ranges studied. No centrality dependence on the correlation peaks, which could have possibly been induced by nuclear-matter effects, is observed within the experimental uncertainties. The limited precision of the results does not provide a further validation of the subtraction technique of the jet-induced correlation peaks, commonly used in analyses searching for positive elliptic flow via two-particle correlations.

Comparison of ALICE results to predictions from Monte Carlo simulations
The azimuthal-correlation functions of D mesons with associated particles, as well as the near-and away-side peak yields and widths measured by ALICE in pp collisions, were compared to expectations from several Monte Carlo event generators.
The PYTHIA event generator [7,8] allows for the generation of high-energy collisions of leptons and/or hadrons. It employs 2 → 2 QCD matrix elements evaluated perturbatively with leading-order precision, with the next-toleading order contributions taken into account during the parton showering stage. The parton showering follows a leading-logarithmic p T ordering, with soft-gluon emission divergences excluded by an additional veto, and the hadronisation is handled with the Lund string-fragmentation model. Two different versions of PYTHIA, with two different parameter tunes, were used in this paper. The PYTHIA 6.4.25 version [7] was employed, incorporating the Perugia 2011 tune [64], which was the first tune considering the data from pp collisions at √ s = 0.9 TeV and √ s = 7 TeV at the LHC. With respect to its predecessor, PYTHIA8 [8] has an improved handling of the multiple-parton interactions and > 1 GeV/c (from top to bottom). Statistical and ϕ-dependent systematic uncertainties are shown as vertical error bars and boxes, respectively, while ϕ-independent uncertainties are written as text. The uncertainties from the subtraction of the baseline are displayed as boxes at ϕ > π the colour reconnection processes. In this paper, it was used with the tune 4C [65].
POWHEG [14,15] is a pQCD generator implementing hard-scattering matrix elements with NLO accuracy, which can be coupled to Monte Carlo generators, like PYTHIA [7,8] or HERWIG [9,10], for the parton showering and hadronisation of the produced partons. In this paper, Monte Carlo simulations were done using the POWHEG-BOX [66] framework coupled to PYTHIA 6.4.25 with the Perugia-2011 tune [64]. A charm-quark mass of m c = 1.5 GeV/c 2 was considered, and the renormalisation and factorisation scales were set to the transverse mass of charm quark, i.e. μ R = μ F = p 2 T + m 2 c . It was verified that simulation  [67]: the variation of the charm-quark mass does not alter the correlation function, while the variation of the renormalisation and factorisation scales produces differences of ±10% (±5%) for the near-side peak yields (widths) and negligible deviations for the away-side peak yields and widths. This can be expected, since the per-trigger correlation function of D mesons with associated particles is scarcely sensitive to the absolute rate of production of D mesons, directly influenced by the aforementioned parameters. An additional set of predictions from POWHEG+PYTHIA was also evaluated (POWHEG LO+PYTHIA6 in the following), by stopping the computation of the hard-scattering matrix elements at leading-order accuracy, before passing the generated partons to PYTHIA for the showering and hadronisation.
The HERWIG 7 [10,11] event generator allows one to perform Monte Carlo simulations at NLO accuracy for most of the Standard Model processes, including heavy-quark production. The parton showering is performed with an angular ordering of the fragments, which correctly takes the coherence effects for soft-gluon emissions into account. In addition, the hadronisation is handled via the cluster hadronisation model, differently from the Lund string fragmentation model employed by PYTHIA.
EPOS 3 [12,13] is a Monte Carlo generator which considers flux tube initial conditions for the collision, generated in the Gribov-Regge multiple-scattering framework, and applies a 3+1D viscous hydrodynamical evolution on the dense core of the collision. Individual scatterings, referred to as Pomerons, are identified with parton ladders, each composed of a pQCD hard process, plus initial-and final-state radiations. The hadronisation is then performed with a string fragmentation procedure. Non-linear effects are considered by means of a saturation scale. An evaluation within the EPOS 3 model shows that the energy density reached in pp collisions at the LHC energies is already sufficient for applying such a hydrodynamic evolution [68]. In the following, due to the limited precision of the available predictions, the comparison between EPOS 3 expectations and data results will < 1 GeV/c, 1 < p assoc T < 2 GeV/c, and 2 < p assoc T < 3 GeV/c (from top to bottom). Statistical and ϕ-dependent systematic uncertainties are shown as vertical error bars and boxes, respectively, while the ϕ-independent uncertainties are written as text. The uncertainties from the subtraction of the baseline are displayed as boxes at ϕ > π be restricted to the kinematic interval 3 < p D T < 16 GeV/c, and will not include the away-side peak observables.
In Fig. 8 the azimuthal-correlation functions of D mesons with associated particles obtained from the aforementioned event generators are compared to the measurements from this analysis, for all the p D T and p assoc T ranges studied, in pp collisions at √ s = 5.02 TeV, after the baseline subtraction. For the models, for which the statistical fluctuations are generally negligible, the baseline was estimated as the minimum of the azimuthal-correlation function, and a systematic uncertainty on the fit parameters was assessed by repeating the fits after fixing the baseline as the weighted average of the two lowest points of the correlation function. Most of the models provide a fair description of the two correlation peaks in the various kinematic ranges studied, though some tensions are visible from this qualitative comparison. In par-ticular, HERWIG underestimates the near-side peak height for p assoc T > 1 GeV/c, especially for low D-meson transverse momentum, while EPOS 3 tends to overestimate the height of the near-side peak and gives a flatter away-side peak. In addition, some systematic hierarchies among the models appear throughout the whole p T ranges analysed, with POWHEG+PYTHIA6 providing the highest near-side peak, and in most of the cases the smallest away-side peak. The overestimation of the near-side peak yield by EPOS 3 is a relevant feature also for the understanding of the dependence of heavy-flavour production on the charged-particle multiplicity measured in the same rapidity window of the heavy-flavour signals [69]. Disentangling the role of jet-biases from effects related to genuine global event properties is fundamental for properly interpreting the measured trends, especially their p T dependence [70].
A more detailed investigation can be performed by quantifying the peak yields and widths extracted from the fit to the correlation functions. In Fig. 9, the comparison of near-side peak yields and widths from data and simulation is shown, as a function of the D-meson p T , for p assoc T > 0.3 GeV/c and for the three p assoc T sub-ranges analysed, in pp collisions at √ s = 5.02 TeV. In the top row (third row down), the absolute value of the yields (widths) are displayed, while the second (fourth) row down reports the ratios of the yields (widths) to those obtained with POWHEG+PYTHIA6, which reduces the visual impact of the statistical fluctuations of the data points. As already visible from Fig. 8, EPOS 3 predicts the largest values of the near-side yields, followed by POWHEG+PYTHIA6, while POWHEG LO+PYTHIA6 shows about 10% lower yields with respect to the version with NLO accuracy. The latter difference could be explained by a different relative contribution of the NLO production mechanisms, in particular the gluon splitting, present already at the level of the hard scattering for POWHEG+PYTHIA6. PYTHIA8 provides near-side yield values comparable to those of POWHEG LO+PYTHIA6, while PYTHIA6 yields are slightly lower. HERWIG expectations for near-side yields are the lowest, except for the 0.3 < p assoc T < 1 GeV/c range, where they are comparable to PYTHIA8 expectations. POWHEG+PYTHIA6 and POWHEG LO+PYTHIA6 provide the best description of the near-side yields, with data points lying between the two predictions. PYTHIA8 also gives a good description of data, especially for p assoc T > 1 GeV/c, while PYTHIA6 predictions are generally lower than data, though in agreement within the uncertainties. The HERWIG expectations for near-side yield describe the data well for the 0.3 < p assoc T < 1 GeV/c range, while they severely underpredict the measurements for the other p assoc T ranges, especially for the lower intervals of the D-meson p T . In particular, for the integrated range p assoc T > 0.3 GeV/c, a discrepancy of 3.3σ (2.9σ ) is found for 3 < p D T < 5 GeV/c (5 < p D T < 8 GeV/c), increasing to 3.4σ (3.6σ ) for the highest associated particle transverse-momentum range 2 < p assoc T < 3 GeV/c. The EPOS 3 model largely overestimates the near-side associated yields, especially at low D-meson p T . For p assoc T > 0.3 GeV/c, the discrepancy between data and predictions ranges between 4.0σ and 5.2σ . Except for EPOS 3, a similar hierarchy among the models also characterises the near-side widths. POWHEG+PYTHIA6 give the broadest peaks, followed by POWHEG LO+PYTHIA6, with increasing difference between the two model expectations with increasing p assoc T . PYTHIA8 gives similar widths as POWHEG LO+PYTHIA6, while PYTHIA6 widths are generally lower. The lowest predictions are provided by EPOS 3. HERWIG predictions are consistent with PYTHIA6 for p assoc T < 1 GeV/c, and are generally lower for p assoc T > 1 GeV/c. POWHEG+PYTHIA6 provides systematically larger widths than data, though still being compatible pointby-point. EPOS 3 predictions tend to underestimate the nearside widths, despite being consistent with data point-bypoint. All the other models provide values of the near-side width closer to data.
The same comparison of model expectations to data is shown for the away-side peak yields and widths in Fig. 10. POWHEG+PYTHIA6 gives the smallest away-side yields, with about 5%-10% smaller values than POWHEG LO+PYTHIA6 predictions. As for the near-side peak yields, this difference could be ascribed to a different contribution from back-to-back topologies of charm-quark pair production. PYTHIA8 and PYTHIA6 yields are rather similar, and systematically larger than POWHEG LO+PYTHIA6 expectations. HERWIG predicts similar yields as POWHEG LO+PYTHIA6 for the integrated p assoc T range (with larger values for 0.3 < p assoc T < 1 GeV/c and smaller values for p assoc T > 1 GeV/c). The best description of the awayside yields is provided, as in the near-side peak case, by POWHEG+PYTHIA6 and POWHEG LO+PYTHIA6 over the whole kinematic range, as well as by HERWIG for p assoc T > 1 GeV/c. As observed for the near-side peak case, the PYTHIA8 and PYTHIA6 expectations tend to overpredict away-side yields in the majority of the transversemomentum intervals studied. For p assoc T > 0.3 GeV/c, about a 2σ difference with respect to the data is present, over the whole p D T interval studied. The largest values of the awayside peak width, in particular for large values of p assoc T , are given by the PYTHIA6 event generator, which tends to systematically overpredict the data points. The predictions from the other models, all in agreement with data, are very similar, with POWHEG+PYTHIA6 being in general the lowest of them. However, the precision of measurements for this observable prevents from discerning the model that best describes the data. Figure 11 shows the baseline values of the measured azimuthal-correlation functions and compares them to predictions from the event generators. The measured baseline values decrease with increasing p T of the associated particle, which is expected as the transverse-momentum distribution of associated particles in pp collisions at √ s = 5.02 TeV peaks at few hundred MeV/c [71]. From the data it cannot be concluded whether the baseline is flat or slightly increasing as a function of D-meson p T . A mildly increasing trend with p D T is predicted by the event generators. However, POWHEG+PYTHIA6 and POWHEG LO+PYTHIA6 predict a larger increase than HERWIG, EPOS 3 and PYTHIA. The same baseline values are obtained by POWHEG+PYTHIA6 and POWHEG LO+PYTHIA6 for all the kinematic ranges. This is not trivial, due to the differ- ent treatment of next-to-leading order contributions to charm production, which can populate the transverse region of the correlation function and, hence, affect the baseline value. The best description of the results, for low values of the associated particle p T , is provided by PYTHIA6, PYTHIA8, and EPOS, while HERWIG overestimates the values by about 15% over the whole p D T range and POWHEG+PYTHIA6 underpredicts them by 20% at low p D T . For p assoc T > 1 GeV/c HER-WIG gives the closest description of the baseline. PYTHIA6, PYTHIA8, EPOS 3, and POWHEG+PYTHIA6 tend to underpredict data values, with the first three models catching well the p D T dependence, while POWHEG+PYTHIA6 also predicting a different behaviour against p D T . The baseline-subtracted azimuthal-correlation functions of D mesons with associated particles measured in p-Pb collisions were compared to simulations from PYTHIA6, PYTHIA8, and POWHEG+PYTHIA6 event generators. The only modifications of the configuration of these models with respect to that used in pp collisions consisted of a rapidity shift of the centre-of-mass system and, for POWHEG+PYTHIA6, a nuclear correction for the parton distribution functions [72], which induced negligible effects on the model expectations. The comparison between these models and the results from p-Pb collision yielded very similar conclusions as those discussed for pp collisions, not only in terms of an overall agreement between data and models, but also for the differences previously mentioned for specific observables and kinematic ranges. This was expected, given the overall agreement of measurements in the two collision systems as discussed in Sect. 5.1, where additional coldnuclear-matter effects, not included in the models, could also be present.

Summary
Measurements of azimuthal-correlation functions of D 0 , D * + , and D + mesons with charged particles in pp collisions at √ s = 5.02 TeV and p-Pb collisions at √ s NN = 5.02 TeV were reported. The results obtained have statistical and systematic uncertainties smaller by a factor of about 2-3 than those reported in our previous paper [2] for the common p T ranges, and extend the p T coverage both for the trigger and associated particles, allowing for a more differential study of the correlation function and charm-jet properties.
After subtracting the baseline, the correlation functions, along with the values and transverse momentum evolution of the near-and away-side peaks, are found to be consistent in pp and p-Pb collisions, in all the kinematic ranges addressed. This suggests that the fragmentation and hadronisation of charm quarks is not strongly influenced by coldnuclear-matter effects, complementing what was observed in previous measurements [25,26,73] that suggested a small impact from cold-nuclear-matter effects on D-meson production in the p T region covered by our measurement.
The analysis in p-Pb collisions was also performed in the 0-20%, 20-60%, and 60-100% centrality intervals, in order to study the possible modifications of the charm fragmentation as a function of the event centrality. The same correlation pattern, along with similar values and p T evolution of the near-side peak observables were found for the three centrality ranges, within large experimental uncertainties.
The baseline-subtracted correlation functions and the near-and away-side peak yields and widths measured by ALICE in pp collisions were compared to predictions by several event generators, with different modelling of charm production, parton showering, and hadronisation. In general, the models describe well the main features of the correlation functions. POWHEG+PYTHIA6 provides the best description to experimental data of near-and awayside yields. PYTHIA8 tends to overestimate the away-side peak yields, while providing a good description of the nearside peak yields and of the widths of both peaks. Overall, PYTHIA6 is more distant from data than PYTHIA8, although in general it is consistent with the measurements. HERWIG largely underestimates the near-side peak yields for p assoc T > 1 GeV/c, while it describes reasonably well the data at lower p assoc T , and provides a good description of the away-side peak features. Finally, EPOS 3 provides a higher near-side peak and qualitatively underestimates the awayside peak. Similar conclusions were obtained when comparing results in p-Pb collisions to predictions from the models available in this collision system.
The agreement between data and model expectations suggests that charm-quark production, fragmentation and hadronisation processes, as implemented in POWHEG+PYTHIA6 and PYTHIA8, provide an overall satisfactory description of the measured correlation functions. Therefore, in view of future analyses in Pb-Pb collisions, these models constitute a valid theoretical baseline for interpreting possible modifications of charm-jet properties and thus of the near-side correlation peak induced by the interactions of charm quarks with the quark-gluon plasma constituents. The same argument holds for the modifications of the whole correlation function, whose characterisation can provide a deeper understanding of heavy-quark dynamics inside the QGP medium. In addition, with the increased precision compared to previous measurement, and being at the same centre-of-mass energy of the available Pb-Pb collision samples, these results constitute the reference for measurements in that collision system.

Acknowledgements
The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector:

Data Availability Statement
This manuscript has associated data in a data repository. [Authors' comment: Manuscript has associated data in a HEPData repository, Link for the record is: https://www.hepdata.net/ record/95121.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .