Consistent simulation of non-resonant diphoton production in hadron collisions including associated jet production up to two jets

# Consistent simulation of non-resonant diphoton production in hadron collisions including associated jet production up to two jets

## Abstract

An event generator for diphoton () production in hadron collisions that includes associated jet production up to two jets has been developed using a subtraction method based on the LLL subtraction. The parton shower (PS) simulation to restore the subtracted divergent components involves both QED and QCD radiation, and QED radiation at very small are simulated by referring to a fragmentation function (FF). The PS/FF simulation has the ability to enforce the radiation of a given number of energetic photons. The generated events can be fed to PYTHIA to obtain particle (hadron)-level event information, which enables us to perform realistic simulations of photon isolation and hadron-jet reconstruction. The simulated events, in which the loop-mediated process is involved, reasonably reproduce the diphoton kinematics measured at the LHC. Using the developed simulation, we found that the 2-jet processes significantly contribute to diphoton production. A large 2-jet contribution can be considered as a common feature in electroweak-boson production in hadron collisions although the reason is yet to be understood. Discussion concerning the treatment of the underlying events in photon isolation is necessary for future higher precision measurements.

## 1 Introduction

The last missing piece of the Standard Model, i.e., the Higgs boson, was discovered by the ATLAS and CMS experiments at the CERN LHC proton-proton () collider [1, 2]. Diphoton () production played a crucial role in this discovery. The experiments have now entered a new phase that aims to investigate the detailed properties of the discovered Higgs boson. Diphoton production is also an important process in this phase [3, 4].

The Higgs boson was discovered as a resonance in the invariant mass distributions of its decay products. In the diphoton mode, the resonance is observed on a large non-resonant diphoton background from other processes. Hence, a precise understanding of this background is indispensable for Higgs-boson studies. However, our present knowledge is not good enough to predict the background properties with sufficient precision. The background contribution is evaluated using data-driven methods in experimental measurements without relying on theoretical predictions or simulations. Because such evaluations are always based on certain assumptions, better theoretical understanding is desired for reliable estimations.

Although photons are produced via well-known quantum electrodynamic (QED) interactions, it is not straightforward to evaluate the properties of diphoton production in hadron collisions. Experiments at Fermilab Tevatron, a proton-antiproton () collider at a center-of-mass energy () of 1.96 TeV, observed that next-to-leading order (NLO) predictions do not precisely reproduce their measurement results [5, 6, 7]. Most remarkably, the CDF experiment found that the NLO predictions significantly underestimate the production of acoplanar two photons [6, 7]. Similar discrepancies have also been observed in experiments at the LHC [8, 9].

Recently, it was demonstrated [10] that the deficit of the NLO predictions can be recovered by improving the approximation to the next-to-next-to-leading order (NNLO) [11]. A similar improvement was previously found in a simulation with a leading-order (LO) event generator, SHERPA [12], which includes associated jet production up to two jets [13]. Here, a jet represents a parton (light quark or gluon) in the final state. These facts imply that the improvements were accomplished predominantly by newly included two-jet production processes, even in the NNLO prediction. More recently, the ATLAS and CMS experiments published their results concerning several kinematical distributions of diphoton production with significantly improved statistics [14, 15]. In these reports, they confirmed the improvements by the NNLO prediction and SHERPA simulation. However, none of the theoretical predictions nor simulations that they examined could reproduce all measurement results within the measurement errors, which are 10% to 20% in most measurement points. Similar improvements and discrepancies were also observed by the CDF experiment [16].

We reported previously [17] that we have successfully developed a consistent Monte Carlo (MC) event generator for non-resonant diphoton production processes, including those associated with one additional jet. Using this event generator, we confirmed that the contributions of the processes induced by constituent quark-gluon () collisions overwhelm the lowest-order quark-antiquark () collision processes in proton-proton collision experiments, owing to the high gluon density inside protons. This observation further implies that gluon-gluon () collisions may significantly contribute to diphoton production. The collisions necessarily produce at least two jets in the final state in association with diphoton production at the tree level, as shown in Figs. 1(a) and 1(b). In addition, quark-quark scattering processes such as those illustrated in Fig. 1(c) may also have a large contribution because any combination of quarks can contribute. The inclusion of these qualitatively new processes may be the reason for the significant improvements accomplished by the NNLO calculation and SHERPA simulation.

The question is now the reason for the remaining discrepancies between the measurements and predictions. In the NNLO predictions, the dominant reason may be the contributions from multiple quantum chromodynamic (QCD) radiations to be resummed. In contrast, soft radiation effects are included by parton shower (PS) simulations in SHERPA. Hence, along with missing higher-orders and loop corrections, implementation of the radiation-merging method and/or modeling are also suspected to be the reason for the inaccuracy of the SHERPA simulation. Therefore, it must be worth cross-checking the performance of the LO simulations with a different merging method and independent implementation.

In this article, we improve the MC event generator that we developed previously [17] to consistently include diphoton + 2-jet production processes by extending our matching method for jet and photon radiations based on the LLL subtraction [18, 19, 20, 21, 22]. The main difference with the SHERPA simulation is the coverage of PS simulations and treatment of non-divergent components. In the SHERPA simulation, which is based on the idea of the CKKW method [23], radiation simulations based on matrix-element (ME) calculations are applied down to a matching scale () that is significantly smaller than the typical energy scale () of the considered hard interaction for which the MEs are evaluated. In order to interpolate between the two energy scales, the Sudakov suppression deduced from multiple QCD radiation is applied to the simulated events by reinterpreting the radiations in the context of a PS. The application of PS simulations is limited in small regions, . Here, the definition of is model-dependent, but at least at the limit it is equal to for spacelike initial-state radiation (ISR) and the squared mass of the radiator for timelike final-state radiation (FSR).

In contrast, in our simulation, PS simulations that resum divergent logarithmic components to all orders are applied up to . ME-based simulations cover radiation contributions in larger regions and those from finite non-logarithmic contributions at . Thus, no interpolation is required and the non-logarithmic contributions in MEs are preserved in all regions. As a drawback, the boundaries between the ME-based and PS simulations for logarithmic components appear in visible regions. Therefore, it is crucial to adopt an appropriate kinematics model for PS branches. We demonstrated that the kinematics model that follows the simplest massless approximation (-prefixed kinematics) provides good matching between the two simulations [17, 19, 20].

An essentially identical matching method is also applied to QED photon radiation from final-state quarks. In order to accomplish the photon-radiation matching, the PS simulation needs to have the ability to radiate photons. The final-state PS that we previously developed [17] simultaneously supports QCD and QED radiation. In addition, photon radiations with very small that PS simulations do not cover are simulated by referring to a fragmentation function (FF) [24]. In principle, we can produce any number of photons in this PS/FF simulation. However, most of the produced photons are too soft to be detected. The simulation has to be repeated until a hard photon that satisfies the detection requirements is produced, in order to simulate the photon production in the ordinary way of PS simulations. Because the probability of producing a detectable hard photon is very small, such a simulation is very inefficient and is practically unrealistic if we require two hard photons.

In order to solve this problem, the previously developed PS/FF simulation has the ability to enforce the radiation of one energetic photon [17]. This simulation is sufficient for diphoton + 1-jet production processes because at least one of the two photons is necessarily well separated from the other final-state particles. However, in diphoton + 2-jet production processes such as those illustrated in Fig. 1, two photons may be collinearly produced with respect to outgoing quarks. In the present study, we extend the simulation so that we can enforce multiple photon radiation. The developed PS/FF simulation is capable of enforcing the radiation of two photons from one outgoing quark (Figs. 1(b) and (d)), as well as those from the combination of two single radiations (Figs. 1(a) and (c)). In principle, we can enforce the radiation of any number of photons in all possible combinations.

In addition to the tree-level processes, we also include the loop-mediated process in the present study. The QCD PS simulations are fully applied to all generated events, and the simulated parton-level events are passed to PYTHIA [25] in order to obtain the event information at the particle (hadron) level. Kinematical distributions are derived from the obtained particle-level events. Photon-isolation cuts and hadron-jet reconstructions can, therefore, be simulated realistically. After checking the internal consistency of the simulation, the simulation results are compared with recent measurement results at the LHC in order to test the capability of the simulation.

The remainder of this article is organized as follows. Section 2 describes the strategy used in our simulation. Together with the overall strategy, the enforced photon radiation in the PS simulation is described. The internal consistency of the simulation is tested in Section 3. After describing the special treatment of the loop-mediated process in Section 4, predictions from the simulation are compared with recent measurement results in Section 5. Remarkable features of the diphoton production processes are discussed in Section 6 using simulation results. Finally, discussions are concluded in Section 7.

## 2 Simulation strategy

### 2.1 Matching method

To simulate events in hadron collisions while allowing associated jet production, it is necessary to avoid double counting the jet production in ME calculations and PS simulations. In our matching method, the matching is accomplished by numerically subtracting divergent leading-logarithmic (LL) terms from the squared MEs of radiative events [18]. Because only the LL components are subtracted, the finite non-logarithmic components in MEs are preserved in the entire phase space. The PS simulations applied to corresponding non-radiative events restore the subtracted terms without double counting or gaps. The divergent LL contributions are regularized in PS simulations by resumming them to all orders. Thus, the cross sections to be considered in the simulation are all finite. Since PS simulations are limited by a certain energy scale (), the subtraction must also be limited by the same energy scale. Hence, we call this method the limited leading-log (LLL) subtraction [19].

We take the above to be equal to the typical energy scale () of the hard interaction for which the MEs are evaluated. Thus, we do not need to account for the Sudakov suppression due to the multiple QCD radiation in the ME evaluation. Because we take be large, the boundaries between the ME-based simulation and PS simulation for the LL contributions emerge in visible regions. Hence, special care is required in both simulations in order to achieve good matching between them. On the ME side, in the subtraction procedure, a radiative event is separated into a hypothetical PS branch and a non-radiative event in which the assumed radiation is removed by exactly reversing the PS procedure. The energy scale that limits the subtraction is determined by referring to the PS scale that is defined for the reconstructed non-radiative event. The subtraction is applied only when the of the separated PS branch is smaller than the obtained PS scale. Here, the subtraction functions are tailored so that they should be exactly identical to the leading term of the PS simulation [21, 17].

On the PS side, the modeling of PS branch kinematics is important for the matching. We adopt the -prefix kinematics in which the value of each branch is fixed to the value that is determined from the simplest massless approximation; i.e., for ISR [20] and for FSR [21]. In order to realize this kinematics, the identity of is allowed to be violated in ISR [19], and the energy conservation is violated in FSR [22]. The energy conservation is restored by adjusting the initial-state momenta after completing the PS simulations. Accordingly, some corrections are applied to compensate for the adjustment in the initial state [17].

We execute the LLL subtraction and PS simulations in a frame where incoming partons are aligned back to back. In such a setup, we do not need to separately consider the soft-gluon (SG) divergence when the studied process includes only one jet in the final state. However, when two or more jets are included, it is necessary to consider the SG divergences together with the collinear divergences that are subtracted by the LLL subtraction to render the cross sections finite. In our simulation, the SG divergences are subtracted simultaneously with collinear divergences using the combined subtraction method [26]. In order to compensate for this alteration in the subtraction, a correction is applied to the events of corresponding non-radiative processes by referring to the hardest branch in the PS simulation (SG correction) [26].

A similar matching method is also applied to QED photon radiation in photon production processes [17, 27]. The cross section diverges if a photon and a quark are produced collinearly. Although such a photon is likely to be confined in a hadron jet induced by the quark, the photon may be detected as an isolated photon when the quark momentum is small. The LLL subtraction is applied to regularize such a final-state divergence, whereas initial-state divergences and soft-photon divergences are not taken into consideration because we always require the detection of energetic photon(s) at large angles when we simulate photon production processes.

In order to restore the subtracted QED LL components, our final-state PS simulation has the ability to radiate photons in the same manner as QCD parton radiation [17]. PS simulations are necessarily cutoff at a small to avoid divergences. We terminate the PS simulation at . Although QCD effects at smaller can be simulated down to the particle level with the help of general-purpose event generators such as PYTHIA, similar QED simulations are not available in these generators. In our PS simulation, photon radiations at are simulated by referring to a fragmentation function (FF) [24]. In addition, in order to improve the simulation efficiency, the PS/FF simulation is capable of enforcing the photon radiation. Although the original version [17] could radiate only one photon in each event, the simulation has been extended to multiple photon radiations, as described later.

### 2.2 Matched diphoton + 2-jet production

Programs for calculating the MEs of diphoton () + 2-jet production processes were produced using the GRACE system [28, 29] and installed in the GR@PPA event generator [21]. The + 2-jet production () MEs have various divergences to be subtracted. Together with the initial-state divergences, we have to subtract the final-state QCD collinear divergences that emerge when the two jets are produced collinearly. Furthermore, in addition to the collinear divergences, we also have to consider the SG divergences when a gluon is included in the final state. This divergence is subtracted simultaneously with the collinear divergences using the combined subtraction method [26].

The subtracted QCD divergences are restored by combining + 1-jet production () events, to which QCD PS simulations are applied. Here, the SG correction is applied in order to compensate for the alteration in the subtraction from due to the SG divergence. Because also has the initial-state QCD divergence, the LLL subtraction is applied and + 0-jet () events are combined to restore the subtracted components, as in the previous study [17].

The final-state QED divergences are also subtracted to make the cross sections finite. The subtracted QED components are restored by combining + 2-jet () events. The events again have both QCD and QED divergences. The QCD divergences are subtracted by recursively applying the combined subtraction and are restored by combining + 1-jet production () events to which the SG correction is applied. The QED divergences are subtracted using the LLL subtraction and restored by combining QCD 2-jet () events, as in the simulation of direct-photon production [27]. The radiation of one energetic photon is enforced in the PS simulation in the generation of and events, and the radiation of two energetic photons is enforced in the event generation. The QED divergences are also subtracted from the events and restored by the forced photon radiation from events.

When one of the four particles in an event is judged to be a soft radiation () for which subtraction should be considered, the reconstructed non-radiative three-body event may again contain a soft radiation. The subtraction is not applied to such doubly-soft events. Instead, they are rejected. Namely, the non-divergent components in doubly-soft contributions are ignored [26]. In addition to the LL components, we further need to subtract next-to-leading logarithmic (NLL) components in order to evaluate these non-divergent components. The evaluation of NLL components is beyond the scope of the present study. Nevertheless, dominant LL components in the doubly-soft configuration are restored by the PS simulations applied to , , and events. Moreover, some higher-order components are restored by PS simulations applied to and events that contain a non-divergent soft radiation.

In total, the events of six processes, , , , , , and , are generated according to their tree-level MEs. The combined subtraction is applied to and , and the LLL subtraction is applied to . One-photon radiation is enforced in the generation of and , and two-photon radiation is enforced in . In addition, the SG correction is applied to and . All these events are combined to compose a consistent simulation of diphoton production allowing associated jet production up to two jets. We take the gluon and light quarks up to the bottom quark as the initial-state partons and final-state jets in the event generation.

In order to carry out the event generation, we need to explicitly define the energy scale of hard interactions. Although any definition is in principle acceptable in our matching method, we adopt the following definition as in the previous study [27];

 μ=max{QT,i}. (1)

Here, quantity is defined for each final-state particle in terms of its mass and , such that

 Q2T,i=m2i+p2T,i, (2)

and the largest value is taken as the energy scale of the event. This definition reduces to the value of the particle that has the largest , as all particles are nearly massless in the present study.

In the event generation, we are required to specify several energy scales, such as the renormalization scale () to determine the coupling parameters in ME calculations and the factorization scale () to resum the initial-state QCD radiations in parton-distribution functions (PDFs). Furthermore, boundary for the LL components described above may be independently defined for ISR and FSR ( and ). In the present study, all these energy scales are taken to be equal to , which is defined in Eq. (1).

The parton-level events to which the PS simulations are fully applied in GR@PPA are passed to PYTHIA 6.425 [25] to simulate the hadronization and particle decays. Small- QCD effects at that are not covered by the GR@PPA PS are also simulated in PYTHIA. The PYTHIA simulation is applied with its default setting, except for the settings of PARP(67) = 1.0 and PARP(71) = 1.0, as in all our previous studies [20, 21]. The event selection is applied to the obtained particle-level events to derive kinematical distributions, although the diphoton kinematics are not significantly affected by this PYTHIA simulation.

The kinematical conditions in the hard-interaction generation and the cuts to select the events to be passed to PYTHIA are sufficiently relaxed in order not to bias the final particle-level selection. Among them, the generation conditions are markedly relaxed because the PS simulations significantly alter the event kinematics. In order to ensure efficient event generation under such conditions, we adopt the LabCut framework [17] for the generation of all processes. In this framework, the PS simulations and event selection are applied before passing the differential cross-section values to the MC integration and event generation utility BASES/SPRING [30, 31] used in GR@PPA. Hence, the distribution of random numbers for the hard-interaction generation is automatically optimized accounting for the event selection conditions after the PS application. In addition, event weights from the SG correction and forced photon radiation can be involved in the differential cross-section values using this framework.

Usually, kinematical cuts are merely applied to the photons in diphoton studies. Accordingly, we impose practically no constraint on the jets in the hard-interaction events. Even when some requirements are imposed on hadron jets in a study, it is dangerous to apply corresponding cuts to the jets (partons) in hard-interaction events because additional hadron jets may be produced by PS simulations. Exceptions are small cutoffs in the values and separation to the other particles, i.e., GeV/ and , where is measured with respect to the beam axis and is defined by the separations in the azimuthal angle () and pseudorapidity () as . These cutoffs are applied to ensure numerical stability of the subtraction. This very loose setting is allowed because all divergences are precisely subtracted from the MEs of radiative processes.

### 2.3 Forced multiple photon radiation

In the present study, the forced photon radiation in the final-state PS simulation developed in a previous study [17] is extended to multiple photon radiation. The procedure is essentially the same as that for single radiation. The difference is mainly in the determination of the values of photon radiation. In the PS simulation, we first determine which quark should radiate what number of photons. This is randomly determined by assuming an equal probability for all combinations. Here, we do not consider any photon radiation from gluons. The number of possible combinations of this assignment can be given by the homogeneous product. Provided there are quarks in the final state, the number of ways to assign photons to them is given as

 Ncomb=H(n,m)=(n+m−1)!(n−1)!m!. (3)

Once the number of photons to radiate is determined for each quark, we can determine the values of the photon radiations in decreasing order using the QED Sudakov form factor  [17]. The values are determined by solving

 SQED(Q2,Q2pre)=ξ, (4)

where is the of the previous radiation and is equal to for the first radiation. In ordinary PS simulations, is a uniform random number between 0 and 1, whereas in the forced radiation, is restricted in the range

 SQED(Q20,Q2pre)<ξ<1, (5)

in order to ensure that we always obtain a solution in the range . Here, is the lower cutoff of the PS simulation. This constraint leads to an event weight of

 w=1−SQED(Q20,Q2pre). (6)

Because the Sudakov form factor represents the no-radiation probability, the event weight in Eq. (6) corresponds to the probability of obtaining a radiation in the relevant range.

The photon radiations at are covered by the FF-based simulation. The FF-based radiation is considered only for the last radiation from each quark. Therefore, the event weight in Eq. (6) is modified to

 wlast={1−SQED(Q20,Q2pre)}+{1−e−PFF(Q20)} (7)

for the last radiation, where is the integration of the radiation probability density given by the FF. We determine whether the photon is radiated by the PS or FF-based simulation according to the ratio of the radiation probabilities described by the first and second terms in Eq. (7). If PS is selected, the value is determined by solving the equation in Eq. (4) with the constraint in Eq. (5).

The base event-weight is determined by the product of the number of possible combinations given in Eq. (3) and the weights in Eq. (6) or Eq. (7) evaluated for all photon radiations. The subsequent procedure is the same as that for the single forced-radiation [17]. Photon radiations are inserted during the QCD evolution when the becomes smaller than the predetermined photon-radiation . A multiple insertion happens if multiple photon-radiations are assigned to a quark. An FF radiation is added after completing the QCD evolution if it is selected for the last radiation. The splitting parameter , which determines the momenta of the branch products, is determined during this procedure. The restrictions on the allowed range to ensure sufficient energy for the radiated photons determine the additional event weights to be multiplied. The event weight is set to zero if the quark does not have sufficient energy.

The differential cross section values to be passed to BASES/SPRING are multiplied with the thus-determined event weight using the LabCut framework of GR@PPA, for appropriate cross-section integration and the generation of unweighted events. The number of forced photon-radiations () and minimum energy () for the photon radiation can be specified by users in the PS simulation developed in this study. Although the cut is applied in the PS simulation that is executed in the center-of-mass frame of the hard interaction, it does not bias the generated events if its value is set to the minimum requirement for the photons.

Diphoton events can be generated by applying this PS simulation to events by setting . Figure 2 shows the distribution of the azimuthal angle separation () between the radiated two photons. The events were simulated under the condition imposed in the comparison with the ATLAS measurement, which is described later. We can clearly observe two different components in this simulation; the two single-radiation contributions such as those corresponding to the collinear components of the processes in Figs. 1(a) and (c) concentrate in a coplanar region, , while the double-radiation contributions such as those corresponding to the processes in Figs. 1(b) and (d) emerge in an acoplanar region, . Events at are suppressed by the required separation between the two photons.

## 3 Matching test

The internal consistency of the radiation matching can be tested by investigating the stability of kinematical distributions of the simulated events against the variation of the energy scale. Here, we take the value given by the definition in Eq. (1) as the standard value and compare the distributions for the settings of , , and . Although usually the larger scale is set to in scale-dependence studies, the choice of very large value is troublesome in our simulation because large PS activities may boost small- photons into detectable regions. Such small- photons cannot be properly simulated because the initial-state QED divergences and soft-photon divergences are not subtracted. This problem is severe in large multiplicity processes such as + 2-jet production. In any case, is large enough to test the stability around . The kinematical distributions are examined in terms of the normalized differential cross section , where represents the tested quantity, because the overall normalization is not relevant to the matching properties.

The event simulation was carried out for collisions at = 7 TeV using the built-in CTEQ6L1 PDF [32]. We deactivated the underlying-event (UE) simulation in PYTHIA by setting MSTP(81) = 0 because it does not affect the properties of energetic objects to be examined in the following procedure. From the simulated particle-level events, we selected those events in which the two generated photons satisfied the conditions

 pT(γ)>100 GeV/c,  |η(γ)|<2.4,  and  ΔR(γγ)>0.4, (8)

where is measured with respect to the beam direction and represents the separation between the two photons. In addition, both photons were required to be isolated from other activities with the condition of

 EconeT(ΔR<0.4)<7 GeV (9)

in order to simulate a realistic detection condition. Cone , , is defined by the sum of the transverse energies of all stable particles other than muons and neutrinos that are contained inside a cone of around the photon.

We examine the properties of + 2-jet events. Here, jet does not refer to a parton but an observable hadron jet. The hadron jets were reconstructed using FastJet 3.0.3 [33], with the application of the anti- algorithm with . All stable particles within , including neutrinos, were used for the reconstruction. The reconstructed jets that satisfied the conditions

 pT(jet)>30 GeV/c,  |η(jet)|<4.4,  and  ΔR(γ,jet)>0.5 (10)

were taken as the detected jets, where is the separation in between the photon and jet. We required the detection of two jets in each event. The large threshold for the photons in Eq. (8) was adopted so that hadron jets sufficiently softer than the energy scale of the event are always allowed to be produced. Such soft jets are predominantly produced by PS simulations.

Figure 3 shows the distribution of the separation between the photon and second jet, i.e., the detected hadron jet having the second largest . The smaller of the two values evaluated for both photons is histogramed in the figure. The sum of the results from the 0- () and 1- ( + ) processes is presented in the upper-left panel, and that from the 2- processes ( + + ) is presented in the upper-right panel. The lower panel shows the total sum of the two results, i.e., the sum of the results from all six processes. The dashed, solid, and dotted histograms illustrate the results for the settings of , , and , respectively.

We expected this distribution to be sensitive to the photon radiation from the final-state quark, as in the + 2-jet production [27]. Indeed, an increase of large-angle photon radiations in the 0- + 1- result and a corresponding decrease in the 2- result are clearly observed when we increase the energy scale from to . However, a similar effect is not clear when the scale is increased from to . This is because defines the boundaries for both QED and QCD radiations, and its effect is complicated for large values. Substantially, the summed distribution is very stable against the variation, implying good matching between the ME-based and PS simulations.

The results for the separation between the leading (largest-) jet and second jet are presented in Fig. 4, and those for the second-jet are presented in Fig 5. In both figures, the sum of the results from the 0-jet () and 1-jet ( + ) processes is shown in the upper-left panel, that from the 2-jet ( + + ) processes in the upper-right panel, and the total sum is shown in the lower panel. The other notations are same as those in Fig. 3. We expected the distribution in Fig. 4 to be sensitive to the final-state QCD radiation and the distribution in Fig. 5 to be sensitive to the initial-state QCD radiation. Although the observations are less clear than those in + 2-jet production, we can see a general tendency for the contribution from smaller jet-multiplicity processes to increase as the energy scale increases, and, accordingly, the contribution from the 2-jet processes decreases. The resultant large energy-scale dependences in the separate results compensate for each other to result in stable distributions in the summed results. These observations again confirm the good matching property of the simulation.

## 4 Loop-mediated gg→γγ process

It is known that the contribution from the fermion-loop mediated process to the diphoton production in collisions is unignorable because of the high gluon density inside protons. We implement this process in GR@PPA by hand-coding the differential cross section formula in a literature [34]. Thus, the lowest-order contribution can be included in the simulation. However, some of the kinematical distributions that are compared with the simulation in later sections are sensitive to extra jet radiations. We simulate such radiation effects by adopting a large energy scale in the event generation. We can choose a very large value because this process is independent of the tree-level matched processes described in Section 2. In the present study, we set

 μ=2mγγ (11)

for the generation of events, where is the invariant mass of the produced two photons. The event generation is carried out in the same manner as the tree-level matched processes using the LabCut framework, and the generated events are processed by PYTHIA to obtain the particle-level event information.

Because defines the upper limit of the PS simulation, the QCD radiations from initial-state partons are simulated up to according to the LL approximation in PS with the above setting. This simulation should not be far away from reality because has a good perturbative property, i.e., the NLO correction is reasonably small [35]. Indeed, our simulation reproduces the distribution of the diphoton system evaluated by a resummed NLO calculation [35] reasonably well.

The dashed histograms in Fig. 6 show the () contribution to the kinematical distributions of the two photons measured by the ATLAS experiment [14]. The simulation was carried out with the above setting under the condition for simulating the measurement that is described later. The simulation results are compared with those of the tree-level matched simulations illustrated with solid histograms. The contribution is approximately 16% of the total yield. Although this contribution is substantial in small and small regions, it provides only minor corrections to the distributions in most regions.

## 5 Comparison with LHC measurements

### 5.1 ATLAS measurement

The ATLAS experiment published their measurement results regarding the kinematical distributions of two photons produced in collisions [14]. The measurement is based on 4.9 data at = 7 TeV provided by the LHC. We simulated this measurement following the signal definition described in their report.

The simulation was carried out using the MRST2007LO* PDF [36] because this PDF is widely used for LO simulations in experimental studies. The two generated photons were required to be observed in angular ranges of

 |η(γ)|<1.37  or  1.52<|η(γ)|<2.37 (12)

after completing the simulation down to the particle level. An asymmetric requirement of

 pT(γ1)>25 GeV/c  and  pT(γ2)>22 GeV/c (13)

and an angular separation requirement of

 ΔR(γγ)>0.4 (14)

were imposed on the observed two photons. In addition, an isolation condition of

 EconeT(ΔR<0.4)<4 GeV (15)

was required for both photons. The definition of is same as that in Eq. (9). The UE simulation in PYTHIA was deactivated as in Section 3 because its effect is subtracted from in the measurement together with the pile-up effects. Contrary to the simulation in Section 3, no constraint was imposed on the hadron jets that may be produced in association with the two photons.

The event simulation was carried out for the tree-level matched diphoton processes in Section 2 with the standard setting of and the process described in Section 4. The simulation results were obtained by adding all these simulations. The resultant contribution of each process is summarized in Table 1. We can see that no process dominates the result. Although the contribution of seems to be relatively small, its contribution to the differential cross section is substantial because the negative-weight fraction is large. Here, is the integration of the absolute value of the differential cross section, i.e.,

 σabs=∫∣∣∣dσdΦ∣∣∣dΦ, (16)

where represents the phase space. Events are generated in proportion to and an event weight of or is assigned according to the sign of . The negative-weight fraction is the fraction of events that have the event weight of . Thus, the true cross section can be obtained as

 σ=(1−2Rneg)σabs. (17)

The same calculation is performed in each measurement bin in order to obtain the kinematical distributions.

The simulated kinematical distributions are compared with the measurement results in Fig. 7. Because LO simulations are not capable of predicting absolute values, the total yield of the simulation is normalized to the total cross section of pb in the ATLAS measurement; namely, the simulation results are multiplied by a factor of 0.910. The results in Fig. 7 show that the simulation is in good agreement with the measurement within the measurement errors. Contrary to the SHERPA simulation examined in the ATLAS report [14], we observe no systematic discrepancy significantly larger than the measurement errors. However, although they are contained within the measurement errors, some systematic tendencies can be observed in the and distributions.

We pursued possible improvements by changing tunable parameters in the simulation. Our simulation includes only a few parameters that we can optimize. The overall normalization has already been optimized. Another tunable parameter is the energy scale . Although the simulation results are very stable against the variation of , small dependencies remain as a result of multiple radiation in PS [20]. As a test, we repeated the simulation of the tree-level matched diphoton processes with the settings of and as in the matching test. The simulation of was unchanged. The total cross section increased by 37% for the setting of and decreased by 11% for . However, after renormalizing the total yield, these changes resulted in alterations of only a few percent in most of the measurement bins, although slightly larger alterations were observed in some edge bins of the and distributions where measurement errors are large. Thus, they do not significantly change the observed tendencies.

The PYTHIA simulation includes many tunable parameters. However, most of them are not relevant to the kinematics of high-energy objects. As far as we know, the simulation in PYTHIA that most significantly affects the kinematics is the primordial- simulation [37], which simulates the motion of partons inside hadrons that cannot be reproduced by PS simulations. The default value of the average is set to 2.0 GeV. We examined the effect of this simulation by changing the average value to 0 and 4.0 GeV. The former was tested by deactivating the primordial- simulation. We applied these PYTHIA simulations to all processes including . As a result, we observed a 20%-level alteration in the smallest bin and nearly a 10%-level alteration in the largest bin. However, we did not observe any changes significantly larger than the statistical errors of the simulation in the other bins. Substantially, the simulation results are very stable against the variation of model parameters in the simulation relative to the measurement errors.

### 5.2 CMS measurement

The CMS experiment also performed a measurement regarding diphoton kinematical distributions [15]. The measurement was based on 5.0 fb data of collisions at = 7 TeV provided by the LHC. Although the same quantities as those in the ATLAS measurement were measured, the signal definition was slightly different. The angular coverage was defined as

 |η(γ)|<1.44  or  1.57<|η(γ)|<2.5. (18)

The two photons in this angular range were required to satisfy the following conditions,

 pT(γ1)>40 GeV/c,  pT(γ2)>25 GeV/c,  and  ΔR(γγ)>0.45. (19)

The threshold values are larger than in the ATLAS measurement. Although the photon-isolation condition is not explicitly defined in their report, a condition of

 EconeT(ΔR<0.4)<5 GeV (20)

was required in their simulation studies.

We simulated the measurement following the above conditions, again using the MRST2007LO* PDF. The UE simulation in PYTHIA was deactivated, and muons and neutrinos were excluded from the calculation of , as in the simulation for the ATLAS measurement. The obtained simulation results are compared with the measurement in Fig. 8. The total cross section was measured to be pb by CMS, while the cross section that we obtained from the simulation is 19.48 pb. Hence, the simulation results are multiplied by a factor of 0.883. This normalization factor is consistent with that for the ATLAS measurement. The difference is well within the measurement errors.

In Fig. 8, we again observe overall good agreement between the simulation and measurement. However, the agreement is marginal in the distribution. A significant deficit of the simulation can be observed around = 10 GeV, and a wavy structure is observed in the ratio. This structure is similar to the tendency observed in the comparison with the ATLAS measurement. However, the details of the structure, such as the minimum/maximum positions, are different. The different thresholds may have resulted in this difference.

A similar deficit of the simulation was also observed at GeV in the SHERPA simulation examined by CMS [15], where the deficit continuously remains down to . Note that the simulation/measurement ratio is inverted and the SHERPA simulation is not normalized to the measurement in the CMS study. The SHERPA simulation in the CMS paper shows better agreement with the measurement in the distribution than ours in Fig. 8. The SHERPA simulation examined by CMS includes 3-jet contributions. This must be the reason for the better performance in the distribution, whereas it does not result in a better performance in .

The CMS measurement applies a very asymmetric cut to the photons. The subleading photon cannot have a value smaller than the threshold for the leading photon if there is no additional transverse activity. Additional QCD activities allow the subleading photon to have values in the range between the two thresholds, and . As a result, the distribution at is reduced. Thus, the observed deficit in this range implies that the QCD activity in the simulation is too strong and the activity in SHERPA is stronger than ours. This discussion seems to be inconsistent with the observation in the distribution, where SHERPA provides a better simulation. However, the relevant range of the QCD activity may be different for the two quantities. The observations can be consistently understood if relatively soft activities that are simulated by PS are effective for small values, while the distribution is sensitive to harder activities that are simulated according to multi-jet MEs.

## 6 Discussion

### 6.1 Acoplanar diphoton

In the Tevatron studies, the most significant discrepancy between the measurement and NLO predictions was observed in the distribution at small values. The distribution of our simulation for the ATLAS measurement is shown in Fig. 9. In the figure, the contributions from , , and processes, which dominate the distribution in acoplanar (small ) regions, are separately presented together with the total sum. We can see that these three processes have nearly equal importance in acoplanar regions.

Among these three processes, the contribution decreases as decreases because there is no enhancement of collinear production. This contribution must have been properly included in the NLO predictions because is a process at NLO. In contrast, as discussed in Section 2.3, shows an enhancement of collinear production due to the double radiation. Formally, this process first appears in the NNLO approximation. In addition, events with a hard radiation should formally be categorized as an NNNLO correction. Therefore, the absence or an inappropriate evaluation of the and/or contributions must have been the reason for the deficit in NLO predictions.

Additional hard radiation effects are necessary to produce acoplanar events. Such effects are simulated by increasing the PS energy scale in our simulation, as described in Section 4. Namely, the hard radiation effects are added according to the LL approximation. Usually, the LL approximation becomes less accurate as the radiation becomes comparable with or higher than the typical energy scale of the hard interaction. An enhancement around that we can observe in the result in Figs. 6 and 9 may have been caused by this inaccuracy. It would also be necessary to apply the radiation matching method to if we need more accurate predictions in acoplanar regions.

The hatched areas in Fig. 10 show the contribution of events in an acoplanar region, , to the other distributions. We can see that the low-mass distribution and bump structure in the distribution are totally determined by these acoplanar events. In addition, the acoplanar events predominantly determine the distribution close to , where all theoretical predictions and simulations failed to reproduce the measurements [14, 15]. Our simulation reproduces the measurements very well in these regions, as we can see in Figs. 7 and 8. However, the scattering angle in the Collins-Soper frame has been introduced to probe the properties of underlying hard interactions. We are not sure whether the of acoplanar events is helpful for such studies. It may be better to apply a certain coplanarity-angle cut in the study of the distribution.

### 6.2 Photon isolation

We deactivated the UE simulation in the comparison with the LHC measurements because ATLAS and CMS state in their reports that the UE contributions are subtracted from the cone for the photon isolation and because the UE effects to the examined kinematical distributions are negligible. However, it is not trivial that the quantities simulated by the UE simulation in PYTHIA are identical to the quantities that the experiments subtracted. As a test, we repeated the simulation for the ATLAS measurement by activating the default UE simulation in PYTHIA 6.425.

Figure 11 shows the rejection rate of the isolation cut, defined as

 EconeT(ΔR<0.4)

The results are presented for three typical processes, and those without and with the UE simulation are compared as a function of . The rejection rate is very small in the events because the produced photons do not have any correlation with additional hadronic activities. In contrast, in the events, photons are produced as PS radiations from quarks. Hence, they are likely to be surrounded by hadronic activities induced by the quarks. As a result, the rejection rate becomes high. The results are observed in approximately the middle of them because one of the two photons originates from the PS simulation.

We observe in these results that the UE effect is process dependent. Because the UE particles alter the value, the effect is larger if the spectrum is steeper. The effect would be negligible if the spectrum was flat. In order to examine such a UE effect, we applied the UE simulation to all processes and applied the isolation cut in Eq. (15) to the simulated events. In spite of a clear process dependence in Fig. 11, the simulated normalized kinematical distributions did not show any significant difference to those in Fig. 7. The differences were well contained within the statistical error of the simulations, although the total cross section decreased by 10% with respect to the value obtained from the simulation in Section 5.1.

The above is a discussion of extreme cases. The ATLAS and CMS experiments describe that they subtract the UE contribution, together with the pile-up effects, from event-by-event using the jet-area/median method [38] implemented in the FastJet package [33]. Although the experiments must have applied the subtraction to the detector-level data, we tried to apply it to particle-level simulation data in order to examine the effect of this method. As in the jet clustering in Section 3, we used all stable particles in the angular range of in the simulated events to which the UE simulation was applied, in order to evaluate the average diffuse-background (DB) density () on the - plane. The recommended Cambridge/Aachen algorithm with = 0.6 was used for the jet clustering. However, the evaluation frequently failed and returned = 0. This was because the DB particles were too sparse with the UE simulation only. Indeed, although the quantity is required to be sufficiently smaller than unity in this method [38], its value was frequently close to unity even when a non-zero was returned. Here, and denote the internal fluctuation parameter and average jet-area, respectively, which are returned by the FastJet program.

We found that this problem can be solved by adding fake particles that imitate the pile-up particles. Here, we call them pseudo-pile-up (PPU) particles. We added the PPU particles uniformly and randomly on the - plane within the range of . The values were determined according to a Gaussian distribution with an average of 1.0 GeV/ and root-mean-square of 0.2 GeV/. A Gaussian distribution was chosen because the average density () is trivial and it is easy to tune the fluctuation. We also tried to use an exponential distribution that simulates the low- spectrum of particles in the simulated events. Although the obtained results were similar, we observed a slightly better performance with the Gaussian distribution in the -parameter distribution.

We added 500 such PPU particles to each event. They were assumed to be massless. From the value returned from the FastJet program, we evaluated the UE density as . The UE contribution to the isolation-cone was estimated to be , where with = 0.4. The estimated UE contribution was subtracted from the measured before applying the isolation cut in Eq. (21). The dotted lines in Fig. 11 show the rejection rates for the isolation condition thus defined. We can see that the subtraction cannot totally remove the activities produced by the UE simulation of PYTHIA. In PYTHIA, UE is simulated by low- QCD 2-jet production. The results in Fig. 11 are reasonable because the thus-simulated UE particles may partly form non-diffusive jet-like structures.

Although we observed no significant difference in the normalized kinematical distributions, the simulation with the DB subtraction and the isolation condition in Eq. (15) gives a total cross section approximately 4% smaller than the value that was obtained in Section 5.1. Hence, looking at the UE-on/off and DB-subtraction results, we have to consider that the cross section result that the experiment presented has an additional 5%-level uncertainty because the corresponding signal definition has an ambiguity of this level. Although this uncertainty is smaller than the current measurement precision, it would be a serious concern in future higher-precision measurements. This uncertainty substantially arises because the experiments do not explicitly define the treatment of the UE effect. At least when numerical results are presented, the corresponding signal definition has to be unambiguously presented by experiments. The current definitions of the LHC experiments are ambiguous because the UE subtraction is not defined using detector-independent physical quantities. Furthermore, because UE is a part of the interaction, it may not be reasonable to exclude its contribution from the definition of the isolation condition. We need further discussion to reach a reasonable consensus on this point.

### 6.3 Two-jet contribution

As discussed in the introduction, we suspected that the contribution of 2-jet production processes would be large in diphoton production based on the observation in direct-photon production [27]. From the results in Table 1, we can evaluate the 2-jet contribution + + to be approximately 20% of the total yield in our standard simulation for the ATLAS measurement. This value is not small, but it is also not very large in the aspect of the QCD correction. Furthermore, the contribution of new processes, which we suspected to be large, is relatively small, as already pointed out in a study of the NNLO approximation [11]. The contribution is approximately 3% from and 5% from () in our simulation. The remainder is dominated by the ordinary gluon-radiation correction to . Its contribution amounts to approximately 10%. The contribution of initial-state processes is again very small (1.8%). Among them, the contribution is only 0.3%. The remainder comes from other new processes: 1.0% from and 0.5% from ().

Simulation results show that the limited 2-jet contribution is partly due to the requirement regarding the photon isolation. As we can expect from interaction diagrams, 2-jet processes, particularly the new processes, have larger rejection probabilities than the other processes. If we remove the isolation requirement, the and contributions increase to 6% and 11%, respectively, and the total 2-jet contribution increases to 34%. This is comparable with the 1-jet () contribution of 41%.

Although the impact of 2-jet processes on the total cross section is not very large under the actual measurement condition, they have a remarkable contribution to kinematical distributions of the two photons. The simulation results in which the 2-jet processes are excluded are compared with the ATLAS measurement results in Fig. 12. These results should be compared with those in Fig. 7. The total yield of the simulation is again normalized to the measured total cross section. From these results, we observe that the contribution of the 2-jet processes depends on the examined quantity.

The distribution of this up-to-1-jet simulation is already in good agreement with the measurement, except for the lowest-mass bin. In contrast, significant discrepancies are observed in the other distributions. This is reasonable because is determined by hard interactions to produce the two photons, while the and distributions are sensitive to additional QCD radiations. The obvious deficit of the simulation in small regions is caused by the absence of events, as discussed in Section 6.1. This deficit results in the discrepancy near in the distribution. The up-to-1-jet simulation is already in good agreement with the measurement in central regions because is predominantly determined by hard interactions. The remarkable improvement in the distribution at high , GeV/, indicates a significant contribution of hard radiation effects in the 2-jet processes. The reason for the remaining deficit around = 10 GeV/ must be different from that for the deficit at small that we observed in Section 5.2, because the difference between the two thresholds is only 3 GeV/ in the ATLAS measurement.

As we discussed above, the 2-jet production processes have a large contribution to diphoton production, at least when we remove the photon-isolation requirement. However, contrary to the initial conjecture, the large contribution is not predominantly caused by the emergence of new processes. The largest contribution comes from the gluon-radiation correction to . A large 2-jet contribution was also observed in direct-photon production in our previous study [27]. Reexamining the simulation results, we found that the underlying properties are almost the same as for diphoton production. While the 2-jet contribution amounts to 42% of the total cross section, the and contributions are only 3% and 7%, respectively, under the condition applied in Section 4.1 of the previous paper [27]. The contribution amounts to approximately 30%. This corresponds to a more than 50% correction to the dominant lowest-order process .

An incomprehensible behavior has also been observed in ordinary perturbative calculations for diphoton production [11, 39]. It was found that the energy-scale dependence increases when the approximation is improved from NLO to NNLO. Namely, the naive perturbative nature seems to be violated at NNLO, which is the approximation order including associated 2-jet production. This strange behavior may be related to the large gluon-radiation contribution in 2-jet processes.

Although the large gluon-radiation contribution is still a puzzle, the large 2-jet contribution can be considered to be a common property of photon-production processes in hadron collisions. Furthermore, this property must also be common to other electroweak-boson production processes, such as single and production and diboson (, , and ) production, because the underlying QCD structure is identical. Indeed, we observed a large 2-jet contribution to -boson production at high in a previous study [26]. The 2-jet contribution is not significant and thus the NLO approximations provide good predictions in low- weak-boson production, probably because the amplitudes of final-state radiation diagrams, such as those corresponding to the diagrams in Fig. 1, are suppressed by the weak-boson mass.

If the above discussion is also valid for high-mass diboson production, we need to take special care regarding the collinear jet production in association with weak-boson production, especially when the weak bosons are detected in hadronic-decay modes as highly-boosted massive objects. It must be difficult to separate the associated jet from the decay products and the weak-boson momenta may be overestimated in such studies. The probability of collinear jet production may be large because approximately 35% of the events are rejected by the isolation requirement in diphoton production under the ATLAS measurement condition. It would hence be necessary to carry out simulation studies with appropriate event generators including associated multi-jet production.

## 7 Conclusion

We have developed an event generator for diphoton () production in hadron collisions that includes associated jet production up to two jets within the framework of the GR@PPA event generator. Processes having different jet and photon multiplicities are combined using a subtraction method based on the LLL subtraction. The subtracted divergent components in radiation-including processes are consistently restored by combining lower multiplicity processes to which PS simulations are applied.

The PS simulation involves QED photon radiations from final-state quarks to restore subtracted final-state QED collinear components. Photon radiations in very small regions that the PS simulation cannot cover are simulated by employing a fragmentation function (FF). This PS/FF simulation has the ability to enforce energetic photon radiations for efficient event generation. In principle, the radiation of any number of photons can be enforced in this simulation. The generated parton-level events can be fed to the PYTHIA event generator to obtain particle-level event information. We can perform realistic photon isolation and hadron-jet reconstruction simulations using the obtained events.

The simulated events, in which the loop-mediated process is also involved, reasonably reproduce the diphoton production kinematics measured at the LHC. The remaining small discrepancies in the azimuthal opening angle of the diphoton system indicate the necessity of further higher-order processes. A question still remains regarding the small discrepancies in the distribution.

We found that the contribution of 2-jet processes is significant in diphoton production, even in the production kinematics of the two photons. Contrary to the initial conjecture, the contribution of new processes and that first emerge in 2-jet production is not very large. The largest contribution comes from gluon-radiation corrections to . Requirements on the photon isolation reduce the 2-jet contribution, especially those from the new processes.

The significant 2-jet contribution can be considered as a common property of photon-production processes in hadron collisions. It may also be common to high- weak-boson production. We need to be careful about the contamination of collinearly produced hadron jets when weak bosons are identified as boosted massive jets.

The diphoton production associated with two jets is not yet well understood. We cannot reasonably explain why the gluon-radiation corrections are large, and the naive perturbative nature does not seem to hold in ordinary perturbative calculations at NNLO. These two observations may be related to each other. We still need further studies in order to improve our understanding.

Improvements in the measurements are also necessary for better understanding. The treatment of the underlying events in the photon-isolation condition may become an obstacle to further improvement. The contribution from the underlying events is subtracted in the LHC measurements, but the subtraction is not defined in a detector-independent form. The ambiguity in the definition may lead to a cross section uncertainty at the level of %. We need to obtain a reasonable consensus concerning the unambiguous signal definition.

## Acknowledgments

This work has been carried out as an activity of the NLO Working Group, a collaboration between the Japanese ATLAS group and the numerical analysis group (Minami-Tateya group) at KEK. The authors wish to acknowledge useful discussions with the members. The programs for the process was coded by Y. Komori of the University of Tokyo.

### Footnotes

1. E-mail: shigeru.odaka@kek.jp.

### References

1. ATLAS Collaboration, G. Aad et al., Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B 716 (2012) 1; arXiv:1207.7214.
2. CMS Collaboration, S. Chatrchyan et al., Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Phys. Lett. B 716 (2012) 30; arXiv:1207.7235.
3. CMS Collaboration, V. Khachatryan et al., Observation of the diphoton decay of the Higgs boson and measurement of its properties, Eur. Phys. J. C 74, 3076 (2014); arXiv:1407.0558.
4. ATLAS Collaboration, G. Aad et al., Measurements of fiducial and differential cross sections for Higgs boson production in the diphoton decay channel at TeV with ATLAS, JHEP 09 (2014) 112; arXiv:1407.4222.
5. D0 Collaboration, V. Abazov et al., Measurement of direct photon pair production cross sections in collisions at TeV, Phys. Lett. B 690 (2010) 108; arXiv:1002.4917.
6. CDF Collaboration, T. Aaltonen et al., Measurement of the Cross Section for Prompt Isolated Diphoton Production in Collisions at TeV, Phys. Rev. Lett. 107, 102003 (2011); arXiv:1106.5123.
7. CDF Collaboration, T. Aaltonen et al., Measurement of the Cross Section for Prompt Isolated Diphoton Production in Collisions at TeV, Phys. Rev. D 84, 052006 (2011); arXiv:1106.5131.
8. ATLAS Collaboration, G. Aad et al., Measurement of the isolated diphoton cross-section in pp collisions at TeV with the ATLAS detector, Phys. Rev. D 85, 012003 (2012); arXiv:1107.0581.
9. CMS Collaboration, S. Chatrchyan et al., Measurement of the Production Cross Section for Pairs of Isolated Photons in pp collisions at TeV, JHEP 01 (2012) 133; arXiv:1110.6461.
10. L. Cieri, Diphoton production at Next-to-Next-to-Leading-Order, Nucl. Phys. B (Proc. Suppl.) 234 (2013) 29; arXiv:1209.3143.
11. S. Catani, L. Cieri, D. de Florian, G. Ferrera, and M. Grazzini, Diphoton production at hadron colliders: a fully-differential QCD calculation at NNLO, Phys. Rev. Lett. 108, 072001 (2012); arXiv:1110.2375.
12. T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert and J. Winter, Event generation with SHERPA 1.1, JHEP 02 (2009) 007; arXiv:0811.4622.
13. S. Hoeche, S. Schumann, and F. Siegert, Hard photon production and matrix-element parton-shower merging, Phys. Rev. D 81, 034026 (2010); arXiv:0912.3501.
14. ATLAS Collaboration, G. Aad et al., Measurement of isolated-photon pair production in collisions at TeV with the ATLAS detector, JHEP 01 (2013) 086; arXiv:1211.1913.
15. CMS Collaboration, S. Chatrchyan et al., Measurement of differential cross sections for the production of a pair of isolated photons in pp collisions at , Eur. Phys. J. C 74, 3129 (2014); arXiv:1405.7225.
16. CDF Collaboration, T. Aaltonen et al., Measurement of the cross section for prompt isolated diphoton production using the full CDF Run II data sample, Phys. Rev. Lett. 110, 101801 (2013); arXiv:1212.4204.
17. S. Odaka and Y. Kurihara, Consistent simulation of non-resonant diphoton production at hadron collisions with a custom-made parton shower, Phys. Rev. D 85, 114022 (2012); arXiv:1203.4038.
18. Y. Kurihara et al., QCD event generators with next-to-leading order matrix-elements and parton showers, Nucl. Phys. B 654 (2003) 301; hep-ph/0212216.
19. S. Odaka and Y. Kurihara, Initial-state parton shower kinematics for NLO event generators, Eur. Phys. J. C 51, 867 (2007); hep-ph/0702138.
20. S. Odaka, Simulation of boson spectrum at Tevatron by leading-order event generators, Mod. Phys. Lett. A 25, 3047 (2010); arXiv:0907.5056.
21. S. Odaka and Y. Kurihara, GR@PPA 2.8: Initial-state jet matching for weak-boson production processes at hadron collisions, Comput. Phys. Commun. 183 (2012) 1014; arXiv:1107.4467.
22. S. Odaka, GR@PPA 2.8.3 update, arXiv:1201.5702.
23. S. Catani, F. Krauss, R. Kuhn, and B. R. Webber, QCD Matrix Elements + Parton Showers, JHEP 11 (2001) 063; hep-ph/0109231.
24. L. Bourhis, M. Fontannaz, and J. Guillet, Quarks and gluon fragmentation functions into photons, Eur. Phys. J. C 2, 529 (1998); hep-ph/9704447.
25. T. Sjostrand, S. Mrenna, and P. Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05 (2006) 026; hep-ph/0603175.
26. S. Odaka, N. Watanabe, and Y. Kurihara, ME-PS matching in the simulation of multi-jet production in hadron collisions using a subtraction method, PTEP 2015 (2015) 053B04; arXiv:1409.3334.
27. S. Odaka and Y. Kurihara, Consistent simulation of direct-photon production in hadron collisions including associated two-jet production, Mod. Phys. Lett. A 31, 1650099 (2016); arXiv:1509.04813.
28. MINAMI-TATEYA Group, T. Ishikawa et al., GRACE manual: Automatic generation of tree amplitudes in Standard Models: Version 1.0, KEK Report, KEK-92-19 (1993).
29. F. Yuasa et al., Automatic computation of cross sections in HEP: Status of GRACE system, Prog. Theor. Phys. Suppl. 138 (2000) 18; hep-ph/0007053.
30. S. Kawabata, A New Monte Carlo Event Generator for High-Energy Physics, Comput. Phys. Commun. 41 (1986) 127.
31. S. Kawabata, A new version of the multidimensional integration and event generation package BASES/SPRING, Comp. Phys. Commun. 88 (1995) 309.
32. J. Pumplin et al., New generation of parton distributions with uncertainties from global QCD analysis, JHEP 07 (2002) 012; hep-ph/0201195.
33. M. Cacciari, G. P. Salam, and G. Soyez, FastJet User Manual, Eur. Phys. J. C 72, 1896 (2012); arXiv:1111.6097.
34. E. L. Berger, E. Braaten, and R. D. Field, Large- production of single and double photons in proton-proton and pion-proton collisions, Nucl. Phys. B 239 (1984) 52.
35. P. M. Nadolsky, C. Balazs, E. L. Berger, and C. P. Yuan, Gluon-gluon contributions to the production of continuum diphoton pairs at hadron colliders, Phys. Rev. D 76, 013008 (2007); hep-ph/0702003.
36. A. Sherstnev and R. S. Thorne, Parton Distributions for LO Generators, Eur. Phys. J. C 55, 553 (2008); arXiv:0711.2473.
37. S. Odaka, Precise simulation of the initial-state QCD activity associated with -boson production in hadron collisions, Mod. Phys. Lett. A 28, 1350098 (2013); arXiv:1301.5082.
38. M. Cacciari and G. P. Salam, Pileup subtraction using jet areas, Phys. Lett. B 659 (2008) 119; arXiv:0707.1378.
39. J. M. Campbell, R. K. Ellis, Y. Li and C. Williams, Predictions for diphoton production at the LHC through NNLO in QCD, JHEP 07 (2016) 148; arXiv:1603.02663.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters