A study of final-state radiation in decays of Z bosons produced in pp collisions at 7 TeV

The differential cross sections for the production of photons in Z to mu+ mu- gamma decays are presented as a function of the transverse energy of the photon and its separation from the nearest muon. The data for these measurements were collected with the CMS detector and correspond to an integrated luminosity of 4.7 inverse femtobarns of pp collisions at sqrt(s) = 7 TeV delivered by the CERN LHC. The cross sections are compared to simulations with POWHEG and PYTHIA, where PYTHIA is used to simulate parton showers and final-state photons. These simulations match the data to better than 5%.


Introduction
We present a study and differential cross section measurements of photons emitted in decays of Z bosons produced at a hadron collider. Such radiative decays of the Z boson were noted in the very first Z boson publications of UA1 and UA2 [1,2], but subsequently have not been given a detailed study in hadron colliders. In 2011, the CERN LHC delivered pp collisions at √ s = 7 TeV, and data corresponding to an integrated luminosity of 4.7 fb −1 were collected with the CMS detector. From these data, we select a sample of events in which a Z boson decays to a µ + µ − pair and an energetic photon. We measure the differential cross sections dσ/dE T with respect to the photon transverse energy E T and dσ/d∆R µγ with respect to the separation of the photon from the nearest muon. Here, ∆R µγ = (φ µ − φ γ ) 2 + (η µ − η γ ) 2 , where φ is the azimuthal angle (in radians) around the beam axis and η is the pseudorapidity. The cross sections include contributions from the Z resonance, virtual photon exchange, and their interference, collectively referred to as Drell-Yan (DY) production.
Photons emitted in Z boson decays, which we call final state radiation (FSR) photons, can be energetic (tens of GeV) and well separated from the leptons (by more than a radian). Quantum electrodynamics (QED) corrections that describe FSR production are well understood. Quantum chromodynamics (QCD) corrections modify the kinematic distributions of the Z boson; in particular, the Z boson aquires a nonzero component of momentum transverse to the beam: q T > 0. The FEWZ program calculates both QCD and QED corrections for the DY process [3]. However, it does not include mixed QCD and QED corrections; the required two-loop integrals are technically very challenging, and progress has been made only recently [4]. In practice, event generators employing matrix element calculations matched to parton showers must be used [5,6]. It is the goal of this analysis to establish the quality of the modeling of FSR over a wide kinematic and angular range. The results will support future measurements of the W mass, the study of Z +γ production, and searches for new particles in final states with photons.
In an attempt to compare photons emitted close to a muon (a process that is modeled primarily by a partonic photon shower) and far from the muons (which requires a matrix element calculation), we measure dσ/dE T for 0.05 < ∆R µγ ≤ 0.5 and 0.5 < ∆R µγ ≤ 3. Furthermore, since the size of the QCD corrections varies with the transverse momentum of the Z boson, we measure dσ/dE T and dσ/d∆R µγ for q T < 10 GeV and q T > 50 GeV, where q T is defined as the vector sum of the transverse momenta of the two muons and the photon. These cross sections are defined with respect to the fiducial and kinematic requirements detailed below; no acceptance corrections are applied. Nonetheless, we do correct for detector resolution and efficiencies.
This article is structured as follows. We briefly describe the CMS detector and the event samples we use in Section 2. The details of the event selection are given in Section 3. Background estimation and the way we unfold the data distributions are discussed in Sections 4 and 5. We discuss the systematic uncertainties in Section 6 and report our results and summarize the work in Sections 7 and 8.

The CMS detector and event samples
A full description of the CMS detector can be found in Ref. [7]; here we briefly describe the components most important for this analysis. The central feature of the CMS experiment is a superconducting solenoid that provides an axial magnetic field of 3.8 T. A tracking system formed to the two muon tracks, and the χ 2 probability of the fit must be at least 0.02. We define the difference between π and the opening angle of the two muons as the acollinearity α, and remove a very small region of phase space where α is less than 5 mrad to reduce contamination by cosmic rays to a negligible level.
Photons are reconstructed using the particle flow (PF) algorithm [18,19] that uses clustered energy deposits in ECAL. The PF algorithm allows us to reconstruct photons at relatively low E T and to maintain coherence with the calculation of the isolation observables described below. Photons that convert to electron-positron pairs are included in this reconstruction. Events selected for this analysis must have at least one photon with E T > 5 GeV, and the separation of this photon with respect to the closest muon must satisfy 0.05 < ∆R µγ ≤ 3. Studies using the simulation show that photons with ∆R µγ < 0.05 are difficult to reconstruct reliably due to the energy deposition left by the muon, and no signal photons with ∆R µγ > 3 are expected. If an event has more than one photon satisfying this selection criteria, we select the one with the highest E T . In events in which one photon is selected, a second photon is present 15% of the time; however, these extra photons are expected to be mostly background, since the fraction of events with a second FSR photon in simulation is approximately 0.5%. More details about these background photons are given in Section 4. All three particles emitted in the Z boson decay-the two muons and the photon-are usually isolated from other particles produced in the same bunch crossing. We can reduce backgrounds substantially by imposing appropriate isolation requirements. The isolation quantities, I µ for the muons and I γ for the photon, are based on the scalar p T sums of reconstructed PF particles within a cone around the muon or photon direction. The cone size for both muons and photons is ∆R = 0.3. The muon p T is not included in the sum for I µ , and the photon E T is not included in the sum for I γ ; these isolation quantities are meant to represent the energy carried by particles originating from the main primary vertex close to the given muon or photon. For a wellisolated muon or photon, I µ or I γ should be small. Special care is taken to avoid inefficiencies and biases occurring when the FSR photon falls close to the muon; in such cases the muon and the photon may appear, superficially, to be nonisolated. To avoid this effect, we exclude any PF photon from the muon isolation sum. Furthermore, since the photon can convert and produce charged particle tracks that cannot always be unambiguously identified as an e + e − pair, we exclude from the muon isolation sum charged tracks that lack hits in the pixel detector or that have p T < 0.5 GeV. Finally, any particle that lands in a cone of ∆R < 0.2 around a PF photon is excluded from the muon isolation sum. After these modifications to the muon isolation variable, the efficiency of the isolation requirement is flat (98%) for all ∆R µγ and is higher than the efficiency of the unmodified isolation requirement by about 0.5%. Adding these modifications does not significantly increase the backgrounds.

Background estimation
The instantaneous luminosity of the LHC was sufficiently high that each bunch crossing resulted in multiple pp collisions (8.2 on average). The extraneous pp collisions are referred to as "pileup" and must be taken into consideration when defining and calculating the muon and photon isolation variables. Charged hadrons, electrons, and muons coming from pileup can be identified by checking their point of origin along the beam line, which will typically not coincide with the primary vertex from which the muons originate. When summing the contributions of charged hadrons, electrons, and muons to the isolation variable, those coming from pileup are excluded. This distinction is not possible for photons and neutral hadrons, however. Instead, an estimate I p of the contribution of photons and neutral hadrons to the sum is made: we use one-half of the (already excluded) contribution from charged hadrons within the isolation cone. This estimate is subtracted from the sum of contributions from photons and neutral hadrons; if the result is negative, we then use a net contribution of zero.
We designate by I h ± the sum of p T for charged particles that are not excluded from the isolation sum. We let I em and I h 0 stand for the sums over the p T of all photons and neutral hadrons, and I p for the estimate of the pileup contribution to I em and I h 0 . The muon isolation variable is, then, Note that the sum is normalized to the p T of the muon. We require I µ < 0.2 for both muons.
The photon isolation variable is calculated as above, except that the muons are not included in the sum, and there is no special exclusion of charged tracks near the photon: We require I γ < 6 GeV.
The emission of FSR photons in Z boson decays reduces the momenta of the muons. Consequently, the dimuon mass M µµ tends to be lower than M Z , the nominal mass of the Z boson. Simulations indicate that, for most of the signal, M µµ < 87 GeV, due to the requirement E T > 5 GeV for the photon. They also show that the M µµ distribution for radiative decays Z → µ + µ − γ ends at M µµ ≈ 30 GeV. A requirement M µµ > 30 GeV also helps to avoid a kinematic region in which the acceptance is difficult to model. Therefore, our signal region is defined by 30 < M µµ < 87 GeV. We also define a control region by 89 < M µµ < 100 GeV, where the contribution of FSR photons is quite small (below 0.5%). The numbers of events we select are 56 005 in the signal region and 45 277 in the control region.

Background estimation
Nearly all selected events have two prompt muons from the DY process. Backgrounds come mainly from "nonprompt" photons, which may be genuine photons produced in the decays of light mesons (such as π 0 and η), a pair of overlapping photons that cannot be distinguished from a single photon, and photons from pileup. We study these backgrounds with simulated DY events and apply corrections so that the simulation reproduces the data distributions, as described in detail below.
Some events come from processes other than DY, such as tt and diboson production. These backgrounds are very small and are estimated using the simulation. Similarly, a small background from the DY production of τ + τ − is also estimated from simulation. The background from multijet events, including events with a W ± → µ ± ν decay, is estimated using events with same-sign muons. Backgrounds from simultaneous nonprompt muon and nonprompt photon sources are negligible. The composition of the signal sample is given in Table 2.
The control region is dominated by nonprompt photons whose kinematic distributions (E T , η, ∆R µγ ) are nearly identical to nonprompt photons in the signal region. Quantitative comparisons of data and simulation revealed significant discrepancies in the control region that prompted corrections to the simulation, which we now explain.
The POWHEG+PYTHIA sample does not reproduce the number of jets per event well, so we apply weights to the simulated events as a function of the number of reconstructed jets in each event. For this purpose, jets are reconstructed from PF objects using the anti-k T algorithm [20] with a size parameter R = 0.5. We consider jets with p T > 20 GeV and |η| < 2.4 that do not overlap with the muons or the photon.
Studies of I γ for events in the control region reveal small discrepancies in the multiplicity and p T spectra of charged hadrons included in the sum. We apply weights to the simulated events to bring the multiplicities into agreement, and we impose p T > 0.5 GeV on charged hadrons. The simulated I γ distributions match those in data very well after applying these weights.
Finally, the E T and η distributions of nonprompt photons in the simulation deviate from those in the data. We fit simple analytic functions to the ratios of data to simulated E T and η distributions and define a weight as the product of those functions. We check that this factorization is valid (i.e., that the E T correction is the same for different narrow ranges of η, and vice versa).
After these three corrections (for jet multiplicity, for the spectrum of charged hadrons in the isolation sums, and for the E T and η of the nonprompt photon), the simulation matches the data in all kinematic distributions in the control region, an example of which is shown in Fig. 1, left. The total change in the background estimate due to these corrections is approximately 5 to 10%. We use the simulation with these weights to model the small background in the signal region ( Fig. 1, right).
Given the definition of the signal region, the contribution of photons emitted in the initial state is very small (on the order of 4 × 10 −4 as determined from the POWHEG+PYTHIA sample) and is counted as signal.

Correcting for detector effects
Our goal is to measure differential cross sections in a form that is optimal for testing FSR calculations. Therefore we are obliged to remove the effects of detector resolution and efficiency (including reconstruction, isolation, and trigger efficiency). The corrections for the muons follow the techniques developed for the DY cross section measurements [17]. The corrections for photons are applied using an unfolding technique, as discussed in this section.
We apply small corrections to the muon momentum scale as a function of muon p T , η µ , and φ µ [21]; they have almost no impact on our measurements. The muon reconstruction efficiency, (taken from simulation and corrected to match the data as a function of p T and η µ ) is taken into account by applying weights on a per-event basis. We do not correct for the approximately 0.5% increase in the isolation efficiency coming from the way we handle FSR photons falling within the muon isolation cone.
The energy scale and efficiencies for photons are more central to our task. Most PF photon energies correspond to the true energies within a few percent. However, in about 13% of the cases the photon energy is significantly underestimated. The simulation reproduces this effect very well. We construct a "response" matrix that relates the PF energy to the true energy as a function of η γ and ∆R µγ . The angular quantities η γ and ∆R µγ are themselves well measured. We use the iterative D'Agostini method of unfolding [22] as implemented in the ROOUNFOLD package [23]. By default, we unfold in the three quantities E T , η γ , and ∆R µγ simultaneously after subtracting backgrounds; as a cross check we also unfold the E T and ∆R µγ distributions one at a time, and we also use a single-value decomposition method [24]. All results are consistent with each other. To verify the independence of the unfolded result on the assumed spectra, we distort the FSR model in an arbitrary manner when reconstructing the response matrix. The unfolded result is no different than the original one we obtained. A closure test in which the simulation is treated as data and undergoes the same unfolding procedure indicates no deviation greater than 1.5%.
The unfolding procedure corrects for the photon reconstruction and isolation efficiency. It also corrects for bias in the PF photon energy assuming that such a bias is reproduced in the simulation. Verification of the photon efficiencies and energy scale in data with respect to the simulation are discussed in Section 6.

Systematic uncertainties
Systematic uncertainties are assigned to each step of the analysis procedure using methods detailed in this section. Tables 3 and 4 present a summary of these uncertainties, which are similar in magnitude to, or somewhat larger than the statistical uncertainties, depending on the photon E T .
The muon efficiency, taken from simulation, is corrected as a function of p T and η using a method derived from data and described in Ref. [17]. The statistical uncertainties for these corrections constitute a systematic uncertainty, which we also take from Ref. [17]. In addition, we assign a 0.5% uncertainty to account for the modifications of the standard isolation variable. We propagate the uncertainty in the muon efficiency by shifting the per-event weights up and down by one unit of systematic uncertainty.
The photon E T scale is potentially an important source of systematic uncertainty although simulations indicate that the bias in PF photon energy is negligibly small. We verify the fidelity of the simulation by introducing an extra requirement, 0.05 < ∆R µγ ≤ 0.9, which gives us a high-purity subset of signal events in which the energy of the photon can be estimated from just the muon kinematics. We refer to this estimate as E kin γ . The quantity Gaussian with a mean close to zero. We conducted detailed quantitative studies of the s distribution in bins of E PF Tγ , separately in the barrel and endcaps. We fit the distributions to a Gaussian-like function in which the width parameter is itself a function of s, namely, σ(s) = c(1 + e bs ), with b and c as free parameters. Examples are given in Fig. 2. Overall, the simulation reproduces the s distributions in data very well. We derive some small corrections from the differences in data and simulation as a function of E PF Tγ and construct an alternate response matrix. The unfolded spectrum we obtained with this alternate response matrix differs from the original by less than 0.2% for E T < 40 GeV, by less than 1% for E T < 75 GeV, and by less than 3% in the highest E T bin. We assign respective systematic uncertainties of 0.5%, 1%, and 3% for these three E T ranges to account for the photon energy scale uncertainty.
fit with a skewed Gaussian as described in the text. The left and right plots pertain to photons in the ECAL barrel with 5 < E T < 10 GeV and in the ECAL endcaps with 20 < E T < 40 GeV, respectively. The circle points and solid curve represent the data and the triangle points and dotted curve represent the simulation.
The photon energy resolution uncertainty is well constrained by studies with electrons and FSR photons [9]. To assess the impact of the uncertainty in the resolution, we degrade the photon energy resolution in simulated events by adding in quadrature a 1% term to the nominal resolution and construct a new response matrix. The differences in the unfolded spectrum relative to the default response matrix are small, and we take these differences as the systematic  uncertainty due to photon energy resolution.
Efficiency corrections for photons are applied as part of the unfolding procedure described in Section 5 and are derived from the simulation. We verify these corrections using the data in the following way. An isolated FSR photon with E T > 5 GeV nearly always produces a cluster in the ECAL. We define an efficiency to reconstruct and select PF photons given such isolated clusters. This efficiency rises from 60% for E T between 5-10 GeV to approximately 90% for E T > 50 GeV and is nearly the same in data and simulation. We take the difference, added in quadrature to the statistical uncertainties of the efficiencies, as the systematic uncertainty.
As described briefly in Section 5, the unfolding procedure has been cross-checked in several ways. To assess a systematic uncertainty due to unfolding, we use the small discrepancies observed in the closure test.
The uncertainty in the background estimate is dominated by the uncertainties associated with the corrections that we obtained from the control region (Section 4). The statistical uncertainty in the weights for jet multiplicity has a negligible impact, as does the correction for charged hadrons in the photon isolation cone. The parameterized functions to correct the photon distributions in E T and η carry statistical uncertainties that we propagate to the measured cross sections through simplified MC models. Since the nonprompt photon E T , η, and ∆R µγ distributions in the control and signal regions are indistinguishable, we do not assess any uncertainty in the modeling.
The uncertainties in the non-DY backgrounds (tt and diboson production) are obtained from the uncertainties in the theoretical cross sections, the luminosity, and the statistical uncertainty in the simulated event samples. We assign 50% uncertainty to the W+jets and multijet background estimates, which are quite small.
The systematic uncertainty from the simulation of pileup depends primarily on the assumed cross section for additional pp collisions (roughly the same as the minimum-bias cross section) [25]. We vary the value of this cross section by 5% and evaluate the impact on the unfolded spectra.
The uncertainty in the integrated luminosity is 2.2% [26].
Theoretical uncertainties have been calculated and pertain to the reported theoretical prediction only. We propagated the uncertainty due to parton distribution functions (PDFs) using the prescription of Ref. [27]. We vary the factorization/renormalization scale parameters by a factor of 2 to estimate associated scale uncertainties introduced due to neglected higher-order quantum corrections. Finally, we include the MC statistical uncertainty.

Results
The differential cross sections are obtained by subtracting the estimated backgrounds from the observed distributions, unfolding the result, and dividing by the bin width and the integrated luminosity, L = 4.7 fb −1 . No acceptance correction is applied, so these cross sections are defined relative to the kinematic and fiducial requirements listed in Table 1.
The measured differential cross sections dσ/dE T and dσ/d∆R µγ are displayed in Fig. 3 and listed in Tables 5 and 6. A bin-centering correction is applied following the method of Ref. [28]; the abscissa of each point is based on the integral of the simulation across the bin and on the bin width. The shaded region represents the prediction and uncertainty from POWHEG+PYTHIA, obtained at the parton level: only the requirements in Table 1 have been applied to the generator-level muons and photons. The agreement with data is good.
Energy spectra for photons closer to (0.05 < ∆R µγ ≤ 0.5) and farther from the muon (0.5 < ∆R µγ ≤ 3) are shown in Fig. 4. The rates for photons with large ∆R µγ and E T are also well reproduced. The number of events with 30 < M µµ < 87 GeV is about 18% of the number with 60 < M µµ < 120 GeV. Of the events with 30 < M µµ < 87 GeV, the fraction of events with at least one photon with E T > 5 GeV and 0.05 < ∆R µγ ≤ 0.5 is 8.7 ± 0.1 (stat) ± 0.2 (syst)%, and with 0.5 < ∆R µγ ≤ 3 is 5.6 ± 0.1 (stat) ± 0.2 (syst)%. Photons with ∆R µγ > 1.2 and E T > 40 GeV constitute a small fraction (1.3 ± 0.5 (stat) ± 0.6 (syst)) × 10 −4 .  In the upper panels, the dots with error bars represent the data, and the shaded bands represent the POWHEG+PYTHIA calculation including theoretical uncertainties. The central panels display the ratio of data to the MC expectation. The lower panels show the standard deviations of the measurements with respect to the calculation. A bin-centering procedure has been applied.
We define two subsamples of signal events, one with the Z boson transverse momentum q T < 10 GeV, and the other with q T > 50 GeV. The measured cross sections shown in Fig. 5 demonstrate rather different energy spectra for these two cases, though dσ/d∆R µγ is basically the same.
As a final illustration of the nature of this event sample, we present distributions of dimuon mass (M µµ ) and the three-body mass (M µµγ ) in Fig. 6. The small increase in the ratio of data to theory for M µµ < 40 GeV reflects the insufficient next-to-leading-order accuracy of the simulation; the kinematic requirements on the muons induce a loss of acceptance that require higherorder QCD corrections, as discussed in Ref. [17]. Although the masses of the dimuon pairs populate the tail of the Z resonance (in fact they were selected this way), the three-body mass distribution displays a nearly-symmetric resonance peak at the mass of the Z boson, thereby confirming the identity of these events as radiative decays Z → µ + µ − γ.     Figure 4: Measured differential cross sections dσ/dE T for photons close to the muon (0.05 < ∆R µγ ≤ 0.5, left) and far from the muon (0.5 < ∆R µγ ≤ 3, right). The dots with error bars represent the data, and the shaded bands represent the POWHEG+PYTHIA calculation including theoretical uncertainties. The central panels display the ratio of data to the MC expectation. The lower panels show the standard deviations of the measurements with respect to the calculation. A bin-centering procedure has been applied.

Summary
A study of final-state radiation in Z boson decays was presented. This study serves to test the simulation of events where mixed QED and QCD corrections are important. The analysis was performed on a sample of pp collision data at √ s = 7 TeV recorded in 2011 with the CMS detector and corresponding to an integrated luminosity of 4.7 fb −1 . Events with two oppositely charged muons and an energetic, isolated photon were selected with only modest backgrounds. The differential cross sections dσ/dE T and dσ/d∆R µγ were measured for photons within the fiducial and kinematic requirements specified in Table 1, and comparisons of dσ/dE T for photons close to a muon and far from both muons were made. In addition, the differential cross sections dσ/dE T and dσ/d∆R µγ were compared for events with large and small transverse momentum of the Z boson, as computed from the two muons and the photon. Simulations based on POWHEG+PYTHIA reproduce the CMS data well, with discrepancies below 5% for 5 < E T < 50 GeV and 0.05 < ∆R µγ ≤ 2 as quantified in Tables 5 and 6.   Figure 5: Measured differential cross sections dσ/dE T and dσ/d∆R µγ for q T < 10 GeV (top row) and q T > 50 GeV (bottom row). The dots with error bars represent the data, and the shaded bands represent the POWHEG+PYTHIA calculation including theoretical uncertainties. The central panels display the ratio of data to the MC expectation. The lower panels show the standard deviations of the measurements with respect to the calculation. A bin-centering procedure has been applied. Distributions of the dimuon mass M µµ (left) and the three-body mass M µµγ (right). The dots with error bars represent the data, and the shaded bands represent the POWHEG+PYTHIA prediction. The central panels display the ratio of data to the MC expectation. The lower panels show the standard deviations of the measurements with respect to the calculation. A bin-centering procedure has been applied.