1 Introduction

The study of the angular correlation of D mesons with charged particles, i.e. the distribution of the differences in azimuthal angles, \(\Delta {\varphi }=\varphi _\mathrm{ch}-\varphi _\mathrm{D}\), and pseudorapidities, \(\Delta \eta =\eta _\mathrm{ch}-\eta _\mathrm{D}\), allows for the characterisation of charm production and fragmentation processes in proton–proton (pp) collisions and of their possible modifications due to nuclear effects in proton–Pb and Pb–Pb collisions [1, 2]. For leading-order (LO) Quantum-ChromoDynamic (QCD) processes, charm quark–antiquark pairs are produced back-to-back in azimuth: the angular correlation of D mesons with charged particles features a “near-side” peak at \((\Delta {\varphi },\Delta {\eta })=(0,0)\) and an “away-side” peak at \(\Delta {\varphi }=\pi \). The former originates from the jet containing the “trigger” D meson, the latter is induced by the recoil jet, which can also include the decay products of the other charmed hadron produced in the collision. The away-side peak extends over a wide range in \(\Delta \eta \). The two peaks lie on top of an approximately flat distribution arising from the correlation of D mesons with charged particles from the underlying event. Next-to-leading order (NLO) production processes can give rise to significantly different correlation patterns [3, 4]. For example, the radiation of a hard gluon from a charm quark smears the back-to-back topology of LO production and broadens both the near- and the away-side peak. In addition, quark–antiquark charm pairs originating from the splitting of a gluon can be rather collimated and, especially at high transverse momentum (\({p}_{\mathrm{T}}\)), can generate sprays of hadrons contributing to a unique and broader “near-side” peak. In such cases, the away-side peak stems from the particles coming from the fragmentation of the recoil parton (typically a gluon or a light quark), which is not aligned with the trigger D meson. Finally, in hard-scattering topologies classified as “flavour excitation” (see e.g. [4]) a charm quark (antiquark) from an initial splitting \(g\rightarrow c\bar{c}\) undergoes a hard interaction. The hadrons originating from the antiquark (quark) can be significantly separated in rapidity with respect to the trigger D meson and contribute with a rather flat term to the \(\Delta {\varphi }\) distribution.

Since the first measurement performed by STAR in Au–Au collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=200~\mathrm {GeV}\) [5], two-particle azimuthal correlations have been exploited at both RHIC and the LHC [6,7,8] to investigate the possible modifications of jet and dijet properties that can be caused by the interaction of high-energy partons with the constituents of the Quark Gluon Plasma (QGP) formed in ultra-relativistic heavy-ion collisions. The most evident effect is the suppression of the away-side correlation peak, commonly attributed to in-medium partonic energy loss. The results allow one to constrain the dependence of the energy loss on the distance covered by partons in the QGP as well as the initial gluon density [9, 10]. The correlation pattern of hadron–hadron pairs primarily arises from the back-to-back production of gluons or light-quarks produced in hard-scattering processes, and their subsequent fragmentation. PHENIX measured the azimuthal correlation of electrons from heavy-flavour hadron decays with charged particles in Au–Au collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=200~\mathrm {GeV}\) [11]. The near- and away-side peaks are suppressed by factors compatible, within uncertainties, to those observed for hadron–hadron correlations, if a similar \({p}_{\mathrm{T}}\) is considered for the trigger hadron and the electron parent hadron. The proper interpretation of nucleus–nucleus results and the connection of the modifications of the correlation peak properties to the parton dynamics in the QGP requires the comparison of data with model predictions. It is crucial that the models reproduce the correlation pattern measured in pp collisions, where nuclear effects are absent, as well as the production spectra in both pp and nucleus–nucleus collisions. Therefore, the measurement of azimuthal correlations of D mesons with charged particles in pp and p–Pb collisions serves not only as a reference for future measurements in Pb–Pb collisions but it also allows for validation of Monte-Carlo generator expectations, which is fundamental for the understanding of the results in all collision systems.

Perturbative QCD calculations relying on the collinear-factorisation approach, like FONLL [12] and GM-VFNS [13], or based on the \(k_\mathrm{T}\)-factorisation approach [14] describe reasonably well the \({p}_{\mathrm{T}}\)-differential production cross sections of D mesons from charm-quark fragmentation measured at central rapidity (\(|y|<0.5\)) in pp collisions at \(\sqrt{s}=7\) and \(2.76~\mathrm {TeV}\) using the ALICE detector [15, 16]. These calculations represent the state of the art for the computation of (\({p}_{\mathrm{T}},y\))-differential cross sections of charm quarks and charmed hadrons. However, the kinematic relationship between D mesons and particles from charm fragmentation and the underlying event is accessible only with event generators coupled with parton-shower Monte-Carlo programs like those provided by PYTHIA [17] and HERWIG [18]. The order of hard-scattering matrix elements used, the specific implementation of parton shower and hadronisation, as well as the modeling of the underlying event have an influence on the angular correlations of D mesons with charged particles produced in the event. For heavy quarks with mass M and energy \(E_\mathrm{Q}\), the suppression of gluon radiation off the quark inside the forward cone with opening angle \(\Theta =M/E_\mathrm{Q}\) (the so-called “dead-cone” effect) reduces the phase space for primary gluon radiation [19]. This implies a harder fragmentation of the quarks into the heavy hadrons and leads to essential differences in the profiles of gluon-, light-quark- and heavy-quark-initiated jets resulting in shape differences of \({p}_{\mathrm{T}}\)-spectra and multiplicity distributions of primary hadrons in the jets [20, 21].

Correlations between D mesons were measured at the LHC in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) with the LHCb experiment [22], providing information on charm production mechanisms and on the properties of events containing heavy flavours. ATLAS measured the production of \({\mathrm{D}^{*+}}\) mesons in jets in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) for jets with \(25<{p}_{\mathrm{T}}<70~\mathrm {GeV}/c\) and \({\mathrm{D}^{*+}}\) carrying a jet momentum fraction (z) in the range \(0.3<z<1\). The results indicate that the production of charm-quark jets or charm-quark fragmentation into \({\mathrm{D}^{*+}}\) mesons is not properly modeled in state-of-the-art Monte-Carlo generators [23]. Azimuthal correlations of electrons from heavy-flavour hadron decays with charged particles were also exploited to study the relative beauty contribution to the population of electrons from heavy-flavour hadron decays in pp collisions at RHIC and at the LHC [24, 25].

The angular distribution of particles produced in an event is sensitive to collective effects that correlate particle production over wide phase-space regions. This is particularly relevant in Pb–Pb collisions with non-zero collision impact parameter, where the azimuthal asymmetry of the overlapping region of the colliding nuclei gives rise to anisotropic pressure gradients inducing an anisotropy in the azimuthal distribution of particle momenta [26, 27]. The main component of the Fourier decomposition used to describe the resulting \(\Delta {\varphi }\) distribution of two particle correlations is the 2nd order term, proportional to \(\cos (2\Delta {\varphi })\), called elliptic flow or \({v}_{2}\). Given that correlations induced by the collective motion of the system extend over large pseudorapidity ranges, the elliptic-flow term manifests itself with the presence of two long-range ridge-like structures in the near and away sides of two-particle angular correlations. Unexpectedly, similar long-range correlation structures were observed in high-multiplicity pp and p–Pb collisions at the LHC [28,29,30,31,32,33]. Also in central d–Au collisions at RHIC [34, 35] similar results were obtained, although contributions from jet-like correlations due to biases on the event selection could not be excluded [36]. The origin of such \({v}_{2}\)-like structures is still debated. Positive \({v}_{2}\) values in high-multiplicity pp collisions and p–Pb (d–Au) collisions at LHC (RHIC) are expected in models that include final-state effects [37,38,39,40,41], as well as initial-state effects related to the Color Glass Condensate [42] or to gluon bremsstrahlung by a quark–antiquark string [43]. A modification of the azimuthal correlations of D mesons with charged particles in p–Pb with respect to pp collisions could be a signal of the presence of long-range \({v}_{2}\)-like correlations for particles originating from hard-scattering processes. This would yield complementary information to that obtained from correlations of light-flavour particles, which at low \({p}_{\mathrm{T}}\) are primarily produced in soft processes. The D-meson \({p}_{\mathrm{T}}\)-differential production cross section in p–Pb collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\) was measured with ALICE in the interval of rapidity in the nucleon-nucleon centre-of-mass system \(-0.96<y_\mathrm{cms}<0.04\) [44]. The data are compatible, within uncertainties, with a Glauber-model-based geometrical scaling of a pp collision reference obtained from the cross sections measured at \({\sqrt{s}}=7~\mathrm {TeV}\) and \({\sqrt{s}}=2.76~\mathrm {TeV}\). This suggests that nuclear effects are rather small for D mesons in the range \(1<{p}_{\mathrm{T}}<24~\mathrm {GeV}/c\). However, they could still affect angular correlations as observed at RHIC for azimuthally-correlated pairs of electrons and muons from decays of heavy-flavour hadrons in d–Au collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=200~\mathrm {GeV}\) [45]. A modification of the azimuthal correlation of heavy-flavour particles in p–Pb collisions could occur at the LHC due to gluon saturation in the heavy nucleus [46]. Moreover, transport models based on the Langevin equation [2, 47] describe, within uncertainties, the nuclear modification factor of D mesons measured in p–Pb collisions at the LHC and that of electrons from heavy-flavour hadron decays measured in d–Au collisions at RHIC [48]. These models assume the formation of a small-size QGP in p–Pb and d–Au collisions and include the possibility of heavy-flavour hadron formation via coalescence of heavy quarks with thermalised light quarks from the medium. These transport calculations predict a positive D-meson \({v}_{2}\) in central p–Pb collisions. As an example, in the case of the POWLANG model [2] the maximum expectation for the 20% most central p–Pb collisions is \({v}_{2}\sim 5\%\) at \({p}_{\mathrm{T}}=4~\mathrm {GeV}/c\). A finite \({v}_{2}\) of muons from heavy-flavour hadron decays in high-multiplicity p–Pb collisions was also suggested in [31] as one of the possibilities for reconciling the measured values of \({v}_{2}\) of inclusive muons with the expectations based on the multi-phase transport model AMPT [49].

In this paper we report the first measurements of azimuthal correlations of D mesons with charged primary particles in pp and p–Pb collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\), respectively. Unless differently specified we always refer to “prompt” D mesons from charm-quark fragmentation. In what follows, primary particles are defined as particles originated at the collision point, including those deriving from strong and electromagnetic decays of unstable particles, and those from decays of hadrons with charm or beauty. The paper is organised as follows. In Sect. 2 the data samples used and the details of the ALICE experimental apparatus relevant for this analysis are described. The analysis strategy, the D-meson signal extraction, the associated-track selection criteria, and the corrections applied to measure the correlations between D mesons and charged primary particles are reported in Sect. 3. In the same section, the fit procedure adopted to quantify the correlation peak properties is described. Section 4 reports the systematic uncertainties affecting the measurement. The results are discussed in Sect. 5. The paper is then summarised in Sect. 6.

2 Experimental apparatus and data samples

2.1 The ALICE detector and event selection

The ALICE apparatus [50, 51] consists of a central barrel embedded in a 0.5 T solenoidal magnetic field, a forward muon spectrometer, and a set of detectors located in the forward- and backward-rapidity regions dedicated to trigger and event characterisation. The analysis reported in this paper is performed using the central barrel detectors. Charged particle tracks are reconstructed using the Inner Tracking System (ITS), consisting of six layers of silicon detectors, and the Time Projection Chamber (TPC). Particle identification (PID) is based on the specific energy loss dE/dx in the TPC gas and on the time of flight from the interaction vertex to the Time-Of-Flight (TOF) detector. The ITS, TPC and TOF detectors provide full azimuthal coverage in the pseudorapidity interval \(|\eta |<0.9\).

The pp data sample consists of about \(3\times 10^8\) minimum-bias events, corresponding to an integrated luminosity of \(L_\mathrm{int} = 5~{\mathrm{nb}^{-1}}\). These collisions are triggered by the presence of at least one hit in one of the V0 scintillator arrays, covering the ranges \(-3.7<\eta <-1.7\) and \(2.8<\eta <5.1\), or in the Silicon Pixel Detector (SPD), constituting the two innermost layers of the ITS, with an acceptance of \(|\eta |<2\) (inner layer) and \(|\eta |<1.4\) (outer layer). The p–Pb data sample consists of about \(10^8\) minimum-bias events, corresponding to an integrated luminosity of about \(L_\mathrm{int}=50~{\upmu \mathrm{b}^{-1}}\). In this case the minimum-bias trigger requires signals in both the V0 detectors.

Only events with a reconstructed primary vertex within \(\pm 10\) cm from the centre of the detector along the beam line are considered for both pp and p–Pb collisions. This choice maximises the detector coverage of the selected events, considering the longitudinal size of the interaction region, and the detector pseudorapidity acceptances (for more details see [51]). For p–Pb collisions, the center-of-mass reference frame of the nucleon-nucleon collision is shifted in rapidity by \(\Delta y_\mathrm{{NN}} = 0.465\) in the proton direction with respect to the laboratory frame, due to the different per-nucleon energies of the proton and the lead beams.

Beam-gas events are removed by offline selections based on the timing information provided by the V0 and the Zero Degree Calorimeters (two sets of neutron and proton calorimeters located around 110 m from the interaction point along the beam direction), and the correlation between the number of hits and track segments in the SPD detector.

The minimum-bias trigger efficiency is 100\(\%\) for events with D mesons with \({p}_{\mathrm{T}}> 1~\mathrm {GeV}/c\) for both pp and p–Pb data sets. For the analyzed data samples, the probability of pile-up from collisions in the same bunch crossing is below 4\(\%\) per triggered pp event and below the percent level per triggered p–Pb event. Events in which more than one primary interaction vertex is reconstructed with the SPD detector are rejected, which effectively removes the impact of in-bunch pile-up events on the analysis. The contribution of particles from pile-up of pp collisions in different bunch crossings is also negligible due to the selections applied to the tracks used in this analysis and the large interval between subsequent bunch crossings in the data samples used.

2.2 Monte-Carlo simulations

Monte-Carlo simulations including a complete description of the ALICE detector are used to calculate the corrections for the azimuthal-correlation distributions evaluated from data. The distribution of the collision vertex along the beam line, the conditions of all the ALICE detectors, and their evolution with time during the pp and p–Pb collision runs are taken into account in the simulations. Proton-proton collisions are simulated with the \(\mathrm{PYTHIA}\) 6.4.21 event generator [17] with the Perugia-0 tune (tune number 320) [52] while p–Pb collisions are simulated using the HIJING v1.36 event generator [53]. For the calculation of D-meson reconstruction efficiencies PYTHIA simulations of pp collisions are used, requiring that in each event a \(\mathrm{c}{\bar{\mathrm{c}}}\) or \(\mathrm{b}{\bar{\mathrm{b}}}\) pair is present. In the simulation used for the analysis of p–Pb data, a p–Pb collision simulated with HIJING is added on top of the PYTHIA event. The generated particles are transported through the ALICE apparatus using the GEANT3 package [54].

The measured angular-correlation distributions are compared to simulation results obtained with the event generators \(\mathrm{PYTHIA}\) 6.4.25 [17] (tunes number 320, 327, and 350, corresponding to the reference versions of the Perugia-0, Perugia-2010, and Perugia-2011 sets [52], respectively), \(\mathrm{PYTHIA}~8.1\) (tune 4C) [55], \(\mathrm{POWHEG}\) [56, 57] coupled to PYTHIA (Perugia-2011 tune), and EPOS 3.117 [58,59,60] (referred to as EPOS 3 hereafter). PYTHIA simulations utilise LO-pQCD matrix elements for \(2\rightarrow 2\) processes, along with a leading-logarithmic \({p}_{\mathrm{T}}\)-ordered parton shower, the Lund string model for hadronisation, and an underlying-event simulation including Multiple-Parton Interactions (MPI). With respect to older tunes, the Perugia tunes use different initial-state radiation and final-state radiation models. One of the main differences is that the parton shower algorithm is based on a \({p}_{\mathrm{T}}\)-ordered evolution rather than a virtuality-ordered one. Significant differences in the treatment of colour reconnection, MPI, and the underlying event were also introduced. Perugia 0 is the first of the series. The Perugia-2010 tunes differ from those of Perugia-0 in the amount of final-state radiation and by a modification of the high-z fragmentation (inducing a slight hardening of the spectra). They are expected to better reproduce observables related to the jet shape. The first LHC data, mainly from multiplicity and underlying-event related measurements, were considered for the Perugia-2011 tunes. PYTHIA 8.1 also includes several improvements in the treatment of MPI and colour reconnection [55]. In the simulations done with \(\sqrt{s}=5.02~\mathrm {TeV}\), the centre-of-mass frame is boosted in rapidity by \(\Delta y_\mathrm{{NN}} = 0.465\) in order to reproduce the rapidity shift of the reference frame of the nucleon-nucleon collision in the p–Pb collision system.

\(\mathrm{POWHEG}\) is a NLO-pQCD generator [56, 57] that, coupled to parton shower programs (e.g. from PYTHIA or HERWIG [18]), can provide exclusive final-state particles, maintaining the next-to-leading order accuracy for inclusive observables. The charm-production cross sections obtained with POWHEG+PYTHIA are consistent with FONLL [12] and GM-VFNS [13] calculations within the respective uncertainties, and are in agreement with measured D-meson production cross sections within the model and experimental uncertainties [61, 62]. The POWHEG+PYTHIA simulations presented in this paper are obtained with the POWHEG BOX framework [63, 64] and the tune Perugia 2011 of PYTHIA 6.4.25. For the comparison with the measured p–Pb collision data, parton distribution functions (PDFs) corrected for nuclear effects (CT10nlo [65] with EPS09 [66]) are used. In addition, a boost in rapidity by \(\Delta y_\mathrm{{NN}} = 0.465\) is applied to the partons generated with POWHEG before the PYTHIA parton shower process.

EPOS 3 [58,59,60] is a Monte-Carlo event generator based on a 3+1D viscous hydrodynamical evolution starting from flux tube initial conditions, which are generated in the Gribov-Regge multiple-scattering framework. Individual scatterings are referred to as Pomerons, and are identified with parton ladders. Each parton ladder is composed of a pQCD hard process with initial and final state radiation. Non-linear effects are considered by means of a saturation scale. The hadronisation is performed with a string fragmentation procedure. Based on these initial conditions, the hydrodynamical evolution can be applied on the dense core of the collision. An evaluation within the EPOS 3 model shows that the energy density reached in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) is high enough to apply such hydrodynamic evolution [60].

3 Data analysis

The analysis procedure consists of three main parts, which are described in the following subsections: D-meson reconstruction and selection of primary particles to be used in the correlation analysis (Sect. 3.1), construction of azimuthal-correlation distribution and corrections, including the subtraction of combinatorial background and beauty feed-down contributions (Sect. 3.2), extraction of correlation properties via fits to the azimuthal distributions (Sect. 3.3).

3.1 D-meson and associated-particle reconstruction

The correlation analysis is performed by associating D mesons (\(\mathrm{D}^{0}\), \({\mathrm{D}^{+}}\), \({\mathrm{D}^{*+}}\) mesons and their antiparticles), defined as “trigger” particles, with charged primary particles in the same event, and excluding those coming from the decay of the trigger D mesons themselves. The \(\mathrm{D}^{0}\), \({\mathrm{D}^{+}}\), \({\mathrm{D}^{*+}}\) mesons and their charge conjugates are reconstructed via their hadronic decay channels \(\mathrm{D}^{0} \rightarrow \mathrm{K}^{-}\pi ^{+}\), with Branching Ratio (BR) of (3.88±0.05)\(\%\), \(\mathrm{D}^{+} \rightarrow \mathrm{K}^{-}\pi ^{+}\pi ^{+}\), BR of (9.13±0.19)\(\%\), and \(\mathrm{D}^{*+} \rightarrow \mathrm{D}^{0} \pi ^{+}\), BR of (67.7±0.5)\(\%\) [67]. The D-meson signal extraction is based on the reconstruction of decay vertices displaced from the primary vertex by a few hundred microns and on the identification of the decay-particle species. The same selection procedures used for the measurements of D-meson production in pp and p–Pb collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\), respectively, are adopted [15, 44]. For both the pp and p–Pb data sets, \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) candidates are formed by combining two or three tracks, respectively, with each track satisfying \(|\eta |<0.8\) and \({p}_{\mathrm{T}}>0.3~\mathrm {GeV}/c\). Additionally, \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) daughter tracks are required to have at least 70 out of a maximum of 159 possible associated space points in the TPC, a \(\chi ^{2}/\)NDF of the momentum fit in the TPC smaller than 2, and at least 2 out of 6 associated hits in the ITS. \({\mathrm{D}^{*+}}\) candidates are formed combining \(\mathrm{D}^{0}\) candidates with tracks with one point in the SPD, \(| {\eta } | <0.8\) and \({p}_{\mathrm{T}}>0.1~\mathrm {GeV}/c\). The main variables used to reject the combinatorial background are the separation between primary and secondary vertices, the distance of closest approach (DCA) of the decay tracks to the primary vertex, and the angle between the reconstructed D-meson momentum and the flight line defined by the primary and secondary vertices. A tighter selection is applied for p–Pb collisions with respect to pp collisions to reduce the larger combinatorial background. Charged kaons and pions are identified using the TPC and TOF detectors. A \(\pm 3\sigma \) cut around the expected value for pions and kaons is applied on both TPC and TOF signals. The D mesons are selected in a fiducial rapidity range varying from \(|y_\mathrm{lab}|<0.5\) at low \({p}_{\mathrm{T}}\) to \(|y_\mathrm{lab}|<0.8\) for D mesons with \({p}_{\mathrm{T}}> 5~\mathrm {GeV}/c\) in order to avoid cases in which the decay tracks are close to the edge of the detector, where the acceptance decreases steeply. The \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) raw yields are extracted using fits to the distributions of invariant mass \(M(\mathrm{K^-}\pi ^\mathrm{+})\) and \(M(\mathrm{K^-}\pi ^\mathrm{+}\pi ^\mathrm{+})\), respectively, with a function composed of a Gaussian term for the signal and an exponential term that models the combinatorial background. In the case of the \({\mathrm{D}^{*+}}\), the raw yield is obtained by fitting the invariant-mass difference \(\Delta M=M(\mathrm{K^-}\pi ^\mathrm{+}\pi ^\mathrm{+}) - M(\mathrm{K^-}\pi ^\mathrm{+})\), using a Gaussian function for the signal and a threshold function multiplied by an exponential (\(a\sqrt{\Delta M - M_\pi } \cdot e^{b(\Delta M - M_\pi )}\)) to describe the background. Relatively wide D-meson \({p}_{\mathrm{T}}\) intervals (\(3<{p}_{\mathrm{T}}<5~\mathrm {GeV}/c\), \(5<{p}_{\mathrm{T}}<8~\mathrm {GeV}/c\), \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\) for pp collisions and \(5<{p}_{\mathrm{T}}<8~\mathrm {GeV}/c\), \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\) for p–Pb collisions) are chosen to reduce the statistical fluctuations in the azimuthal-correlation distributions. Figure 1 shows the \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) invariant mass, and \({\mathrm{D}^{*+}}\) invariant-mass difference distributions in the \(3<{p}_{\mathrm{T}}<5~\mathrm {GeV}/c\) interval for pp collisions and in the \(5<{p}_{\mathrm{T}}<8~\mathrm {GeV}/c\), \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\) intervals for p–Pb collisions. The fits used to evaluate the raw yields are also shown.

The statistical uncertainty of the D-meson raw yields in the \({p}_{\mathrm{T}}\) intervals analyzed varies from about 5 to 8% (3 to 5%) in pp (p–Pb) collisions for the \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) mesons and from about 5 to 6% (5 to 10%) for the \({\mathrm{D}^{*+}}\) mesons, depending on \({p}_{\mathrm{T}}\). For both collision systems, the signal over background ratio of the signal peaks is between 0.2 and 1 for the \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) mesons, and up to 2.6 for the \({\mathrm{D}^{*+}}\) meson. In the interval \(3<{p}_{\mathrm{T}}<5~\mathrm {GeV}/c\) the D-meson yield can be extracted from the invariant mass distribution with statistical uncertainty smaller than 3% in both pp and p–Pb collisions. However, in the latter case, the near- and away-side peaks of the azimuthal-correlation distribution, that have a small amplitude at low D-meson \({p}_{\mathrm{T}}\), cannot be disentangled from the statistical fluctuations of the baseline, which is related to the multiplicity of the event and thus higher in p–Pb than in pp collisions. Therefore, for this \({p}_{\mathrm{T}}\) interval, the results are shown only for pp collisions.

Fig. 1
figure 1

Distributions of \(\mathrm{D}^{0}\) (left column) and \({\mathrm{D}^{+}}\) (middle column) candidate invariant mass and of the \({\mathrm{D}^{*+}}\) candidate invariant-mass difference (right column). The distributions are shown for pp collisions in the \(3<{p}_{\mathrm{T}}<5~\mathrm {GeV}/c\) range (top row) and for p–Pb collisions in the \(5<{p}_{\mathrm{T}}<8~\mathrm {GeV}/c\) (middle row) and \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\) (bottom row) ranges. The fits to the invariant mass distributions and the Gaussian mean and sigma values are also shown

Associated particles are defined as all charged primary particles with \({p}_{\mathrm{T}}^\mathrm{assoc}>0.3\) \(\mathrm {GeV}/c\) and with pseudorapidity \(|\eta |<0.8\), except for the decay products of the trigger D meson. Particles coming from other weak decays or originating from interactions with the detector material are defined as secondary particles and are discarded. Reconstructed tracks with at least 70 points in the TPC and 3 in the ITS, and a \(\chi ^{2}/\)NDF of the momentum fit in the TPC smaller than 2 are associated to D-meson candidates. Using Monte Carlo simulations (see Sect. 2.2), these selection criteria yield an average track reconstruction efficiency for charged primary particles of about 85% in the pseudorapidity range \(|\eta |<0.8\) and in the interval \(0.3<{p}_{\mathrm{T}}<24~\mathrm {GeV}/c\), with variations contained within \(\approx \)5% for \({p}_{\mathrm{T}}<1.5~\mathrm {GeV}/c\). Negligible variations are observed at higher \({p}_{\mathrm{T}}\). The contamination of secondary particles is removed by requiring the DCA of the associated tracks to the primary vertex to be less than \(2.5~\mathrm{mm}\) in the transverse (xy) plane and less than \(1~\mathrm{cm}\) along the beam line (z direction). This selection identifies primary particles with a purity (\(p_\mathrm{prim}\)) of approximately \(96\%\) and an efficiency higher than \(99\%\), also for particles originating from decays of charm or beauty hadrons, which can be displaced by several hundred micrometers from the primary vertex. The purity is independent of \({p}_{\mathrm{T}}\) in the measured \({p}_{\mathrm{T}}\) range. For the \(\mathrm{D}^{0}\)-meson case, the low-\({p}_{\mathrm{T}}\) pion produced from the \({\mathrm{D}^{*+}}\rightarrow \mathrm{D}^{0}\pi ^+\) decay is removed from the sample of associated particles by rejecting tracks that yield a \(\Delta M\) compatible within \(3\sigma \) with the value expected for \({\mathrm{D}^{*+}}\) mesons. It was verified with Monte Carlo simulations that this selection rejects more than 99% of the pions from \({\mathrm{D}^{*+}}\) decays in all D-meson \({p}_{\mathrm{T}}\) intervals considered and has an efficiency larger than \(99\%\) for primary particles with \({p}_{\mathrm{T}}>0.3~\mathrm {GeV}/c\).

3.2 Azimuthal-correlation distributions and corrections

D-meson candidates with invariant mass (M) in the range \(|M-\mu |<2\sigma \) (peak region), where \(\mu \) and \(\sigma \) denote the mean and width of the Gaussian term of the invariant-mass fit function, are correlated to tracks selected with the criteria described above, and the difference in the azimuthal angle (\(\Delta \varphi \)) and in pseudorapidity (\(\Delta \eta \)) of each pair is computed. In order to correct for the acceptance and reconstruction efficiency (\(\mathrm {Acc} \times \epsilon \)) of the associated tracks and for the variation of (\(\mathrm {Acc} \times \epsilon \)) of prompt D mesons inside a given \({p}_{\mathrm{T}}\) interval, a weight equal to the inverse of the product of both (\(\mathrm {Acc} \times \epsilon \)) is assigned to each pair. The dependence of the associated-track efficiency on transverse momentum, pseudorapidity, and position of the primary vertex along the beam axis is taken into account. The dependence of the track reconstruction efficiency on the event multiplicity is negligible and therefore neglected. The reconstruction efficiency of prompt D mesons is calculated as a function of \({p}_{\mathrm{T}}\) and event multiplicity. It is on the order of few percent in the lowest D-meson \({p}_{\mathrm{T}}\) interval, about 20\(\%\) at high \({p}_{\mathrm{T}}\) [15, 44], and varies within each \({p}_{\mathrm{T}}\) interval by up to a factor 2–3 (1.5–2) at low (high) \({p}_{\mathrm{T}}\), depending on the D-meson species and collision system. The D-meson (\(\mathrm{{Acc}} \times \epsilon \)) factor also accounts for the \({p}_{\mathrm{T}}\)-dependent fiducial rapidity range of the selected D mesons (Sect. 3.1) in order to normalise the results to one unit of rapidity.

The obtained distribution, \(C(\Delta \varphi ,\Delta \eta )_\mathrm{peak}\), also includes the angular correlation of combinatorial D-meson candidates in the peak range, which is a source of background and needs to be subtracted. This contribution is estimated via the per-trigger correlation distribution of background candidates in the sideband invariant-mass range, \(1/B_\mathrm{sidebands} \times C(\Delta \varphi ,\Delta \eta )_\mathrm{sidebands}\), where \(B_\mathrm{sidebands}\) is the amount of background in the sideband region \(4\sigma<|M-\mu |<8\sigma \) (right side only, \(4\sigma<M-\mu <15\sigma \), in the case of \({\mathrm{D}^{*+}}\) mesons). The term \(C(\Delta \varphi ,\Delta \eta )_\mathrm{sidebands}\) represents the correlation distribution obtained as described above, but selecting trigger D-meson candidates with invariant mass in the sidebands. The background contribution is then subtracted from \(C(\Delta \varphi ,\Delta \eta )_\mathrm{peak}\) after being normalised to the amount of combinatorial background in the peak region, \(B_\mathrm{peak}\). The latter is obtained from the counts in the invariant-mass distribution in the peak region, after subtracting the signal, \(S_\mathrm{peak}\), estimated from the invariant-mass fit. Note that \(S_\mathrm{peak}\), \(B_\mathrm{peak}\) and \(B_\mathrm{sidebands}\) are calculated from the invariant-mass distributions weighted by the inverse of the prompt D-meson reconstruction efficiency.

The correlation distributions \(C(\Delta \varphi ,\Delta \eta )_\mathrm{peak}\) and \(C(\Delta \varphi ,\Delta \eta )_\mathrm{sidebands}\) are corrected for the limited detector acceptance and spatial inhomogeneities using the event mixing technique. In this approach, D-meson candidates found in a given event are correlated with charged tracks from other events with similar multiplicity and primary-vertex position along the beam axis. The distribution obtained from the mixed events, \(\mathrm{ME}(\Delta \varphi ,\Delta \eta )\), shows a typical triangular shape as a function of \(\Delta \eta \), due to the limited \(\eta \) coverage of the detector, and is approximately flat as a function of \(\Delta \varphi \). The event-mixing distribution is rescaled by its average value in the range (\(-0.2<\Delta \varphi <0.2\),\(-0.2<\Delta \eta <0.2\)) and its inverse is used as a map to weight the distributions \(C(\Delta \varphi ,\Delta \eta )_\mathrm{peak}\) and \(C(\Delta \varphi ,\Delta \eta )_\mathrm{sidebands}\). A correction for the purity of the primary-particle sample (\(p_\mathrm{prim}\), see Sect. 3.1) is applied and the per-trigger normalisation is obtained dividing by \(S_\mathrm{peak}\). The above procedure is summarised in Eq. 1, where the notation \(\tilde{C}\) refers to angular-correlation distributions normalised by the number of trigger particles:

$$\begin{aligned}&\tilde{C}_\mathrm{inclusive}(\Delta \varphi ,\Delta \eta )\nonumber \\ \quad&= \frac{p_\mathrm{prim}}{S_\mathrm{peak}}\left( \left. \frac{C(\Delta \varphi ,\Delta \eta )}{\mathrm{ME}{(\Delta \varphi ,\Delta \eta )}}\right| _\mathrm{peak} -\frac{B_\mathrm{peak}}{B_\mathrm{sidebands}}\left. \frac{C(\Delta \varphi ,\Delta \eta )}{\mathrm{ME}{(\Delta \varphi ,\Delta \eta )}}\right| _\mathrm{sidebands}\right) ,\nonumber \\&\mathrm{ME}(\Delta \varphi ,\Delta \eta )=\left( \frac{C(\Delta \varphi ,\Delta \eta )}{\langle C(\Delta \varphi ,\Delta \eta )\rangle _{|\Delta \varphi |,|\Delta \eta |<0.2}}\right) _\mathrm{Mixed~Events}. \end{aligned}$$
(1)

Finally, the per-trigger azimuthal distribution \(\tilde{C}_\mathrm{inclusive}(\Delta \varphi )\) is obtained by integrating \(\tilde{C}_\mathrm{inclusive}(\Delta \varphi ,\Delta \eta )\) in the range \(|{\Delta \eta }| <1\).

It was verified using Monte-Carlo simulations based on PYTHIA (Perugia-2011 tune) that the per-trigger azimuthal correlation of D mesons and secondary particles not rejected by the track selection has a \(\Delta \varphi \)-dependent modulation with a maximum variation of \(7\%\) with respect to the azimuthal correlation of D mesons and primary particles. This \(\Delta \varphi \)-dependent contamination has a negligible impact on the final results, considering the \(4\%\) level of contamination of secondary particles in the sample of associated tracks, hence, it was neglected.

A fraction of the reconstructed D mesons consists of secondary D mesons coming from B-meson decays. The topological cuts, applied to reject combinatorial background, preferentially select displaced vertices, yielding a larger (by about a factor 2 for \(\mathrm{D}^{0}\) mesons in the measured \({p}_{\mathrm{T}}\) range) efficiency for secondary D mesons than for prompt D mesons. Therefore, the fraction \(f_\mathrm{prompt}\) of reconstructed prompt D mesons does not coincide with the natural fraction and depends on the analysis details. The different fragmentation, as well as the contribution of B-meson decay particles and a possible different contribution of gluon splitting to charm- and beauty-quark production, imply a different angular-correlation distribution of prompt and secondary D mesons with charged particles, as it was verified with the Monte-Carlo simulations described in Sect. 2.2. The contribution of feed-down D mesons to the measured angular correlation is subtracted as follows:

$$\begin{aligned} \tilde{C}_\mathrm{prompt}(\Delta \varphi )= & {} \frac{1}{f_\mathrm{prompt}}(\tilde{C}_\mathrm{inclusive}(\Delta \varphi )\nonumber \\&-(1-f_\mathrm{prompt})\tilde{C}_\mathrm{{\text {feed-down}}}^\mathrm{MC~templ}(\Delta \varphi ) ). \end{aligned}$$
(2)

In Eq. 2, \(\tilde{C}_\mathrm{prompt}(\Delta \varphi )\) is the per-trigger azimuthal-correlation distribution after the subtraction of the feed-down contribution, \(f_\mathrm{prompt}\) is the fraction of prompt D mesons and \(\tilde{C}_\mathrm{feed-down}^\mathrm{MC~templ}(\Delta \varphi )\) is a template for the azimuthal-correlation distribution of the feed-down component. Using the same method described in [15], \(f_\mathrm{prompt}\) was evaluated on the basis of FONLL calculations of charm and beauty \({p}_{\mathrm{T}}\)-differential production cross sections [12] and of the reconstruction efficiencies of prompt and secondary D mesons, calculated using Monte-Carlo simulations. The value of \(f_\mathrm{prompt}\), which depends on the D-meson species and varies as a function of \({p}_{\mathrm{T}}\), is estimated to be larger than \(75\%\). The azimuthal correlation of feed-down D mesons, \(\tilde{C}_\mathrm{feed-down}^\mathrm{MC~templ}\), was obtained from PYTHIA (tune Perugia 2011 [52]) simulations of pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and \({\sqrt{s}}=5.02~\mathrm {TeV}\) for the analysis of pp and p–Pb data, respectively. In order to avoid biases related to the different event multiplicity in real and simulated events, the correlation distribution was shifted to have its minimum coinciding with the baseline of the data azimuthal-correlation distribution before feed-down subtraction. A difference smaller than \(8\%\) was observed in the simulation between the baseline values of the azimuthal-correlation distributions for prompt and feed-down D mesons. Considering the typical values of \(f_\mathrm{prompt}\), this difference results in a shift of the baseline of \(\tilde{C}_\mathrm{prompt}(\Delta \varphi )\) smaller than \(2\%\), negligible with respect to the other uncertainties affecting the measurement.

3.3 Characterization of azimuthal-correlation distributions

In order to quantify the properties of the measured azimuthal correlations, the following fit function is used:

$$\begin{aligned} f(\Delta \varphi )=b+\frac{A_{\mathrm{NS}}}{\sqrt{2\pi }\sigma _\mathrm{fit,NS}}e^{-\frac{(\Delta \varphi )^{2}}{2\sigma ^{2}_{\mathrm{fit,NS}}}}+\frac{A_{\mathrm{AS}}}{\sqrt{2\pi }\sigma _{\mathrm{fit,AS}}}e^{-\frac{(\Delta \varphi -\pi )^{2}}{2\sigma ^{2}_{\mathrm{fit,AS}}}}. \end{aligned}$$
(3)

It is composed of two Gaussian terms describing the near- and away-side peaks and a constant term describing the baseline. A periodicity condition is also imposed to the function, requiring \(f(0) = f(2\pi )\).

The integrals of the Gaussian terms, \(A_{\mathrm{NS}}\) and \(A_{\mathrm{AS}}\), correspond to the associated-particle yields for the near (NS)- and away (AS)-side peaks, respectively, while \(\sigma _{\mathrm{fit,NS}}\) and \(\sigma _{\mathrm{fit,AS}}\) quantify the widths of the correlation peaks. By symmetry considerations, the mean of the Gaussian functions are fixed to \(\Delta \varphi = 0\) and \(\Delta \varphi = \pi \). The baseline b represents the physical minimum of the \(\Delta {\varphi }\) distribution. To limit the effect of statistical fluctuations on the estimate of the associated yields, b is fixed to the weighted average of the points in the transverse region, defined as \(\pi /4< |\Delta \varphi | < \pi /2\), using the inverse of the square of the point statistical uncertainty as weights. Given the symmetry of the correlation distributions around \(\Delta \varphi =0\) and \(\Delta \varphi =\pi \), the azimuthal distributions are reported in the range \(0<\Delta \varphi <\pi \) to reduce statistical fluctuations. The effect of a \({v}_{2}\)-like modulation in the \(\Delta \varphi \) distribution, which could be present in p–Pb collisions, was estimated and assessed in Sect. 5.

In the case of the simulations, for which statistical fluctuations are negligible, the baseline is estimated as the minimum of the azimuthal-correlation distribution. An alternative fitting procedure based on a convolution of two Gaussian functions for the description of the NS peak was performed for Monte Carlo simulations. The resulting NS yields were found to be compatible with those obtained with the standard procedure, with a maximum variation of 7% (10%) in pp (p–Pb) collisions in case of EPOS 3 simulations.

4 Systematic uncertainties

The fit of the D-meson invariant-mass distribution introduces systematic uncertainties on \(S_\mathrm{peak}\) and \(B_\mathrm{peak}\) (Sect. 3.2, Eq. 1). The uncertainty on the correlation distribution was estimated by calculating \(B_\mathrm{peak}\) from the integral of the background term of the invariant-mass fit function in the range \(|M-\mu |<2\sigma \) and by varying the fit procedure. In particular, the fit was repeated modeling the background distribution with a linear function and a parabola instead of an exponential function (for \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) mesons only), considering a different histogram binning, and varying the fit range. A \(10\%\) systematic uncertainty was estimated from the corresponding variation of the azimuthal-correlation distribution. No significant trend was observed as a function of \(\Delta {\varphi }\) and the same uncertainty was estimated for all D-meson species in all \({p}_{\mathrm{T}}\)-intervals and in both pp and p–Pb collision systems.

A 5% uncertainty (10% for \({\mathrm{D}^{+}}\) mesons in p–Pb collisions) arises from the possible dependence of the shape of \(\tilde{C}(\Delta \varphi ,\Delta \eta )_\mathrm{sidebands}\) on the sideband range. This source of uncertainty was determined by restricting the invariant-mass sideband window to the intervals \(4\sigma < |{M-\mu }| < 6\sigma \) or to \(6\sigma<\langle {M-\mu }\rangle <8\sigma \) for all the D mesons, and also by considering, for \(\mathrm{D}^{0}\) and \({\mathrm{D}^{+}}\) mesons, only the left or only the right sideband.

The uncertainty on the correction for the associated-particle reconstruction efficiency was assessed by varying the selection criteria applied to the reconstructed tracks, removing the request of at least three associated clusters in the ITS, or demanding a hit on at least one of the two SPD layers. A \(\pm 4\%\) uncertainty was estimated for \(\text{ p--Pb }\) collisions, while a \(^{+10\%}_{-5\%}\) contribution was obtained for the pp analysis, with the +10% contribution arising from the request of hits in the SPD. No significant trend in \(\Delta {\varphi }\) was observed.

The uncertainty on the residual contamination from secondary tracks was evaluated by repeating the analysis varying the cut on the DCA in the (xy) plane from \(0.1~\mathrm{cm}\) to \(1~\mathrm{cm}\), and re-evaluating the purity of charged primary particles for each variation. This resulted in a 5% (3.5%) systematic uncertainty in pp (\(\text{ p--Pb }\)) collisions, independent of \(\Delta {\varphi }\) and \({p}_\mathrm{T}^\mathrm{assoc}\).

A 5% systematic effect originating from the correction of the D-meson reconstruction efficiency was evaluated by applying tighter and looser topological selections on the D-meson candidates. No significant dependence on \(\Delta {\varphi }\) was observed and the same uncertainty was estimated for the three D-meson \({p}_{\mathrm{T}}\) intervals, apart from \({\mathrm{D}^{+}}\) meson in \(\text{ p--Pb }\) collisions, for which a 10% uncertainty was assigned.

The uncertainty on the subtraction of the beauty feed-down contribution was quantified by generating the templates of feed-down azimuthal-correlation distributions, \(\tilde{C}_\mathrm{feed-down}^\mathrm{MC~templ}(\Delta \varphi )\) in Eq. 2, with different PYTHIA 6 tunes (Perugia 0, Perugia 2010, see Sect. 2.2), and by considering the range of \(f_\mathrm{{prompt}}\) values obtained by varying the prompt and feed-down D-meson \({p}_{\mathrm{T}}\)-differential production cross sections within FONLL uncertainty band, as described in [15]. The effect on the azimuthal-correlation distributions is \(\Delta {\varphi }\) dependent and contained within \(8\%\) and is more pronounced in the near side, in particular in the low and mid D-meson \({p}_{\mathrm{T}}\) intervals.

The consistency of the whole correction procedure, prior to the feed-down subtraction, was verified by performing the analysis on simulated events (“Monte-Carlo closure test”) separately for prompt and feed-down D mesons. For prompt D mesons, no effect was found for both pp and \(\text{ p--Pb }\) collision systems. Conversely, for feed-down D mesons, an overestimate by about 20% in the near side was found for both collision systems. It was verified that the source of this excess is related to a bias induced by the topological selection applied to D mesons, that tends to favour cases with a small angular opening between the products of the beauty-hadron decay, thus between the D meson and the other decay particles. This effect results in a \(\Delta {\varphi }\)-dependent overestimate of the feed-down subtracted correlation distribution in the near side, contained within 2%.

The systematic uncertainties affecting the \(\Delta {\varphi }\)-correlation distributions are summarised in Table 1 for both pp and \(\text{ p--Pb }\) collision systems. The \(\Delta {\varphi }\)-dependent parts of the uncertainties arising from the feed-down subtraction and the Monte-Carlo closure test define the \(\Delta {\varphi }\)-uncorrelated systematic uncertainties. All the other contributions, correlated in \(\Delta {\varphi }\), act as a scale uncertainty. No significant dependence on the transverse momentum of D mesons and associated particles was observed for both \(\Delta {\varphi }\)-correlated and uncorrelated uncertainties, except for the feed-down systematic uncertainty.

Table 1 List of systematic uncertainties for the \(\Delta {\varphi }\)-correlation distributions in pp and \(\text{ p--Pb }\) collisions. See text for details

Different approaches were applied to estimate the systematic uncertainty on the near-side peak associated yield and peak width and on the baseline, obtained from the \(A_\mathrm{NS}\), \(\sigma _{\mathrm{fit,NS}}\), and b parameters of the fit of the azimuthal-correlation distribution, as described in Sect. 3.3. The main source of uncertainty originates from the definition of the baseline itself, which is connected to the assumption that the observed variation of the azimuthal-correlation distribution in the transverse region is determined mainly by statistical fluctuations rather than by the true physical trend. The variation of \(A_\mathrm{NS}\), \(\sigma _{\mathrm{fit,NS}}\), and b values obtained when considering a \(\pm \pi /4\) variation of the \(\Delta {\varphi }\) range defining the transverse region is interpreted as the systematic uncertainty due to the baseline definition. In addition, the fits were repeated by moving upwards and downwards the data points by the corresponding value of the \(\Delta {\varphi }\)-uncorrelated systematic uncertainty. The final systematic uncertainty was calculated by summing in quadrature the aforementioned contributions and, for the associated yields and baseline, also the systematic uncertainty correlated in \(\Delta \varphi \). The values of the total systematic uncertainties on the near-side peak yield, width, and baseline are reported in Table 2, for two intervals of transverse momentum of D mesons and associated particles. Considering all the measured kinematic ranges, the uncertainties vary from \(\pm 12\) to \(\pm 25\%\) for the near-side peak yield, from \(\pm 2\) to \(\pm 13\%\) for the near-side peak width and from \(\pm 11\) to \(\pm 16\%\) for the baseline. Typically, lower uncertainties are obtained for p–Pb collisions, where the larger available statistics of the correlation distributions allow for a more precise estimate of the baseline height, which constitutes the main source of uncertainty also on the evaluation of the near-side peak associated yield and width.

Table 2 List of systematic uncertainties for near-side (NS) peak associated yield, near-side peak width, and baseline in pp and \(\text{ p--Pb }\) collisions, for two different kinematic ranges of D mesons and associated particles. See text for details

5 Results

Fig. 2
figure 2

Comparison of the azimuthal-correlation distributions of D mesons with charged particles obtained for \(\mathrm{D}^{0}\), \({\mathrm{D}^{+}}\) and \({\mathrm{D}^{*+}}\) mesons for \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\), \({p}_\mathrm{T}^\mathrm{assoc}> 1~\mathrm {GeV}/c\) in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) (left panel) and for \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\), \({p}_\mathrm{T}^\mathrm{assoc}> 1~\mathrm {GeV}/c\) in \(\text{ p--Pb }\) collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\) (right panel). The statistical uncertainties are shown as error bars, the \(\Delta {\varphi }\)-uncorrelated systematic uncertainties as boxes, while the part of systematic uncertainty correlated in \(\Delta {\varphi }\) is reported as text (scale uncertainty). The latter is largely uncorrelated among the D-meson species

Fig. 3
figure 3

Average of the azimuthal-correlation distributions of \(\mathrm{D}^{0}\), \({\mathrm{D}^{+}}\) and \({\mathrm{D}^{*+}}\) mesons with \(3<{p}_\mathrm{T}^\mathrm{D}<5~\mathrm {GeV}/c\) (left column), \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\) (middle column), and \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\) (right column), with charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>0.3~\mathrm {GeV}/c\) (top row), \(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\) (middle row), and \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) (bottom row), measured in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and in p–Pb collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\). The statistical uncertainties are shown as error bars, the \(\Delta {\varphi }\)-uncorrelated systematic uncertainties as boxes, while the part of systematic uncertainty correlated in \(\Delta {\varphi }\) is reported as text (scale uncertainty)

Fig. 4
figure 4

Comparison of the azimuthal-correlation distributions of D mesons with \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\) (left column) and \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\) (right column) with charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>0.3~\mathrm {GeV}/c\) (top row), \(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\) (middle row), and \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) (bottom row) in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and in p–Pb collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\), after baseline subtraction. The statistical uncertainties are shown as error bars, the \(\Delta {\varphi }\)-uncorrelated systematic uncertainties as boxes around the data points, the part of systematic uncertainty correlated in \(\Delta {\varphi }\) is reported as text (scale uncertainty), the uncertainties deriving from the subtraction of the baselines are represented by the boxes at \(\Delta {\varphi }>\pi \)

Fig. 5
figure 5

Examples of the fit to the azimuthal-correlation distribution, for D mesons with \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\) with charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) (left), and for D mesons with \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\) with charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) in p–Pb collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\) (right). The statistical uncertainties are shown as error bars, the \(\Delta {\varphi }\)-uncorrelated systematic uncertainties as boxes, while the part of systematic uncertainty correlated in \(\Delta {\varphi }\) is reported as text (scale uncertainty). The terms of the fit function described in Sect. 3.3 are also shown separately: near-side Gaussian function (blue dashed line), away-side Gaussian function (green dashed–dotted line) and baseline constant term (magenta dotted line)

Fig. 6
figure 6

Comparison of the near-side peak associated yield (top row) and peak width (bottom row) in pp and \(\text{ p--Pb }\) collisions as a function of \({p}_\mathrm{T}^\mathrm{D}\), for \({p}_\mathrm{T}^\mathrm{assoc}>0.3~\mathrm {GeV}/c\) (left column), \(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\) (middle column), and \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) (right column). The points and error boxes for p–Pb collisions are shifted by \(\Delta {p}_{\mathrm{T}}= +0.3~\mathrm {GeV}/c\). Statistical and systematic uncertainties are shown as error bars and boxes, respectively

Fig. 7
figure 7

Comparison of \(\Delta {\varphi }\)-correlation distributions of D mesons with charged particles measured in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and Monte-Carlo simulations performed with different event generators, after the subtraction of the baseline. The statistical and systematic uncertainties of the measured distributions are displayed as in Fig. 4

The azimuthal-correlation distributions of \(\mathrm{D}^{0}\), \({\mathrm{D}^{+}}\), \({\mathrm{D}^{*+}}\) mesons with charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) are compared in Fig. 2 for \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\) in pp collisions (left panel) and for \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\) in \(\text{ p--Pb }\) collisions (right panel). The distributions obtained with the three D-meson species are compatible within the quadratic sum (\(w_{i}\), \(i=\mathrm{D}^{0}\), \({\mathrm{D}^{+}}\), \({\mathrm{D}^{*+}}\)) of the statistical uncertainty and of the systematic uncertainties on the signal, background normalisation, and on the background shape (see Table 1), that are uncorrelated among the three meson species. The \(\mathrm{D}^{0}\)-, \({\mathrm{D}^{+}}\)-, \({\mathrm{D}^{*+}}\)-meson data are averaged using \(1/w_{i}^{2}\) as weights. The averages of the distributions are shown, for all the considered kinematic ranges, in Fig. 3 for pp and \(\text{ p--Pb }\) collisions. A rising trend of the height of the near-side peak with increasing D-meson \({p}_{\mathrm{T}}\) is observed for both collision systems. A similar trend is present for hadron–hadron correlations measured at Tevatron and LHC energies [68,69,70,71]: an increase of hadron multiplicity in jets with increasing jet energy is expected from the evolution of parton cascade with the parton energy for both light and heavy quarks [19]. A decrease of the baseline level with increasing \({p}_{\mathrm{T}}\) of the associated particles can also be noticed.

Figure 4 shows the \(\Delta {\varphi }\) distributions after the subtraction of the baseline, calculated as described in Sect. 3.3. The distributions show a near-side peak along with a wider and lower peak in the away-side region. The results obtained for the two collision systems are compatible within the total uncertainties. According to simulations of pp collisions performed using PYTHIA 6 (Perugia-0, -2010, and -2011 tunes), the different centre-of-mass energy and the slightly different D-meson rapidity range of the two measurements should induce variations in the baseline-subtracted azimuthal-correlation distributions smaller than 7% in the near- and away-side regions. The same estimate is obtained with POWHEG+PYTHIA simulations including the EPS09 parametrisation of nuclear PDFs (see Sect. 2.2). Such differences are well below the current level of uncertainties.

A further comparison of the results from pp and p–Pb collisions is done by quantifying the integrals and the widths of the near-side correlation peaks by fitting the measured distributions as described in Sect. 3.3. The fit results are reported only for the near-side peak parameters and the baseline because of the poor statistical precision on the fit parameters of the away-side peaks. Figure 5 shows an exemplary fit to the azimuthal-correlation distributions of D mesons with charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\), for \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\) in pp collisions (left panel) and for \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\) in p–Pb collisions (right panel). The curves superimposed to the data represent the three terms of the function defined in Eq. 3.

Within the uncertainties, the fit function describes the measured distributions in all kinematic cases considered, yielding \(\chi ^{2}/\mathrm{NDF}\) values close to unity. The evolution of the near-side peak associated yield as a function of the D-meson \({p}_{\mathrm{T}}\) is reported in Fig. 6 (top row), for pp and \(\text{ p--Pb }\) collisions, for \({p}_\mathrm{T}^\mathrm{assoc}>0.3~\mathrm {GeV}/c\) (left panel) and for the two sub-intervals \(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\) (middle panel) and \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\) (right panel). The near-side peak associated yield exhibits an increasing trend with D-meson \({p}_{\mathrm{T}}\) and has similar values, within uncertainties, for the softer (\(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\)) and the harder (\({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\)) sub-ranges of \({p}_\mathrm{T}^\mathrm{assoc}\) used, in each D-meson \({p}_{\mathrm{T}}\) interval considered. The values obtained for pp and \(\text{ p--Pb }\) collision data are compatible within statistical uncertainties. In the bottom row of the same figure the width of the near-side Gaussian term (\(\sigma _\mathrm{fit,NS}\)) is shown. Although the case with \({p}_\mathrm{T}^\mathrm{assoc}>0.3~\mathrm {GeV}/c\) seems to suggest that \(\sigma _\mathrm{fit,NS}\) does not strongly depend on D-meson \({p}_{\mathrm{T}}\) in the range of the measurement, the current level of uncertainty does not allow for quantification of the dependence of \(\sigma _\mathrm{fit,NS}\) on D-meson and associated charged particle \({p}_{\mathrm{T}}\), as well as any potential difference between the values extracted using pp and p–Pb data. In particular, our approach for baseline calculation (Sect. 3.3) guarantees a robust estimate of the minimum, but the baseline uncertainty and its impact on the associated-yield uncertainty are rather large (Sect. 4). This systematic uncertainty is expected to be significantly reduced in future measurements with larger data samples, where a smaller \(\Delta {\varphi }\) range for the baseline calculation could be used.

A \({v}_{2}\)-like modulation of the baseline would introduce a bias in the measurement of the associated yield and peak width and that needs to be taken into account while interpreting the measured quantities in terms of charm-jet properties. In order to get an estimate of this possible effect, for the p–Pb case the fit was repeated by subtracting from the correlation distribution a \({v}_{2}\)-like modulation assuming \({v}_{2}=0.05\) for D mesons and \({v}_{2}=0.05~(0.1)\) for associated charged particles with \({p}_{\mathrm{T}}>0.3~(1)~\mathrm {GeV}/c\). These values were chosen on the basis of charged-particle measurements in high-multiplicity p–Pb collisions [30] and assuming for D mesons the maximum value predicted in [2] for the 20% most central p–Pb collisions as a test case. With such assumptions, rather extreme also considering that this measurement is performed without any selection on event multiplicity, \(A_\mathrm{NS}\) varies by \(-10\%\) (\(-6\%\)) for D mesons with \(5<{p}_{\mathrm{T}}<8~\mathrm {GeV}/c\) and for \(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\) (\({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\)). The variations on \(\sigma _{\mathrm{fit,NS}}\) and on the baseline are below 4 and 1%, respectively. Significantly smaller modifications result for D mesons with \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\). With the available statistics, the precision of the measurement is not sufficient to observe or exclude these modifications.

Figure 7 shows the comparison of the averaged azimuthal-correlation distributions measured in pp collisions with expectations from simulations performed with PYTHIA, POWHEG+PYTHIA, and EPOS 3 (see Sect. 2.2), after the baseline subtraction. The average of the two lowest values of the azimuthal-correlation distribution is used to define the uncertainty related to the baseline definition in Monte-Carlo simulations (see Sect. 3.3). This uncertainty is negligible and not displayed in the figures. The distributions obtained with the different generators and tunes do not show significant differences in the near side, except from EPOS 3 which tends to have higher and wider distributions. In the away side, the PYTHIA 6 tunes Perugia 0 and Perugia 2010 tend to have higher correlation values, especially for \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\), compared to the other simulation results. Similar considerations hold for EPOS 3 in the case of D mesons with \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\). The considered Monte-Carlo simulations describe, within the uncertainties, the data in the whole \(\Delta {\varphi }\) range. The comparison of the associated yield in the near-side peak in data and in simulations is displayed in the top row of Figs. 8 and 9, for pp and p–Pb collisions, respectively. The simulations obtained with EPOS 3 provide a better description of the near-side yields for D mesons with \(8<{p}_{\mathrm{T}}<16~\mathrm {GeV}/c\) in both pp and p–Pb collisions. At lower D-meson \({p}_{\mathrm{T}}\) a better agreement is obtained with PYTHIA and POWHEG+PYTHIA simulations. The width of the near-side peaks, shown in the second row of the same figures, is better reproduced by the simulations in the case of p–Pb than of pp results. The evolution of the baseline value as a function of the D-meson \({p}_{\mathrm{T}}\) is compared for pp-collision data to expectations from PYTHIA simulations in the bottom row of Fig. 8 for the three ranges of \({p}_\mathrm{T}^\mathrm{assoc}\) considered in the analysis. The value of the baseline, mainly determined by the event multiplicity, does not show substantial variations as a function of D-meson \({p}_{\mathrm{T}}\), as expected also from PYTHIA and EPOS 3 simulations, which reproduce the observed values within the uncertainties.

Fig. 8
figure 8

Comparison of near-side peak associated yield (top row), near-side peak width (middle row), and baseline (bottom row) values measured in pp collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) with the expectations from simulations performed with different Monte-Carlo event generators. Statistical and systematic uncertainties are shown as error bars and boxes, respectively

Fig. 9
figure 9

Comparison of near-side peak associated yield (top row) and near-side peak width (bottom row) values measured in p–Pb collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\) with the expectations from simulations performed with different Monte-Carlo event generators. Statistical and systematic uncertainties are shown as error bars and boxes, respectively

6 Summary

The first measurements of the azimuthal correlations between D mesons with charged particles in pp and p–Pb collisions at \({\sqrt{s}}=7~\mathrm {TeV}\) and \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\), respectively, performed with the ALICE apparatus at the LHC were presented. The \(\Delta {\varphi }\) distributions were studied in pp collisions in three different D-meson transverse-momentum intervals, \(3<{p}_\mathrm{T}^\mathrm{D}<5~\mathrm {GeV}/c\), \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\), and \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\), for associated charged particles with \({p}_\mathrm{T}^\mathrm{assoc}>0.3~\mathrm {GeV}/c\), and in the two sub-ranges \(0.3<{p}_\mathrm{T}^\mathrm{assoc}<1~\mathrm {GeV}/c\) and \({p}_\mathrm{T}^\mathrm{assoc}>1~\mathrm {GeV}/c\). For p–Pb collisions, the results were reported in two D-meson \({p}_{\mathrm{T}}\) ranges, \(5<{p}_\mathrm{T}^\mathrm{D}<8~\mathrm {GeV}/c\), and \(8<{p}_\mathrm{T}^\mathrm{D}<16~\mathrm {GeV}/c\). The baseline-subtracted azimuthal-correlation distributions observed in the two collision systems are compatible within uncertainties. The variations expected from the lower nucleon-nucleon centre-of-mass energy of p–Pb collisions and from the slightly different D-meson rapidity ranges used for the p–Pb analysis were studied with simulated pp collisions at the two centre-of-mass energies and are well below the sensitivity of the measurements.

The properties of the near-side correlation peak, sensitive to the characteristics of the jet containing the D meson, were described in terms of the yield of associated charged particles and peak width, obtained by fitting the \(\Delta {\varphi }\) distributions with a function composed of a constant term, representing the physical minimum of the distribution, and two Gaussian terms modeling the near- and away-side peaks. The values measured in the two collision systems are compatible within uncertainties.

The measured azimuthal distributions, as well as the properties of the correlation peaks, were compared to expectations from simulations performed with different Monte-Carlo generators. The simulations reproduce the correlation distributions within uncertainties.

Considering that the overall uncertainty is dominated by the statistical component, the data collected from pp collisions at \({\sqrt{s}}=13~\mathrm {TeV}\) in the ongoing Run 2 at the LHC will allow for a more precise measurement. In particular, the predicted increase of the cross section for charm production by more than a factor 2 at \({p}_{\mathrm{T}}=10~\mathrm {GeV}/c\) at the higher collision energy [12], along with the foreseen larger integrated luminosity, will allow for a significant reduction of the statistical uncertainty, providing a more quantitative and constraining comparison of the data with expectations from Monte-Carlo generators. As mentioned in Sect. 5, with larger data samples a different determination of the baseline of the azimuthal-correlation distribution will become possible, bringing to a significant reduction of the systematic uncertainty on the measurement of the associated yields. The data that will be collected in next p–Pb collision runs at the LHC may also allow for a study of the evolution of the azimuthal-correlation distribution as a function of the event multiplicity, searching for possible long-range ridge-like structures already observed with angular correlation of light particles.

The results reported in this paper represent a first step towards the measurement of possible modifications concerning the azimuthal correlation of D mesons with charged particles in Pb–Pb collisions, which has the potential to provide important information on the charm-quark energy-loss mechanisms in the presence of the medium formed in heavy-ion collisions at LHC energies. Given the same collision energy, the p–Pb results presented in this paper could serve as a reference to study medium effects in Pb–Pb collisions at \({\sqrt{{{s}}_{\scriptscriptstyle {\mathrm{NN}}}}}=5.02~\mathrm {TeV}\) collected during the LHC Run 2.