UvA-DARE Measurement of the distributions of event-by-event flow harmonics in lead-lead collisions at √sNN = 2.76 TeV with the ATLAS detector at the LHC

: The distributions of event-by-event harmonic ﬂow coeﬃcients v n for n = 2– 4 are measured in √ s NN = 2 . 76 TeV Pb+Pb collisions using the ATLAS detector at the LHC. The measurements are performed using charged particles with transverse momentum p T > 0 . 5 GeV and in the pseudorapidity range | η | < 2 . 5 in a dataset of approximately 7 µ b − 1 recorded in 2010. The shapes of the v n distributions suggest that the associated ﬂow vectors are described by a two-dimensional Gaussian function in central collisions for v 2 and over most of the measured centrality range for v 3 and v 4 . Signiﬁcant deviations from this function are observed for v 2 in mid-central and peripheral collisions, and a small deviation is observed for v 3 in mid-central collisions. In order to be sensitive to these deviations, it is shown that the commonly used multi-particle cumulants, involving four particles or more, need to be measured with a precision better than a few percent. The v n distributions are also measured independently for charged particles with 0 . 5 < p T < 1 GeV and p T > 1 GeV. When these distributions are rescaled to the same mean values, the adjusted shapes are found to be nearly the same for these two p T ranges. The v n distributions are compared with the eccentricity distributions from two models for the initial collision geometry: a Glauber model and a model that includes corrections to the initial geometry due to gluon saturation eﬀects. Both models fail to describe the experimental data consistently over most of the measured centrality range.


Introduction
Heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) create hot, dense matter that is thought to be composed of strongly interacting quarks and gluons. A useful tool to study the properties of this matter is the azimuthal anisotropy of particle emission in the transverse plane [1,2]. This anisotropy has been interpreted as a result of pressure-driven anisotropic expansion (referred to as "flow") of the created matter, and is described by a Fourier expansion of the particle distribution in azimuthal angle φ, around the beam direction: where v n and Φ n represent the magnitude and phase of the n th -order anisotropy of a given event in the momentum space. These quantities can also be conveniently represented by the per-particle "flow vector" [2]: ⇀ v n = (v n cos nΦ n , v n sin nΦ n ). The angles Φ n are commonly referred to as the event plane (EP) angles.
In typical non-central [2] heavy ion collisions, the large and dominating v 2 coefficient is associated mainly with the "elliptic" shape of the nuclear overlap. However, v 2 in -1 -

JHEP11(2013)183
central (head-on) collisions and the other v n coefficients in general are related to various shape components of the initial state arising from fluctuations of the nucleon positions in the overlap region [3]. The amplitudes of these shape components, characterized by eccentricities ǫ n , can be estimated via a simple Glauber model from the transverse positions (r, φ) of the participating nucleons relative to their centre of mass [4,5]: ǫ n = r n cos nφ 2 + r n sin nφ 2 r n . (1. 2) The large pressure gradients and ensuing hydrodynamic evolution can convert these shape components into v n coefficients in momentum space. Calculations based on viscous hydrodynamics suggest that v n scales nearly linearly with ǫ n , for n < 4 [6]. The proportionality constant is found to be sensitive to properties of the matter such as the equation of state and the ratio of shear viscosity to entropy density [7,8]. In particular, the proportionality constant is predicted to decrease quickly with increasing shear viscosity [9]. Hence detailed measurements of v n coefficients and comparisons with ǫ n may shed light on the collision geometry of the initial state and transport properties of the created matter [10,11].
Significant v n coefficients have been observed for n ≤ 6 at RHIC and the LHC [12][13][14][15][16][17][18]. These observations are consistent with small values for the ratio of shear viscosity to entropy density, and the existence of sizable fluctuations in the initial state. Most of these measurements estimate v n from the distribution of particles relative to the event plane, accumulated over many events. This event-averaged v n mainly reflects the hydrodynamic response of the created matter to the average collision geometry in the initial state. More information, however, can be obtained by measuring ⇀ v n or v n on an event-by-event (EbyE) basis.
Some properties of the v n distributions, such as the mean v n , the standard deviation (hereafter referred to as "width") σ vn , the relative fluctuation σ vn / v n , and the root-meansquare v 2 n ≡ v n 2 + σ 2 vn , were previously estimated from a Monte Carlo template fit method [19], or two-and four-particle cumulant methods [20][21][22]. The value of σ v 2 / v 2 was measured to be 0.3-0.7 in different centrality bins. However, these methods are reliable only for σ vn ≪ v n , and are subject to significant systematic uncertainties. In contrast, v n , σ vn and higher-order moments can be calculated directly from the full v n distributions.
The EbyE distributions of ⇀ v n or v n also provide direct insight into the fluctuations in the initial geometry [23]. If fluctuations of ⇀ v n relative to the underlying flow vector associated with the average geometry, ⇀ v RP n , in the reaction plane 1 (RP) [23,24] are described by a two-dimensional (2D) Gaussian function in the transverse plane, then the probability density of ⇀ v n can be expressed as: Model calculations show that this approximation works well for central and mid-central collisions [23,25]. Integration of this 2D Gaussian over the azimuthal angle gives the 1 The reaction plane is defined by the impact parameter vector and the beam axis.

JHEP11(2013)183
one-dimensional (1D) probability density of v n = | ⇀ v n | in the form of the Bessel-Gaussian function [7,26]: (1.4) where I 0 is the modified Bessel function of the first kind. Additional smearing to eq. (1.3) also arises from effects of the finite number of particles produced in the collision. If it is Gaussian, this smearing is expected to increase the observed δ vn value, but the value of v RP n should be stable. The parameters v RP n and δ vn in eq. (1.4) are related to v n and σ vn , and can be estimated directly from a fit of the measured p(v n ) distribution with eq. (1.4). For small fluctuations δ vn ≪ v RP n [23]: which is equivalent to the "fluctuation-only" scenario, i.e. v RP n = 0. In this case, both the mean and the width are controlled by δ vn [27]: and hence: σ vn v n = 4 π − 1 = 0.523, v 2 n = 2 √ π v n = 1.13 v n .
(1. 8) In the intermediate case, δ vn ≈ v RP n , a more general approximation to eq. (1.4) can be used via a Taylor expansion of the Bessel function, I 0 (x) = e x 2 /4 1 − x 4 /64 + O(x 6 ) : Defining α ≡ δ vn /v RP n , eqs. (1.9) and (1.10) imply that for v n ≪ 2 √ 2δ vn α, the shape of the distribution is very close to that of eq. 1.6, except for a redefinition of the width. For example, the deviation from the fluctuation-only scenario is expected to be less than 10% over the range v n < 1.6δ vn α. Hence the reliable extraction of v RP n requires precise determination of the tails of the v n distributions, especially when v RP n is smaller than δ vn .
-3 -selection requires a time difference |∆t| < 3 ns between the MBTS trigger counters on either side of the interaction point to suppress non-collision backgrounds. A coincidence between the ZDCs at forward and backward pseudorapidity is required to reject a variety of background processes, while maintaining high efficiency for non-Coulomb processes. Events satisfying these conditions are required to have a reconstructed primary vertex with z vtx within 150 mm of the nominal centre of the ATLAS detector. About 48 million events pass the requirements for the analysis. The Pb+Pb event centrality is characterized using the total transverse energy (ΣE T ) deposited in the FCal over the pseudorapidity range 3.2 < |η| < 4.9 measured at the electromagnetic energy scale [33]. A larger ΣE T value corresponds to a more central collision. From an analysis of the ΣE T distribution after all trigger and event selections, the sampled fraction of the total inelastic cross section is estimated to be (98±2)% [34]. The uncertainty associated with the centrality definition is evaluated by varying the effect of the trigger, event selection and background rejection requirements in the most peripheral FCal ΣE T interval [34]. The FCal ΣE T distribution is divided into a set of 5%-wide percentile bins, together with five 1%-wide bins for the most central 5% of the events. A centrality interval refers to a percentile range, starting at 0% for the most central collisions. Thus the 0-1% centrality interval corresponds to the most central 1% of the events; the 95-100% centrality interval corresponds to the least central (i.e. most peripheral) 5% of the events. A standard Glauber model Monte Carlo analysis is used to estimate the average number of participating nucleons, N part , for each centrality interval [34,35]. These numbers are summarized in table 1.
The v n coefficients are measured using tracks reconstructed in the ID that have p T > 0.5 GeV and |η| < 2.5. To improve the robustness of track reconstruction in the highmultiplicity environment of heavy ion collisions, more stringent requirements on track -5 - quality, compared to those defined for proton-proton collisions [36], are used. At least 9 hits in the silicon detectors (compared to a typical value of 11) are required for each track, with no missing Pixel hits and not more than 1 missing SCT hit, after taking into account the known non-operational modules. In addition, at its point of closest approach the track is required to be within 1 mm of the primary vertex in both the transverse and longitudinal directions [15]. The efficiency, ǫ(p T , η), of the track reconstruction and track selection cuts is evaluated using Pb+Pb Monte Carlo events produced with the HIJING event generator [37]. The generated particles in each event are rotated in azimuthal angle according to the procedure described in ref. [38] to give harmonic flow consistent with previous ATLAS measurements [15,16]. The response of the detector is simulated using GEANT4 [39] and the resulting events are reconstructed with the same algorithms as applied to the data. The absolute efficiency increases with p T by 7% between 0.5 GeV and 0.8 GeV, and varies only weakly for p T > 0.8 GeV. However, the efficiency varies more strongly with η and event multiplicity [40]. For p T > 0.8 GeV, it ranges from 72% at η = 0 to 51% for |η| > 2 in peripheral collisions, while it ranges from 72% at η = 0 to about 42% for |η| > 2 in central collisions. The fractional change of efficiency from most central to most peripheral collisions, when integrated over η and p T , is about 13%. Contributions of fake tracks from random combinations of hits are generally negligible, reaching only 0.1% for |η| < 1 for the highest multiplicity events. This rate increases slightly at large |η|.

Method and data analysis
To illustrate the level of EbyE fluctuations in the data, the top panels of figure 1 show the azimuthal distribution of charged particle tracks with p T > 0.5 GeV for three typical events in the 0-5% centrality interval. The corresponding track-pair ∆φ distributions from the same events are shown in the bottom panels. For each pair of particles two ∆φ entries,  Figure 1. Single-track φ (top) and track-pair ∆φ (bottom) distributions for three typical events (from left to right) in the 0-5% centrality interval. The pair distributions are folded into the [−0.5π, 1.5π] interval. The bars indicate the statistical uncertainties of the foreground distributions, the solid curves indicate a Fourier parameterization including the first six harmonics: dN/dφ = A(1 + 2 6 i=1 c n cos n(φ − Ψ n )) for single-track distributions and dN/d∆φ = A(1 + 2 6 i=1 c n cos n(∆φ)) for track-pair distributions, and the solid points indicate the eventaveraged distributions (arbitrary normalization). distributions shown by the solid points (arbitrary normalization), are observed. These EbyE distributions are the inputs to the EbyE v n analyses.
The azimuthal distribution of charged particles in an event is written as a Fourier series, as in eq. (1.1): v obs n,x cos nφ + v obs n,y sin nφ , , v obs n,x = v obs n cos nΨ obs n = cos nφ , v obs n,y = v obs n sin nΨ obs n = sin nφ , where the averages are over all particles in the event for the required η range. The v obs n is the magnitude of the observed EbyE per-particle flow vector: ⇀ v obs n = (v obs n,x , v obs n,y ). In the limit of very large multiplicity and in the absence of non-flow effects, it approaches the true flow signal: v obs n → v n . The key issue in measuring the EbyE v n is to determine the response function p(v obs n |v n ), which can be used to unfold the smearing effect due to the finite number of detected particles. Possible non-flow effects from short-range correlations, such as resonance decays, Bose-Einstein correlations and jets, also need to be suppressed.
The rest of this section sets out the steps to obtain the unfolded v n distribution. Since the data-driven unfolding technique has rarely been used in the study of flow phenomena, details are provided to facilitate the understanding of the methods and systematic uncertainties. Section 4.1 explains how v obs n and the associated response function can be obtained from the EbyE single-particle distributions, such as those shown in the top panels of figure 1. Section 4.2 describes how v obs n and the response function can be obtained from EbyE two-particle correlations (2PC), similar to those shown in the lower panels of figure 1. In this paper the 2PC approach is used primarily as a consistency check. The Bayesian unfolding procedure, applicable to either the single-particle or 2PC data, is described in section 4.3. The performance of the unfolding is described in section 4.4, while the systematic uncertainties are discussed in section 4.5.

Single-particle method
The azimuthal distribution of particles in figure 1 needs to be corrected for non-uniform detector acceptance. This is achieved by dividing the foreground distribution (S) of a given event by the acceptance function (B) obtained as the φ distribution of all tracks in all events (see top panels of figure 1): where the index i runs over all tracks in an event, ǫ(η, p T ) is the tracking efficiency for a given centrality interval, and v det n,x and v det n,y are Fourier coefficients of the acceptance function in azimuth, also weighted by the inverse of the tracking efficiency. The influence of the structures in the acceptance function can be accounted for by taking the leading-order term of the Taylor expansion of 1/B(φ) in terms of cos nφ and sin nφ: where the values of v det n,x or y are less than 0.007 for n = 2-4. The higher-order corrections to eq. (4.5) involve products of v raw n,x or y and v det n,x or y . They have been estimated and found to have negligible impact on the final v n distributions for n = 2-4. Figure 2 shows the distribution of the EbyE per-particle flow vector ⇀ v obs 2 = (v obs 2,x , v obs 2,y ) and v obs 2 obtained for charged particles with p T > 0.5 GeV in the 20-25% centrality interval. The azimuthal symmetry in the left panel reflects the random orientation of ⇀ v obs 2 because of the random orientation of the impact parameter. Due to the finite track multiplicity, the measured flow vector is expected to be smeared randomly around the true flow vector by a 2D response function p( ⇀ v obs n | ⇀ v n ). In order to determine p( ⇀ v obs n | ⇀ v n ), the tracks in the entire inner detector (referred to as full-ID) for a given event are divided into two subevents with symmetric η range, η > 0 and η < 0 (labelled by a and b and referred to as half-ID). The two half-IDs have the same average track multiplicity to within 1%. The distribution of flow vector differences between  Figure 3. Left: the distribution of the difference between the EbyE per-particle flow vectors of the two half-IDs for events in the 20-25% centrality interval for n = 2. Middle: the x-projection overlaid with a fit to a Gaussian. Right: the y-projection overlaid with a fit to a Gaussian. The width from the fit, δ 2SE , and the quality of the fit, χ 2 /DOF, are also shown. the two subevents, p sub ( ⇀ v obs n ) a − ( ⇀ v obs n ) b , is then obtained and is shown in the left panel of figure 3. The physical flow signal cancels in this distribution such that it contains mainly the effects of statistical smearing and non-flow. The middle and right panels of figure 3 show the x-and y-projections of the distribution, together with fits to a Gaussian function. The fits describe the data very well (χ 2 /DOF≈ 1) across five orders of magnitude with the same widths in both directions, implying that the smearing effects and any effects due to non-flow short-range correlations are purely statistical. This would be the case if either the non-flow effects are small and the smearing is mostly driven by finite particle multiplicity, or the number of sources responsible for non-flow is proportional to the multiplicity and they are not correlated between the subevents [31]. The latter could be true for resonance decays, Bose-Einstein correlations, and jets. Similar behaviour is observed for all harmonics up to centrality interval 50-55%. Beyond that the distributions are found to be described better by the Student's t-distribution, which is a general probability density function for the difference between two estimates of the mean from independent samples. The t-distribution approaches a Gaussian distribution when the number of tracks is large.

JHEP11(2013)183
Denoting the width of these 1D distributions by δ 2SE , the widths of the response functions for the half-ID and the full-ID are δ 2SE / √ 2 and δ 2SE /2, respectively. The response functions themselves can be obtained by rescaling the left panel of figure 3 in both dimensions by constant factors of 2 and √ 2 for the full-ID and half-ID, respectively [31]: This scaling behaviour was found to be valid in a Monte-Carlo study based on the HIJING event generator [31]. Integrated over azimuth, eq. 4.6 reduces to a Bessel-Gaussian function in 1D: The difference between the observed and the true signal, denoted by s = v obs n −v n , accounts for the statistical smearing. The similarity between eq. (4.7) and eq. (1.4) is a direct consequence of the 2D Gaussian smearing. However, the smearing leading to eq. (4.7) is due to the finite charge-particle multiplicity, while the smearing leading to eq. (1.4) is due to the intrinsic flow fluctuations associated with the initial geometry. Hence the smearing in eq. (4.7) is expected to increase the observed δ vn value but the value of v RP n should be relatively stable.
The analytical expression eq. (4.7) can be used to unfold the v obs n distribution, such as that shown in the right panel of figure 2. Alternatively, the measured distribution p(v obs n |v n ), taking into account sample statistics, can be used directly in the unfolding. This measured distribution is obtained by integrating out the azimuthal angle in the 2D response function, and the response function is obtained by rescaling the left panel of figure 3 as described earlier. This approach is more suitable for peripheral collisions where the analytical description using eq. (4.7) is not precise enough.

Two-particle correlation method
The EbyE two-particle correlation (2PC) method starts from the ∆φ information in each event, where ∆φ is calculated for each pair of charged tracks as described at the start of section 4. In order to reduce the effect of short-range correlations in η, the tracks in each pair are taken from different half-IDs. This procedure corresponds to convolving the azimuthal distributions of single particles in the two half-IDs: where A n = cos n∆φ and B n = sin n∆φ . The parameters A n and B n are calculated by averaging over the pairs in each event, with each track weighted by the tracking efficiency, as in eq. (4.3). Due to a large rapidity gap on average between the two particles in each pair, the non-flow effects in eq. 4.8 are naturally suppressed compared with the single particle distribution of eq. 4.1.
An EbyE track-pair variable v obs n,n is subsequently calculated for each event: v obs n,n ≡ A 2 The observed flow signal from the two-particle correlation analysis is then calculated as: where s a = v obsa n − v n and s b = v obs b n − v n are independent variables described by the probability distribution in eq. (4.7) with δ = δ 2SE / √ 2. The response function for v obs,2PC n is very different from that for the single-particle method, but nevertheless can be either calculated analytically via eq. (4.7) or generated from the measured distribution such as that shown in figure 3. For small v n values, the s a s b term dominates eq. (4.10) and the distribution of v obs,2PC n is broader than v obs n . For large v n values, the distributions of s a and s b are approximately described by Gaussian functions and hence: where the fact that the average of two Gaussian random variables is equivalent to a Gaussian with a narrower width has been used, and the smearing of the flow vector for the half-IDs (s a and s b ) is taken to be a factor of √ 2 broader than that for the full-ID (s). Hence the distribution of v obs,2PC n is expected to approach the v obs n distribution of the full-ID when v n ≫ δ 2SE / √ 2.

Unfolding procedure
In this analysis, the standard Bayesian unfolding procedure [41], as implemented in the RooUnfold framework [42], is used to obtain the v n distribution. In this procedure, the true v n distribution ("cause"ĉ) is obtained from the measured v obs n or v obs,2PC n distribution ("effect"ê) and the response function A ji ≡ p(e j |c i ) (i,j are bins) as: where the unfolding matrixM 0 is determined from the response function and some initial estimate of the true distributionĉ 0 (referred to as the prior). The matrixM 0 is used to obtain the unfolded distributionĉ 1 andM 1 , and the process is then iterated. More iterations reduce the dependence on the prior and give results closer to the true distributions but with increased statistical fluctuations. Therefore the number of iterations N iter is -11 -

JHEP11(2013)183
adjusted according to the sample size and binning. The prior can be chosen to be the v obs n distribution from the full-ID for the single-particle unfolding, or the v obs,2PC n distribution obtained by convolving the two half-IDs (eq. (4.8)) for the 2PC unfolding. However, a more realistic prior can be obtained by rescaling the v obs n in each event by a constant factor R to obtain a distribution with a mean that is closer to that of the true distribution: where v obs n and σ obs vn are the mean and the standard deviation of the v obs n distribution, respectively, and v EP n is measured using the FCal event plane method from ref.
[16] with the same dataset and the same track selection criteria. The EP method is known to measure a value between the mean and the root-mean-square of the true v n [25,28] (see figure 13): (4.14) The lower limit is reached when the resolution factor [16] used in the EP method approaches one, and the upper limit is reached when the resolution factor is close to zero. Thus f = 0 (default choice) gives a prior that is slightly broader than the true distribution, f = 1 gives a distribution that has a mean close to the true distribution, and f > 1 typically gives a distribution that is narrower than the true distribution. The unfolding procedure in this analysis has several beneficial features: 1. The response function is obtained entirely from the data using the subevent method described above (eq. (4.6)).
2. Each event provides one entry for the v obs n distribution and the response function (no efficiency loss), and these distributions can be determined with high precision (from about 2.4 million events for each 5% centrality interval).

Unfolding performance
This section describes the unfolding based on the single-particle method and summarizes a series of checks used to verify the robustness of the results: a) the number of iterations used, b) comparison with results obtained from a smaller η range, c) variation of the priors, d) comparison with the results obtained using the 2PC method, and e) estimation of possible biases due to short-range correlations by varying the η gap between the two half-IDs. Only a small subset of the available figures is presented here; a complementary selection can be found in appendix A.
The left and middle panels of figure 4 show the convergence behaviour of the unfolding based on the single-particle method for v 2 in the 20-25% centrality interval measured with the full-ID. Around the peak of the distribution, the results converge to within a few percent of the final answer by N iter = 4, but the convergence is slower in the tails and there are small, systematic improvements at the level of a few percent for N iter ≥ 8. The f=0 as prior   f=0 as prior  refolded distributions (right panel), obtained by convolving the unfolded distributions with the response function, indicate that convergence is reached for N iter ≥ 8. Figures 5 and 6 show similar distributions for v 3 and v 4 . The performance of the unfolding is similar to that shown in figure 4, except that the tails of the unfolded distributions require more iterations to converge. For example, figure 6 suggests that the bulk region of the v 4 distributions has converged by N iter = 32, but the tails still exhibit some small changes up to N iter = 64. The slower convergence for higher-order harmonics reflects the fact that the physical v n signal is smaller for larger n, while the values of δ 2SE change only weakly with n. These studies are repeated for all centrality intervals. In general, more iterations are needed for peripheral collisions due to the increase in δ 2SE .
The statistical uncertainties in the unfolding procedure are verified via a resampling technique [43]. For small N iter , the statistical uncertainties as given by the diagonal entries of the covariance matrix are much smaller than √ N , where N is the number of entries in each bin, indicating the presence of statistical bias in the prior. However, these uncertainties f=0 as prior  increase with N iter , and generally approach √ N for 64 ≤ N iter ≤ 128. In this analysis, the centrality range for each harmonic is chosen such that the difference between N iter = 32 and N iter = 128 is less than 10%. The centrality ranges are 0-70% for v 2 , 0-60% for v 3 and 0-45% for v 4 .
The robustness of the unfolding procedure is checked by comparing the results measured independently for the half-ID and the full-ID. The results are shown in figure 7. Despite the large differences between their initial distributions, the final unfolded results agree to within a few percent in the bulk region of the unfolded distribution, and they are nearly indistinguishable on a linear scale. This agreement also implies that the influence due to the slight difference (up to 1%) in the average track multiplicity between the two subevents is small. Systematic differences are observed in the tails of the distributions for v 4 , especially in peripheral collisions, where the half-ID results are slightly broader. This behaviour reflects mainly the deviation from the expected truth (residual non-convergence) for the half-ID unfolding, since the response function is a factor of √ 2 broader than that for the full-ID.
A wide range of priors has been tried in this analysis, consisting of the measured v obs n distribution and the six rescaled distributions defined by eq. (4.13). Figure 8 compares the convergence behaviour of these priors for v 3 in the 20-25% centrality interval. Despite the significantly different initial distributions, the unfolded distributions converge to the same answer, to within a few percent, for N iter ≥ 16. A prior that is narrower than the unfolded distribution leads to convergence in one direction, and a broader prior leads to convergence from the other direction. Taken together, these checks allow a quantitative evaluation of the residual non-convergence in the unfolded distributions.  non-convergence in the 2PC method, which has a broader response function than the single-particle method.
One important issue in the EbyE v n study is the extent to which the unfolded results are biased by non-flow short-range correlations, which may influence both the v obs n distributions and the response functions. This influence contributes to both the v obs , and hence are expected largely to cancel out in the unfolding procedure. This conclusion is supported by a detailed Monte-Carlo model study based on the HIJING event generator with a realistic flow afterburner [31], where the unfolding performance was evaluated. It is also  Figure 9. Comparison of the input distributions (solid symbols) and unfolded distributions for N iter = 128 (open symbols) between the single-particle unfolding and 2PC unfolding in the 20-25% centrality interval for v 2 (left panels), v 3 (middle panels) and v 4 (right panels). The ratios of 2PC to single-particle unfolded results are shown in the bottom panels.
supported by the consistency between the single-particle and 2PC methods (figure 9), which have different sensitivities to the non-flow effects. Furthermore, both unfolding methods have been repeated requiring a minimum η gap between the two subevents used to obtain the input distributions and the response functions. Six additional cases, requiring η gap = 0.2, 0.4, 0.6, 0.8, 1.0, 2.0, have been studied and the results have been compared (see figure 24). The unfolded v n distributions are observed to narrow slightly for larger η gap , reflecting the fact that the true v n decreases slowly at large |η| [16] and a larger η gap on average selects particles at large |η|. However, the results are always consistent between the two methods independent of the η gap value used. This consistency supports further the conclusion that the influence of the short-range non-flow correlations on the final unfolded results is not significant.
The dependence of the EbyE v n on the p T of the charged particles has also been checked: the particles are divided into those with 0.5 < p T < 1 GeV and those with p T > 1 GeV, and the EbyE v n measurements are repeated independently for each group of particles. About 60% of detected particles have 0.5 < p T < 1 GeV, and this fraction varies weakly with centrality. The unfolding performance is found to be slightly worse for charged particles with 0.5 < p T < 1 GeV than for those with p T > 1 GeV, due to their much smaller v n signal. Hence the v n range of the unfolded distribution for the final results is chosen separately for each p T range.
The final v n distributions are obtained using the single-particle unfolding with the full-ID and N iter = 128, separately for charged particles in the two aforementioned p T ranges and the combined p T range. The prior is obtained by rescaling the v obs n distribution according to eq. (4.13) with f = 0, and the response function is measured from correlations of the two half-IDs with no η gap in between.

Systematic uncertainties
The systematic uncertainties associated with the unfolding procedure include contributions from the residual non-convergence, dependence on the prior, uncertainty in the response function, the difference between the single-particle method and the 2PC method, and the tracking efficiency. The residual non-convergence is estimated from the difference between the results for N iter = 32 and N iter = 128 or between the results for the half-ID and full-ID. These two estimates are strongly correlated, so in each bin of the unfolded distribution the larger of the two is used. The dependence on the prior is taken as the difference between the results for f = 0 and f = 2.5. Note that the prior for f = 0 (f = 2.5) is broader (narrower) than the final unfolded distributions. The uncertainty of the response function is estimated from the difference between results obtained using the analytical formula eq. (4.7) and results obtained using the measured distribution, as well as the change in the results when the small dependence of δ 2SE on the observed v obs n is taken into account. Results for the p T -dependence of the v n distributions (see section 5 and figure 11) show that the mean values vary with p T , but, after they are rescaled to a common mean, the resulting shapes are almost identical. Motivated by this finding, every source of systematic uncertainty is decomposed into two components: the uncertainty associated with the v n or the v n -scale, and the uncertainty in the shape after adjustment to the same v n or the adjusted v n -shape. The uncertainties are then combined separately for the v n -scale and the adjusted v n -shape. Most shape uncertainties can be attributed to σ vn , such that the remaining uncertainties on the adjusted v n -shape are generally smaller.
To estimate the uncertainty due to the tracking efficiency, the measurement is repeated without applying the efficiency re-weighting. The final distributions are found to have almost identical shape, while the values of v n and σ vn increase by a few percent. This increase can be ascribed mainly to the smaller fraction of low-p T particles, which have smaller v n , so this increase should not be considered as a systematic uncertainty on the v n -scale. Instead, the scale uncertainty is more appropriately estimated from the change in the v EP n when varying the efficiency correction within its uncertainty range. On the other hand, small changes are observed for σ vn / v n and the adjusted v n -shape. Since these changes are small, they are conservatively included in the total systematic uncertainty in the v n -shape.
The variation of the efficiency with the detector occupancy may reduce the v n coefficients in eq. 1.1. This influence has been studied by comparing the v n values reconstructed via the two-particle correlation method with the generated v n signal in HJING simulation with flow imposed on the generated particles. The influence is found to be 1% or less, consistent with the findings in [15], and it is included in the uncertainty due to tracking efficiency. It should be pointed out that the influence of detector occupancy is expected to be proportional to the magnitude of the v n signal, and hence it mainly affects the v n scale, not the adjusted v n shape.
Additional systematic uncertainties include those from the track selection, dependence on the running period, and trigger and event selections. These account for the influence of fake tracks, the instability of the detector acceptance and efficiency, variation of the -17 -

JHEP11(2013)183
Uncertainty in v2 or σv 2 Uncertainty in σv 2 / v2 Centrality 0-10% 10-30% 30-50% 50-70% 0-10% 10-30% 30-50% 50-70%  Table 2. Summary of systematic uncertainties as percentages of v n , σ vn and σ vn / v n (n = 2-4) obtained using charged particles with p T > 0.5 GeV. The uncertainties for v n and σ vn are similar so the larger of the two is quoted. The uncertainties associated with track selection, the trigger and stability are taken from ref. [16]. For completeness, the model dependent estimates of the residual non-flow effects derived from ref. [31] for unfolding method with the default setup are also listed. Most uncertainties are asymmetric; the quoted numbers refer to the maximum uncertainty range spanned by various centrality intervals in each group. centrality definition, respectively. All three sources of systematic uncertainty are expected to change only the v n -scale but not the v n -shape, and they are taken directly from the published v EP n measurement [16]. Table 2 summarizes the various systematic uncertainties as percentages of v n , σ vn and σ vn / v n , for charged particles with p T > 0.5 GeV. The uncertainties on v n and σ vn are strongly correlated, which in many cases leads to a smaller and asymmetric uncertainty on σ vn / v n . In most cases, the uncertainties are dominated by tracking efficiency, as well as the track selection, trigger, and run stability assessed in ref. [16]. The uncertainties -18 -

JHEP11(2013)183
associated with the unfolding procedure are usually significant only in peripheral collisions, except for v 4 , where they are important across the full centrality range. The relative uncertainties in table 2 have also been evaluated separately for charged particles with 0.5 < p T < 1 GeV and p T > 1 GeV. In general, the systematic uncertainties are larger for the group of particles with 0.5 < p T < 1 GeV, due mainly to the increased contributions from residual non-convergence and the choice of priors.
The final v n distributions are shown over a v n range that is chosen such that the statistical uncertainties in all bins are less than 15%, and the results obtained with the default setups between N iter = 32 and N iter = 128 are consistent within 10%. The systematic uncertainties on the adjusted shape from the sources discussed above are then combined to give the final uncertainty for the v n -shape. The total systematic uncertainties are typically a few percent of v n in the main region of the v n distributions and increase to 15%-30% in the tails, depending on the value of n and centrality interval (see figure 10). Within the chosen v n ranges, the statistical uncertainties are found to be always smaller than the systematic uncertainties for the v n -shape, and the integrals of the v n distributions outside these ranges are typically < 0.5% for the 5%-wide centrality intervals and < 1% for the 1%-wide centrality intervals.
This analysis relies on the data-driven unfolding method to suppress the non-flow effects. In fact, many of the cross-checks presented in section 4.4 are sensitive to the residual non-flow, where "residual non-flow" refers to that component of the non-flow effects that is not removed by the unfolding method. Hence the total experimental systematic uncertainties quoted for the v n -scale in table 2 and the systematic uncertainties on the v n -shape discussed above already include an estimate of the residual non-flow effects. An alternative, albeit model-dependent approach, is to rely on simulations. One such study is carried out in ref. [31] based on HIJING with EbyE v n imposed on the generated particles. This study demonstrates that most non-flow effects are indeed suppressed by the datadriven unfolding method used in this analysis. This study also shows that the residual nonflow effects for the unfolding method with the default setup have no appreciable impact on the v 3 distributions, but broaden slightly the v 2 and v 4 distributions. Furthermore, ref. [31] also shows that most of these changes can be absorbed into simultaneous increases of v n and σ vn values by a few percent. For completeness, these model-dependent estimates of the residual non-flow contribution on the v n -scale are quoted in table 2: they are found to be smaller than the total systematic uncertainties of the measurement, especially for the higher-order harmonics. Ref. [31] also estimates the influence of residual non-flow on the intrinsic v n -shape, obtained by comparing the shapes of the unfolded and the input distributions after both distributions are rescaled to have the same v n . The variations of the intrinsic v n -shape due to residual non-flow are found to reach a maximum of 5%-15% (of p(v n )) in the tails of the distributions, depending on the choice of n and the centrality interval, but these variations are typically much smaller than the total systematic uncertainties on the v n -shape in the data (see figure 10). As noted above, a large fraction of the residual non-flow effects estimated by ref. [31] are expected to be already included in various data-driven cross-checks performed in section 4.4. Hence, to avoid double counting, the total systematic uncertainties applied to the data plots do not include the estimates from ref. [31]. However, the uncertainties from ref. [31] for the v n -scale (table 2) and v n -shape, are well within the total systematic uncertainties derived from the data analysis. Figure 10 shows the probability density distributions of the EbyE v n in several centrality intervals obtained for charged particles with p T > 0.5 GeV. The shaded bands indicate the systematic uncertainties associated with the shape. These uncertainties are strongly correlated in v n : the data points are allowed to change the shape of the distribution within the band while keeping v n unchanged. The v n distributions are found to broaden from central to peripheral collisions (especially for v 2 ), reflecting the gradual increase of the magnitude of v n for more peripheral collisions [15,16]. The shape of these distributions changes quickly with centrality for v 2 , while it changes more slowly for higher-order harmonics. These distributions are compared with the probability density function obtained from eq. (1.6) (v RP n = 0), which represents a fluctuation-only scenario for v n . These functions, indicated by the solid curves, are calculated directly from the measured v n values via eq. (1.7) for each distribution. The fluctuation-only scenario works reasonably well for v 3 and v 4 over the measured centrality range, but fails for v 2 except for the most central 2% of collisions, i.e. for the centrality interval 0-2%. Hence for v 2 the solid curve representing the fluctuation-only scenario is shown only for the 0-1% centrality interval (the data for the 1-2% interval are not shown). However, there is a small systematic difference between the data and the curve in the tails of the v 3 distributions in mid-central collisions, with a maximum difference of two standard deviations. Using eq. (1.9), this difference is compatible with a non-zero v RP 3 similar to the findings reported in ref. [44]. . Top panels: the unfolded distributions for v n in the 20-25% centrality interval for charged particles in the p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV ranges. Bottom panels: same distributions but rescaled horizontally so the v n values match that for the p T > 0.5 GeV range. The shaded bands represent the systematic uncertainties on the v n -shape. Figure 11 compares the unfolded v n distributions for charged particles in three p T ranges: p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV. The v n distributions for p T > 1 GeV are much wider than for 0.5 < p T < 1 GeV, reflecting the fact that the v n values increase strongly with p T in this region [16]. However, once these distributions are rescaled to the same v n as shown in the lower row of figure 11, their shapes are quite similar except in the tails of the distributions for n = 2. This behaviour suggests that the hydrodynamic response of the medium to fluctuations in the initial geometry is nearly independent of p T in the low-p T region; it also demonstrates the robustness of the unfolding performance against the change in the underlying v n distributions and response functions. Figure 12 shows a summary of the quantities derived from the EbyE v n distributions, i.e. v n , σ vn and σ vn / v n , as a function of N part . The shaded bands represent the total systematic uncertainties listed in table 2, which generally are asymmetric. Despite the strong p T dependence of v n and σ vn , the ratio σ vn / v n is relatively stable. For v 2 , the value of σ vn / v n varies strongly with N part , and reaches a minimum of about 0.34 at N part ∼ 200, corresponding to the 20-30% centrality interval. For v 3 and v 4 , the values of σ vn / v n are almost independent of N part , and are consistent with the value expected from the fluctuation-only scenario ( 4/π − 1 via eq. (1.8) as indicated by the dotted lines), except for a small deviation for v 3 in mid-central collisions. This limit is also reached for   Figure 12. The N part dependence of v n (left panels), σ vn (middle panels) and σ vn / v n (right panels) for n = 2 (top row), n = 3 (middle row) and n = 4 (bottom row). Each panel shows the results for three p T ranges together with the total systematic uncertainties. The dotted lines in the right column indicate the value 4/π − 1 expected for the radial projection of a 2D Gaussian distribution centred around origin (see eq. (1.8)). The values of σ vn / v n are compared with the σ ǫn / ǫ n given by the Glauber model [35] and MC-KLN model [45].

Results
peripheral collisions, where the resolution factor used in the EP method is small. The results in figures 10 and 12 imply that the distributions of v 2 in central collisions (centrality interval 0-2%), and of v 3 and v 4 in most of the measured centrality range are described by a 2D Gaussian function of ⇀ v n centred around the origin (eq. 1.6). On the other hand, the deviation of the v 2 distribution from such a description in other centrality intervals suggests that the contribution associated with the average geometry, v RP 2 , becomes important. In order to test this hypothesis, the v 2 distributions have been fitted to the Bessel-Gaussian function eq. (1.4), with v RP 2 not constrained to be zero. The results of the fit are shown in figure 14 for various centrality intervals. The fit works reasonably well up to the 25-30% centrality interval, although systematic deviations in the tails are apparent already in the 15-20% centrality interval. The deviations increase steadily for more peripheral collisions, which may be due to the fact that the fluctuations of ǫ 2 (eq. 1.2) are no longer Gaussian in peripheral collisions where N part is small [25] (see also figure 18).
The values of v RP centrality binning has been checked and corrected as follows. Taking the 20-25% interval as an example, the results obtained using the full centrality range within this interval are compared to the results obtained from the average of the five individual 1% intervals: 20-21%,. . . ,24-25%. This procedure has been carried out for each 5% centrality interval, and the difference is found to be significant only for the 0-5% and 5-10% intervals, and negligible for all the others. For the 0-5% interval, results are reported in the individual 1% bins. For the 5-10% bin, results are averaged over the five individual 1% bins.
As a cross-check, the Bessel-Gaussian fits are also performed on the v obs 2 distributions before the unfolding. Systematic devations are also observed between the fit and the v obs 2 data, but the deviations are smaller than those shown in figure 14. The value of v RP 2 from the v obs 2 distribution is found to agree to within a few percent with that from the unfolded v 2 distribution, while the value of δ v 2 from the v obs 2 distribution is significantly larger. This behaviour is expected since the smearing by the response function (eq. 4.7) increases mainly the width, and the value of v RP 2 should be stable.
Similarly, the value of v RP n can be estimated from σ vn and v n without relying on the fit: Both relations are shown in the left panel of figure 15 for v 2 . Good agreement with the data is observed for 100 < N part < 350, corresponding to 5-45% centrality interval in table 1 description [21]. The second-order coefficients for the 2k th -order cumulants, v 2 {2k}, are generally calculated via 2k-particle correlations [20]. These coefficients can also be calculated analytically from the measured v 2 distributions, labelled as v calc 2 {2k}. The first four terms can be expressed [23] as: The last part of each equation is exact when the v 2 distribution follows the Bessel-Gaussian function eq. (1.4). In this case only the first two cumulants are independent, and the higher-order cumulants do not provide new information. For the same reason, it has been argued that differences between the higher-order cumulants v 2 {2k} (k > 1) are sensitive to non-Gaussian behaviour in the underlying v 2 distribution [23]. Figure 16 shows results for v calc 2 {4}, v calc 2 {6} and v calc 2 {8}, calculated directly from the measured v 2 distributions (eq. (5.3)). They are compared with the results for v RP 2 obtained from the Bessel-Gaussian fits shown in figure 14. The results calculated from the cumulants agree with the results obtained from the fit, except for N part < 100. However, the calculated coefficients from higher-order cumulants agree with each other over the whole centrality range to within 0.5%-2%, despite the poor description of the v 2 distribution by the Bessel-Gaussian function for the centrality interval 25-70% ( N part < 200), shown in figure 14. It therefore follows that similar values of v 2 {2k} for k ≥ 2 observed in previous measurements [21,22] are no guarantee that the distribution is well described by a Bessel-Gaussian over the full range in non-central collisions, since the uncertainties of these measurements are bigger than the differences seen between  In previous analyses based on cumulants [7,22], several other quantities have been used to study the v 2 fluctuations. The definitions of these quantities are given below and their physical meaning can be understood from eq. (5.3): where F 1 , F 2 , and F 3 are calculated from the unfolded distributions, using eq. (5.3). The approximation for F 3 is valid when v RP 2 ≫ δ v 2 . In central collisions where v RP 2 ≪ δ v 2 , the value of F 3 is expected to approach one.  Figure 18 compares the EbyE v 2 distributions with the distributions of the eccentricity ǫ 2 of the initial geometry, calculated via eq. (1.2) from the Glauber model [35] and the MC-KLN model [45]. The MC-KLN model is based on the Glauber model but takes into account the corrections to the initial geometry due to gluon saturation effects. Three million events have been generated and grouped into centrality intervals according to the impact parameter. The ǫ 2 distribution for each centrality interval is rescaled to match the v 2 of the data, and then normalized to form a probability density function. Since v 2 is expected to be proportional to ǫ 2 in most hydrodynamic calculations [6], the deviations between the v 2 distributions and the rescaled ǫ 2 distributions can be used to improve the modeling of the initial geometry. Figure 18 shows that the rescaled ǫ 2 distributions describe the data well for the most central collisions, but not so well for non-central collisions. In peripheral collisions, both the Glauber and MC-KLN models fail to describe the data.   figure 19 and figure 20, respectively. Agreement with the models is better than in the n = 2 case. However, this could simply reflect the fact that these distributions are dominated by Gaussian fluctuations, which have a universal shape. The shape differences between the v n and ǫ n distributions are also quantified in the right panels of figure 12 by comparing the values of σ vn / v n with the values of σ ǫn / ǫ n . Clearly, both the Glauber and MC-KLN models fail to describe the data consistently, in particular for n = 2, across most -28 -

JHEP11(2013)183
of the measured centrality range. It should be noted that in the centrality intervals that correspond to the more peripheral collisions (e.g centrality interval 55-60%), the ǫ n values have been scaled down by a large factor and the sharp cutoffs of the ǫ n distributions are a natural consequence of the kinematic constraint ǫ n ≤ 1 (see eq. (1.2)). However, this behaviour also implies that v n is not proportional to ǫ n for large ǫ n values.

Summary
Measurements of the event-by-event harmonic flow coefficients v n for n = 2, 3 and 4 have been performed using 7 µb −1 of Pb+Pb collision data at √ s N N = 2.76 TeV collected by the ATLAS experiment at the LHC. The observed v n distributions are measured using charged particles in the pseudorapidity range |η| < 2.5 and the transverse momentum range p T > 0.5 GeV, which are then unfolded via a response function to estimate the true v 2 distributions. The response function is constructed via a data-driven method, which maps the true v n distribution to the observed v n distribution. The influence of residual non-flow effects are studied by varying the pseudorapidity gap between the two subevents used to obtain the observed v n distributions and the response functions, as well as by comparing the two different analysis methods based either on the single particle distribution or two-particle correlations. The influence of the residual non-flow is also compared to an estimation based on model simulations [31], and the latter is found to be within the total systematic uncertainty of this measurement.
The v n distributions are obtained in various centrality intervals: over the 0-70% centrality range for v 2 , 0-60% centrality range for v 3 and 0-45% centrality range for v 4 . The measured v 2 distributions are found to approach that of a radial projection of a 2D Gaussian distribution centred around zero in the 0-2% centrality range, which is consistent with a scenario where fluctuations are the primary contribution to the overall shape (fluctuationonly scenario) for these most central collisions. Starting with the centrality interval 5-10%, the v 2 distributions differ significantly from this scenario, suggesting that they have a significant component associated with the average collision geometry in the reaction plane, v RP 2 . In contrast, the v 3 and v 4 distributions are consistent with a pure 2D Gaussian-fluctuation scenario (i.e., v RP n = 0) over most of the measured centrality range. However, a systematic deviation from this fluctuation-only scenario is observed for v 3 in mid-central collisions, the presence of a non-zero v RP 3 can be allowed. Similarly, due to the limited range of the measured v 4 distribution, a non zero v RP 4 on the order of δ v 4 can not be excluded by this analysis.
The v n distributions are also measured separately for charged particles with 0.5 < p T < 1 GeV and p T > 1 GeV. The shape of the unfolded distributions, when rescaled to the same v n , is found to be nearly the same for the two p T ranges. This finding suggests that the hydrodynamic response to the eccentricity of the initial geometry has little variation in this p T region. The ratios of the width to the mean, σ vn / v n , of these distributions are studied as a function of the average number of participating nucleons ( N part ) and p T . The values of σ v 2 / v 2 are observed to reach a minimum of 0.34 for N part ≈ 200, while the values of σ v 3 / v 3 and σ v 4 / v 4 are nearly independent of N part and are close to the value expected for a pure Gaussian-fluctuation scenario.

JHEP11(2013)183
To understand further the role of average geometry and fluctuations for v 2 , the v 2 distributions have been fitted to a Bessel-Gaussian function to estimate the value of v RP 2 and the width of the fluctuation δ v 2 . In central collisions, where v RP 2 values are small and change rapidly with N part , narrow binning in centrality is necessary in order to obtain reliable estimates for these parameters. The values of δ v 2 /v RP 2 are found to decrease with N part and reach a minimum of 0.38 ± 0.02 at N part ≈ 200, but then increase and are greater than one in central collisions. Furthermore, a systematic deviation of the fit from the data is observed for centralities starting in the 15-20% centrality interval, and becoming more pronounced for the more peripheral collisions, suggesting significant non-Gaussian behaviour in the flow fluctuations for collisions with small N part . Multi-particle cumulant methods have been used to study such non-Gaussian behaviour. However, the v 2 coefficients for the four-, six-and eight-particle cumulants calculated from the v 2 distribution are found to agree with each other over the full centrality range to within 0.5%-2%. Hence the precision of experimental measurements of higher-order cumulants needs to be better than a few percent in order to be sensitive to the non-Gaussian behaviour in the v n distributions.
To elucidate the relation between the azimuthal anisotropy and underlying collision geometry, the measured v n distributions are compared with the eccentricity distributions of the initial geometry from the Glauber model and the MC-KLN model. Both models fail to describe the data consistently over most of the measured centrality range.

JHEP11(2013)183
A Comprehensive performance and data plots For reference, a more complete set of plots detailing the unfolding performance and crosschecks is presented. Figure 21 shows the distribution of the difference of the flow vectors ⇀ v obs 2 obtained for two half-IDs in peripheral collisions, calculated either without an η gap or with η gap = 2 between the two half-IDs. As mentioned in section 4.1, due to the small number of charged particles in the events, this distribution is expected to be described by the Student's t-distribution instead of the Gaussian distribution. Figures 22 and 23 show the dependence on the priors for v 2 and v 4 in the 20-25% centrality interval. Despite the very different shapes of the underlying flow distributions, the convergence is robust and independent of the prior. Figure 24 compares the unfolding performance for the singleparticle and 2PC methods for several choices of η gap and demonstrates the consistency of the two methods for all η gap values considered. Figures 25 and 26 Figure 33 shows the probability density distributions of v 2 for the full 0-5% centrality interval and the five individual 1% centrality intervals, together with the fits to the Bessel-Gaussian function. The deviation of the data from the Bessel-Gaussian description for the 0-5% centrality interval is due mainly to the rapid change of v RP 2 with centrality within this interval. This supports the need to use small centrality intervals for the most central collisions when calculating quantities such as v 2 , σ v 2 , v RP 2 and δ v 2 (see figure 15).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. ) b (arbitrary normalization) for the 65-70% centrality interval with η gap = 0 (left panel) and η gap = 2 (right panel), together with fits to a Gaussian function and a Student's t-distribution. Bottom panels: the ratios of the data to the t-distribution fits. The distributions are wider for η gap = 2 since the average number of tracks used is smaller.  Single v unfolded 4 2PC v Figure 24. The unfolded distributions from 2PC method (top row), single-particle method (middle row) and the ratios of the two (bottom row) for different values of η gap in the 20-25% centrality interval and for n = 2 (left panels), n = 3 (middle panels) and n = 4 (right panels). Rescaled to the same mean Figure 25. Top panels: the unfolded distributions for v n in the 0-5% centrality interval for charged particles in the p T > 0.5 GeV, 0.5 < p T < 1 GeV and p T > 1 GeV ranges. Bottom panels: same distributions but rescaled horizontally so the v n values match that for the p T > 0.5 GeV range. The shaded bands represent the systematic uncertainties on the v n -shape.     Figure 29. Comparison in 0.5 < p T < 1 GeV (left panel) and p T > 1 GeV (right panel) of the v RP 2 obtained from the Bessel-Gaussian fit of the v 2 distributions with the values for four-particle (v calc 2 {4}), six-particle (v calc 2 {6}) and eight-particle (v calc 2 {8}) cumulants calculated directly from the v 2 distributions via eq. (5.3). The bottom part of each panel shows the ratios of the cumulants to the fit results, with the error bars representing the total uncertainties. [16] ATLAS collaboration, Measurement of the azimuthal anisotropy for charged particle production in -38 -