Edinburgh Research Explorer Measurement of the production cross section of three isolated photons in $pp$ collisions at $\sqrt{s}$ = 8 TeV using the ATLAS detector

A measurement of the production of three isolated photons in proton–proton collisions at a centre-of- mass energy √ s = 8 TeV is reported. The results are based on an integrated luminosity of 20.2 fb − 1 collected with the ATLAS detector at the LHC. The differential cross sections are measured as functions of the transverse energy of each photon, the difference in azimuthal angle and in pseudorapidity between pairs of photons, the invariant mass of pairs of photons, and the invariant mass of the triphoton system. A measurement of the inclusive ﬁducial cross section is also reported. Next-to-leading-order perturbative QCD predictions are compared to the cross-section measurements. The predictions underestimate the measurement of the inclusive ﬁducial cross section and the differential measurements at low photon transverse energies and invariant masses. They provide adequate descriptions of the measurements at high values of the photon transverse energies, invariant mass of pairs of photons, and invariant mass of the triphoton system. © 2018 The 3 .


Introduction
The production of three prompt photons in proton-proton (pp) collisions, pp → γ γ γ + X , provides a testing ground for perturbative quantum chromodynamics (pQCD). This process is rare in the Standard Model (SM) since the leading-order (LO) contribution to triphoton production is of order α 3 EM . The measurement of triphoton production can be performed in a broader range of kinematic regions than in 2 → 2 reactions such as inclusive-photon [1][2][3][4] and diphoton [5][6][7] production. This provides a complementary test of pQCD in processes with photons in the final state.
Precise measurements of triphoton production can be used to improve the description of this process in Monte Carlo (MC) models. In addition, SM triphoton production provides one of the main irreducible backgrounds for some beyond-the-SM (BSM) searches. Potential BSM processes include the associated production of a photon and an exotic neutral particle decaying into a photon pair (qq → X 0 γ ), where X 0 can be a Kaluza-Klein graviton (G KK ) [8-10] or a pseudoscalar (a) [11]. Moreover, triphoton production is also the main background to the predicted decay of the Z boson into three photons. The current upper limit at 95% confidence level on the branching fraction for Z → 3γ is 2.2 × 10 −6 [12].
Three photons can be produced via two main mechanisms: direct and fragmentation production. In the case of the direct E-mail address: atlas .publications @cern .ch. production process, three photons are produced in the hard interaction via the annihilation of an initial-state quark-antiquark pair (qq → γ γ γ ). In the fragmentation process, at least one of the photons arises from the fragmentation of a high-transversemomentum (high-p T ) parton (qg → γ γ q[γ ]). Direct photons are typically isolated, while those originating from the fragmentation process are usually accompanied by nearby partons. Measurements of final-state photons include an isolation requirement to reduce background contributions from neutral-hadron decays into photons. As a consequence, signal processes with one or more fragmentation photons are also suppressed.
This Letter presents measurements of three-photon production. The analysis is performed using 20.2 ± 0.4 fb −1 of ATLAS data at a centre-of-mass energy of √ s = 8 TeV [13]. The measurements study the topology and kinematics of the individual photons, pairs of photons, and the three-photon system. Differential cross sections are measured as functions of the transverse energy 1 of the leading photon (E γ 1 in azimuthal angle and in pseudorapidity between pairs of photons ( φ γ 1 γ 2 , φ γ 1 γ 3 , φ γ 2 γ 3 , | η γ 1 γ 2 |, | η γ 1 γ 3 |, | η γ 2 γ 3 |); the invariant mass of pairs of photons (m γ 1 γ 2 , m γ 1 γ 3 and m γ 2 γ 3 ); and the invariant mass of the triphoton system (m γ γ γ ). A measurement of the inclusive fiducial cross section is also reported. Photons are required to be isolated based on the amount of transverse energy, excluding the photon contribution, inside a cone of size R ≡ (η − η γ ) 2 + (φ − φ γ ) 2 = 0.4 centred around each photon direction (defined by the photon pseudorapidity η γ and azimuthal angle φ γ ). Finally, the measurements are compared to next-toleading-order (NLO) QCD calculations.

ATLAS detector
The ATLAS detector [14] is a multi-purpose detector with a forward-backward symmetric cylindrical geometry. The most relevant systems for the present measurement are the inner detector, immersed in a 2 T magnetic field produced by a thin superconducting solenoid, and the calorimeters. At small radii, the inner detector is made up of fine-granularity pixel and microstrip detectors.
These silicon-based detectors cover the pseudorapidity range |η| < 2.5. A gas-filled straw-tube transition radiation tracker complements the silicon tracker at larger radii in the range |η| < 2.0 and also provides electron identification capabilities based on transition radiation. The electromagnetic calorimeter is a lead/liquid-argon sampling calorimeter with accordion geometry. The calorimeter is divided into a barrel section covering |η| < 1.475 and two endcap sections covering 1.375 < |η| < 3.2. For |η| < 2.5 it is divided into three layers in depth, which are finely segmented in η and φ.
A thin presampler layer, covering |η| < 1.8, is used to correct for fluctuations in upstream energy losses. The hadronic calorimeter in the region |η| < 1.7 uses steel absorbers with scintillator tiles as the active medium. Liquid-argon with copper absorbers is used in the hadronic end-cap calorimeters, which cover the region 1.5 < |η| < 3.2. Events are selected using a first-level trigger implemented in custom electronics, which reduces the event rate to a value of 75 kHz using a subset of detector information. Software algorithms with access to the full detector information are then used in the high-level trigger to yield a recorded event rate of about 400 Hz [15].

Monte Carlo simulations
The MC samples were generated to study the characteristics of the signal and background events. The MC program MadGraph 5.1.4.4 [16] interfaced with Pythia 8.186 [17] was used to simulate signal events. The partonic subprocess was simulated by MadGraph to include the leading-order matrix element (qq → γ γ γ ), whereas Pythia was added to include the initial-and final-state parton showers and the fragmentation of partons into hadrons. The LO CTEQ6L1 parton distribution functions (PDFs) [18] are used to parameterise the parton momentum distributions in the proton. To study the effect of the contribution of photon fragmentation, a Pythia MC sample supplemented by QED final-state radiation was generated with LO CTEQ6L1 PDFs. This sample includes the LO diphoton, photon+jet and dijet processes with initial-state and final-state radiation modelled by the parton shower (PS).
The MC program Sherpa 1.4.1 [19] was used to estimate the background arising from electrons misreconstructed as photons. Three processes were simulated with at least one high-p T electron and photon in the final state: e + e − γ , e + e − γ γ , and e ± ν e γ γ .
The matrix elements were calculated with up to three final-state partons at LO in pQCD and used the CT10 PDFs at NLO [20]. The matrix elements were merged with the Sherpa parton-shower algorithm [21] following the ME+PS@LO prescription [22].
The generated signal and background event samples were passed through the Geant4-based [23] ATLAS detector and trigger simulation programs [24]. The signal and background samples include a simulation of the underlying event (UE) where Pythia event-generator parameters were set according to the "AU2" tune [25]. The generation of the simulated event samples includes the effect of multiple pp interactions per bunch crossing, as well as the effect of the detector response to interactions from bunch crossings before or after the one containing the hard interaction. These MC events were weighted to reproduce the distribution of the average number of interactions per bunch crossing observed in the data. The generated MC events are reconstructed and analysed with the same program chain as the data.

Next-to-leading-order pQCD predictions
The NLO pQCD predictions presented in this Letter are computed using the programs MCFM [26,27] and MadGraph5_ aMC@NLO 2.3.3 [28]. The strong coupling constant is calculated at two loops with α S (m Z ) = 0.118 and the electromagnetic coupling constant is set to α EM = 1/137. In addition, the number of massless quark flavours is set to five and the CT10 parameterisations of the proton PDFs at NLO are used.
The MCFM program includes NLO pQCD calculations of the direct contribution, whereas the production of a photon via parton fragmentation is estimated from the LO QCD matrix element multiplied by the BFG II parton-to-photon fragmentation functions [29].
The renormalisation scale μ R , factorisation scale μ F and fragmentation scale μ f are chosen to be μ R = μ F = μ f = m γ γ γ . In addition, the MCFM calculations are performed using an isolation criterion which requires the total transverse energy from the partons inside a cone of size R = 0.4 around the photon direction to satisfy E iso T < 10 GeV. The MCFM NLO pQCD predictions refer to the parton level while the measurements are performed at the particle level. Since the E iso T requirement at the particle level is applied after the subtraction of the UE transverse energy, it is expected that parton-to-hadron corrections to the NLO pQCD predictions are small. This is confirmed by computing the ratio of the particle-level cross section for a MadGraph sample interfaced with Pythia with UE effects to the computed cross section without hadronisation and UE effects. The ratio is consistent with unity over the measured range of the variables under study. Therefore, no correction is applied to the MCFM NLO pQCD calculations.
Deviations from unity of O (1%) on the parton-to-hadron correction factors are found when the hadronisation and UE effects are included using Herwig++ 7.0.1 [30]. Predictions based on other proton PDF sets, namely MSTW2008 [31] and NNPDF2.1 [32], are also computed. Differences of +5% and +6% in the calculation of the inclusive fiducial cross section are found using the MSTW2008 and NNPDF2.1 PDF sets, respectively, whereas the dependence of the shape of the differential cross sections on the PDF sets is found to be small.
MadGraph5_aMC@NLO calculations include the NLO pQCD contribution of direct processes and apply a smoothly varying isolation cone to the photons [33]. This isolation requirement regularises the photon collinear divergences which appear in the calculation of the matrix element and removes the contribution of photons resulting from the fragmentation of a parton: E iso is the sum of the transverse energies of the particles around the photon up to R. The MadGraph5_aMC@NLO calculations are interfaced with Pythia 8.212 [34] in the NLO+PS prescription to include the initial-and final-state parton showers and the hadroni-sation [35]. The renormalisation and factorisation scales are chosen to be equal to the transverse mass of the clustered jets from the final state partons and photons defined in the matrix element. This choice follows the recommendations in Ref. [28] when interfacing the MadGraph5_aMC@NLO calculations to Pythia. After the generation, the isolation value of the photon is computed by summing the transverse energy of all final-state particles (excluding muons and neutrinos) inside a cone of size R = 0.4 around the photon candidate. Events with E iso T > 10 GeV for any of the photons are excluded.

Event selection
The data considered in this analysis were taken in stable beam conditions and satisfy detector and data-quality requirements. Events are recorded using a diphoton trigger with a transverse energy threshold of 20 GeV. The trigger efficiency for pairs of isolated photons with E γ T > 22 GeV and |η γ | < 2.37 is higher than 99%. Events are required to have a reconstructed primary vertex with at least two associated tracks with p T > 500 MeV and |η| < 2.5, consistent with originating from the same threedimensional spot within the luminous region of the colliding proton beams. If multiple primary vertices are reconstructed, the one with the highest sum of the p 2 T of the associated tracks is selected as the primary vertex.
Photon and electron candidates are reconstructed from clusters of energy deposited in the electromagnetic calorimeter. Candidates without a matching track or reconstructed conversion vertex in the inner detector are classified as unconverted photons [36]. Those with a matching reconstructed conversion vertex or a matching track consistent with originating from a photon conversion are classified as converted photons. Photons reconstructed within |η γ | < 2.37 are retained. Those in the transition region between the barrel and end-caps (1.37 < |η γ | < 1.56) or regions of the calorimeter affected by read-out or high-voltage failures are not considered in the event reconstruction.
Photon candidates passing loose identification requirements, based on the energy leaking into the hadronic calorimeter and the lateral shower shape in the second layer of the electromagnetic calorimeter, are retained [1,2]. The photon cluster energies are corrected using an in situ calibration based on the Z → e + e − reconstructed mass peak [37]. Once these corrections are applied, the three reconstructed photons with the highest transverse energies E R distance in the η-φ plane above 0.45 between pairs of photons, are selected. Additionally, the invariant mass of the triphoton system m γ γ γ is required to be above 50 GeV. This requirement corresponds to the minimum value of m γ γ γ predicted at particle level by the signal MC sample described in Section 3.
Two further criteria are used to define the signal region and the background-enriched regions used to estimate the jet-to-photon misidentification background. A tight photon-identification selection [36] is applied to reject hadronic jet background, by imposing requirements on nine discriminating variables (referred to as "shower shapes") computed from the energy leaking into the hadronic calorimeter and the lateral and longitudinal shower development in the electromagnetic calorimeter. The efficiency of this selection for one photon is ≈67% (>90%) for E γ T ≈ 15 GeV (>100 GeV). For the MC simulations, the shower-shape variables are shifted to correct for small differences in the average values between data and the simulation. In addition, E γ T -and η γ -dependent factors are applied to correct for the residual mismatch between the photon identification efficiencies in the simulation and the data. The isolation of the photon E iso T is based on the amount of transverse energy inside a cone of size R = 0.4 in the η-φ plane around the photon candidates, excluding an area of size η × φ = 0.125 × 0.175 centred on the photon energy cluster.
The isolation transverse energy is computed from the topological clusters of calorimeter cells [38]. The measured E iso T is corrected for the leakage of the photon's energy into the isolation cone and the estimated contributions from the UE and pile-up. These latter two corrections are computed simultaneously on an eventby-event basis and the combined correction is typically between 1.5 and 2.0 GeV [3]. The E iso T value for isolated photons is required to be lower than E iso . The efficiency of the isolation requirement is typically above 80% and increases as a function of E γ T . The number of data events selected in the signal region is 1085. For background studies, two alternative categories of photons are defined. First, non-tight photon candidates are defined as those passing the loose selection but not satisfying the tight identification criteria for at least one of the shower-shape variables computed from the energy deposits in cells of the first layer of the EM calorimeter. Second, non-isolated photon candidates are defined to have E iso

Background estimation and signal extraction
The background contributions to the signal come from high-p T jets and electrons that are misidentified as isolated photons (referred to as jet and electron backgrounds). The estimation of these backgrounds is explained in the following.

e-γ misidentification
The number of background events due to e-γ misidentification is estimated using the MC samples listed in Section 3.1. The Sherpa MC events were weighted to correct the e-γ misidentification rates to match those found in data (referred to as e-γ scale factors in the following). These weights were estimated from Z → e + e − events where at least either the electron or the positron was reconstructed as a photon. The expected number of electron background events in the signal region is 71 ± 2 (stat), which corresponds to (6.5 ±0.2)% of the selected events. A systematic uncertainty is computed by propagating the uncertainty in the e-γ scale factors to the estimation of the yield (see Section 7).
The normalisation of the MC samples is tested by fitting the signal, e-γ and jet-γ misidentification contributions to the data as a function of m γ γ γ in the region 50 < m γ γ γ < 125 GeV. Since 86% of electron background events come from processes where a photon is emitted by an electron or positron originating from the decay of a Z boson (pp → Z → e + e − γ ), a peak around m γ γ γ ≈ m Z is expected. To enhance the relative contribution of electrons that are misidentified as photons, only events with at least one converted photon are considered. Signal and electron background MC events are used to describe the shape of the m γ γ γ distribution, whereas data events with at least one non-tight identified photon are used to describe the jet background contribution. The fit gives an electron background yield that is consistent with the MC estimation, since it predicts a correction factor equal to 1.0 ± 0.4 (stat). Moreover, the result of the fit is found to be independent of the definition of non-tight identified photons and a change of <2% is found when the isolation requirement is loosened by 1 GeV.

Jet-γ misidentification
A large background from jet-γ misidentification remains in the selected sample, even after imposing the tight identification and isolation requirements on the photons. The jet background originates from multi-jet ( jjj), photon + jets (γ jj), and diphoton + jets (γ γ j) processes in which at least one jet is misidentified as a photon. The two-dimensional-sideband method exploited in Refs. [2,3,5,39-41] to measure the inclusive photon and diphoton differential cross sections is used to perform an in situ statistical subtraction of the background. The method uses the photon isolation energy and photon identification criteria to discriminate prompt photons from jets. It relies on the fact that the correlations between the isolation and identification variables in jet background events are small, and that the signal contamination in the non-tight or non-isolated control region is low.
The two-dimensional-sideband method counts all combinations of photons meeting or failing to meet the tight identification or isolation criteria. Four categories are defined for each photon, resulting in 64 categories of events where 63 of these categories correspond to jjj, γ jj, and γ γ j background-enriched regions. The inputs of the method are the number of events in each category, the correlation between the isolation and identification variables in jet background events (R bg ), the signal leakage fractions in nontight and non-isolated regions, and the expected number of electron background events in each category. The correlation between the isolation and identification variables is taken to be negligible (R bg = 1.0) based on studies in simulated background samples and on data in a background-dominated region [3]. The signal leakage fractions and electron-background events are estimated using the MC samples described in Section 3.1.
The method allows the extraction of the number of true threephoton signal events (N γ γ γ ), the number of events where at least one, two and three candidates are true jets and the tight and isolation efficiencies for fake photon candidates from jets ("fake rates"). The number of events in each category is expressed as a function of the following parameters: signal, electron-and jet-background yields, signal leakage fractions, fake rates, and R bg . Then, the system of 64 independent equations is grouped into 21 dependent linear equations which are solved iteratively using a χ 2 minimisation procedure. The size of each bin of the observables under study is chosen to have a sufficiently large number of events to apply this method bin-by-bin. The statistical uncertainty of the signal and jet background-enriched regions is propagated to the estimation of the three-photon signal yield via pseudo-experiments.
The signal purity, defined as N γ γ γ /N SR , where N SR is the number of selected events in the signal region, is found to be (55 ± 5)% (stat), with a value of ≈45% (≈60%) at low (high) E γ T .

Unfolding to particle level
The production cross section for three isolated photons is measured as functions of E | η γ 1 γ 2 |, | η γ 1 γ 3 |, | η γ 2 γ 3 |, m γ 1 γ 2 , m γ 1 γ 3 , m γ 2 γ 3 and m γ γ γ . The fiducial phase-space region is listed in Table 1. The predictions of the MC generators at particle level are defined using those particles with a lifetime τ longer than 30 ps; these particles are referred to as "stable". The particles associated with the overlaid pp collisions are not considered. The particle-level isolation requirement on the photons is built by summing the transverse energy of all stable particles, except for muons and neutrinos, in a cone of size R = 0.4 around the photon direction. The contribution from the UE is subtracted using the same procedure as applied to the data at the reconstruction level [3]. The data distributions after background subtraction are unfolded to the particle level using Table 1 Fiducial phase-space region defined at particle level.

Requirements on
Phase-space region bin-by-bin correction factors determined using the signal MC sample. The correction factors take into account the efficiency of the event and photon selection criteria and the small migration effects.
Of the signal events reconstructed in a given bin, the fraction that are generated in the same bin is typically found to be >93%. The data distributions are unfolded to the particle level via the formula dσ d A where for a given bin i, (dσ /dA) is the differential cross section as a function of observable A, N sig is the number of backgroundsubtracted data events, C is the correction factor, L is the integrated luminosity and A is the width of the bin. The correction factors are computed using the MC sample of events as

Experimental uncertainties
The sources of experimental systematic uncertainty that affect the measurements are the photon energy scale and resolution, photon identification, jet and electron background subtraction, modelling of the photon isolation, the photon fragmentation contribution, the unfolding procedure and the luminosity.
• Photon energy scale and resolution. The uncertainty due to the photon energy scale is estimated by varying the photon energies in the MC simulation [37]. This uncertainty mostly affects the C (i) correction factor. The effect of this variation on the estimation of the cross section is typically <2%. In addition, the uncertainty in the energy resolution is estimated by smearing photon energies in the MC simulation as described in Ref. [37]. The resulting uncertainty in the cross section is typically < 0.1%. • Photon identification efficiency. The uncertainty in the photon identification efficiency is estimated from the effect of differences between shower-shape variable distributions in data and simulation [36]. This uncertainty affects the estimation of the non-tight signal leakage fractions and the C (i) correction factor and is fully correlated between photons. The correlation between tight and non-tight identification variables is also considered in the propagation of the uncertainty. The resulting uncertainty in the cross section is ≈10% (≈4%) at low (high) E γ T .
• Photon identification and isolation correlation in the background. The photon isolation and identification variables used to define the two-dimensional background sidebands are assumed to be independent in jet background events (R bg = 1.0).
Any correlation between these variables affects the estimation of signal purity and leads to systematic uncertainties in the background-subtraction procedure. The value of R bg is estimated using background MC samples and is found to be consistent with unity within ±10% [3,41]. This value of R bg is verified using background-enriched regions in data. The assumption of R bg = 1.0 is found to hold within ±10% in the kinematic region of the measurements presented here. The resulting uncertainty in the cross section is ≈8% (≈4%) at low • Photon isolation modelling. Differences between data and signal MC events in the modelling of the isolation distribution can lead to systematic uncertainties in the estimation of the non-isolated signal leakage fractions and the C (i) correction factor. Two subsamples are selected from data by applying either the tight or non-tight identification criteria to each photon; the subsample selected with non-tight identification criteria is expected to be enriched in background candidates. The E iso T value for the non-tight candidates is scaled so that the integral for E iso T > 10 GeV, where the contribution from the signal is expected to be negligible, matches that of the tight candidates. The rescaled background distribution is subtracted from that of the tight photon candidates to extract the isolation profile of signal-like candidates. These distributions are used to derive Smirnov transformations [36]. The Smirnov transformation shifts the photon isolation values event-byevent in MC simulation to match the isolation distribution found in data. This Smirnov-transformed MC sample is used to estimate new differential cross sections. Differences from the nominal results are taken as systematic uncertainties. The resulting uncertainty in the cross section is ≈7% (≈4%) at low • Photon fragmentation contribution. The admixture of direct and fragmentation photons affects the estimation of the signal leakage fractions which are used in the jet background subtraction procedure and the C (i) correction factor. A photon originating from the fragmentation of a parton can be modelled in the MC simulation by allowing the radiation of a photon by a parton. A sample of fragmentation photons is selected by applying the event selection to a diphoton MC sample (see Section 3.1). This selects three-photon events where at least one of the final-state photons results from fragmentation. The diphoton MC sample predicts that for more than 98% of the events the sub-sub-leading photon originates from parton bremsstrahlung. Differences in the isolation distributions between direct and fragmentation photons are expected. Therefore, a template fit to the sub-sub-leading photon isolation distribution is performed to determine the optimal admixture of the nominal and diphoton MC samples. The direct and fragmentation isolation templates are given by the nominal and diphoton MC samples respectively, whereas the jet background template is taken from a data control region where the subsub-leading photon candidate satisfies the non-tight selection. The fit estimates that about 40% of the sub-sub-leading photons originate from fragmentation, as modelled by the diphoton MC sample. This value is used to merge the nominal and diphoton MC samples. The new MC sample is used to estimate the signal leakage fractions and the C (i) correction factors.
The deviation of the differential cross section from the value obtained using the Smirnov-transformed MC sample is taken as the systematic uncertainty. This avoids double counting the effect of the photon isolation modelling. The resulting uncertainty in the cross section is ≈ 4%. • e-γ misidentification. The uncertainty in the electron background contamination is estimated by propagating the uncertainty in the e-γ scale factors (see Section 5.1), which affects the prediction of the e-γ misidentification rates, to the estimation of the cross section. The resulting uncertainty is ≈0.1%. • Unfolding procedure. The effect of unfolding is investigated by using smooth functions to re-weight the signal MC simulation to match the data distributions after background subtraction. The data are unfolded using this reweighted MC sample and the resulting cross sections are compared to the nominal measurements. The differential cross sections are found to differ by <1%.  Table 2 shows the breakdown of the systematic uncertainties in the measurement of the inclusive fiducial cross section. The statistical uncertainty in the measured inclusive fiducial cross section is ≈9%.

Theoretical uncertainties
The following sources of uncertainty in the theoretical predictions are considered for the MCFM and MadGraph5_aMC@NLO calculations.
• The uncertainty in the NLO QCD calculations due to terms beyond NLO is estimated by repeating the calculations using values of μ R , μ F and μ f scaled by factors 0.5 and 2. For the MadGraph5_aMC@NLO calculations, only the μ R and μ F scales are varied. In addition, the scales are either varied simultaneously, individually or by fixing one and varying the other two. The final uncertainty is taken as the largest deviation of the possible variations with respect to the nominal value.
• The uncertainty in the NLO QCD calculations due to uncertainties in the proton PDFs is estimated by repeating the calculations using the 52 additional sets from the CT10 error analysis [20].
• The uncertainty in the NLO QCD calculations due to the value of α S (m Z ) = 0.118 is estimated by repeating the calculations The dominant theoretical uncertainty in the predicted cross section arises from the missing terms beyond NLO and amounts to 10-12%. The uncertainty arising from the PDF variations amounts to 2-3% and the uncertainty arising from the value of α S (m Z ) is below 2%. The total theoretical uncertainty is obtained by adding in quadrature the individual uncertainties listed above and amounts to 10-13%.
The NLO QCD calculations underestimate the measured inclusive fiducial cross section by factors of 2.3 and 1.6 for MCFM and MadGraph5_aMC@NLO, respectively. The addition of the parton shower to the MadGraph5_aMC@NLO prediction improves the agreement with the measured value. The NLO electroweak corrections are small and cannot account for the observed differences between NLO QCD and the measurements [42]. Similar discrepancies between the NLO calculations and the measurements are found for the prediction of the inclusive fiducial cross section for γ γ , W γ γ and Z γ γ production [5,43,44]. The NNLO calculations, which are available for the computation of γ γ but not for γ γ γ production, significantly improve the description of the diphoton fiducial cross section [6,45]. Fig. 1 shows the three-isolated-photons differential cross sections as functions of E   Fig. 2. The MCFM calculations underestimate the measurements by 50% in the low invariant mass regions, whereas the differences are 30-40% for m γ 1 γ 2 150 GeV, m γ 1 γ 3 75 GeV, m γ 2 γ 3 75 GeV and m γ γ γ 150 GeV. The MadGraph5_aMC@NLO calculations also underestimate the data by 30-50% in the low invariant mass regions. However, they tend to give a better description of the measurements for m γ 1 γ 2 150 GeV, m γ 1 γ 3 75 GeV, m γ 2 γ 3 75 GeV and m γ γ γ 150 GeV. For such regions, MadGraph5_aMC@NLO predictions are 25-30% higher than the MCFM estimates. Fig. 3 shows the three-isolated-photons differential cross sections as functions of φ γ 1 γ 2 , φ γ 1 γ 3 , φ γ 2 γ 3 , | η γ 1 γ 2 |, | η γ 1 γ 3 | and | η γ 2 γ 3 |. The theoretical calculations underestimate the normalisation of the measurements. This is due to the fact that these distributions are mainly populated by low-E γ T photons. Both NLO QCD calculations give an adequate description of the shape of the differential cross sections as functions of | η γ 1 γ 2 |, | η γ 1 γ 3 | and | η γ 2 γ 3 |. A quantitative comparison of the NLO QCD predictions to the measurements as functions of φ γ 1 γ 2 , φ γ 1 γ 3 and φ γ 2 γ 3 is performed with a χ 2 fit to the cross-section normalisation including both statistical and systematic uncertainties. This tests the description of the shape of the differential cross sections. The total systematic uncertainty is considered to be fully correlated across bins and is included in the χ 2 definition using nuisance parameters. After the χ 2 minimisation, scale factors equal to ≈ 1.6 (MadGraph5_aMC@NLO) and ≈ 2.3 (MCFM) are found for each angular distribution independently. Both theoretical predictions give an adequate description of the shape of dσ /d φ γ 2 γ 3 (χ 2 /ndof = 6/5 and 7/5 for MadGraph5_aMC@NLO and MCFM, respectively, where ndof is the number of degree of freedom). In addition, MadGraph5_aMC@NLO calculations describe adequately the shape of dσ /d φ γ 1 γ 2 and dσ /d φ γ 1 γ 3 (χ 2 /ndof = 6/5 and 7/5, respectively) but not MCFM (χ 2 /ndof = 13/5 and 14/5, respectively). This shows the importance of the addition of the parton shower to improve the description of the shape of dσ /d φ γ 1 γ 2 and dσ /d φ γ 1 γ 3 .

Summary
A measurement of the production cross section of three isolated photons in pp collisions at √ s = 8 TeV with the ATLAS detector at the LHC is presented using a data set with an integrated luminosity of 20.2 fb −1 . Differential cross sections as functions of E T > 15 GeV, m γ γ γ > 50 GeV, and |η γ | < 2.37, excluding the region 1.37 < |η γ | < 1.56. The distance between pairs of photons in the η-φ plane is required to be R > 0.45. The selection of isolated photons is ensured by requiring that the transverse energy in a cone of size R = 0.4 around the photon is smaller than 10 GeV. The inclusive fiducial cross section is measured to be σ meas = 72.6 ± 6.5 (stat.) ± 9.2 (syst.) fb 150 GeV, m γ 1 γ 3 75 GeV, m γ 2 γ 3 75 GeV and m γ γ γ 150 GeV. Both NLO calculations give an adequate description of the shape of the measured cross section as functions of | η γ 1 γ 2 |, | η γ 1 γ 3 | and | η γ 2 γ 3 |, whereas they underestimate the normalisation of the measurements. In addition, both theoretical predictions inadequately describe the normalisation of the measurements as functions of φ γ 1 γ 2 , φ γ 1 γ 3 and φ γ 2 γ 3 . MCFM predictions give an adequate description of the shape of dσ /d φ γ 2 γ 3 and fail to describe the shape of dσ /d φ γ 1 γ 2 and dσ /d φ γ 1 γ 3 , whereas MadGraph5_aMC@NLO predictions give an adequate description of the shape of the measured cross sections as functions of all three angular variables. The measurements provide a test of pQCD for the description of the dynamics of triphoton production and indicate the need for improved modelling of this process in MC models.