Prompt photon production and photon-hadron jet correlations with POWHEG

We present a calculation of direct photon production at next-to-leading order of QCD and a matching of this calculation with parton showers using POWHEG BOX. Based on simulations with POWHEG+PYTHIA, we perform a detailed phenomenological analysis of PHENIX data on prompt photon production and photon-hadron jet correlations in pp collisions at RHIC, considerably improving the description of these data with respect to previous calculations, and we suggest additional interesting analyses.


Introduction
In heavy-ion collisions at RHIC [1,2] or the LHC [3][4][5], a new state of strongly interacting matter can be created, the so-called quark-gluon plasma (QGP). Important probes of this hot medium of deconfined quarks and gluons are thermal photons, which interact only electromagnetically and can thus leave the medium without strong modifications of their thermal spectrum. The exponential falloff of the photon spectrum at low transverse momenta (p T ) can then be related to the temperature of the QGP [6][7][8][9]. The extraction of the true temperature of the QGP at the time of its creation is complicated by the fact that the medium is rapidly expanding and cooling, that photons are radiated at all stages of the collision including the phases before thermalisation and after recombination of the quarks and gluons into charged and neutral hadrons, that neutral pions decay preferably into pairs of photons, and that photons are also produced promptly in partonic collisions, either directly or through fragmentation processes.
In Refs. [6,9], a first observation of a low-p T photon signal after subtraction of the meson decay background in Pb-Pb collisions at the LHC has been reported by the ALICE collaboration, and an inverse slope parameter was extracted from the p T -spectrum for 0-20% central collisions. Using next-to-leading order (NLO) QCD calculations, the relative contributions to prompt-photon production from different initial and final states and the theoretical uncertainties coming from independent variations of the renormalisation and factorisation scales, the nuclear parton densities and the fragmentation functions have been analysed. Based on different fits to the unsubtracted and prompt-photon subtracted ALICE data for 0-40% central collisions, we found effective temperatures of T = 304 ± 58 MeV and 309±64 MeV at p T ∈ [0.8; 2.2] GeV and p T ∈ [1.5; 3.5] GeV as well as a power-law (p −4 T ) behaviour for p T > 4 GeV as predicted by QCD hard scattering [7,8]. In lowerenergy Au-Au collisions at RHIC, a smaller effective temperature of T = 221 ± 27 MeV had previously been measured [10].
Precise calculations of prompt photon production in hadronic collisions are not only imperative for a reliable extraction of the thermal photon spectrum, but also for measurements of photon-hadron and photon-jet correlations, which represent a second important probe of the hot medium due to the p T -imbalance and azimuthal asymmetries induced by jet quenching [11,12]. In both cases, additional parton emission can significantly modify the spectra and thus the physical conclusions. So far, the theoretical description of prompt photon production has relied on NLO calculations [13] and in particular JETPHOX [14] with at most one additional parton for direct and fragmentation processes. The latter dominate at low p T and require a convolution with insufficiently determined non-perturbative fragmentation functions [15], unless one applies photon isolation criteria [16] or departs from real to slightly virtual photons [17][18][19][20].
An alternative approach consists in the combination of NLO calculations with parton showers (PS). There, the photon fragmentation function can be modelled by an interleaved QCD+QED parton shower, leading e.g. to a correct description of the photon fragmentation function at LEP [21,22]. In addition, the exclusive events produced in this Monte Carlo approach allow for a detailed comparison to experiment, in particular realistic isolation cuts and even a combination with detector simulations. Furthermore, the parton shower resums the leading logarithmic contributions from multiple additional parton emissions, thus providing considerably more realistic kinematic distributions. This applies in particular to photon-hadron and photon-jet correlations with their unrealistic δ-functions or divergences predicted at leading order (LO) and NLO in the regions of vanishing photonhadron transverse momentum p γh T → 0 and back-to-back azimuthal angle ∆φ → π. As we will see explicitly, the collinear region ∆φ → 0 is of course closely related to photon fragmentation processes.
The combination of NLO QCD corrections with PS requires a careful treatment of soft/collinear regions in order to avoid double counting. Methods like MC@NLO [23] and POWHEG [24] are now well established for QCD processes. The treatment of photons is more intricate, as it also requires a QED parton shower. It has previously been achieved for di-photon production as a Higgs boson background [25] and as an electroweak correction to single-W production [26,27]. In this paper, we report on a re-calculation and validation of direct photon production at NLO in Sec. 2, a matching of this calculation with PS using POWHEG BOX [28] in Sec. 3, and in Sec. 4 on a successful phenomenological reanalysis of PHENIX data on photon and photon-hadron production in pp collisions at √ s = 200 GeV, which form the baseline for the corresponding analyses in heavy-ion collisions. Our conclusions are presented in Sec. 5.
2 Direct photon production at NLO Direct photon production proceeds at LO through the partonic processes qq → γg and qg → γq. We computed the spin-and colour-averaged cross sections for these processes analytically using FormCalc 8.4 [29] and checked the results against MadGraph 5 [30] and the literature [31][32][33][34]. The same procedure was applied to the real emission processes with an additional parton in the final state. The virtual one-loop corrections were computed with FORM [35] in D = 4−2ε dimensions and reduced from tensor to scalar integrals using the Passarino-Veltman procedure [36]. After renormalisation of the ultraviolet divergences in the MS scheme [37], infrared divergences remained which could be shown to cancel against those from the integrated Catani-Seymour dipoles [38] as computed e.g. with AutoDipole 1.2.3 [39]. For the finite remainders of the one-loop contributions, agreement with those produced by MadLoop [30] was then obtained.
As is well known [13], a consistent calculation of prompt photon production up to NLO requires also the inclusion of fragmentation processes at least in LO in order to cancel the divergences from collinear quarks and photons. For photons with finite transverse momenta, they appear only in the final state and are canceled by the corresponding dipole terms arising from collinear factorisation. The LO direct and purely partonic fragmentation processes scale formally with O(αα s ) and O(α 2 s ), respectively. However, the latter must still be convoluted with fragmentation functions (FFs), which scale as so that both contributions eventually have the same scaling behaviour. For the numerical evaluation of our NLO direct and LO fragmentation results, we computed the scalar loop integrals using LoopTools 2.13 [29]. The partonic cross sections were then convoluted with parton density functions (PDFs) and FFs and compared to JETPHOX [14]. As expected, good numerical agreement was found. Important advantages of NLO over LO calculations are a more reliable (typically larger) normalisation of the total cross section, its stabilisation with respect to variations of the unphysical renormalisation and factorisation scales, and improved descriptions of kinematic distributions. Disadvantages with respect to Monte Carlo generators are the restriction to at most one additional parton and the absence of hadronisation effects.

Prompt photon production with POWHEG
The POWHEG BOX [28] provides a framework to smoothly incorporate NLO corrections in general-purpose event generators such as PYTHIA [40] or HERWIG [22], as long as they allow for a p T -ordered parton shower or have the ability to veto radiation with a p T higher than that of the first radiation. Usually, it suffices to provide the Born amplitudes along with their colour-and spin-correlated counterparts, the Born phase space, a decomposition of the amplitudes in the colour flow basis, the finite part of the virtual corrections, and the real correction amplitudes as FORTRAN routines. These inputs are then linked against the core POWHEG BOX code. In our case, however, it was necessary to modify small parts of the POWHEG BOX code itself, primarily to accommodate a consistent treatment of photon radiation off quarks and furthermore to introduce an artificial enhancement of photon radiation. Since the POWHEG BOX in version 2 is already able to handle photon radiation off leptons as reported in Ref. [26], only minor modifications of the code were necessary. In the following, we use the notation established in Ref. [24].
In POWHEG [24], the real processes are subdivided into parts corresponding to collinear and soft regions, such that for a specific flavour structure (i.e. partonic subprocess) the real process R can be written as a sum R = αr R αr (3.1) with an index α r denoting the different singular regions. The individual contributions R αr are chosen such that, for some region of the real correction phase space where the configuration of two particles produces a collinear or soft singularity, only one R α ′ r becomes singular, while all other R αr with α r ≠ α ′ r vanish. Hence, every region α r with n+1 particles corresponds to an underlying Born flavour structure with n particles, denoted by f b and obtained by replacing the two particles in the singular configuration by the particle from which they emerged in a splitting process. This defines a mapping α r → f b . In addition, there exists for all α r a decomposition of the n+1-particle phase space, denoted by Φ n+1 , into n-particle kinematics Φ αr n and radiation variables Φ rad , giving a mapping of real kinematics to Born kinematics. The decomposition in Eq. (3.1) and the relation between Φ αr n , Φ rad and Φ n+1 are central to the Frixione-Kunszt-Signer (FKS) subtraction method [41], which is employed by the POWHEG BOX to regularise infrared singularities, but also to the formulation of the POWHEG Sudakov form factor and the POWHEG cross section. The latter is defined as Here, B f b is the NLO inclusive n-particle cross section, where all radiative corrections from regions α r with underlying Born structure f b -denoted by the set {α r f b } -have been integrated out, i.e.
with the Born amplitude B, the subtracted virtual corrections V sv , the real counterterms C and the collinear remnants G ⊕ and G ⊖ . Eq. (3.2) is the cross section for events with at most one radiation off the Born flavour structure with p T > p min T , or in parton shower terminology the cross section with the first radiation evolved down to p min T . Compared to the usual parton shower prescription for the first radiation, the POWHEG modifications are given by the replacement B → B and the substitution of the parton shower splitting kernel (usually the Altarelli-Parisi splitting kernel) with the ratio of the real and Born amplitudes R αr (Φ n+1 ) Φ αr n =Φn B f b (Φ n ) for each radiation region α r ∈ {α r f b }. Consequently, the POWHEG Sudakov form factor is given by In our treatment of photon production at NLO we follow the approach of Ref. [25] and include all real amplitudes R with one photon in the final state. Our goal is then to simulate the fragmentation contribution arising from the QED radiation off partons (as described for a tree-level merging approach in Ref. [21]). Thus, we need to include the Born flavour structures f QCD b for LO dijet production along with the photon production flavour structures f QED b as underlying Born processes. In this case the sum over the flavour structures f QCD b in Eq. (3.2) includes the photon fragmentation contribution through the combination of QED non-branching and branching probabilities, i.e. the terms in the brackets.

Born amplitudes and phase space
As described above, we implement both the QED and QCD Born amplitudes, expanded in D = 4 − 2ε dimensions up to O(ε 2 ) for reasons specified in Sec. 3.4. Aside from the rather straight-forward implementation of the Born amplitudes themselves, an assignment of colours to the external legs has to be given according to the large-N c limit. To facilitate this, one identifies the squared diagrams for a specific amplitude which have a planar colour flow and assigns a colour to each of their lines. If then several diagrams have conflicting colour assignments, one is chosen according to the relative weight of its planar diagram with respect to the sum of all planar diagrams. However, in the case of QED Born amplitudes this is not necessary, since there is only one qqg-vertex. It is therefore trivial to assign a colour to each quark. For the more complicated case of the QCD Born amplitudes, we use the colour assignment of the POWHEG BOX implementation of jet pair production [42].
Since we treat all particles as massless, the Born phase space is the same for our process and dijet production, enabling us to also use the phase space routines from that code. As the Born amplitudes diverge for vanishing squared partonic momentum transfert, a phase space cut on the final state transverse momentum k T with respect to the beam axis is mandatory. To regularise the divergence, we make use of the Born suppression factor [42] This feature is activated in the parameter file by setting the option bornsuppfact to a value for k T,supp and has the effect of replacing the functions One thus avoids the divergence of the Born amplitudes. This change of the cross section is corrected for by weighting the generated events with the inverse of Eq-(3.5). By default, we use k T,supp = 100 GeV and i = 3 in Eq. (3.5), but we have checked (with k T,supp = 10 GeV and i = 4) that within statistics and at sufficiently large p T our results are independent of these choices. Alternatively, a simple cut on k T can be activated with the parameter bornktmin. After showering, this cut would, however, lead to a loss of events in the low-p T region, which is precisely the region in heavy-ion collisions, where one wants to extract the thermal photon spectrum. We have instead tested the continuous approximation of the Heaviside function which does, however, not improve significantly on the statistics. Independently of the method used to regulate the LO divergence, it must be noted that the impact of a choice of k T,supp (or a k T cut) on the photon spectrum is not obvious, since k T is not generally the photon transverse momentum, but the momentum of some particle prior to any generation of radiation. In particular, in the original POWHEG scheme, only the splitting q → qγ, but not q → γq, was generated, which could yield high-p T photon events with low statistics, but large weights after showering. Adding the line doublefsr 1 to the file powheg.input allows to generate also the splitting q → γq, which avoids this statistical problem [43].

Colour-correlated Born amplitudes
In the soft limit, a real correction is given by the Born amplitude times an eikonal factor, where the latter depends on the colour correlations between coloured legs of the Born amplitude. To enable the POWHEG BOX to compute these limits, the colour-correlated Born amplitudes, defined by are given as an input. In Eq. (3.7), the Born matrix elements are denoted by M {c k } with the colours of the external particles specified by the index set {c k } and the averaging factors subsumed in N . The colour indices for the legs i and j are contracted with colour charge operators T a . The colour correlations for processes with less than four coloured partons -as in our QED Born processes -are readily reduced to sums of Casimir operators by making use of the colour conservation relations with the sum running over all coloured legs i and an implicit contraction of the colour index c i . Physically, these equations simply state that an infinitesimal global transformation of all initial and final state colours simultaneously has no effect. From this it follows that with the Casimir constants C f j . Eq. (3.9) gives three equations for three coloured legs, which can be solved for all the colour correlations, since B ij is symmetric under the exchange of i and j. Thus for qq → γg with momenta p 1 , p 2 , p 3 and p 4 , respectively, we have The results for all other partonic QED subprocesses are related to this result via crossing symmetry. The colour correlations for the QCD Born flavour structures are not needed, since we do not implement O(α 3 s ) corrections.

Spin-correlated Born amplitudes
Similarly to the colour correlations for the construction of soft limits, spin-correlated amplitudes are required to construct the collinear limits for gluons splitting into two partons. The spin-correlated Born amplitude is defined via where now the index s j represents the spin of the gluon on leg j and ε µ s j is a polarisation vector for leg j. This prescription amounts to replacing the polarisation vectors of leg j in the matrix element M and its Hermitian conjugate M † with the physical polarisation sums, e.g.
with η some light-like vector spanning Minkowski space together with p j and the two physical polarisations ε s j . As is usual, the polarisation vectors are chosen to be space-like, orthogonal, and normalised to unity, leading to Taking again the partonic process qq → γg and performing the described substitutions of polarisation vectors for the gluon on leg 4 in FORM, we find In practice, we choose η µ = g µµ p µ 4 (no sum over µ). Again, the other QED processes are obtained using crossing symmetry. Here, again, the correlations for QCD Born amplitudes are not needed, since only QED radiation is allowed off them. We also do not need the spin correlations for photons, since they are not allowed to split, i.e. we do not include real corrections with an internal photon line.

Virtual corrections
The virtual corrections are provided to the POWHEG BOX in the form of the finite part V fin of the MS-renormalised virtual amplitude V . The relation between V and V fin is given in the conventional dimensional regularisation (CDR) scheme via While the coefficients a and c ij are independent of ε, the amplitudes B, B ij are given in D = 4 − 2ε dimensions and thus depend on ε. Since our results computed in Sec. 2 have the form with all ε-dependence explicit, it was necessary to compute from our finite result V (0) the CDR-finite result V fin . By expanding the Born amplitudes in ε, and comparing Eqs. (3.17) and (3.18), it is easy to see that the relation holds, which we implement in the code. The coefficients a and c ij can be extracted from Ref. [41], giving where i, j run over all coloured legs and µ R is the renormalisation scale. The constants γ f i are given by where T F = 1 2 as usual.
We implement the virtual O(α s ) corrections to the photon production processes using Eq. (3.20), but do not include the virtual O(α s ) and O(α) corrections to the dijet processes, as they lead to higher-order corrections for prompt photon production. We comment on the cancellation of divergences and the inclusion of finite remnants of the subtraction method in the next section.

Real corrections and the cancellation of divergences
In our implementation of the real corrections with a final state photon, we have to make sure that the QED collinear and soft divergences are correctly identified. In order for the POWHEG BOX to find the QED singularities in addition to the QCD singularities, it has to treat the photon as just another massless parton, which is achieved by setting the variable flst lightpart to 3. Then all singular regions α r of the real amplitudes are automatically identified and subtracted. As a check of consistency of the real corrections and the colour-and spin-correlated Born amplitudes from Secs. 3.2 and 3.3, the POWHEG BOX tests numerically if, for each singular limit, the ratio of R αr and the corresponding limiting expression tends to one. Our implementation passes these tests.
The subtraction of singularities is handled internally by the POWHEG BOX with the FKS subtraction method. The soft and collinear counterterms (denoted by C in Eq. (3.3)) are automatically assembled from the soft and collinear limiting expressions and subtracted from the real corrections. On the other hand, the counterterms, integrated in D = 4 − 2ε dimensions over the momentum of the emitted parton, are added to the virtual corrections V in Eq. (3.17), leading to a cancellation of poles apart from collinear initial-state singularities, which are automatically absorbed into parton distribution functions in the MS factorisation scheme. The leftovers from the absorption of the divergences of the initialstate collinear counterterms into the PDFs are the collinear remnants G ⊕ and G ⊖ in Eq. (3.3). Since these are automatically computed by the POWHEG BOX for all Born flavour structures, we implemented if-clauses, which disable the (in our case inconsistent) computation of these terms for the f QCD b amplitudes. However, even apart from the collinear initial-state singularities, the addition of integrated counterterms and virtual corrections does not equal V fin . Rather there are some finite terms, that are computed automatically by the POWHEG BOX giving the soft-virtual contribution that enters Eq. (3.3). The definitions of I and Q are provided in Ref. [24]. Here, we again made sure that these terms are only computed for the f QED b amplitudes, i.e. those amplitudes for which we actually implemented V fin .
Since we do not implement the virtual QED corrections to the f QCD b flavour structures, as would be required in a fully consistent treatment, there is some ambiguity in the choice of the finite remnants of the soft and collinear QED singularities. In a fixed-order calculation, this freedom of choice amounts to the choice of a factorisation scheme for the photon fragmentation function. In analogy to the MS factorisation scheme, we cancel only the poles of the QED singularities, which (by inspection of Eq. The contribution of Eq. (3.25) for QED corrections can be activated by setting the flag flg with em = .TRUE. in the code, whose function we have extended to massless quarks. With the I-terms for QED already included in the POWHEG BOX version 2, we just added the Q-term for photon radiation off massless quarks with (cf. Eq. (2.100) in Ref. [24] with δ o = 2), where the sum is over all charged legs, Q f i is the charge of particle f i , and In summary, our soft-virtual term for the QED singularities is where σ i = 0, if f i is a (initial-) final-state (anti-) particle, and σ i = 1, if it is a (final-) initial-state (anti-) particle.

Enhanced QED radiation
As it stands, the implementation of single-photon production in the POWHEG BOX framework described in the preceding sections leads to a very low photon production rate. For example, a pp-run at the PHENIX energy √ s = 200 GeV contains a photon in only about 2% of the events, while the remaining events are made up of QCD Born configurations. Reasons for this behaviour are a relative suppression of photonic vs. purely partonic processes from the ratio of electromagnetic and strong coupling constants (α α s ), QCD colour factors larger than unity vs. squared fractional quark charges smaller than unity, and a larger multiplicity of contributing processes in QCD. To boost the contribution of photons, we implement a procedure described in Ref. [21]: We multiply the integrand in the exponent of the POWHEG Sudakov form factor in Eq. (3.4), denoted by f (Φ rad ) in the following, for QED radiation with a constant c > 1, thus decreasing the no-branching probability, and compensate for it by reweighting the event.
Usually, the transverse momentum of a radiation is generated by first solving the equation for a uniformly distributed random number r ∈ [0, 1], where ∆ (U ) (p T ) is a lower bound for Eq. (3.4), obtained by replacing f (Φ rad ) by an upper bounding function U (Φ rad ) > f (Φ rad ). Afterwards, the generated p T is accepted with a probability given by the ratio f U or vetoed with a probability 1 − f U . In the latter case Eq. (3.33) is solved again, but this time restricting p T to values below the vetoed value. The whole procedure is then reiterated, until a p T is accepted or when a low p T -cutoff of the order of Λ QCD is reached. Replacing Eq. (3.33) with has in the POWHEG BOX the same effect as multiplying both f and U by c. As described in App. B of Ref. [21], we then compensate for this by weighting the event with where the product runs over all QED radiation vetoes and f i , U i are the values of f and U at the respective vetoed p T . A similar procedure, described in Ref. [44], has also been tested and produces consistent results.

Parton shower with PYTHIA 8
The events generated by the POWHEG BOX have to be passed to a parton shower generator to produce complete events. Every parton shower generator that is p T -ordered or has facilities to veto radiation with a p T higher than the scale of the hardest event is viable. PYTHIA is employing a p T -ordered parton shower, and we use PYTHIA 8 [40] in this work. However, since the definitions of the transverse momentum of a radiation differ for the POWHEG BOX and PYTHIA, it is suggested to use the class PowhegHooks to account for the differences. The preferred mode of usage is to have PYTHIA evolve the shower starting from the kinematical limit rather than the scale passed by the POWHEG BOX, translate the p T of a generated radiation from the PYTHIA definition to the POWHEG BOX definition, and then veto radiation harder than the hard POWHEG scale. In addition, there is in our case the question of how to handle the scales for QED and QCD radiation. The default would be to make no distinction and veto the evolution of QED and QCD radiation above the POWHEG scale independently of the type of event. Instead, we follow the approach presented in Ref. [25], which suggests employing two different hard scales, one for the QED and one for the QCD shower. A discussion of this approach and another approach, that includes a competition between QED and QCD radiation already at the level of the hard process (i.e. an implementation of R ∼ O(α 3 s )), can be found in Ref. [27]. In the nomenclature of that reference, we use the NC-scheme.
To allow for the distinctive treatment of QED and QCD radiation, we modify the POWHEG BOX code to pass the scale of the underlying Born event, in addition to the scale of the first radiation, to PYTHIA. According to our implementation, the Born scale is the hard scale for QED radiation, if the underlying Born event includes a photon, or the hard scale for QCD radiation, if the underlying Born event is a pure QCD event. The scale of the emission that is usually passed by the POWHEG BOX corresponds in those two cases to the hard scale for QCD and QED radiation, respectively. By a modification of the PowhegHooks class, we ensure that the QED and QCD showers are vetoed accordingly.

Comparison with PHENIX data
In this section, we present the numerical results of our implementation of the NLO corrections to prompt photon production into the POWHEG BOX matched to the PYTHIA 8 parton shower. 1 We take the opportunity to reanalyse data taken by the PHENIX collaboration in pp collisions at √ s = 200 GeV on nearly-real virtual and real inclusive photons at low (p T ∈ [1; 6] GeV and p T ∈ [5; 16] GeV, respectively) [48] and higher transverse momenta (p T ∈ [5; 25] GeV) [49]. These "vacuum" data are not only important as a baseline for heavy-ion collisions, but the low-p T region is also interesting for studies of perturbative QCD itself, in particular of the fragmentation contribution, photon-jet and photon-hadron correlations, and soft radiation effects. Besides comparing our new results to data, we also validate them against a pure NLO calculation with JETPHOX, pointing out important similarities and differences, i.e. the common NLO normalisation, but the presence of only one additional parton and of fragmentation contributions parameterised by non-perturbative fragmentation functions in the latter. In our comparison with the stand-alone PYTHIA 8 Monte Carlo generator, we emphasise the common multiple parton emission, which leads to a better description of kinematic distributions, and the different (NLO vs. LO) normalisation.
In all our theoretical calculations, we employ CT14NLO parton densities in the proton [45]. Proton PDFs including photons and photon radiation through LO QED evolution are also available [46]. However, we did not use these for several reasons: they pertain to only initial and initial-state soft and collinear photons, would require the implementation of the full (i.e. also virtual) QED corrections, affect the production of prompt photons with finite p T only beyond LO in QED, and are thus here of little numerical importance. In our comparison with JETPHOX, we employ the BFG II photon fragmentation functions [47], which we have previously shown to be favoured by the PHENIX low-p T data [15]. By default, the renormalisation and factorisation scales are set to the p T of the underlying Born process, which can be a parton for LO QCD processes. As mentioned in the previous section, we have used the Born suppression factor in Eq. (3.5) with k T,supp = 100 GeV and i = 3 by setting bornsuppfact 100. No generation cut on the Born k T is applied. QED radiation is enhanced from ∼ 2 % to ∼ 25 % through Eq. (3.34) and Eq. (3.35) with c = 50, and independence of the precise value of c has been checked by varying it to c = 100.
The PHENIX experiment has detected real photons with two electromagnetic calorimeter (EMCal) arms covering the pseudorapidity range η γ < 0.35. Conversion photons were identified with additional electron hits in the ring imaging Cerenkov detector. During the 2006 RHIC runs, integrated luminosities of L = 4.0 and 8.0 pb −1 were collected in the low [48] and higher p T ranges [49]. After subtracting the decay background, the prompt photon with the hardest p T was selected in each event. In the higher-p T analysis, also the effects of an isolation cut on the photons was analysed. There, the hadronic energy fraction in a cone of radius R = (∆η) 2 + (∆φ) 2 ≤ 0.5 was restricted to be less than 10% of the photon energy. In a previous publication, the PHENIX collaboration have published data on photon-hadron jet correlations [11]. Using integrated luminosities of 3.0 and 10.7 pb −1 collected during 2005 and 2006 RHIC runs, they identified charged-hadron jets with a tracking system composed of a drift chamber and pixel pad chambers. Like the PHENIX collaboration, we take into account in our calculations all charged hadrons associated with a photon and within the PHENIX acceptance. We have checked that modeling the charged-hadron jet in a modern way with the anti-k T cluster algorithm and a distance parameter R = 0.4 and rejecting events with jets that do not contain a charged hadron produces similar results.

Transverse-momentum spectrum of prompt photons
We begin the discussion of our numerical results with the transverse-momentum spectrum of prompt photons in pp collisions with √ s = 200 GeV at RHIC, shown in Fig. 1. As one can see, the LO+PS prediction with PYTHIA alone (green) starts to describe the data only at p γ T > 14 GeV, but falls short of the experimental cross section below this value and by up to two orders of magnitude in the lowest bin. The NLO prediction with JETPHOX, the fragmentation function BFG II and central scale choices (µ R = µ p = µ γ = p T , blue) describes the measured inclusive photon spectrum reasonably well as expected. 2 This is even more true for the NLO+PS prediction with POWHEG+PYTHIA (red), which coincides with the NLO prediction within statistical errors over a wide range of p T > 10 GeV, thus validating the calculation in a region that should be insensitive to multiple soft/collinear parton emissions. At lower p T , where these emissions become relevant, the NLO+PS prediction exhibits a characteristic increase (lower panel of Fig. 1) and follows the data very well, while the pure NLO prediction over-/undershoots the lowest and second-lowest data points.
Strictly speaking, a description of fragmentation photons with partons showers captures only the pointlike fragmentation component. As mentioned above, this leads indeed to a correct description of the photon fragmentation function at LEP [21,22]. It is, however, well-known that the photon fragmentation function also has a non-perturbative component that is traditionally described with the Vector Meson Dominance (VMD) model  [48,49]. The ratio of NLO+PS over pure NLO is shown in the lower panel.
where f V are the vector-meson decay constants [13]. We estimate the possible contributions from these long-range processes in Fig. 2. As expected, the VMD contributions are dominated by the lightest vector meson (ρ) and fall off more rapidly in p T than the pointlike photon contribution. Their contributions can only be substantial at very low p T . Like the vector mesons themselves, their fluctuations into photons might there indeed be sensitive to strong medium effects, contrary to the naive expectation that photons interact only electromagnetically. A consistent combination of the pointlike and VMD contributions is beyond the scope of this paper and is left for future work. It has, however, previously been argued that in specific factorisation schemes such as the DIS γ scheme the VMD component is completely negligible [50]. Since in this scheme an additional soft/collinear term ln[x 2 (1 − x)] is resummed to all orders in the fragmentation function, similarly to the parton shower, we do not expect a large VMD contribution in our NLO+PS approach.

Fraction of isolated photons
The fraction of isolated photons in the higher-p T PHENIX data set is shown in Fig. 3 (black). The NLO calculation (blue) overestimates this fraction considerably in the lowand intermediate-p T region. This remains true for all standard scale choices as already discussed in the experimental publication [49]. There, the difference was temptatively attributed to the underlying event activity or quark fragmentation contributions. A com- parison of the data with a LO+PS calculation from PYTHIA did, however, not show any drop in the low-p T region, either. As one can see in Fig. 3, our new NLO+PS calculation with POWHEG+PYTHIA (red) describes the PHENIX data better, although the statistical error bars are still relatively large. Note also that in this calculation the scale uncertainty cancels completely, as the ratio is constructed from the same event sample in the numerator and denominator and does not include any contributions from the scale-dependent fragmentation function that affects the NLO calculation differently in the numerator and denominator.

Transverse-momentum spectrum of the associated charged hadron
Important observables for the quark-gluon plasma are hadron energy loss and jet quenching, i.e. the energy loss of a hadronic jet induced by the hot medium. The p T -distribution of the charged hadrons produced in association with the photon is therefore shown in Fig.  4. Unfortunately, it has not been measured in the cited PHENIX publications [48,49]. We therefore use the PHENIX detector acceptance for charged hadrons η h < 0.35 quoted in the preceding analysis of photon-hadron jet correlations (see below) [11]. In LO, the leading jet transverse momentum equals that of the photon, and indeed the LO+PS (green) and NLO+PS (red) distributions follow those of the photon in Fig. 1 except for a shift of the p T -scale of roughly 20% due to the missing neutral-hadron contribution. In contrast to the photon spectrum, however, the NLO K-factor remains almost constant at high p T due . Ratio of isolated photons in NLO (blue) and NLO+PS (red) and compared to PHENIX data (black) [49].
to the fact that this region is very much QCD-like in the sense that the observed charged hadron can also be balanced by other partons rather than only a photon.

Transverse-momentum balance of photons and charged hadrons
Going one step further, one can also measure photon-hadron jet correlations in pp or heavyion collisions. They have the advantage that the photon (or more generally electroweak boson) balancing the charged-hadron jet is not strongly influenced by the hot medium and can thus serve as a gauge for the initial jet transverse momentum. An exact balance holds, however, only at LO of perturbative QCD, so that deviations can not be uniquely attributed to medium effects, but must also take into account higher-order QCD corrections.
In Fig. 5 we therefore show the distribution in the combined photon-hadron transverse momentum p γh T . To avoid the non-perturbative region when all transverse momenta vanish, we here apply individual cuts on p γ T > 1 GeV and p h T > 1 GeV. 3 At LO (not shown), one then obtains a δ-distribution as the combined p γh T → 0, while at NLO (also not shown) the differential cross section is sensitive to the incomplete cancellation of infrared divergences from the emission of one additional soft parton. This region is resummed to all orders by the parton shower in PYTHIA already at LO (green), so that this prediction exhibits a finite and physical turnover. The LO normalisation is, however, still incorrect and modified by up to two orders of magnitude in POWHEG+PYTHIA(red). Only at this order can the  p T imbalance be reliably predicted and this observable subsequently applied to heavy-ion collisions in order to extract genuine medium effects.
An important aspect of our NLO calculation is the appearance of new partonic processes with no LO correspondence. In particular, the process qq → γqq with a recoiling quark jet first enters at this order in addition to the process qq → γgg, which has just an additional gluon compared to the LO process qq → γg with a recoiling gluon. The recoiling quarks and gluons are then of course expected to behave differently in the medium due to their different colour charges and infrared behaviour. Another additional process first entering at NLO is gg → γqq with no corresponding process at LO. This gluon-initiated process is expected to be more sensitive to the initial conditions of the heavy-ion collision, in particular due to shadowing and saturation effects.

Azimuthal correlation of photons and charged hadrons
Besides jet quenching and the transverse-momentum imbalance, azimuthal correlations of photons and hadron jets represent another important probe of the quark-gluon plasma. They have therefore indeed been measured by the PHENIX collaboration [11], using unequal transverse-momentum cuts on the photon (p γ T ∈ [5; 7] GeV) and charged-hadron jet (p h T ∈ [3; 5] GeV) as required [51]. In particular the away-side correlation (∆φ → π) has been found to be suppressed in 0-20% central Au-Au collisions for both decay and direct photons, which we can interpret as an indication of decorrelation due to rescattering on the medium.
In Fig. 6 we reproduce the PHENIX data (black) [11], which have been obtained with a statistical subtraction of the decay photon background. In particular, the cross section at the minimum of the correlation function has been subtracted, assuming a Zero-Yield at Minimum (ZYAM), and in addition it has been normalised to the total number of trigger photons. Both the LO and NLO calculations predict unphysical results for this observable (a δ-distribution at π and vanishing results for jets with ∆φ below π 2 [52]) and are therefore not shown. In contrast, our NLO+PS calculation with POWHEG+PYTHIA (red) exhibits a physical behavior with finite results in both limits ∆φ → {0; π} and describes the PHENIX data quite well, while the LO+PS prediction with PYTHIA alone (green) cannot correctly describe the region ∆φ → 0 (see below) and also has a wrong normalisation (not shown).
To better understand the individual contributions to the near-and away-side correlation function, we show in Fig. 7 the isolated (red), non-isolated (blue) and total (black) photon contributions to the azimuthal-angle correlation of the photon-hadron pair. As expected, the non-isolated photons originating from fragmentation processes dominate the cross section in the near-side region, i.e. they are mostly collinear to the parton fragmenting into the observed hadron (∆φ ≃ 0), and these NLO processes are of course missing in the LO+PS calculation. At the level of 35%, they also contribute in the away-side region (∆φ ≃ π), where they originate from fragmentation from the (unobserved) second parton. This region is, however, dominated by back-to-back photons and partons as in LO, so that PYTHIA captures the essence of the physics in this region. Note that our predictions in Fig. 7 are neither subtracted to ZYAM nor normalised for possible future comparisons of absolute cross sections. Remember also that in heavy-ion collisions isolation cuts can not be used due to the large hadronic underlying event in the low-p T region, where one wants to extract the thermal-photon spectrum, so that a correct description of fragmentation processes is very important.

Pseudorapidity correlation of photons and charged hadrons
For completeness, we show in Fig. 8 the correlation in radidity of the photon-hadron pair with the same individual pseudorapidity cuts as in the previous subsection and p γ T , p h T > 1 GeV. While the pseudorapidity range accessible to the RHIC detectors has so far been quite limited, future upgrades at RHIC or measurements at the LHC bear important potential for studies of the low-x region and therefore of the initial conditions of the formation of the QGP e.g. from a colour-glass condensate.

Conclusion
In this paper we have presented a calculation of direct photon production at NLO QCD, a matching of this calculation with parton showers using POWHEG BOX, and a detailed phe- nomenological analysis of PHENIX data on prompt photon production and photon-hadron jet correlations in pp collisions at RHIC energies. Our work was motivated by the facts that the inclusion of photons in parton showers is highly non-trivial, that prompt photons are important probes of the QGP in heavy-ion collisions, for which pp form the indispensable baseline, and that they give furthermore access to photon-hadron jet correlations and studies of jet quenching. To this end, we have described in detail the analytical and numerical validations of our calculations at different stages and then solutions to various encountered difficulties such as the suppression of divergent Born contributions, the symmetrisation of parton splittings involving photons, and the enhancement of QED radiation.
The application of our NLO+PS calculations to PHENIX data taken at RHIC has led to important improvements compared to both LO+PS and pure NLO: in the description of the low-p T inclusive photon spectrum, of the fraction of isolated photons contained in this data sample, and of the azimuthal correlations of photons and charged hadrons. In addition, we have made predictions for the p T spectra of the associated charged hadron and for the p T balance of the photon-hadron pair as well as their pseudorapidity correlation and have decomposed the azimuthal correlation function into fragmentation and non-fragmentation components.
As the next step, our calculation can easily be applied to pA and AA collisions for studies of cold nuclear effects at RHIC or the LHC. Subsequently, one can tackle the implementation of rescattering on the medium along the lines of Ref. [53], but including NLO corrections. A consistent combination of pointlike and VMD contributions to photon fragmentation in Monte Carlo generators is also envisaged in future work.