Measurements of $e^+e^-$ pairs from open heavy flavor in $p$+$p$ and $d$+$A$ collisions at $\sqrt{s_{NN}}=200$ GeV

TL;DR

This study measures electron-positron pairs from heavy-flavor decays in proton-proton collisions at 200 GeV, comparing data with simulations and examining potential nuclear effects.

Contribution

It provides the first detailed separation of $bar{b}$ and $car{c}$ contributions in $p$+$p$ collisions at this energy using multiple event generators.

Findings

All three generators describe the data within detector acceptance.

Significant differences in total cross sections when extrapolating to full phase space.

No observable nuclear modification in $p$+$p$ versus $d$+$A$ collisions.

Abstract

We report a measurement of $e^+e^-$ pairs from semileptonic heavy-flavor decays in $p$+$p$ collisions at $\sqrt{s_{NN}}=200$~GeV. The $e^+e^-$ pair yield from $b\bar{b}$ and $c\bar{c}$ is separated by exploiting a double differential fit done simultaneously in dielectron invariant mass and $p_T$. We used three different event generators, {\sc pythia}, {\sc mc@nlo}, and {\sc powheg}, to simulate the $e^+e^-$ spectra from $c\bar{c}$ and $b\bar{b}$ production. The data can be well described by all three generators within the detector acceptance. However, when using the generators to extrapolate to $4\pi$, significant differences are observed for the total cross section. These difference are less pronounced for $b\bar{b}$ than for $c\bar{c}$. The same model dependence was observed in already published $d$+$A$ data. The $p$+$p$ data are also directly compared with $d$+$A$ data in mass and $p_T$, and within the statistical accuracy no nuclear modification is seen.

Tables

Table 1. Table 1: Fit parameters for p 𝑝 p + + p 𝑝 p collisions derived from a simultaneous fit to the π 0 superscript 𝜋 0 \pi^{0} and charged pions data using the modified Hagedorn function (Eq.  4 ).

Parameter	Value
$dN/dy$	1.139 $\pm$ 0.10
$A$[$mbGe{V}^{-2}{c}^3$]	492 $\pm$ 67
$a$[$(GeV/$c$){}^{-}1$]	0.266 $\pm$ 0.031
$b$[$(GeV/$c$){}^{-}2$]	0.092 $\pm$ 0.021
$p_0$[$GeV/$c]	0.68 $\pm$ 0.02
$n$	8.27 $\pm$ 0.07

Table 2. Table 2: Rapidity density for the mesons extracted from the fits and used in the EXODUS decay generator.

Meson	${dN/dy|}_{y=0}$	Data source
${\pi }^0,{\pi }^{+},{\pi }^{-}$	1.139 $\pm$ 0.10	Adler et al. (2003b); Adare et al. (2007a, 2011b)
$\eta$	0.093 $\pm$ 0.0002	Adler et al. (2007); Adare et al. (2011c)
$\omega$	0.0744 $\pm$ 0.0017	Adare et al. (2011d)
$\varphi$	0.009 $\pm$ 0.0002	Adare et al. (2011e, d)
${\eta }^{\prime }$	0.0123 $\pm$ 0.0008	Adare et al. (2011d)
$J/\psi$	1.74$\times {10}^{-05}\pm$ 5.1$\times {10}^{-7}$	Adare et al. (2007b, 2012b)
${\psi }^{\prime }$	3.1$\times {10}^{-06}\pm$ 6.2$\times {10}^{-7}$	Adare et al. (2012b)

Table 3. Table 3: Summary of the various systematic uncertainties considered in this analysis.

Source	Syst. uncertainty
	(mass $\le 1.0$GeV/$c^2$)	(mass $>1.0$ GeV/$c^2$)
Data systematics
eID	15%	10%
Input model	15%	15%
ERT	10%	5%
Fiducials	10%
$\alpha -$ correction	5%
BBC bias	10%
Cocktail systematics
Hadronic cocktail	20%
$c\overline{c}$ cross section	32%
$b\overline{b}$ cross section	36%

Table 4. Table 4: Summary of c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} and b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} cross section measured in p 𝑝 p + + p 𝑝 p collisions using three different generators pythia , mc@nlo , and powheg . These are obtained by extrapolating to 4 ​ π 4 𝜋 4\pi the fitting results from the measured e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} pairs from heavy flavor.

$p$$+$$p$	pythia ($\mu$b)	mc@nlo ($\mu$b)	powheg ($\mu$b)
$c\overline{c}$	356 $\pm$ 27 (stat) $\pm$ 89(syst)	708 $\pm$ 55 (stat) $\pm$ 175 (syst)	267 $\pm$ 19 (stat) $\pm$ 67 (syst)
$b\overline{b}$	4.81 $\pm$ 0.71 (stat)$\pm$ 1.00 (syst)	3.85 $\pm$ 0.73 (stat)$\pm$ 0.8 (syst)	2.91 $\pm$ 0.63 (stat)$\pm$ 0.61 (syst)

Table 5. Table 5: Step by step extrapolation from the number of e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} pairs for m e + ​ e − subscript 𝑚 superscript 𝑒 superscript 𝑒 m_{e^{+}e^{-}} ≥ \geq 1.16 GeV/ c 2 superscript 𝑐 2 c^{2} from c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} in the PHENIX acceptance to the number of c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} pairs in 4 ​ π 4 𝜋 4\pi for pythia , mc@nlo , and powheg . Numbers are in units of pairs per event using the c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} cross sections determined in this paper. The factors in brackets quantify the increase in number of pairs. We have factored out the effective branching ratio BR=0.094 for decays of c → → 𝑐 absent c\rightarrow e in the step from e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} to c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} pairs. The number of c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} pairs in 4 ​ π 4 𝜋 4\pi is equal to the c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} cross section in table  4 divided by the inelastic p 𝑝 p + + p 𝑝 p cross section σ p ​ p = 42 subscript 𝜎 𝑝 𝑝 42 \sigma_{pp}=42 mb.

$c\overline{c}$	pythia	mc@nlo	powheg
	3.20$\times$10-8	3.55$\times$ 10-8	3.61$\times$10-8
$m_{e^{+}{e}^{-}}$ $\ge$ 1.16 GeV/$c^2$
	1.66$\times$10-7 (5.19)	2.55 $\times$10-7 (7.18)	1.93$\times$10-7 (5.33)
$\left|y_{c\overline{c}}\right|\le$0.5	2.33$\times$10-3 (124/BR2)	5.09 $\times$10-3 (176.6/BR2)	1.80$\times$10-3 (82.5/BR2)
$4\pi$	8.48$\times$10-3 (3.64)	16.9 $\times$10-3 (3.31)	6.36$\times$10-3 (3.53)

Table 6. Table 6: Step by step extrapolation from the number of e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} pairs for m e + ​ e − subscript 𝑚 superscript 𝑒 superscript 𝑒 m_{e^{+}e^{-}} ≥ \geq 1.16 GeV/ c 2 superscript 𝑐 2 c^{2} from b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} in the PHENIX acceptance to the number of b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} pairs in 4 ​ π 4 𝜋 4\pi for pythia , mc@nlo , and powheg . Numbers are in units of pairs per event using the b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} cross sections determined in this paper. The factors in brackets quantify the increase in number of pairs. We have factored out the effective branching ratio BR=0.158 for decays of b → → 𝑏 absent b\rightarrow e in the step from e + ​ e − superscript 𝑒 superscript 𝑒 e^{+}e^{-} to b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} pairs. The number of b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} pairs in 4 ​ π 4 𝜋 4\pi is equal to the b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} cross section in table  4 divided by the inelastic p 𝑝 p + + p 𝑝 p cross section σ p ​ p = 42 subscript 𝜎 𝑝 𝑝 42 \sigma_{pp}=42 mb.

$b\overline{b}$	pythia	mc@nlo	powheg
	10.3$\times$10-9	8.34$\times$10-9	6.99$\times$10-9
$m_{e^{+}{e}^{-}}$ $\ge$ 1.16 GeV/$c^2$
	2.18$\times$10-8 (2.11)	1.83 $\times$10-8 (2.19)	1.46$\times$10-8(2.12)
$\left|y_{b\overline{b}}\right|\le$0.5	4.47$\times$10-5 (51.1/BR2)	3.49 $\times$10-5 (47.6/BR2)	2.61$\times$10-5 (44.6/BR2)
$4\pi$	11.5$\times$10-5 (2.56)	9.17 $\times$10-5 (2.62)	6.93$\times$10-5 (2.66)

Table 7. Table 7: Summary of c ​ c ¯ 𝑐 ¯ 𝑐 c\bar{c} and b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} cross section in d 𝑑 d + + Au collisions expressed as nucleon-nucleon equivalent cross section by dividing the d 𝑑 d + + Au cross section by the average number of binary nucleon-nucleon collisions N coll subscript 𝑁 coll N_{\rm coll} = 7.6 ± 0.4 absent plus-or-minus 7.6 0.4 =7.6\pm 0.4 .

$d$$+$Au/$N_{\mathrm{coll}}$	pythia ($\mu$b)	mc@nlo ($\mu$b)	powheg ($\mu$b)
$c\overline{c}$ (Reanalysis)	385 $\pm$ 34 (stat) $\pm$ 119 (syst)	795 $\pm$ 80 (stat) $\pm$ 275 (syst)	303 $\pm$ 26 (stat) $\pm$ 94 (syst)
$b\overline{b}$ (Reanalysis)	3.40 $\pm$ 0.65 (stat)$\pm$ 1.10 (syst)	2.95 $\pm$ 0.67 (stat)$\pm$ 0.95 (syst)	2.0 $\pm$ 0.6 (stat)$\pm$ 0.65 (syst)

Equations

$$\alpha (m,p_T)$$

$$\mbox{$S_{+-}$}(m,p_{T})=\mbox{$N_{+-}$}(m,p_{T})-\alpha(m,p_{T})\times\mbox{$N_{\pm\pm}$}(m,p_{T})$$

$$\displaystyle\frac{d^{2}N_{\mbox{$e^{+}e^{-}$}}}{dm\ dp_{T}}$$

$$E\frac{d{\sigma }^3}{d{p}^3}=A(e^{-(a{p}_T+b{p}_T^2)}+p_T/p_0{)}^{-n}$$

$${\mu }_0^2=\frac{1}{2}(m_T^2(Q)+m_T^2(\bar{Q})),$$

$${\mu }_0^2=p_T^2+m^2,$$

Full text

I Introduction

Heavy quarks such as charm and bottom are excellent probes to understand the properties of the quark gluon plasma (QGP) created in high energy heavy-ion collisions. Both charm and bottom quarks have masses significantly larger than the quantum chromodynamics (QCD) scale parameter $\Lambda_{QCD}\approx$ 0.2 GeV, and as such, their production is limited to the primordial nucleon-nucleon collisions. Heavy flavor production in the subsequent early, hot stages of heavy-ion collisions is not significant and thus any modification of the primordial heavy flavor phase space distributions in heavy ion collisions will be the result of the quarks traversing the QGP and later phases in the space time evolution.

Prior to the studies of heavy flavor production done at the Relativistic Heavy Ion Collider (RHIC), the high $p_{T}$ suppression Adcox et al. (2002); Adler et al. (2003a); Adams et al. (2003) of light flavor hadrons was primarily associated to radiative energy loss via medium-induced gluon radiation. This predicted a distinctive mass hierarchy of high $p_{T}$ suppression as measured via the nuclear modification factor $R_{AA}$, implying that hadrons with heavy flavor will have a smaller suppression: $R_{AA}^{\pi^{0}}<R_{AA}^{c}<R_{AA}^{b}$. $R_{AA}^{\pi^{0}}$ denotes the nuclear modification of $\pi^{0}$, defined as the ratio of yield measured in $AA$ collisions to the yield measured in p$$+$$p collisions scaled by the number of binary collisions for $AA$ system, and $R_{AA}^{c}$(or $R_{AA}^{b}$) denotes the same for charm (or bottom) quarks. However, the measurements showed similar suppression for light and heavy flavor hadrons. Including collisional energy loss via elastic scattering, which is more important for heavy flavor than for the light quarks, leads to a qualitative explanation of the large energy loss for heavy flavor  He et al. (2012); Sharma and Vitev (2013). But other approaches are similarly successful, including Langevin-based transport models Moore and Teaney (2005); van Hees et al. (2006) and AdS/CFT (anti de-Sitter-space/conformal field theory) string drag energy loss models Horowitz and Gyulassy (2008). Despite significant effort, a full quantitative understanding of the energy loss has not been achieved yet.

To test different theoretical approaches, it is crucial to understand primordial heavy flavor production, and any modifications there in the presence of nuclei. Primordial heavy flavor production can be studied in p$$+$$p collisions. When nuclei are involved in a collision, one might expect modifications to the initial state, which can be described as shadowing or anti shadowing of the parton distribution functions. Also modifications in the final state that can be expressed as changes of the fragmentation process are possible, for example, via energy loss or re-scattering in cold nuclear matter. It is commonly accepted that these effects are observable in p(d)$$+A collisions, where QGP formation is not expected. Differences between the single electron spectra from heavy flavor decays from d$$+Au data and p$$+$$p data have been interpreted as cold nuclear matter effects Adare et al. (2012a).

Recently hints of collectivity have been found in high multiplicity events from collisions of small nucleus with a large nuclei, which suggests that hot matter might even be formed in small systems. However, one would not expect sizable collective effects on the heavy flavor phase space distributions even if hot matter is created due to the small reaction volume in these collisions.

The primordial heavy flavor production can be calculated in the framework of perturbative QCD (pQCD). Therefore, measurement of heavy flavor in p$$+$$p serves as a test for these calculations and can be used to improve Monte-Carlo (MC) generators. Results from MC generators can be scaled to A$$+$$A or p(d)$$+A collision systems with the number of binary collisions and serve as a reference for observables in the absence of p$$+$$p data.

At RHIC, open heavy flavor production has been measured by both the PHENIX and STAR experiments in different collision systems, spanning p$$+$$p, d$$+Au, Cu$+$Cu and Au$+$Au systems, and by exploiting various techniques such as single electrons/muons via semi-leptonic decays Adare et al. (2011a, 2012a), electron-hadron correlations Aggarwal et al. (2010), $e-\mu$ Adare et al. (2014), $e^{+}e^{-}$ Adare et al. (2015) and also via reconstruction of $D$-mesons Adamczyk et al. (2014). This paper reports the measurement of heavy flavor production via dielectrons in p$$+$$p collisions at midrapidity. The $e^{+}e^{-}$ pairs coming from the semi-leptonic decays of charm and bottom dominate different regions in mass and $p_{T}$ allowing to disentangle the two contributions. Studying the $e^{+}e^{-}$ pairs from heavy flavor may also provide sensitivity to the heavy quark correlations which is important to constrain the MC models. The results from the p$$+$$p data from this paper can be directly compared to the previously published d$$+Au data Adare et al. (2015) that exploited the same technique.

The paper is organized as follows: Section II describes the experimental apparatus and trigger. Section III details the data analysis including electron identification, background subtraction, and efficiency corrections. A description of the hadronic cocktail and heavy flavor generators is outlined in Section IV, followed by studies of systematic uncertainties in Section V. The data are presented as double differential spectra in mass and $p_{T}$ in Section VI. The final results and the comparison of p$$+$$p and d$$+Au, as well as the comparison to several models of charm and bottom production are discussed in Section VII. Section VIII gives our summary and conclusions.

II Experiment

A detailed description of the PHENIX detector is available in Adcox et al. (2003a). We focus here on the components of the two central arm spectrometers and the beam-beam counters (BBCs) that are critical for the analysis of $e^{+}e^{-}$ pairs. Each of the two central arms cover a pseudorapidity range of $|\eta|<$0.35 (70${}^{\circ}<\theta<$110∘) and 90∘ in azimuthal angle $\phi$. They are located almost back-to-back, with an angular gap of 67.5∘ between them at the top. They span a range from about 220 cm to 500 cm radially from the beam axis. The location of collision vertex in the beam direction, the collision time, and the minimum bias (MB) trigger are provided by a system of two beam-beam counters (BBC) that are located at a distance of 144 cm from the nominal interaction point on either side. Each BBC covers the full azimuth and a rapidity range of 3.1 $<|\eta|<$ 3.9. The collision vertex resolution in the beam direction is approximately 2 cm for p$$+$$p collisions. The MB trigger requires a coincidence between both sides with at least one hit on each side, and accepts the events if the BBC vertex is within 38 cm of the nominal interaction point. The BBC cross section in p$$+$$p collisions was determined via the van der Meer scan technique Dre (2001) and was found to be \sigma_{\rm BBC}^{\mbox{p$+$p}}~{}=~{}23.0~{}\pm~{}2.2 mb or 0.545 $\pm$ 0.06 of the total inelastic p$$+$$p cross section of \sigma_{inel}^{\mbox{p$+$p}}~{}=~{}42~{}\pm~{}3 mb.

There are two primary charged particle tracking subsystems in PHENIX: drift chambers (DC) and pad chambers (PC) Adcox et al. (2003b). The DC along with first layer of PC (PC1) form the inner tracking system used here. The DC measures the trajectories of charged particles in the plane perpendicular to the beams and allows one to determine their charge and transverse momentum $p_{T}$. The PC1 provides a space point along the trajectory of charged particles, which is used to determine the polar angle $\theta$ and $z$-coordinate of the track. The momentum resolution for this data set is $\delta p/p$ = 0.011 $\oplus$ 0.0116$\,p$[GeV/$c$].

Each central arm is equipped with a ring imaging Čerenkov (RICH) detector that serves as the primary device for electron identification. With $CO_{2}$ as a radiator gas, an $e/\pi$ rejection of better than one part in $\sim$103 is achieved, for the tracks with momenta below the pion Čerenkov threshold of $\sim$ 4.87 GeV/$c$. For each electron on average 10 Čerenkov photons are reconstructed on a ring of 11.8 cm diameter with an array of photo multiplier tubes. Further electron identification is provided by the electromagnetic calorimeters (EMCal) that measure the spatial position and energy of the electrons. This is achieved by placing a cut on the ratio of the energy measured by EMCal and momentum given by the DC Adare et al. (2011a).

To select potentially interesting events containing electrons, PHENIX uses a hardware trigger known as ERT (EMCal-RICH) trigger. The trigger is based on the online sum of the energy signals in a tile of $2\times 2$ EMCal towers Aphecetche et al. (2003). For all EMCal trigger tiles above a predetermined threshold value, the location of the EMCal tile is matched with hits in the corresponding RICH tile ($4\times 5$ PMTs). The location of the RICH tile depends on the energy of the trigger particle and is determined from a look-up table, assuming that the trigger particle is an electron. If a spatial match is found, an ERT trigger is issued.

III Data Analysis

The data reported in this paper were collected during the 2006 RHIC p$$+$$p run. The data were recorded with the PHENIX detector using a MB trigger and the ERT trigger. The ERT energy threshold $E_{th}$ was set to 400 MeV for majority of the run, but was raised to 600 MeV towards the end of the run. A total of 855 million ERT triggered events corresponding to 143 billion sampled MB events were analyzed. The corresponding integrated luminosity is $6.6~{}p$b*-1*.

III.1 Event selection and electron identification cuts

The p$$+$$p analysis described here is very similar to the analysis of $e^{+}e^{-}$ pairs from d$$+Au collisions published in Adare et al. (2015). A detailed description of electron identification as well as pair cuts can be found in Adare et al. (2015); Sharma (2010). Events selected were required to have a reconstructed $z$-vertex within 30 cm of the nominal interaction point. Charged tracks reconstructed using the DC and PC1 must pass stringent quality cuts and an explicit cut of $p_{T}$ $>$ 0.2 GeV/$c$. The track is then selected as an electron if at least two photomultiplier tubes registered Čerenkov photons on the expected ring. Additionally, electron tracks are required to have a good match to a cluster in EMCal, and the energy of the cluster must satisfy the requirement $E/p~{}>$ 0.5, where $p$ is the momentum measured by the DCs.

III.2 Combining tracks to electron pairs

All electron tracks in a given event are combined to form pairs. We apply a minimum cut on the transverse mass of the pair, $m_{T}=\sqrt{(mc)^{2}+p_{T}^{2}}>$ 650 MeV/$c$. For the data taken using an ERT trigger, we require that one of the tracks of the pair has a $p_{T}$ of at least 500 (700) MeV/c exceeding the nominal energy threshold 400 (or 600) MeV of the trigger.

These pairs can be subdivided into three groups: (i) Signal pairs that we want to extract. In p$$+$$p collisions these are mostly from the decays of pseudoscalar mesons, vector mesons, heavy flavor mesons. (ii) Combinatorial pairs, which are an undesired background. These result from the combinations of unrelated tracks in any given event, such as combining tracks that originate from two different decays. (iii) Correlated background pairs, which are also undesired, but these pairs do not result from random combinations of tracks. The combinatorial and correlated background pairs should be removed to extract the signal pairs. Most of this is done via a statistical subtraction discussed in detail in Section III.3. However, some of the correlated background can be removed through cuts on the pairs referred to as pair cuts.

There are several sources of correlated pairs which are treated separately. One type of correlated pairs result from detector problems or ambiguities in the pattern recognition. The most important contributor are hadron tracks that are parallel to electron tracks in the RICH. Both tracks share the same ring and are identified as electrons. These pairs can be removed by placing a cut on the distance between the projections of both pairs to the RICH focal plane. Similar cuts to avoid detector overlaps are placed on all detector systems.

Another type of correlated pairs are the ones that originate from the photons that convert to $e^{+}e^{-}$ pairs in the detector material in front of the tracking detectors, e.g. in the beryllium beam pipe (0.3% of a radiation length ($X_{0}$) for the year 2006). The tracks from these pairs get reconstructed with an incorrect momentum, because the tracking algorithm assumes that all tracks originate from the vertex and hence traverse the full magnetic field. This leads to an artificial opening angle of the pairs that is always oriented perpendicular to the axial magnetic field. A cut on the orientation of the opening angle removes these pairs from the sample. See Adare et al. (2015, 2010) for a full description of the pair cuts.

There are also correlated pairs that are from the same p$$+$$p interaction, these are two tracks that share the same ancestry. These pairs can arise if there are two $e^{+}e^{-}$ pairs in the event from the same parent meson, e.g. from a double Dalitz decay of $\pi^{0}$ or $\eta^{0}$ or from a $\gamma\gamma$ decay where both photons convert in the detector material. In this case, the cross-combination of electrons that do not result from the same real or virtual photon are possible. Another source of these correlated pairs are hadrons from jet fragmentation, either within the same jet or in back-to-back jets, that decay into electron pairs. These pairs are part of the statistical subtraction discussed in the next section.

III.3 $e^{+}e^{-}$ pair spectrum

Because the source of any electron or positron is unknown, we combine all the electrons and positrons in a given event into like-sign ($N_{\pm\pm}$) foreground pairs, which is defined as sum of pairs of electrons and pairs of positrons, and $e^{+}e^{-}$ pairs referred to as unlike-sign ($N_{+-}$) foreground pairs. The unlike-sign foreground spectrum $N_{+-}$ measures the sum of signal, combinatorial and correlated background. For this analysis we use the like-sign pairs to determine the backgrounds. The like-sign subtraction method compared to the event-mixing technique has the advantage that it also accounts for the the correlated pair background that exists in the unlike-sign pairs. However, one first needs to correct the like-sign spectrum for the relative acceptance difference between $N_{\pm\pm}$ and $N_{+-}$ pairs.

The relative acceptance correction $\alpha$ which is purely due to the detector geometry is determined via an event mixing technique and is given as the ratio of unlike-sign ($B^{\rm comb}_{+-}$) to like-sign ($B^{\rm comb}_{\pm\pm}$) pair spectrum from the mixed events. The mixed events are generated from MB events and are subject to the same requirement as the ERT data, i.e. each pair is required to have at least one track above 500 (or 700) MeV and this track should have fired the ERT trigger. $\alpha$ is given by the following equation:

[TABLE]

Fig. 1 shows the $p_{T}$-integrated $\alpha$ correction as a function of mass. The acceptance difference is largest around 500 MeV/$c^{2}$. For larger masses, the acceptance difference becomes smaller, and consequently $\alpha$ approaches unity as the mass increases. In the analysis we apply the $\alpha$-correction double differentially in mass and $p_{T}$. The errors on the $\alpha$-correction are propagated to the final spectrum. For systematics, the analysis was checked for $p_{T}$ dependent fixed $\alpha$-values at high masses and results obtained were consistent within 5%. Fig. 2 shows the $p_{T}$ integrated $N_{+-}$ and relative acceptance corrected like-sign mass spectrum ($\alpha~{}\times$ $N_{\pm\pm}$). The acceptance corrected like-sign spectrum is subtracted from unlike-sign $N_{+-}$ spectrum to extract the signal spectrum, $S_{+-}$, as defined by Eq. 2.

[TABLE]

III.4 Efficiency corrections

Eq. 3 gives the invariant yield corresponding to a $S_{+-}$ pair with mass $m$ and transverse momentum $p_{T}$ into the PHENIX aperture:

[TABLE]

Here $\Delta m$ and \Delta\mbox{p_{T}} are the bin width in mass and $p_{T}$, respectively. There are two efficiency corrections that are applied in order to obtain the invariant $e^{+}e^{-}$ yield. These are the inverse of the pair reconstruction efficiency $\varepsilon_{\rm rec}(m,p_{T})$ and pair trigger efficiency $\varepsilon_{\rm ERT}(m,p_{T})$. The $\varepsilon_{\rm rec}(m,p_{T})$ accounts for losses due to track reconstruction, electron identification, pair cuts and detector dead areas. The $\varepsilon_{\rm ERT}(m,p_{T})$ describes the bias introduced by the trigger requirements. Here the BBC efficiency of $\varepsilon_{\rm BBC}=0.545~{}\pm~{}0.06$ is the fraction of inelastic p$$+$$p collisions recorded by the BBC. The BBC trigger bias $\varepsilon_{\rm bias}$ factor takes into account the fact that for the events with tracks in the central arms, the BBC trigger requirement is fulfilled only by $0.79~{}\pm~{}0.02$ of the events.

The pair reconstruction efficiency $\varepsilon_{\rm rec}(m,p_{T})$, as well as pair trigger efficiency $\varepsilon_{\rm ERT}(m,p_{T})$ are determined using a GEANT based simulation of the PHENIX detector. The GEANT simulation is tuned to describe the performance of each detector subsystem and includes all necessary detector characteristics (dead and hot channel maps, gains, noise, etc.).

We simulate $e^{+}e^{-}$ pairs with a constant yield in $m,~{}p_{T},~{}\phi,~{}|y|<1$, and in the mass range 0<m_{\mbox{e^{+}e^{-}}}<16 GeV/$c$2 with $p_{T}$ in the range from 0 to 10 GeV/$c$. These simulated pairs are processed through the PHENIX GEANT framework, and are then weighted with the expected yield from hadron decays for a given pair $[m,p_{T}]$. A detailed description about pair efficiency and trigger efficiency determination can be found in Adare et al. (2015, 2010, 2009). The efficiency corrections are applied double differentially in mass and $p_{T}$, and similar to the previously published PHENIX dielectron analyses, the data are presented in the PHENIX acceptance. The mass spectrum with all corrections is shown in Figure 4, together with the expected sources discussed in the next section.

IV Expected pair sources

The expected yield of $e^{+}e^{-}$ pairs from various sources needs to be simulated in order to interpret the experimental data. This so called cocktail of sources includes the contributions from pseudoscalar and vector meson decays, semileptonic decay of heavy flavor, and $e^{+}e^{-}$ pairs originated via Drell-Yan mechanism.

IV.1 Hadron decays to $e^{+}e^{-}$ pairs

To model the yield of the pseudoscalar mesons $\pi^{0}$, $\eta$, $\eta^{\prime}$, and vector mesons, $\rho$, $\omega$, $\phi$, $J/\psi$, $\psi^{\prime}$, $\Upsilon$, we use a detailed fast Monte Carlo software package called EXODUS developed within the PHENIX framework Adare et al. (2010). EXODUS is a phenomenological event generator that simulates the particle distributions and their decays. EXODUS applies the branching ratios Yao et al. (2006) and decay kinematics according to Adare et al. (2006). External bremsstrahlung in the PHENIX detector material is approximated by placing all the detector material to be traversed by the electron at the radius of the beampipe. The pair mass distribution from Dalitz decays ($\pi^{0},\eta,\eta^{\prime}\rightarrow ee\gamma$) and $\omega\rightarrow ee\pi^{0}$ follows the Kroll-Wada expression Kroll and Wada (1955) multiplied by the electromagnetic form factors measured by the Lepton-G collaboration Landsberg (1985); Dzhelyadin et al. (1981). The vector mesons (\rho,\omega,\phi,J/\psi,\psi^{\prime}\rightarrow\mbox{e^{+}e^{-}}) are assumed to be unpolarized and for their decay the Gounaris/Sakurai expression is used Gounaris and Sakurai (1968). For the Dalitz decays in which the third body is a photon, the angular distribution is sampled according to $1+\lambda cos^{2}\theta_{CS}$ distribution. $\theta_{CS}$ is the polar angle of the electrons in the Collins-Soper frame and $\lambda$ is an angular parameter.

The hadrons are generated with a uniform rapidity density $dN/dy$ within $|\eta|\leq$0.35 and a homogeneous azimuthal distribution in 2$\pi$. Once generated, these hadrons are filtered through the ideal PHENIX acceptance while applying the measured momentum resolution from the data. The key input is the parameterization of the $p_{T}$ dependence of the invariant cross section of neutral pions. To obtain this reference we fit the $p_{T}$ distribution of $\pi^{0}$ and $\pi^{\pm}$ data, as reported by PHENIX Adler et al. (2003b); Adare et al. (2007a, 2011b), to a modified Hagedorn function (Eq. 4):

[TABLE]

The fit parameters and resulting $dN/dy$ values for p$$+$$p collisions are tabulated in Table 1. These values supersede those published in Adare et al. (2010, 2009) as they are based on new and/or more precise data from larger data sets. The pion parameterization determined here deviates by about 3% from the one used in earlier publications.

The $p_{T}$ distribution of other mesons is parameterized by fixing all but the normalization parameter ($A$) from the pion spectrum, and assuming scaling with $m_{T}$, i.e. replacing $p_{T}$ by \sqrt{(\mbox{p_{T}}^{2}-(m_{\pi 0}c)^{2}+(m_{h}c)^{2})}, where $m_{h}$ is the mass of the hadron. The normalization parameter $A$ relates the total $dN/dy$ of a given hadron to the $dN/dy$ of the pions. The successful description of $m_{T}$ scaling is apparent in Fig. 3 which shows measured $p_{T}$ spectra of various mesons as published by PHENIX. In order to extract the meson yield the fits were integrated over all the $p_{T}$. For the $\rho$ meson, we assume $\sigma_{\rho}/\sigma_{\omega}=1.15~{}\pm~{}0.15$ consistent with the values found in the jet fragmentation Yao et al. (2006).

A compilation of the $dN/dy$ values for the various mesons extracted from the fits and the references for the data used are shown in Table 2. These values agree with those from Adare et al. (2010, 2009) within the systematic uncertainties. The differences reflect that more precise data for the pion and other mesons are available today.

IV.2 $e^{+}e^{-}$ pairs from Drell Yan

We used pythia event generator with same settings as mentioned in Adare et al. (2015) to simulate $e^{+}e^{-}$ pairs from the Drell-Yan mechanism. For the normalization we used a cross section of 42 nb as was used in Adare et al. (2015, 2009). We also performed a study where the DY contribution was left as a free parameter. This affected the $b\bar{b}$ cross section by 20% and we assigned that as a systematic uncertainty on the cross section determination.

IV.3 Heavy flavor contribution to $e^{+}e^{-}$ pairs

The $e^{+}e^{-}$ pairs that originate from the semileptonic decays of $c\bar{c}$ and $b\bar{b}$ are collectively referred to as heavy flavor pairs. The heavy flavor yield was simulated using three different event generators. The details of these event generators are described below.

IV.3.1 pythia

pythia Sjostrand et al. is a multi-purpose leading order event generator. It generates heavy quark pairs with massive matrix elements and fragmentation and hadronization is based on the Lund string model. Additional transverse momentum is generated in pythia by virtue of the assumed intrinsic (primordial) transverse momentum $k_{T}$. We used pythia in forced $c\bar{c}$ or $b\bar{b}$ production mode, and CTEQ5L was used as the input parton distribution function. The same settings as published in the d$$+Au paper Adare et al. (2015) are hereby used.

IV.3.2 mc@nlo

The mc@nlo (Monte Carlo at next-to-leading order) formalism is described in detail in Frixione and Webber ; Frixione et al. (a), and is a method for matching next-to-leading order (NLO) QCD calculations to parton shower Monte Carlo (pSMC) simulations. Parton showers will generate terms that are already present in the NLO calculations. To avoid double counting, the mc@nlo scheme removes such terms from the NLO expression. As a result, mc@nlo output contains events with negative weight.

In this work, mc@nlo v4.10 (interfaced with herwigv6.521 Corcella et al. ) was used. The default package was altered to enable charm production by changing the process code from -1705 (H_{1}H_{2}\rightarrow\mbox{b\bar{b}}+X) to -1704 (H_{1}H_{2}\rightarrow\mbox{c\bar{c}}+X) and the heavy quark mass was adjusted to the charm quark mass i.e. 1.29 GeV/$c^{2}$. $H_{1,2}$ represent hadrons (in practice, nucleons or antinucleons). The bottom quark mass was set to 4.1 GeV/$c^{2}$. The default scale choice was used:

[TABLE]

where $m_{T}^{2}=p_{T}^{2}+m^{2}$ and $p_{T}$ is the transverse momentum of the heavy flavor in the underlying Born configuration. $Q$ and $\bar{Q}$ correspond to the heavy quark and antiquark. No other parameters were modified. CTEQ6M Pumplin et al. was used to provide the input parton-distribution function.

IV.3.3 powheg

The powheg (Positive Weight Hardest Emission Generator) formalism is described in detail in Frixione et al. (b). Compared to mc@nlo, powheg generates positive weighted events only, and can be interfaced to any shower MC that is either $p_{T}$-ordered (e.g. pythia), or allows the implementation of a $p_{T}$ veto (e.g. herwig ++), while avoiding any double counting when matching NLO calculations and parton shower Monte Carlo. In this work, powheg v1.0 was interfaced with pythia v8.100 Sjostrand et al. (2008). Parton showering in pythia is $p_{T}$ ordered and merges naturally with powheg. CTEQ6M Pumplin et al. was used to provide the input parton distribution function. Similar to the other two frameworks, the charm and bottom masses were set to 1.29 GeV/$c^{2}$ and 4.1 GeV/$c^{2}$ respectively. The default scale choice was used:

[TABLE]

where $p_{T}$ is the transverse momentum of the heavy flavor in the underlying Born configuration. No other parameters were modified.

The electrons and positrons from all the above mentioned generators are filtered through the PHENIX acceptance Adare et al. (2010) and are folded with the experimental momentum resolution as well as with the energy loss due to bremsstrahlung. The $e^{+}e^{-}$ pair acceptance depends on the production process, which determines the correlation between the electron and positron. More detailed description about the $e^{+}e^{-}$ pair acceptance on (i) the QCD production of the $q\bar{q}$ pair and (ii) the decay kinematics of the two independent semi-leptonic decays has been discussed in Adare et al. (2015). Because the heavy flavor generators discussed above treat the $q\bar{q}$ correlations differently, the number of $e^{+}e^{-}$ pairs that fall into PHENIX acceptance varies from one generator to the other.

V Systematic uncertainties

In this section we summarize the systematic uncertainties on the data and expected sources. Systematic uncertainties on the data are due to limitations in the determination of the relative acceptance correction, the electron identification efficiency, model input used to evaluate the efficiency, and the ERT trigger efficiency. These uncertainties are evaluated by varying all the electron identification cuts and pair cuts, by varying the ERT trigger efficiency within its statistical accuracy and by using different cuts and sub-samples of the data to determine the relative acceptance correction. For all the variations the final result was determined and found stable within the quoted systematic uncertainties.

The main systematic uncertainties on the hadron cocktail comes from the measured uncertainty on the $dN/dy$ of pions. For the heavy flavor part of the cocktail, the assigned uncertainty to $c\bar{c}$ and $b\bar{b}$ normalization comes from this analysis.

Table 3 gives a summary of the systematic errors. The total systematic error on data are added in quadrature and the same is done for the expected sources.

VI Results

VI.1 Heavy-flavor $e^{+}e^{-}$ pairs from p$$+$$p collisions

Figure 4 shows the measured double differential $e^{+}e^{-}$ pair yield in the PHENIX acceptance projected onto the mass axis. The figure also shows the distributions of $e^{+}e^{-}$ pairs from charm, bottom and Drell-Yan obtained using the pythia event generator. The mass region below 1.0 GeV/$c$2 is comprised of resonances and a continuum dominated by three body decays of pseudoscalar and vector mesons. In this mass region, all cocktail contribution, with exception of the heavy flavor meson decay contributions, are absolutely normalized as discussed previously. The contributions of various hadronic decay sources to the cocktail are shown in the inset that highlights the mass spectrum up to 4.5 GeV/$c^{2}$. The mass spectrum above 1.0 GeV/$c^{2}$ is dominated by the $e^{+}e^{-}$ pairs from decays of heavy flavor mesons. The heavy flavor contributions to the dilepton continuum above 1.0 GeV/$c$2 are normalized to the data. Good agreement between data and cocktail over the entire mass range is evident from the ratio of data to the cocktail shown in the lower panel of Fig. 4. We note that below 0.6 MeV/$c^{2}$ there are large systematic uncertainties resulting from the ERT trigger efficiency correction. In this mass region, the results published in Adare et al. (2009) are more accurate due to large sample of MB data available for that analysis. The bulk of the 2006 data used here was taken with the ERT trigger. Our current heavy flavor analysis is based on the mass region above 1.16 GeV/$c^{2}$ and thus not affected by systematic uncertainties around 0.5 GeV/$c^{2}$.

The $e^{+}e^{-}$ pair spectrum from heavy flavor decays is determined using the technique developed for d$$+Au collisions Adare et al. (2015). The expected yield of $e^{+}e^{-}$ pairs from pseudoscalar and vector meson decays as well as Drell-Yan pairs is subtracted from the $e^{+}e^{-}$ pair spectra. The subtraction is done double differentially in mass and $p_{T}$. The resulting mass spectra of $e^{+}e^{-}$ pairs from heavy flavor decays are shown in Fig. 5 for different pair $p_{T}$ ranges. Below 1.0 GeV/$c$2, the yield of $e^{+}e^{-}$ pairs is dominated by hadronic decay contributions and after the subtraction the $e^{+}e^{-}$ pair yield from heavy flavor decays cannot be extracted with sufficient accuracy. Therefore Fig. 5 is truncated just below 1 GeV/$c^{2}$. For those mass regions above 1 GeV/$c^{2}$ where the inclusive $e^{+}e^{-}$ yield is dominated by vector meson decays to $e^{+}e^{-}$ the subtracted yield can not be determined accurately, and hence upper limits are quoted for the subtracted spectra. We use $p_{T}$ bins of width of 500 MeV/$c$ up to $p_{T}$ = 2.5 GeV/$c$. For pair \mbox{p_{T}}>3.0 GeV/$c$, statistical limitations dictate the use of broader $p_{T}$ bins.

The $e^{+}e^{-}$ pair distributions from heavy flavor decays were simulated using three Monte Carlo generators, pythia, mc@nlo, and powheg with parameter settings as discussed above. The results are shown in Fig. 9. The three generators are compared using the normalization

obtained from fitting the data to the respective event generators as described below. As seen in Fig. 9 and already described in detail in Adare et al. (2015), the separation of $e^{+}e^{-}$ pairs from $c\bar{c}$ and $b\bar{b}$ is more evident when one simultaneously analyzes mass and $p_{T}$ of the pairs. The yield from $c\bar{c}$ is dominant for masses below 3 GeV/$c^{2}$ and pair $p_{T}$ less than 2 GeV/$c$, whereas $b\bar{b}$ is dominant across all mass region for higher $p_{T}$. For the pairs with $p_{T}$ $>$ 3.5 GeV/$c$, the largest contribution to the $e^{+}e^{-}$ yield comes from single $b$ decay chains with a semileptonic decay of the parent B meson followed by a semileptonic decay of the daughter D meson.

The generated distributions are fitted simultaneously to all data in $p_{T}$ and mass in the mass regions between 1.15 <\mbox{m_{e^{+}e^{-}}}< 2.4 GeV/$c^{2}$ and 4.1 <\mbox{m_{e^{+}e^{-}}}< 8.0 GeV/$c^{2}$. The mass region from 2.4 to 4.15 GeV/$c^{2}$ is excluded to avoid any remnant contributions to the $e^{+}e^{-}$ yield from $J/\psi$ and $\psi^{\prime}$ decays after the subtraction. Such remnant yield could result from an imperfect description of the line shapes, in particular of the low mass tail due to bremsstrahlung. For each MC generator there are two independent parameters that are fitted, which are the $c\bar{c}$ and $b\bar{b}$ cross sections in $4\pi$. Figs. 9, 9, and 9 show the comparison of fitted distributions to the data for pythia, mc@nlo, and powheg, respectively. The $\chi^{2}/NDF$ values are 1.2, 1.5, and 1.4 for pythia, mc@nlo, and powheg, respectively, with an $NDF$ equal to 65. Here, only statistical errors are used in the fit. Because the $c\bar{c}$ simulated pairs have smaller statistics at high masses for $p_{T}$ $>$5 GeV/$c$, we include the errors from simulations into the fitting routine. Any improvement from additional statistics is expected to be minimal unless significant computing resources are allocated.

The fitted cross sections are tabulated in Table 4. For the $c\bar{c}$ cross section we find 356, 708, and 267 $\mu$b for pythia, mc@nlo and powheg respectively. For each the statistical uncertainty is about 8%, while the systematic uncertainty due to the data is approximately 25%. The values cover a range of $\sim\pm 220~{}\mu$b around the average value, indicating large model dependencies that are further discussed in the following. The $c\bar{c}$ cross section values are consistent with earlier measurements from single electron spectra that gave $\sigma_{c\bar{c}}=567\pm 57(stat)\pm 244(syst)~{}\mu$b Adare et al. (2006) and from $e^{+}e^{-}$ pairs that resulted in $\sigma_{c\bar{c}}=544{\pm}39(stat){\pm}244(syst){\pm}200^{\rm model}~{}\mu$b Adare et al. (2009). For the $b\bar{b}$ cross section we find values of 4.81, 3.85, and 2.91, again for pythia, mc@nlo, and powheg, respectively. The statistical uncertainties are 15–22% and the systematic uncertainties are 21%. The observed model dependence is about $\sim 0.85~{}\mu$b around the average, which is significantly smaller than for $c\bar{c}$ cross section.

Despite the differences between the MC generators, each one achieves an adequate description of the data within the uncertainties. This may be more easily seen in the projections of the $e^{+}e^{-}$ yield from heavy-flavor decays onto the mass and $p_{T}$ axes as shown in Fig. 10.

As a consistency check and to see if more discrimination power between the models can be achieved in terms of different projections of the data, we also looked at the $\Delta\phi$ distribution for $e^{+}e^{-}$ pairs. Because the analysis was done in 2 dimensions, mass and $p_{T}$, some extra steps were necessary. We first generated $\Delta\phi$ distributions for foreground and mixed unlike-sign and like-sign pairs for the mass region between 1.15 <\mbox{m_{e^{+}e^{-}}}< 2.4 GeV/$c^{2}$ and 4.1 <\mbox{m_{e^{+}e^{-}}}< 8.0 GeV/$c^{2}$. The relative-acceptance corrected like-sign foreground $\Delta\phi$ distribution is subtracted from the unlike-sign pairs, which results in the $\Delta\phi$ distribution for heavy flavor pairs. These $\Delta\phi$ distributions were then efficiency corrected.

The data are compared to $\Delta\phi$ distributions from simulated $e^{+}e^{-}$ pairs from $c\bar{c}$, $b\bar{b}$, and Drell-Yan. For each generator, the $c\bar{c}$ and $b\bar{b}$ contributions were normalized using the cross section values from Table 4. For the $b\bar{b}$ contribution the like-sign pairs were subtracted to match the procedure used in the data. The $\Delta\phi$ distributions from pythia, mc@nlo, and powheg are shown for different pair $p_{T}$ ranges and compared to the data in Fig. 11. Note that these distributions are for $e^{+}e^{-}$ pairs within the PHENIX acceptance. Again, all three generators describe the data reasonably well. The conclusions are consistent with those drawn from the comparison in $p_{T}$ and mass. At lower $p_{T}$ the yield is dominated by $c\bar{c}$ production. The yield peaks at large opening angle $\Delta\phi$, which is characteristic for back-to-back production. At the same $p_{T}$, the pairs from $b\bar{b}$ production show no pronounced back-to-back structure. This is consistent with the $e^{+}e^{-}$ pair opening angle being less correlated with the $b\bar{b}$ opening angle due to the decay kinematics of the much heavier B mesons. For larger $p_{T}$ $b\bar{b}$ production dominates, and the $e^{+}e^{-}$ pair opening angle $\Delta\phi$ distribution peaks for opening angles smaller than 90 degrees. This is due to the fact that these pairs result from the decay products of a single B-meson rather than from the $b\bar{b}$ pair.

Only moderate differences are observed between the generators. While there are differences in the shape of the $\Delta\phi$ distributions for $c\bar{c}$ and $b\bar{b}$, the main structure seen in Fig. 11 is given by the two arm detector acceptance. We find that the statistical significance of our data is insufficient to add more discriminating power between the generators by looking at the $\Delta\phi$ projections.

While the data are well described by all three generators within the PHENIX central arm acceptance and over the range they were fitted to the data, the obtained cross section values, tabulated in Table 4, indicate that there are large systematic differences when extrapolated beyond the range where the models were fitted to the data.

The $c\bar{c}$ cross sections found using pythia and powheg differ by about 30%, while for mc@nlo a much larger $c\bar{c}$ cross section is determined. This may be due to the fact that our powheg simulation uses the pythia fragmentation scheme. Such differences can have important consequences if the generators are used to estimate yields from $c\bar{c}$ outside the fit range, even within the PHENIX acceptance. This was first pointed out in Adare et al. (2016) and is apparent when one looks at the $e^{+}e^{-}$ pair mass distributions below 1 GeV/$c^{2}$, depicted in Fig. 12. For pythia and powheg there is very little difference going from mass larger than 1 GeV/$c^{2}$ to zero mass, while for mc@nlo the $e^{+}e^{-}$ pair yield is much larger. This is an important contribution to the larger $c\bar{c}$ cross section determined with mc@nlo.

To get a better quantitative understanding, we divided the extrapolation into following three steps: the first step is the extrapolation from the fitted $e^{+}e^{-}$ pairs in the PHENIX acceptance to $e^{+}e^{-}$ pairs at all masses, then to the $q\bar{q}$ rapidity density, and finally to $4\pi$. These factors are tabulated in Table 5 and Table 6 for $c\bar{c}$ and $b\bar{b}$, respectively. For $c\bar{c}$ production the number of $e^{+}e^{-}$ pairs in the fit range is similar for pythia, mc@nlo, and powheg. This is expected, because the normalization is essentially fitted in the range from 1 to 2 GeV/$c^{2}$ where $c\bar{c}$ dominates. The extrapolation to zero mass is different only for mc@nlo, and is responsible for about 50% of the larger cross section for mc@nlo. Going from $e^{+}e^{-}$ pairs in the PHENIX acceptance to the rapidity density dN_{\mbox{c\bar{c}}}/dy at $y=0$ has the largest variations between models. The final step from $c\bar{c}$ rapidity density to $4\pi$ has little model dependence indicating that the underlying rapidity distribution for $c\bar{c}$ is similar in all the generators.

The situation is however different for $b\bar{b}$ production. From Table 6 it is evident that every step of the extrapolation from the fit range to $4\pi$ is very similar for all three generators. Again this is expected because the $e^{+}e^{-}$ pair distributions from $b\bar{b}$ production are dominated by decay kinematics Adare et al. (2015). However, the number of $e^{+}e^{-}$ pairs in the fit range is different, which leads to different $b\bar{b}$ cross section values. The extracted $b\bar{b}$ cross section value using pythia is larger as compared to the one derived from mc@nlo, with the latter being larger than powheg. From Figs. 9 and 10, one can see that the shape of the $e^{+}e^{-}$ pair distributions from $b\bar{b}$ production are very similar among the three generators. However, this is not the case for $e^{+}e^{-}$ pairs from $c\bar{c}$ production, in particular for powheg, the $e^{+}e^{-}$ pair momentum distribution is much harder as compared to other generators. Because the $c\bar{c}$ contribution is essentially fixed in the mass region between 1.0 to 2.0 GeV/$c^{2}$ at low pair $p_{T}$, a harder distribution can only be accommodated in the overall fit by reducing $b\bar{b}$ production, which we expect to account for all the seen variation between the three generators. Additional differences in the rapidity and momentum distribution also contribute to the very model dependent extrapolations of the $c\bar{c}$ cross section in $4\pi$.

VI.2 Comparison of p$$+$$p and d$$+Au results

The results of the analysis of p$$+$$p data presented here can be directly compared to the already published d$$+Au results Adare et al. (2015). Because we are now including powheg and are using a newer version of mc@nlo for the p$$+$$p analysis, we refitted the data published in Adare et al. (2015) with the generator versions used for p$$+$$p. We scaled down the d$$+Au data by the average number of binary nucleon-nucleon collisions of $N_{\rm coll}$ ($=7.6\pm 0.4$). Therefore the resulting normalization constants represent the equivalent nucleon-nucleon cross section, and can be directly compared to the p$$+$$p results.

Table 7 summarizes the $c\bar{c}$ and $b\bar{b}$ nucleon-nucleon equivalent cross sections extracted from the d$$+Au data. We note that the numbers quoted here for the mc@nlo simulation are 17% and 12% smaller for $c\bar{c}$ and $b\bar{b}$, respectively, compared to the numbers quoted in Adare et al. (2015). This is potentially due to using a newer mc@nlo version, which needed to be modified to generate charm, or a previous inaccuracy in how the negative weights should be used to avoid double counting in the herwig fragmentation Adare et al. (2015). In either case the difference is small enough to change the conclusions neither here nor in the original paper Adare et al. (2015).

The comparison of the numbers in Table 4 and Table 7 is shown graphically in Fig. 13. We see the same model dependence for d$$+Au as was seen for p$$+$$p. For a given model, the obtained $c\bar{c}$ cross sections are consistent within the given uncertainties in p$$+$$p and d$$+Au. We also looked at the ratio (or nuclear modification) of cross sections of $c\bar{c}$ and $b\bar{b}$ in d$$+Au and p$$+$$p and this is plotted in Fig. 16. This ratio is similar for all the event generators and no deviation from unity is observed.

Fig. 15 and Fig. 15 show a direct comparison of the measured mass and $p_{T}$ spectra of $e^{+}e^{-}$ pairs from heavy flavor decays between p$$+$$p and d$$+Au systems. The top panels show an overlay of mass and $p_{T}$ spectra in p$$+$$p and MB d$$+Au collisions, where we scaled the p$$+$$p yield by $N_{\rm coll}$ ($=7.6\pm 0.4$), corresponding to MB d$$+Au collisions. Within the statistical precision of the data, the mass and $p_{T}$ spectra in p$$+$$p and d$$+Au agree with each other. The bottom panel in these figures show the ratio of d$$+Au to p$$+$$p data. Given the uncertainties, the ratios are consistent with 1. While the $e^{+}e^{-}$ pair data shows no evidence for any nuclear modification to the $c\bar{c}$ and $b\bar{b}$ production, due to the large statistical and systematic uncertainty, they would not be sensitive to effects smaller than 30%. For example, the observed modification of single electron spectra seen in d$$+Au collisions Adare et al. (2012a) could result in a change of 30% in the $e^{+}e^{-}$ pair mass and $p_{T}$ distributions, but that might not be seen here due to the large uncertainties.

VII Summary and conclusions

We present $e^{+}e^{-}$ pair measurements from heavy flavor decays in p$$+$$p collisions at $\sqrt{s}=200$ GeV. The data are shown multi-differential as a function of pair mass,  $p_{T}$, and $\Delta\phi$. By comparing the $e^{+}e^{-}$ pair data to pQCD calculations, the $c\bar{c}$ and $b\bar{b}$ production cross sections can be constrained. Three different pQCD based Monte-Carlo models are used: pythia, mc@nlo, and powheg. We find that the $c\bar{c}$ production cross section ranges from 267 to 708 $\mu$b with a statistical (systematic) uncertainty of about 8% (25%). The $b\bar{b}$ production cross section ranges from 2.9 to 4.8 $\mu$b with a statistical (systematic) uncertainty of 15–22% (21%).

The $e^{+}e^{-}$ pair distributions obtained from pythia, mc@nlo, and powheg within the PHENIX acceptance, once normalized to data, were found to be consistent in mass, $p_{T}$ and $\Delta\phi$. In case of $c\bar{c}$, the extrapolation beyond the measured range shows substantial model dependence. This is evident by more than 400 $\mu$b difference between the obtained $c\bar{c}$ cross sections, which is more than 100% compared to the average value.

We find a smaller variation for $b\bar{b}$, which is less than 50% of the average $b\bar{b}$ cross section value. This variation is entirely due to the model dependence of $c\bar{c}$ production. The extrapolation of $b\bar{b}$ from the measured range shows little model dependence, because in our acceptance the decay kinematics dominate the $e^{+}e^{-}$ pair distributions from $b\bar{b}$.

We compare our p$$+$$p results directly to $e^{+}e^{-}$ pair measurements from MB d$$+Au collisions. The $c\bar{c}$ and $b\bar{b}$ cross sections are determined in the same way for both the systems. Although there is significant model dependence in extracting the cross sections, within a given model, there is no difference between the cross sections determined from p$$+$$p and the equivalent nucleon-nucleon cross section obtained from d$$+Au. Furthermore, we compare directly the measured $e^{+}e^{-}$ pair mass and $p_{T}$ distributions from p$$+$$p and d$$+Au. After scaling with the number of binary collisions, we observe no evidence for nuclear modifications of heavy flavor production in the d$$+Au system within our experimental uncertainties.

ACKNOWLEDGMENTS

We thank the staff of the Collider-Accelerator and Physics Departments at Brookhaven National Laboratory and the staff of the other PHENIX participating institutions for their vital contributions. We acknowledge support from the Office of Nuclear Physics in the Office of Science of the Department of Energy, the National Science Foundation, a sponsored research grant from Renaissance Technologies LLC, Abilene Christian University Research Council, Research Foundation of SUNY, and Dean of the College of Arts and Sciences, Vanderbilt University (U.S.A), Ministry of Education, Culture, Sports, Science, and Technology and the Japan Society for the Promotion of Science (Japan), Conselho Nacional de Desenvolvimento Científico e Tecnológico and Fundação de Amparo à Pesquisa do Estado de São Paulo (Brazil), Natural Science Foundation of China (People’s Republic of China), Croatian Science Foundation and Ministry of Science and Education (Croatia), Ministry of Education, Youth and Sports (Czech Republic), Centre National de la Recherche Scientifique, Commissariat à l’Énergie Atomique, and Institut National de Physique Nucléaire et de Physique des Particules (France), Bundesministerium für Bildung und Forschung, Deutscher Akademischer Austausch Dienst, and Alexander von Humboldt Stiftung (Germany), National Science Fund, OTKA, EFOP, and the Ch. Simonyi Fund (Hungary), Department of Atomic Energy and Department of Science and Technology (India), Israel Science Foundation (Israel), Basic Science Research Program through NRF of the Ministry of Education (Korea), Physics Department, Lahore University of Management Sciences (Pakistan), Ministry of Education and Science, Russian Academy of Sciences, Federal Agency of Atomic Energy (Russia), VR and Wallenberg Foundation (Sweden), the U.S. Civilian Research and Development Foundation for the Independent States of the Former Soviet Union, the Hungarian American Enterprise Scholarship Fund, and the US-Israel Binational Science Foundation.