M Monasor
email: monasor@uchicago.edu

J R Vázquez

D García-Pinto

F Arqueros

The impact of the fluorescence yield on the reconstructed shower parameters of ultra-high energy cosmic rays

An accurate knowledge of the fluorescence yield and its dependence on atmospheric properties such as pressure, temperature or humidity is essential to obtain a reliable measurement of the primary energy of cosmic rays in experiments using the fluorescence technique. In this work, several sets of fluorescence yield data (i.e. absolute value and quenching parameters) are described and compared. A simple procedure to study the effect of the assumed fluorescence yield on the reconstructed shower parameters (energy and shower maximum depth) as a function of the primary features has been developed. As an application, the effect of water vapor and temperature dependence of the collisional cross section on the fluorescence yield and its impact on the reconstruction of primary energy and shower maximum depth has been studied.

Introduction

Fluorescence telescopes record the longitudinal profile of extensive air showers induced by very energetic cosmic rays through the detection of the fluorescence light generated by secondary charged particles. This technique allows an accurate determination of the shower maximum depth X max . In addition, since the fluorescence intensity is proportional to the deposited energy, the integration in depth of the fluorescence profile allows a calorimetric determination of the shower energy. The total primary energy is later obtained by applying a correction accounting for the so-called missing energy. A key parameter for the reconstruction of the calorimetric energy is the fluorescence yield (FY), defined as the number of fluorescence photons emitted per unit deposited energy. The FY which depends on the atmospheric parameters (i.e. pressure P , temperature T , humidity) is measured in dedicated laboratory experiments. Presently one of the largest contributions to the total systematic error in the primary energy comes from the absolute value of this parameter. Also the uncertainties in the dependence of the FY on atmospheric parameters might have a non-negligible effect on both primary energy and X max measurements.

In this paper a simple procedure to study the effect of the assumed fluorescence yield on the reconstructed shower parameters will be shown. Preliminary results were published in [START_REF] Monasor | The impact of the fluorescence yield on the reconstructed shower parameters of ultra-high energy cosmic rays[END_REF]. In section 2 a brief summary of the FY properties will be presented. In the next section the features of some data sets of FY commonly used will be described and compared. A simple method to evaluate the effect of the FY selection on shower reconstruction is shown in section 4. The shower development will be described by a Gaisser-Hillas profile. The effect of a variation in the FY (including its atmospheric dependence) is a change in the reconstructed longitudinal development of the deposited energy which can be easily evaluated. Results for the dependence of both reconstructed X max and primary energy on the FY selection will be presented in section 5.

Keilhauer and Unger [START_REF] Keilhauer | Fluorescence emission induced by extensive air showers in dependence on atmospheric conditions[END_REF] have carried out a similar study focused on estimating the influence of the atmospheric conditions on the reconstructed shower parameters. However our approach is somewhat different since we only need the Gaisser-Hillas profile, while in [START_REF] Keilhauer | Fluorescence emission induced by extensive air showers in dependence on atmospheric conditions[END_REF] simulated showers are used for the analysis.

Air-Fluorescence Yield

Air-fluorescence in the near UV range is basically produced by the de-excitation of atmospheric nitrogen molecules excited by the shower electrons. The spectrum of fluorescence consists of a set of molecular bands represented by their wavelengths λ and therefore the total fluorescence yield Y in a given wavelength interval ∆λ can be obtained by adding up the contributions of all the molecular bands Y λ .

Y (X) = ∆λ Y λ (X) . (1)
Excited molecules can also decay by collisions with an environmental molecule. Because of this effect, the FY in the absence of quenching Y 0 λ (that is, at null pressure) is reduced by the so-called Stern-Volmer factor.

Y λ (P, T ) = Y 0 λ 1 + P/P (λ, T ) . (2)
The dependence of Y λ on atmospheric conditions can be described by a single parameter, the characteristic pressure P which is defined as the pressure for which both radiative and collisional de-excitation have the same probability. In general P contains a contribution of all possible quenchers i (i.e., N 2 , O 2 , H 2 O).

1 P = i f i P i , P i = kT τ 1 σ Ni vNi , vNi = 8kT πµ Ni . (3)
In the above expressions f i is the fraction of molecules of type i in the mixture, σ Ni is the collisional cross section which depends on the particular band, and v N i and µ N i are the relative velocity and reduced mass of the two body system N-i respectively; k is the Boltzman constant and τ the radiative lifetime of the corresponding level.

A number of experimental results of the FY in dry air are available (a review is presented in [START_REF]Proc. 5th Fluorescence Workshop[END_REF] and a critical comparison of the absolute measurements can be found in [START_REF] Rosado | [END_REF]). Usually the absolute values of Y λ are measured at atmospheric pressure and room temperature and the P parameters are evaluated at the same temperature. Equation [START_REF] Keilhauer | Fluorescence emission induced by extensive air showers in dependence on atmospheric conditions[END_REF] allows the determination of Y 0 λ and the dependence of the FY with pressure. Unfortunately measurements of the P parameters show large discrepancies and therefore the Y 0 λ parameter is not suitable for comparison between different experiments.

According to expression (3), the P value for dry air can be expressed as

1 P air (λ, T ) = f O P O (λ, T ) + f N P N (λ, T ) , (4)
where f O = 0.21 and f N = 0.79 are the fraction of oxygen and nitrogen molecules in air respectively. Measurements of P O (λ) and P N (λ) are also available [5].

Humidity dependence

Water molecules are also effective quenchers for excited N 2 and therefore fluorescence emission is partly suppressed in humid air. The characteristic pressure for humid air P hum (λ, T ) can be calculated from equation ( 3) by including the water contribution. The relationship between the characteristic pressures of dry P dry (λ, T ) and humid air

P hum (λ, T ) is 1 P hum (λ, T ) = 1 P dry (λ, T ) 1 - P w P + P w P 1 P w (λ, T ) , (5)
where P w is the partial pressure of water and P w is the characteristic pressure for fluorescence quenching by collisions with water molecules. Measurements of P w values have been performed by several authors. The AIRFLY collaboration has published results for four bands at 293 K [START_REF] Ave | Proc. 5th Fluorescence Workshop[END_REF]. In ref. [5] experimental values of P w for three relevant bands are also reported. More details can be found in [START_REF]Proc. 5th Fluorescence Workshop[END_REF].

Temperature dependence

The T dependence of the fluorescence yield is given by equations ( 2) and (3). Neglecting the dependence of the collisional cross section on the kinetic energy of the colliders, P is proportional to √ T . Assuming a measurement of the characteristic pressure at a reference temperature T 0 , the T dependence of P i (λ, T ) is given by

P i (λ, T ) = P i (λ, T 0 ) T T 0 . (6)
However, as it is well known, the average collisional cross section follows a behavior which can be described by a power law (∼ T α ) [START_REF] Holtermann | [END_REF]. Therefore equation ( 6) becomes

P i (λ, T ) = P i (λ, T 0 ) T T 0 • T α i λ 0 T α i λ = P i (λ, T 0 ) T 0 T α i λ -1/2 . ( 7
)
The dependence of the FY with both T and P can be easily predicted by combining equations (2), (3) and [START_REF] Holtermann | [END_REF]. The temperature dependence of the collisional cross section for dry air has been measured recently by the AIRFLY collaboration and the corresponding values of α λ parameter for four molecular bands have been published [START_REF] Ave | Proc. 5th Fluorescence Workshop[END_REF]. In principle, this parameter depends on the quencher (nitrogen, oxygen, water). However independent measurements of α i λ for each air component is difficult and usually it is assumed that α i λ takes the same value for all components. On the other hand α λ values for water are not presently available. However, as we will show below all these uncertainties and approximations have a marginal impact on the shower reconstruction.

Air-fluorescence yield data

Required information

The reconstruction of the shower parameters requires the following information on the fluorescence yield.

1) The absolute value in dry air at a given pressure and temperature for all bands within the spectral range of the telescope. As an alternative, the absolute value for a reference transition (e.g. 337 nm) and the relative intensities of the bands at a given pressure will also provide the total fluorescence yield.

2) The values of P in dry air at a reference temperature for all bands. Measurements of P N (λ, T 0 ) and P O (λ, T 0 ) could be also used for this purpose.

3) The values of P w (λ, T 0 ) for all wavelengths.

4) The T dependence of σ Ni for all wavelengths. In principle α λ values for each quencher would be needed, however this information is, at present, rather limited.

A data set containing 1) and 2) allows us to determine the total FY for dry air as a function of P at the reference temperature T 0 . Neglecting the T dependence of σ N i , the total FY can be calculated at any P and T condition. Adding 3), the air-fluorescence yield can be evaluated for humid air. Finally, 4) provides a more accurate extrapolation at temperatures far from the reference one (i.e. at high altitude). Apart from pioneering works [START_REF]Proc. 5th Fluorescence Workshop[END_REF], in the last 10 years several measurements of the FY have been carried out in laboratory experiments injecting accelerated electrons into air targets [START_REF]Proc. 5th Fluorescence Workshop[END_REF]. Nowadays three data sets combining some of these measurements are mainly being used in cosmic ray experiments using fluorescence telescopes. They are those of Kakimoto-Bunner, Nagano and the combination Nagano-AIRFLY. These data sets are described below. In all of them the humidity effect is neglected and the T dependence of quenching is calculated assuming a constant collisional cross section. The effect of these approximations on the fluorescence yield is not negligible and will be also evaluated in this section.

Kakimoto-Bunner (K-B)

This FY data set which was used by the HiRes collaboration in 2001 [START_REF] Dawson | Private Communication[END_REF] is a combination of fluorescence yield measurements of Kakimoto et al. [START_REF] Kakimoto | [END_REF] and the relative intensities of the molecular bands reported by Bunner [START_REF] Bunner | Cosmic Ray Detection by Atmospheric Fluorescence[END_REF].

Kakimoto et al. measured the absolute value of the air-fluorescence yield expressed in photons per meter λ for the three primary molecular bands of nitrogen (337, 357 and 391 nm) at 1000 hPa and 288 K. The remaining fluorescence bands are found to contribute with about a 30% of the total fluorescence photons when the filter employed in the HiRes telescopes (bandpass filter between 295 and 405 nm) is considered. As pointed out by [START_REF] Dawson | Private Communication[END_REF], these data translate into an absolute FY contribution of 0.46 ph/m for these bands. The remaining fluorescence spectrum (13 bands) is distributed according to the relative intensities collected by Bunner [START_REF] Bunner | Cosmic Ray Detection by Atmospheric Fluorescence[END_REF].

Kakimoto et al. calculated the fluorescence efficiency of the primary bands (defined as the fraction of deposited energy emitted as fluorescence radiation) φ λ as

φ λ = λ ρ (dE/dX) dep • hc λ , ( 8
)
where h is the Planck constant, c is the speed of light, (dE/dX) dep is the energy deposited by the electron per unit of atmospheric depth, and ρ is the air density. Experimental results of pressure dependence for the primary bands are also reported by Kakimoto et al. In that work the FY is parameterized as

λ (P, T ) = ρA λ 1 + ρB λ √ T . (9)
In this expression A λ is the number of photons emitted per electron and unit column density (g/cm 2 ) in the absence of quenching (i.e. P =0) and the B λ parameter is related with the characteristic pressure through ( 2) and ( 3) by

B λ = R gas √ T P air (λ, T ) , (10)
where R gas is the specific gas constant. The implementation of φ λ and B λ in equation ( 3) can be easily performed. Firstly, the fluorescence yield is related with the fluorescence efficiency by

Y λ = φ λ • λ hc . ( 11
)
Therefore the Y 0 λ values can be obtained from the fluorescence efficiency at P 0 =1000 hPa and the B λ parameters at T 0 = 288 K from

Y 0 λ = φ λ (P 0 ) λ hc 1 + P 0 B λ R gas √ T 0 . (12)
Finally, the values of P λ inferred from the reported B λ parameters of Kakimoto et al.

for the 337 nm band are used for all 2P components while that of the 391 nm is used for the 1N system.

Nagano

Nagano et al. [START_REF] Nagano | [END_REF] provided experimental results of λ in dry air for 15 nitrogen bands at 1013 hPa and 293 K as well as the corresponding P λ values at the same temperature. Using these measurements, the authors calculate the fluorescence efficiency at null pressure φ 0 λ . The corresponding Y 0 λ values can be easily obtained applying expression [START_REF] Bunner | Cosmic Ray Detection by Atmospheric Fluorescence[END_REF]. Notice that both Kakimoto et al. and Nagano et al. assume for this calculation that the energy lost by the electron is fully deposited in the field of view of the experimental set-up. Implication of this approximation has been discussed in [START_REF] Rosado | [END_REF], [START_REF]43rd Rencontres de Moriond: Electroweak Interactions and Unified Theories[END_REF] and [14].

Nagano-AIRFLY (N-A)

In this data set the absolute scale is given by the φ 0 337 value reported by Nagano et al. for the 337 nm band which translates into a Y 337 value of 6.38 ph/MeV at P 0 =800 hPa and T 0 =293 K. Relative intensities I λ at T 0 and P 0 for 34 wavelengths and the corresponding P (λ) values in dry air at 293K reported by AIRFLY [15] are used to get a complete set. The corresponding Y 0 λ for any wavelength is thus obtained as

Y 0 λ = Y 337 (P 0 , T 0 ) • I λ (P 0 , T 0 ) • 1 + P 0 P air (λ, T 0 ) . ( 13
)
This data set is presently used by the Pierre Auger Observatory [16].

Comparison of data sets

The FY for the above described data sets has been compared. In figure 1 the ratios Y Nagano /Y N-A and Y K-B /Y N-A against atmospheric depth are shown. The large deviations observed at low depth are due to discrepancies in the P values in the data sets. Obviously these discrepancies have no relevant impact in shower reconstruction since in our energy range they take place before the development of the shower. On the other hand, at larger depths differences around 20% between N-A and K-B data sets are found while a relative small discrepancy of around 2% is observed between N-A and Nagano descriptions.

In practice the comparison between FY data sets should include other ingredients. Fluorescence photons generated by cosmic ray showers are detected with ground telescopes located at a large distance from the emission point and thus suffering a nonnegligible attenuation which is wavelength dependent. In addition, the optical systems of the telescope, including the filter, the mirror and the PMT photocathode have a response which is also wavelength dependent. Therefore the relative differences between two FY data sets are not directly translated to the energy reconstruction since the relative intensities of the molecular bands are modified by the detector efficiency. For instance, discrepancies between two data sets which only affect wavelength regions for which the filter has a low transmittance might have a marginal impact on energy reconstruction although the discrepancy in total FY is large. The parameter more closely related with possible deviations in the shower reconstruction is the product Y λ (X)ε λ T λ (X, X 0 ) where ε λ is the optical efficiency of the telescopes including all components and T λ (X, X 0 ) represents the atmospheric transmittance between the emission point and the telescope location.

The effect of the optical filter is particularly relevant in the comparison between K-B and N-A data sets. As shown in figure 1 the large disagreement between them nearly disappears when the efficiency of the optical system of the Auger telescopes [17] is included. Figure 2 shows the fluorescence yield for both data sets (left) and the effect of the optical efficiency (right). Note that the large discrepancy in the bin centered at 317 nm is smoothed by the strong absorption of the filter at this wavelength. However the discrepancies at wavelengths with good transmission, and thus less sensitive to the filter, nearly compensate each other.

In a similar way, it has been shown [START_REF] Nagano | [END_REF] that the large disagreement between the Nagano and K-B data sets are significantly reduced when the optical efficiency of HiRes is included.

Effect of temperature and humidity on the FY

The above described data sets do not include the contribution of quenching by water molecules in humid air. In addition, the temperature dependence of the characteristic pressure is assumed to follow a √ T law (i.e. α λ is assumed to be zero). The effect of these approximations on the FY profile has been evaluated using the monthly average description for the atmosphere measured at the Auger site [18]. As an example, in figure 3 the effect of including the humidity and temperature corrections (December month) in the N-A data set is shown against atmospheric depth. For these calculations the quenching parameters of [START_REF] Ave | Proc. 5th Fluorescence Workshop[END_REF] have been used. As expected, humidity decreases the fluorescence yield at large depth (i.e. near ground) while the effect of the T term is only relevant at high altitudes. Neglecting the possible dependence of the collisional cross section N * 2 -H 2 O with T (i.e. assuming α λ is zero for water), the total effect of humidity and temperature can be easily computed. The result is also shown in figure 3. A similar behavior has been found in a study carried out in [START_REF] Keilhauer | Proc. 5th Fluorescence Workshop[END_REF] using several sets of fluorescence yield data including the temperature and humidity dependence of quenching from various authors.

As expected, neither the optical efficiency of the detector nor the atmospheric transmission introduce any significant effect in these results since atmospheric conditions (i.e. temperature and humidity) has a negligible impact on the fluorescence spectrum.

Method

In this note we propose a simple procedure to evaluate the effect of the fluorescence yield selection on the reconstruction of shower parameters. Assuming a certain longitudinal profile of deposited energy dE(X)/dX, the number of fluorescence photons generated per unit atmospheric depth is determined by the fluorescence yield in the wavelength in-terval of the telescope.

dn gen γ dX = dE dX Y (X) . (14)
The total calorimetric energy E can be obtained from the integral of the longitudinal profile and therefore

E = ∞ 0 1 Y (X) dn gen γ dX dX . ( 15
)
If the FY assumption is now changed to a new data set Y (X), the profile of deposited energy consistent with the fluorescence profile is

dE dX = dE dX Y (X) Y (X) , (16)
and thus the calorimetric energy associated to the new FY selection will be

E = ∞ 0 dE dX Y (X) Y (X) dX . (17)
Therefore, in principle, the effect of changing the FY selection on the primary energy can be evaluated by comparing (15) and (17). On the other hand, the effect on the X max value can be obtained by comparing the shape of the dE/dX and dE /dX profiles. Only a small fraction of the fluorescence photons generated by the shower reaches the PMT camera of the fluorescence detector. In the first place, the number of photons is strongly reduced by a factor A/(4πR 2 (X)) where A is the area of the telescope and R is the distance from the emission point to the telescope location. In addition, as previously pointed out, the atmosphere scatters a non-negligible fraction of photons on their way to the telescope. Finally the optical elements of the telescope also absorb a certain fraction of those photons reaching the telescope window. Both atmospheric transmission and optical efficiency are wavelength dependent. Therefore the profile of deposited energy is calculated from the profile of observed photons using the expression

dE dX = dn obs γ dX • 4πR(X) 2 A • 1 ∆λ ε λ • T λ (X) • Y λ (X) , (18)
and thus, relationship (16) becomes

dE dX = dE dX Y λ (X) • ε λ • T λ (X) Y λ (X) • ε λ • T λ (X) . ( 19
)
Therefore the shower maximum depth has to be evaluated using [START_REF] Keilhauer | Proc. 5th Fluorescence Workshop[END_REF] instead of ( 16) and the calorimetric energy E inferred from the new FY selection is given by

E = ∞ 0 dE dX Y λ (X) • ε λ • T λ (X) Y λ (X) • ε λ • T λ (X) dX . (20)
The variation of both X max and E due to a change in the FY selection is not sensitive to the fine details of the longitudinal development of dE/dX and thus instead of using real data or simulated showers it is easier and more convenient for this study to describe the longitudinal development of deposited energy by a Gaisser-Hillas (GH) profile [START_REF] Gaisser | Proc. 15th International Cosmic Ray Conference[END_REF] given by

dE dX = dE dX Xmax X -X 0 X max -X 0 (Xmax-X0)/λ e Xmax-X λ , (21)
where X 0 and λ are shape parameters and X max is the shower maximum depth.

We have studied the impact of the FY selection for proton and iron showers with zenith angle θ of 30 o and 60 o and fixed primary energies E 0 of 10 19 and 10 20 eV. The corresponding calorimetric energies have been calculated assuming an invisible energy as parameterized by [START_REF] Song | [END_REF]. The values of X 0 λ and X max have been obtained by fitting the average longitudinal development of a sample of CORSIKA showers [START_REF] Perrone | Private Communication[END_REF] to eq. ( 21). These longitudinal developments are fully compatible with those shown in [START_REF] Song | [END_REF]. The value of (dE/dX) Xmax is obtained as a normalization constant to account for the total deposited energy. For a given energy, the shower reaches its maximum development at an altitude which grows with θ and therefore, since the FY varies with altitude (through pressure, temperature and humidity), the impact of the FY selection depends on the shower geometry. The effect of atmospheric transmission increases with the distance between the telescope and the shower axis. Values of 10 and 30 km have been tried for this parameter. As already mentioned the atmosphere has been modeled using monthly average profiles measured at the Auger site [18].

As an example we show in figure 4 the effect in the profile of deposited energy of a 10 19 eV proton shower when changing from Y N-A 16), the modified profile (i.e., the profile of deposited energy which would be obtained using the K-B fluorescence yield), has been calculated (red dotted line). The integral of this profile results in a total deposited energy significantly larger (by about a 20%). However when taking into account the effect of the optical efficiency, the ratio

Y K-B λ (X) • ε λ / Y N-A λ (X)
• ε λ is close to unity except at very low atmospheric depth (dashed blue line measured in the blue axis on the right) and therefore the reconstructed profile from expression (19) (dotted blue) is nearly unaffected giving a deviation in the deposited energy around 2% in full agreement with previous studies [24]. The effect of the atmospheric transmission does not introduce significant differences between FY data sets in this case. This example has been calculated using the atmospheric profile of December. Similar results have been found for other seasons. The effect of humidity and that of a T -dependent collisional cross section has also been studied with this method. Some typical examples are shown in figures 5 and 6. In figure 5 the effect of water vapor (December month) in the N-A data set for a Fe shower of 10 20 eV and θ angles of 30 o and 60 o is studied. The development of the shower depends on the slant depth and therefore the smaller is θ, the closer to the ground the maximum shower development is reached. Since humidity is stronger at low altitudes, its effect is more important for vertical showers. In this case the calorimetric energy has to be increased by about a 4% at 30 o while at 60 o the correction is negligible (< 1%). On the contrary the effect of neglecting the T dependence of the collisional cross section is stronger at low temperatures (i.e. at high altitudes) since the T 0 value at which the FY is measured are close to that at ground. This feature can be observed in figure 6 for the same kind of shower assuming the same atmospheric profile. In this case the effect of temperature can be corrected by increasing E by about 2% (6%) at 30 o (60 o ).

The corresponding effect of the correction on the shower maximum depth is also displayed in the above figures. As expected the X max value increases (decreases) if the slope of the FY ratio versus depth function at X = X max is negative (positive).

The longitudinal profiles have been calculated until the showers reaches ground level that in the case of the Auger site corresponds to a vertical atmospheric depth of about 880 gcm -2 .

Results

The method described above has been used to predict the impact of various FY assumptions on the shower reconstruction. The relative difference in reconstructed primary energy between several FY selections and the N-A one δ E = (E 0 -E N-A 0 )/E N-A 0 has been calculated as well as the absolute difference of the shower maximum depth ∆X max = X max -X N-A max . The effect of the optical efficiency of the telescope [17] (in particular that of the filter) and the atmospheric absorption to these deviations has also been evaluated.

As already mentioned the contribution of the invisible energy has been calculated following the parameterization of [START_REF] Song | [END_REF]. Since typical variations in deposited energy due to the FY selection are smaller than 10% and the invisible energy contributes an amount ranging from 8% and 12%, the corresponding deviations in E 0 are smaller than the one of the calorimetric energy by as much as one percent unit.

In the first place the FY data sets of K-B and Nagano have been compared with the N-A one. The results for total energy are displayed in tables 1 (10 20 eV) and 2 (10 19 eV) while those of X max are shown in tables 3 (10 20 eV) and 4 (10 19 eV). As can be seen in tables 1 and 2, the Nagano FY would lead to a small decrease in energy (around or smaller than 2%) without any additional relevant effect when the optical efficiency of the fluorescence telescope or the atmospheric transmission is included in the calculation. On the contrary, using the K-B FY would give rise to a significant energy increase of nearly 20%. However, as anticipated in the previous section, this deviation is reduced down to 12

K-B/N-A Nagano/N-A 30 o 60 o 30 o 60 o p Fe p Fe p Fe p Fe Y λ (X) 19 19 19 20 -1 -1 -1 -1 Y λ (X) • F λ 7 7 7 8 0 0 0 0 Y λ (X) • ε λ 2 2 2 3 -2 -2 -1 -1 Y λ (X) • ε λ • T λ (X) (10km) 1 1 1 2 -2 -2 -1 -1 Y λ (X) • ε λ • T λ (X) (30km) 0 0 0 0 -2 -2 -2 -2
Table 1: Percentage differences (100δ E ) between the primary energy reconstructed using either K-B or Nagano instead of the N-A FY data set, for typical GH profiles of p and iron 10 20 eV showers. These results have been obtained using the atmospheric profile of December. The effect of the filter F λ , the total optical efficiency ε λ and the atmospheric transmission T λ assuming shower-telescope distances of 10 and 30 km is shown. Uncertainties in 100δ E are estimated in about ±0. about 2% when the optical efficiency of other components are included as pointed out by [24]. In particular the filter reduces this nearly 20% energy increase down to about 7%, with the optical efficiency of the remaining components of the telescope responsible for the further decrease down to that final 2%. The effect of the atmospheric transmission is basically negligible In regard to X max reconstruction, using the FY data sets of Nagano or K-B would decrease the shower maximum depth by about 1 gcm -2 (see tables 3 and4). The effect of the atmospheric transmission is nearly negligible. These results both of reconstructed energy and maximum depth are basically independent on primary energy and shower angle. As expected, they are also nearly independent of season.

-1 -1 -1 Y λ (X) • F λ 7 7 8 8 0 0 0 0 Y λ (X) • ε λ 2 2 3 3 -2 -2 -1 -1 Y λ (X) • ε λ • T λ (X) (10km) 1 1 2 2 -2 -2 -1 -1 Y λ (X) • ε λ • T λ (X) (30km) 0 0 0 0 -2 -2 -2 -2
The effect of the humidity and the T -dependent collisional cross section on the N-A data set has also been studied. As already mentioned these effects vary with geometry and season due to the different atmospheric conditions undergone by the shower track, as well as with the nature of the primary (mass and energy) due to the different longitudinal development.

When the water vapour quenching is taken into account, large deviations in the calorimetric energy (up to 7%) are found for humid months (e.g. February) and shower maximum close to ground (i.e., either high energy, light primary, vertical incidence or a combination of them). This deviation lowers down to about 1% at 60 o . The effect of a collisional cross-section dependent on temperature ranges from 6.5% for inclined showers to about 2% for vertical ones. The total effect in the primary energy is displayed for all months in figure 7. Excluding the February -May period, it ranges typically between 5.5 and 6.5%. Notice the significant variation of δ E with season for vertical showers, which is due to the humidity effect. The correction in shower maximum depth due to the temperature effect is ∆X max ≈-5 gcm -2 , nearly independent on the shower and atmospheric features while the correction for the humidity effect ranges from about 8 gcm -2 for vertical showers in austral summer down to about 1 gcm -2 at 60 o . The total X max deviations including both effects are displayed for all months in figure 8. As can be seen in this figure, the corrections for inclined showers are negative and nearly season independent. This behaviour was expected since the slope of the function FY ratio versus depth is positive for the T effect (see figure 3) and, as already mentioned, the temperature profile is nearly constant along the year. On the contrary the correction for vertical showers is strongly dependent on the month again due to the humidity effect.

K-B/N-A Nagano/N-A 30 o 60 o 30 o 60 o p Fe p Fe p Fe p Fe Y λ (X) -1 -1 -2 -2 -1 -1 -1 -1 Y λ (X) • F λ -1 -1 -2 -2 -1 -1 -1 -1 Y λ (X) • ε λ -1 -1 -2 -1 0 0 -1 -1 Y λ (X) • ε λ • T λ (X) (10km) 0 0 -1 -1 0 0 -1 -1 Y λ (X) • ε λ • T λ (X) (30km) 1 0 -2 -2 0 0 -1 -1
(X) -1 -1 -2 -2 -1 -1 -1 1 Y λ (X) • F λ -1 -1 -2 -2 -1 -1 -1 -1 Y λ (X) • ε λ -1 -1 -2 -2 0 -1 -1 -1 Y λ (X) • ε λ • T λ (X) (10km) 0 -1 -1 -1 0 -1 -1 -1 Y λ (X) • ε λ • T λ (X) (30km) 0 0 -2 -2 0 0 -1 -1
The experimental analysis of the fluorescence quenching due to water vapor as well as that of the temperature dependence of the quenching cross section is rather difficult. As a consequence these measurements are very scarce and subject to large uncertainties. The procedure presented here can be used for the evaluation of the effects of the errors in the P w and α parameters on the total uncertainty of the reconstructed values of the shower energy and X max . For this purpose we have recalculated δE and ∆X max after increasing (decreasing) α λ by a factor 1.5 (0.5) and increasing (decreasing) P w by a factor 1.2 (0.8) for all molecular bands. In general the effect of such large eventual deviations in these parameters is very small. The largest effect is found for high humidity months and/or showers more sensitive to T and humidity effects.

For instance for 30 o protons of 10 20 eV the 8% increase in energy for the February month would move within the 6% to 10% interval while the corresponding shift in maximum depth of ∆X max = 9 gcm -2 remains unchanged. For the same kind of shower (p, 30 o , 10 20 eV) the 5% energy increase in the August month would be within the 3 to 6% interval and the negligible X max effect reported in fig. 8 would be in between -2 and -1 gcm -2 . Finally for Fe showers of 60 o and 10 20 eV the 6% energy increase in August would be in between 3 and 9% and the corresponding ∆X max = -4 gcm -2 value would be within the -5 to -1 gcm -2 interval. Note that the above examples represent the cases with the largest dependence on possible uncertainties in the P w and α parameters.

Conclusions

A simple analytical method has been proposed to quantify the influence of the assumed FY on the shower reconstruction, in particular on primary energy and shower maximum depth. The longitudinal development of the shower is described by a Gaisser-Hillas profile which is slightly modified by a variation in the FY assumption. Several data sets of FY including absolute values, wavelength spectra, as well as pressure, temperature and humidity dependencies have been used for this study.

Using this simple procedure we have confirmed that the Nagano-AIRFLY, Nagano and Kakimoto-Bunner data sets lead to close values (within around a 2%) of the reconstructed energy, as long as the optical efficiency of the Auger telescopes is taken into account. The corresponding effects on X max are typically smaller than 1 gcm -2 .

The dependence of the FY with atmospheric properties (P , T and humidity) and its effect on reconstructed parameters of the shower have also been analyzed. These effects, when combined, introduce typical deviations at the level of 6% in the reconstructed primary energy. However significantly larger deviations can be found for vertical showers (30 o ) in humid months. The deviations induced in X max by temperature and humidity are opposite and they nearly cancel when combined. Deviations of about -4 gcm -2 are found for inclined showers (60 o ) while for vertical ones the effect is strongly dependent on humidity with large positive deviations in some months.

These deviations in primary energy and X max values are not negligible when compared with typical uncertainties (both systematic and statistical) in the reconstruction of the shower parameters. Notice, however, that they can be easily corrected as far as the temperature and humidity effects are included in the reconstruction algorithms. On the other hand, we have found that even relatively large uncertainties in the quenching parameters associated to humidity or the T-dependence of the collisional cross section have not a significant impact on the reconstructed parameters. In fact the derived systematic uncertainties are below 2% for primary energy and below 2gcm -2 for X max , and therefore they represent a small contribution to the total uncertainty in the reconstruction of the shower parameters.

The results shown here have been obtained from the whole longitudinal development of the fluorescence light. In a real case the reconstruction of the longitudinal development is restricted to an interval given by the field of view of the telescope. In addition, the contribution of Cherenkov light is a very important ingredient in the reconstruction of real data [25] not used in our simple procedure. These simplifications could give rise to small discrepancies between the predictions of our simple algorithm and those from the analysis of real data. Nevertheless our results are basically in agreement with those of [START_REF] Keilhauer | Fluorescence emission induced by extensive air showers in dependence on atmospheric conditions[END_REF] obtained from a detailed reconstruction of a sample of MC events.

Figure 1 :

 1 Figure 1: Ratio of Nagano and Kakimoto-Bunner fluorescence yields to that of Nagano-AIRFLY versus atmospheric depth. The ratio K-B/N-A of ΣY λ (X)ε λ , also displayed, indicates that the effect of the optical efficiency (filter and telescope) reduces significantly the large disagreement in the total FY of both data sets.

Figure 2 :

 2 Figure 2: Fluorescence yield in absolute units versus wavelength for the N-A and K-B data sets (left). Since K-B intensities are given in 9 nm bins [9], the N-A fluorescence yields have been recalculated for the same wavelength bins for a better comparison. The effect of the Auger optical efficiency (black line) on the FY of both data sets is illustrated (right).

Figure 3 :

 3 Figure 3: Effect of the humidity and a T-dependent collisional cross section on the N-A data set. In this plot the ratio of fluorescence yields for the month of December versus atmospheric depth is shown.

Figure 4 :

 4 Figure 4: Comparison of FY N-A and K-B descriptions and their impact in the shower reconstruction for a E = 10 19 eV proton shower with 30 o zenith angle (black profile). The total FY of K-B and N-A (red axis on the right) are significantly different and therefore with a large impact on the profile of deposited energy (red profile). When including the optical efficiency of the telescopes the FY ratio (blue axis on the right) is very close to one and their impact on deposited energy in nearly negligible (blue profile).

λ

  (X) to Y K-B λ (X) fluorescence yield. The black solid line is the GH profile for a shower with an E 0 value of 10 19 eV. The red dashed line represents the ratio Y K-B λ /Y N-A λ (red axis on the right) as a function of depth. Applying (

Figure 5 :

 5 Figure 5: Effect of including the humidity contribution in the profile of energy deposited for 30 o (left) and 60 o (right) zenith angles for a Fe shower of 10 20 eV. Notice that the dependence of the FY ratio with slant depth varies with the incoming angle. See text for details.

Figure 6 :

 6 Figure 6: Effect of including the temperature dependence of the cross-section in the shower reconstruction for 30 o (left) and 60 o (right) zenith angles for the same case of figure 5. See text for details.

Figure 7 :

 7 Figure 7: Percentage increase in reconstructed E 0 value when including humidity and a collisional cross section dependent on temperature to the N-A data set, versus month. Results are shown for both protons and Fe nuclei of 10 19 eV (triangles) and 10 20 eV (circles) energies.

Figure 8 :

 8 Figure 8: Increase of shower maximum depth due to humidity and a collisional cross section dependent on temperature to the N-A data set, versus month. See figure 7 for more details.

Table 2 :

 2 Same as table 1 for 10 19 eV.

Table 3 :

 3 ∆Xmax differences (gcm -2 ) between the reconstructed values using K-B or Nagano instead the N-A data set, for typical GH profiles of p and iron 10 20 eV showers. Uncertainties in ∆Xmax are estimated in about ±0.5 gcm -2 . Null results represent absolute deviations below 0.5gcm -2 . More details on first column are given in table 1.

	K-B/N-A	Nagano/N-A
	30 o	60 o	30 o	60 o
	p Fe p Fe p Fe p Fe
	Y λ

Table 4 :

 4 Same as table 3 for 10 19 eV.

Acknowledgements

This work has been supported by the Spanish Ministerio de Ciencia e Innovacion (FPA2009-07772 and CONSOLIDER CPAN CSD2007-42) and "Comunidad de Madrid" (Ref.: 910600). M. Monasor acknowledges the "Consejería de Educación y Ciencia de Castilla-La Mancha" and the "Fondo Social Europeo" for a postdoctoral fellowship. Very fruitful discussion with our colleagues of the Auger collaboration are acknowledged, in particular with B. Keilhauer, J. Matthews, V. H. Ponce and M. Unger.


Q Avril

V Gouranton

B Arnaldi

Fast Collision Culling in Large-Scale Environments Using GPU Mapping Function

This paper presents a novel and efficient GPU-based parallel algorithm to cull non-colliding object pairs in very large-scale dynamic simulations. It allows to cull objects in less than 25ms with more than 100K objects. It is designed for many-core GPU and fully exploits multi-threaded capabilities and data-parallelism. In order to take advantage of the high number of cores, a new mapping function is defined that enables GPU threads to determine the objects pair to compute without any global memory access. These new optimized GPU kernel functions use the thread indexes and turn them into a unique pair of objects to test. A square root approximation technique is used based on Newton's estimation, enabling the threads to only perform a few atomic operations. A first characterization of the approximation errors is presented, enabling the fixing of incorrect computations. The I/O GPU streams are optimized using binary masks. The implementation and evaluation is made on largescale dynamic rigid body simulations. The increase in speed is highlighted over other recently proposed CPU and GPU-based techniques. The comparison shows that our system is, in most cases, faster than previous approaches. evolution of graphics hardware and multi-GPU platforms have provided programmers with better processing tools. With these new tools it is now essential to take into account

Introduction

Collision detection is a well-studied and active research field where the main problem is to determine how and if one or more objects collide in a 3D virtual environment. Collision detection is an issue affecting many different fields of study, including physical-based simulation, computer animation, robotics, haptic applications and video games. In these applications, real-time performance, efficiency and robustness are key issues. In the field of Virtual Reality, physical virtual environments in digital mock-ups and industrial applications are now commonplace, and are of increasing complexity. The expected level of real time performance is becoming harder to ensure in such large-scale virtual environments. Unsurprisingly, collision detection has been an integral part of virtual reality bottlenecks for over thirty years.

The use of parallel processing has become necessary to take advantage of recent gains in Moore's Law. Over the past several years, processor specialists have provided clock frequency increases and parallelism improvements in instruction sets. Now, to better manage power consumption, they promote multi-core architectures. It is no longer possible to rely on the evolution of processing power to overcome the problem of real-time collision detection. The impressive run-time architectures to improve collision detection performance in large-scale virtual environments.

Main Results: A fast and novel approach to cull noncolliding object pairs in large-scale environments is proposed. The approach is based on the use of new semi bruteforce GPU mapping function used during the broad phase step. New access functions are proposed to optimize execution time of the GPU kernel reducing global GPU memory accesses. The kernel time is also reduced using square root approximation based on Newton's approximation where a first way to characterize and fix errors is presented. Input and output data transfers are optimized using binary masks to fully use the entire binary stream. This new broad phase algorithm enables to get an "interactive" computation time (<25ms) with more than 100K objects.

Organization: The rest of the paper is organized as follows: in Section 2 related work on parallel collision detection and broad phase techniques is reported. Section 3 describes the new brute force broad phase algorithm and Section 4 gives details on the dynamic spatial subdivision technique adapted to a GPU use. The implementation and performance measurements are described in Section 5. The comparison and analysis, over major CPU and GPU-based techniques, are done in Section 6. Then, the conclusion and the opening on future work end the paper.

Related Work

A brief overview is presented on parallel solutions for collision detection focusing on the broad phase step. For more details on the collision detection field, the reader may refer to surveys on the topic [LG98, TKH * 05].

Broad Phase Collision Detection

The collision detection process is organized and built as a pipeline [START_REF] Hubbard | Collision detection for interactive graphics applications[END_REF] composed by two main parts: broad and narrow phase. The broad phase step is in charge of a quick and efficient removal of object pairs that are not in collision. In the bounding volume family many models have been proposed such as spheres [START_REF] Hubbard | Collision detection for interactive graphics applications[END_REF], Axis-Aligned-Bounding-Box (AABB) [Ber97], Oriented-Bounding-Box (OBB) [GLM96], discrete oriented polytopes (k-DOP) and many others. Other methods have been proposed to use spatial partitioning and divide space into unit cells: regular grid [START_REF] Overmars | Point location in fat subdivisions[END_REF] or non-regular [Mir97, EG07], quadtree, octree [BT95], Binary Space Partitioning (BSP) or k-d tree structure. Topological approach is based on the positions of objects (Sweep and Prune (SaP) [CLMP95]). There are two related but different concepts on the way the SaP operates internally: starting from scratch each time (brute force) or updating internal structures (persistent). The persistent algorithm has been enhanced by using a segmented interval list combined with subdivision [START_REF] Tracy | Efficient large-scale sweep and prune methods with AABB insertion and removal[END_REF]. It provides faster sequential execution in large-scale environments. The brute-force SaP has been paralyzed on multi-core architecture [START_REF] Gouranton | Dynamic adaptation of broad phase collision detection algorithms[END_REF] and it was shown that, according to the number of objects and the number of CPU cores, the parallelization of the brute force method may become more efficient than the persistent one. Luque et al. [LCF05] proposed a new structure called "Semi-Adjusting BSP-tree" to represent scenes consisting of thousands of objects. They showed that the tree does not require a complete restructuring even for highly dynamic scenes.

Parallel Solution for Collision Detection

The parallel solution of collision detection algorithms is a recently explored and active field in high performance computing. There are three different families of algorithms, namely CPU, GPU and hybrid-based solutions.

A new task splitting approach for implicit time integration and collision handling on multi-core architecture has been proposed [TPB08]. Tang et al. [TMT09] propose to use a hierarchical representation to accelerate collision detection queries and an incremental algorithm exploiting temporal coherence. The overall is distributed among multiple cores. They obtained a 4X-6X speed-up on an 8-core processor based on several deformable models. Kim et al. [KHeY08] propose to use a feature-based bounding volume hierarchy (BVH) to improve performances of continuous collision detection. They also propose a new method of task decomposition for their BVH-based collision detection and dynamic task assignment methods. They obtained a 7X-8X speed-up using an 8-core architecture compared to a single-core. Hermann et al. [HRF09] propose a parallelization of interactive physical simulations. They obtain a 14X-16X speed-up on a 16-core architecture compared to a single-core. Tang et al. [TLW11] proposed a new parallel approach based on a new concept of filters into subspaces. They propose a new phase between the broad and narrow phase where they apply two successive filters (linear and planar filters). This technique enables to widely prune elementary tests and significantly reduce the computation time.

Cinder [START_REF] Knott | Cinder: Collision and interference detection in real-time using graphics hardware[END_REF] is an algorithm exploiting the GPU to implement a ray-casting method to detect collision. GPU-based algorithms for self-collision and cloth animation have also been introduced by Govindaraju et al. [GLM05]. A solution using image-space visibility queries has been proposed for the broad phase [START_REF] Govindaraju | Cullide: Interactive collision detection between complex models in large environments using graphics hardware[END_REF]. Le Grand [LG07] has presented a GPU implementation of the SaP broad phase algorithm based on CUDA faster than the CPU implementation. He described how to implement a spatial subdivision method for GPU between particles. Liu et al. [START_REF] Liu | Real-time collision culling of a million bodies on graphics processing units[END_REF] recently presented a new GPU-based SaP. They propose to reduce the density of intervals along the sweep axis by using principal component analysis to choose the best sweep direction. They couple it with spatial subdivisions to further reduce the number of false positive overlaps. Recent works use thread and data parallelism on a single GPU to perform fast hierarc The Eurographics Association 2012. chy construction, updating, and traversal using tight-fitting bounding volumes such as oriented bounding boxes (OBB) and rectangular swept spheres (RSS) [START_REF] Lauterbach | gproximity: Hierarchical gpu-based operations for collision and distance queries[END_REF]. Tang The total number of objects in recent virtual environments and games increase and become more dynamic. The performance of the most-used algorithms for the broad phase step is becoming a serious issue. The Sweep and Prune algorithm does not scale very well and its performance decreases in very large scale environments. More precisely, the insertion and sorting of new objects as well as the interaction with far away objects is very time-consuming.

Optimized Broad Phase

In a n-body simulation, a semi-brute force algorithm performs n 2 -n 2 computations at each loop (cf. Table 1). The pairs that look like a (x, x) form are not tested (Deletion of the matrix diagonal). The symmetrical pairs are only computed once , i.e. (a, b) is tested but (b, a) is not. It is therefore a right upper triangular matrix. The computation made on each pair consists in determining if two bounding volumes overlap or not. As AABBs are used, each object has 6 values: its minimal and maximal positions on the 3 axis. To detect overlaps, an algorithm is used to determine if one of the maximum borders of an object is inferior to the minimum border of the other object on the same axis. If so, it means that there is a maxmin space between both bounding volumes and there is no collision.

The GPU Mapping Function

In order to parallelize computations on GPU, each thread is in charge of one object pair. For a n-body simulation, n 2 -n 2 threads will be created. It is necessary to supply in every thread the pair of objects which it has to compute. Each object consists in 2 vectors (min and max) and each of them consist in three values (x, y and z) so 12 values must be transmitted to each of the threads. The accuracy of these values is essential to ensure a consistent result for the presence or absence of a collision, it is thus essential that each of these 12 values must be coded on 4 bytes (floating point number). 12 * 4 = 48 bytes have to be transmitted to each thread. With one or ten thousand of objects, the transfer size would be respectively 23MB or 2.3GB. It is quite obvious that the transfer size quickly grows not linearly.

Each GPU thread has its own unique ID, the thread 0 deals with the pair 0 until the thread n 2 -n 2 dealing with the ( n 2 -n 2 ) th pair. As the CPU-GPU transfer occurs at each time step, the data should not be as big as shown previously. The main goal is to minimize the transfer size while preserving the quality of the information. In the previous example, each object was transmitted n -1 times to the GPU at each time step, it is possible to consider sending each object only once per step that would significantly reduce the transfer size. For a simulation with 1K objects, the size would be: 23KB. But it is no longer possible for a thread to know which pair it has to compute. The first idea would be to allocate n 2 -n 2 locations in the GPU memory. They would correspond to the pairs' ID that each thread has to compute. Working with simulations with a fixed number of objects, it is possible to create and to store the list of the pairs' identifiers during an initialization phase. For example, let pairs_memory be the allocated memory containing pairs of objects' ID. The thread k accesses on the k th location of pairs_memory and reads two values ob ject1 ID and ob ject2 ID . These two IDs are then used to localize the objects in the global memory in order to get their bounding volume's coordinates. In this case, each thread does two global memory accesses to determine which objects it has to compute, followed by two other global memory accesses to get bounding volume's coordinates. Each thread performs four memory accesses before performing its computations.

To reduce the number of memory accesses, a novel approach is proposed to determine the pair (a, b) that a thread has to compute without any memory access only with the number of objects n and the identifier of the thread k.

Example: In a simulation with 1000 objects, the pairs to handle are: (1,2),(1,3),(1,4), ... , (998,1000),(999,1000). The goal is to find a method to determine, without memory access, that the thread 71452 has to handle the pair (75,303) and the thread 46108, the pair (48,285).

1 2 3 ... n-1 n 1 X 0 1 ... n -3 n -2 2 X X n-1 ... 2n -5 2n -4 3 X X X ... 3n -8 3n -7 ... ... ... ... ... ... ... n-1 X X X ... X n 2 -n 2 n X X X ... X X
Table 1: The pairs matrix in a simulation of n objects.

Linear Function of Access

As per Table 1, it is possible to determine the objects identifiers for a given pair number, using computation loops based on mathematical series. But more the identifier of the pair is high and more the number of computation loops will be important. If each GPU kernel function has a different computation time due to an evolutionary number of loops, the c The Eurographics Association 2012.

global computation time will be impacted by the slowest thread. The computation time in GPU kernels has to be homogeneous regardless of the pair identifier to handle.

Turn the Index Into Unique Pair

The numbering of the matrix rows and columns is done from 1 to n (from the top left to the bottom left). Let a be the index of the row and b the column one. In this case, the first line contains (n -1) numbers, the second one (n -2), ..., the k th (nk),... and the n th 0. From line 1 to line a -1 there is:

Sa = (n -1) + (n -2) + ... + (n -(a -1)) elements
As the pairs numbering starts at 0, Sa is the first number of the line a (0, n -1, ..., n 2 -n 2 on Table 1). Sa is a series, the operation can be simplified:

Sa = (a-1)(2n-a) 2
As this first number is in the column a + 1, the number that is in the column b is thus Sa + b -(a + 1). With the identifier k of the thread and the number of objects n, the a and b values that correspond to objects to compute have to be determined. The inequality to solve is the following:

Sa ≤ k ≤ S a+1 (a-1)(2n-a) 2 ≤ k ≤ ((a+1)-1)(2n-(a+1)) 2
But, to solve Sa ≤ k ≤ S a+1 , it is not necessary to use both inequalities because it is enough to say that:

• a is the larger integer such as Sa ≤ k • or that a is the smaller integer such as k ≤ S a+1 Which, in both cases, leads to seek for the real solution of the equation and then to take the integer part. Starting from the special case where the number of the thread corresponds to the beginning of one line of the matrix (k = Sn):

k = (a-1)(2n-a) 2 → -a 2 + a(2n + 1) -2n -2k = 0
This second degree polynomial equation has two roots:

X 1 = -(2n+1)+ √ ∆ -2 X 2 = -(2n+1)- √ ∆ -2 with ∆ = 4n 2 -4n -8k + 1.
During the roots resolution, only X 1 is located between the borders 1 and n. X 2 is out of bounds of the solution and leads to a wrong coordinate in the matrix (up to the number of objects), which has no meaning. X 1 needs to be solved in order to get the y coordinate of the matrix. The integer part of the root X 1 enables to get the matrix coordinate.

a = -(2n + 1) + √ 4n 2 -4n -8k + 1 -2 (1)
To get b, the following equation needs to be solved:

b = (a + 1) + k -Sa = (a + 1) + k - (a -1)(2n -a) 2 (2)
Figure 2: General overview of our approach that reduces the GPU global memory accesses done by the threads.

where Sa is the first value of line a. With this polynomial resolution, it is now possible for each thread to determine (using n (number of objects) and its index k) the two objects it has to deal. It is no longer necessary to allocate the pairs in the GPU memory, which saves four memory accesses per thread and a fairly substantial area of memory. The only data to transfer at each time step is the two Min and Max vectors of each object. Which, for a n objects simulation is n * 2 * 3 * 4 = 24 * n bytes. With 10K objects, the size of data to transfer is 234Kb and the transfer can be performed in less than 50 µs (0.050 ms) on recent GPU architectures. This time is almost negligible compared to our problem.

Comparison with Lazy Way

It could have been possible to use a simpler way to retrieve the pairs with the IDs that consists in generating pairs without computing any square roots and checking if they are valid. The algorithm is for a given ID k and N objects, a = int(k/N) and b = k%N. Then the result is tested before handing this thread ID number back: if a ≥ b, then this is an invalid ID that has generated an already existing pair. The calculations are stopped for this thread. While if a < b, then the pair is valid and the calculation on the potential overlap of objects is pursued. The point is, for a large N, about half of the IDs are thrown away in the 0 to n 2 -n 2 list of numbers and pass the other half on. It is therefore necessary to test n 2n IDs at the cost of an additional "if", integer divide and modulo operation. When comparing this approach on sequential CPU with the previously detailed method, it quickly appears that, in a sequential CPU context, it is better suited than the previous approach. In spite of the twice higher number of loops, the atomic functions are simpler and without computing the square root, the algorithm is four times faster when running both on sequential CPU. But with 10K objects, it takes 4.51s to compute on CPU... Thus the parallel optimization becomes necessary. When comparing both techniques on GPU, the results are inversed. The Table 2 shows results obtained on a Nvidia Quadro FX 3600M. Performing twice as more computations in parallel on GPU, even if the atomic operations are simpler, may be less efficient. That shows the inconsistency of the lazy way in a parallel context. The comparison of our approach and the use of no square root to determine the objects pair with an ID.

Square Root Optimization

Computing the square root of a floating point number is essential in many 3D applications. The new formulas used by the threads that have been proposed use square roots. A brief overview on the study on different possible optimizations for the square root computation is presented. Square root functions are very time consuming on most CPU even with particular instructions. That is why the GPU have dedicated hardware to perform such functions faster. But it still remains an important time consuming operation even on GPU. In 2005, id Software(www.idsoftware.com) released the source code of the Quake III engine under GPL license.

In this code, two methods were found by engineers and assigned to the programmer John Carmack (Few investigations showed that the code had older roots in the hardware and software side of the computer graphics community). One of these functions is designed to compute the square root of a floating point number without any loop, only using linear elementary computations. This function is based on a very rapid Newton's-method-based approximation to the square root. The C code was essentially (without original comments): This method does not contain any loop thanks to the use of the "magic number" 0x5f3759df. This hexadecimal constant is used as a first approximation and significantly reduces the iterations number needed to obtain a satisfactory approximation. This code is based on the fact that the IEEE 754-2008 floating point format approximates base-2 logarithm. The Carmack's constant was then studied and improved, first by Chris Lomont [START_REF] Lomont | Fast Inverse Square Root[END_REF] followed by Charles McEniry [McE07]. Lomont has proposed a new constant 0x5f375a86 which appears to perform slightly better than the previous one. McEniry obtained the same constant that Lomont determined but he has proved that a constant R can be used to approximate the integer value of the inverse square root of a floating point number but he has not proved that this constant R assumes the value used in the code itself.

Lomont and McEniry reported that this square root approximation can lead to errors. As shown previously, in both formulas 1 and 2, the square root is only used with integers and the decimal part of the value is useless. In preliminary tests it has been found that for very high values of n (number of objects), errors persist and lead in our case to a wrong result. An error means that two threads with a consecutive identifier will compute the same pair which means that one pair will not be processed. To avoid getting errors in the computation of the pairs of objects, we analyzed them and successfully manage to characterize them. The square root in the GPU kernel function is √ 4n 2 -4n + 8k + 1 with n, the number of objects in the simulation and k, the thread identifier. Errors can only occur in only two different situations of bad approximation of the square root. These two situations give a wrong pair (a, b) that looks like:

• b > Nb_ob j • or a > b
In both cases, there is a way to fix it. In the first case, the rule to apply to get the exact pairs coordinates is:

1 if(b > Nb_obj){ 2 a++; 3 b = a+(b-Nb_obj);
Whereas in the second case, the rule to apply is:

1 if(a >= b){ 2 b = nb_obj -(a -b); 3 a--;}
The next Table 3 shows the % of errors of the Carmack's and Lomont's approximation without correcting errors. Two different tests were performed with respectively two and three accuracy steps. An error occurs when the lower integral part of the computed square root and the approximated square root are different. The percentage of error with 4 accuracy steps is exactly the same than the 3 step's one. From 1K to 100M tests, a 0% of error is obtained, when applying the conditions (cf. above) after the computation of the square root approximation to fix the errors leading to bad pairs. It is now possible to use the method of the approximation of the square root in the kernel function being sure of the consistency of the result. The Carmack's and Lomont's square root approximations outperform the standard CUDA sqrt. The difference is not as significant as expected but is still a time reduction. On average, the percentage of obtained computation time is 13.7%.

accuracy steps 3 accuracy steps

Output Transfer Optimization

It is also essential to optimize as much as possible the output transfer time (from the GPU to the CPU). This step is made at each time step so shorter it is and better it is. Each thread is in charge of determining if two bounding volumes of objects are in collision or not. The answer "yes" or "no" can be represented with only one bit (1 or 0). For n objects in the simulation, there are n 2 -n 2 pairs to compute, as many bits have just to be allocated in the memory. If the bit number X is 1, it means that the bounding volumes of the objects of the pair number X are in collision. In the CUDA/C++ programming context, it is not possible to directly manipulate 1 bit parameters or attributes. Even a character or a Boolean value is coded on 1 byte (8 bits). Allocate a character per pair would be 8 times more important than just one bit and the transfer time would also be longer.

The 8 bits of one character are used to store the collision test result of 8 pairs. The GPU memory allocation can only be made in bytes so the number of needed bytes is nb_Bytes = N/8 . Thus, (N/8 + (1 or 0)) bytes are allocated to store results. When the collision test is positive, the thread creates a character equal to 1 and then proceeds to a binary shift of (ID%8) bits into the character (ID is the thread identifier). Then it has to write this character in the (ID/8) th memory location (previously allocated). But 8 threads will write to the same memory location so the result information has to be maintained in the same character. To do so, before writing on the memory, each thread has to read the (ID/8) th memory location and apply an OR logical function of the two binary strings. This elementary operation enables to do the union of the results from the 8 threads in the same character.

The read of this binary string transmitted to the CPU is But if 8 threads write at the same memory location, they cannot write it in parallel. These accesses are serialized. The differences between 8 parallel writings in the global GPU memory (using a 1-byte data type) and 8 serialized writings (using a 1-bit data type) were evaluated. Table 5 presents the results of these performance measurements from 100 to 7500 objects in the simulation. The kernel execution time and the output transfer time were measured. The input transfer was not controlled because the choice between 1-bit or 1-byte does not impact it. The Figure 3 shows the difference between both configurations. The kernel execution time is exactly the same between using 1-byte or 1-bit so between parallel or serialized accesses on the global memory of the GPU. While the output transfer time is significantly longer when using 1-byte per pair. Table 5: Kernel and output with 1-byte or 1-bit per pair.

Intermediate Results

Intermediate results obtained by the new broad phase algorithm based on a semi-brute force approach are presented. The graph shown on Figure 4 summarizes the evolution of the GPU kernel time according to the number of objects and the changes in the input/output data size. These results show c The Eurographics Association 2012.

Intermediate Comparison

For the simulations with less than 3000 objects, the new model outperforms the latest approaches. If we compare the simulation results from 2500 objects with Liu et al.

[LHLK10], our approach is 20% faster. Their approach is based on a standard GPU SaP coupled with an algorithm to choose the best projection axis. Compared with the approach of Pabst et al.

[PKS10] and Kim et al. [KHH * 09] that test their hybrid model with the n-body benchmark of the UNC (http ://gamma.cs.unc.edu/), our broad phase algorithm is, respectively, 5 and 3 times faster. Kim's approach is based on a tree parallel traversal and the Pabst's one is based on a uniform spatial subdivision coupled with a non-optimized for triangles brute force SaP.

Intermediate Synthesis

The initial objective was to show that a very simple algorithm (brute force) properly optimized, could be very powerful. It is partly achieved. The new proposed algorithm is more efficient than existing approaches for scenes less than 3000 objects but as expected, given the combinatorial nature of the algorithm with a quadratic complexity in O(n 2 ), the results collapse when the number of objects increases. The objective is clearly the simulation of large-scale virtual environments, so it is necessary to revise this model to add conditions to prune the number of pairs to test. To do so this combinatorial effect is removed by coupling this approach with a spatial subdivision technique. The optimizations on the size and transfer time and the GPU mapping function that enables the threads to compute their current pair without any memory access are still used with the spatial subdivision.

Spatial Subdivision

To associate the broad phase algorithm with a spatial subdivision has become necessary to continue the reduction of the computation time by avoiding useless calculations of overlapping [START_REF] Liu | Real-time collision culling of a million bodies on graphics processing units[END_REF]. To further reduce the computation time and break the quadratic complexity, the new semi-brute force approach is associated to spatial subdivision. Indeed, a 3D grid is proposed to reduce the number of pairs to compute and to locally perform our brute force approach. This approach can be compared to the Multi-SAP concept, first proposed and thought by Erwin Coumans and Pierre Terdiman [START_REF] Terdiman | Sweep-and-prune: Multi-sap[END_REF] but still marked as experimental. It has been previously shown that the strategy which consists only in applying the brute-force method was adapted well to "small environments" but totally unsuitable for larger environments. It is therefore useful to widely prune the number of pairs. The grid-type structure is perfectly suited to this type of problem. This grid is set at the beginning of the simulation and then updated at each time step of the simulation. The number of cells and their size depend on the number of objects and their size. In the following, the establishment and the update of the grid structure are presented. The construction of the grid is done in two steps: (1) determining number and size of the cells and (2) putting objects in the grid cells.

Cells Characteristics

To ensure uniformity and consistency, the density of objects within the environment is considered in order to determine the characteristics of the grid cells. The dimensions of the environment and the number of objects are used in the construction. First, the objects density is evaluated in the environment, dividing the number of objects by the cubic volume of the environment. Then, a X value is set in order to determine the ideal volume that a cell should have to contain X objects. X is divided by the density of the environment for the optimal volume of cell. Then the cubic root of this volume is computed in order to determine the size of the side of a cell. In this approach the grid has cubic cells. After obtaining the size of a cell, the number of cells per axis can be determined allowing refining the cell size for a perfect fit to the size of the environment. The value X varies according to simulations.

Grid Fill in

The goal is then to link the objects to their respective cells.

To do it the GPU is used to massively parallelize the computations. Each calculation has to determine in which grid cell is an object. The computation for each object is totally independent of the other computations, they can be performed in parallel on the GPU. The cells to which belong objects are determined. The cells' identifiers of one object are computed using the position of the object and the maximum and minimum bounds of it in order to loop through all relevant cells c The Eurographics Association 2012.

and mark these as occupied by that object. At the end of this step, each grid cell is composed of a list of objects' identifiers. The semi-brute force algorithm presented previously is then locally applied according to the cells of the grid. The occupancy of the GPU is maximal to simultaneously run a maximum of cells in parallel.

Grid Update

This update step is essential when using a grid-like structure.

The goal of the grid being roughly to prune candidate pairs to compute during the brute force algorithm, the update of its cells must therefore be fast. This step is also performed on the GPU. The GPU kernel function computes the cell identifier of each object in parallel. It thus takes into account their new positions and assigns them to the corresponding cells.

Implementation and Performance

In this section, the performance measurements of the new parallel broad phase algorithm is presented. In all tests, all objects are moving in order to evaluate the performance in the most extreme situations. The algorithm was developed an tested on a Intel Core-Duo CPU X7900 of 2,8 GHz with 3 Gb of memory with a Nvidia Quadro FX 3600M. The operating system is Windows XP 32 bits. The second column of the Table 6 corresponds to the step of the update of the bounding volumes (here AABBs). The third column is the time spent to update the grid on the GPU. The fourth column gives the computation time of the brute force Sweep and Prune on GPU. This time includes the time of the input and output transfers. The Table 7 shows a comparison of our approach with or without the use of the 3D grid. The use of the grid is useless in "small" environments, the execution time is better when not using it. These small environments are those composed of less than 1000 objects. Beyond these 1000 objects, using a 3D grid leads to obtain more efficient results.

Comparison and Analysis

In the following, a comparison is made between our approach and existing and most used algorithms from the recent literature.

The approach is compared with the last releases of two other broad phase algorithms of the Bullet library [Bul], the 3D incremental sweep and prune (SAP) [CLMP95] (a 64 bits version has been developed to handle more objects) and the btDbvt algorithm which is based on dynamic AABB trees. Two different tests were performed. The first one is to determine the impact of the objects' movement on the computation time because many broad phase algorithms are optimized when many objects are static. Thus several performance measurements with varying percentage of moving objects are run. Results (cf. Figure 6) show that the approach is not impacted by the fact that the objects are static or dynamic. Unlike other approaches, whatever the percentage of moving object, our computation time remains constant. The second test is to assess the robustness of the algorithms according to an important number of objects. The algorithms is evaluated from one thousand to one million of objects, 100% moving. These tests were made with the same scenario where an object column of the same size falls on to a plan (cf. Figure 1) with the same Nvidia Quadro FX 3600M and a Intel Core2 CPU X7900 @ 2.8GHz. Results are illustrated on the Figure 5. The algorithm significantly outperforms the two existing methods with large-scale environments.

It is a bit difficult to compare results with the literature because there are only few approaches that have tried to compute the collision detection with more than a million of objects only using one GPU. It is also difficult to compare with other approaches because the GPU are different. There are several works studying the particles-based simulations with millions of elements but it goes away from our domain. The works of Tracy and al. [START_REF] Tracy | Efficient large-scale sweep and prune methods with AABB insertion and removal[END_REF] focused on the improvement of the SaP on the CPU and their results present only environments with a low percentage of moving c The Eurographics Association 2012. object. The unique result they indicate with 100% of moving objects consists of 3000 Objects. Their algorithm takes 4ms whereas ours takes 0,46ms (including transfer time). Liu and al. [START_REF] Liu | Real-time collision culling of a million bodies on graphics processing units[END_REF] have tested their approach on an environment consists of 960k spheres of three different sizes. However, they do not give the exact numerical values of the computation time of their algorithm, which makes the comparison difficult. The only comparable data are those from three measurements done on three different numbers of spheres: 64K, 128K and 256K. Their approach takes respectively 13.71ms, 48.16ms and 241ms to apply the algorithm during the broad phase stage. Our algorithm takes respectively 15.67ms, 30.68ms and 50.33ms for the same conditions of simulation. We are thus slightly below the performances for the environment of 64K spheres, 1.5 times faster for the one with 128K spheres and 5 times faster for the last one. Other approaches such as Pabst and al. [PKS10] and Lauterbach and al. [START_REF] Lauterbach | gproximity: Hierarchical gpu-based operations for collision and distance queries[END_REF] are only tested on small n-body environments (300 Objects). The computation times of their algorithms of broad and narrow phase are not differentiated, only the total time is given. For a simulation of 300 spheres and 4 cones, the approach of Pabst takes 128.1ms to make the calculations in discrete mode and 184.8ms in continuous mode. It is always very difficult to estimate time made by the broad phase in an approach. Taking a broad phase which would represent 10%, 5%, 2% or 1 % of the total execution time, theirs would be respectively 12.81ms, 6.40ms, 2.56ms and 1.28ms. If we cross our results, our algorithm of SaP-Grid on GPU for the same type of environment takes exactly 0.41ms. Even for a broad phase which would represent only 1% of their total time, our algorithm would be three times faster.

Conclusion and Future Work

A novel and efficient GPU solution has been presented to cull non-colliding object pairs. This new algorithm is based on a grid-based spatial technique and a new GPU mapping function that enables GPU threads to turn their ID into a unique objects pair without any global memory access. This paper brings new directions for the collision detection optimization. It has presented a different direction from recent collision detection works by going back to a semi brute-force approach. The novel approach enables to reduce the global memory accesses done by threads and leads to a faster execution time. To break definitively with the combinatorial of the algorithm, the approach is coupled with spatial subdivision. The kernel execution time is further reduced with the use of the square root approximation where a first characterization of the errors is presented. It enables the fixing of wrong computations leading to incorrect object pairs. The logical optimization with binary masks of the output data transfer has minimized the size and the time of the GPU-CPU transfer. Numerical and comparative results showed that this new algorithm enables to cull objects in less than 25ms with more than 100K moving objects on a Quadro FX 3600M. The comparisons made with the last releases of the well-known and most-used algorithms showed that this new model enables to further reduce the collision detection time. Compared to the literature, this model is, in most cases, as or faster than other ones.

There are many ways for future work. The GPU mapping function is fully scalable on multi-GPU platform so it could be interesting to evaluate its performance on bigger platform. The next step of our work will focus on the second stage (narrow phase) of the collision detection pipeline in order to know if it is possible to apply such a parallel spatial brute force approach. Using this approach not only on bounding volumes but also on triangles or polygons could be an interesting way. The grid establishment and update can be improved and better mapped on to the GPU architecture. A study on the optimal number of cells into the grid is crucial and also on the cell size. It is now essential to take into account the run-time architecture to speed-up collision detection algorithms. We think that the key challenge is in the establishment of mapping models between the algorithms and the architectures.

Figure 1 :

 1 Figure 1: 250K cubes, 100K spheres or 75K rectangles. Using the GPU mapping function, culling non-colliding pairs is performed in less than 25ms on a Nvidia Quadro FX3600M in the spheres environment.

  increase accuracy(6) 8 y = y * (1.5f -(nb_half * y * y)); 9 return number * y;}

Figure 3 :

 3 Figure 3: The output transfer and kernel times between using 1-byte or 1-bit per objects pair. Time was measured in milliseconds from 100 to 7500 objects.

Figure 4 :

 4 Figure 4: GPU Kernel time, input and output data size of our model from 100 to 5000 objects.

Figure 6 :

 6 Figure 6: Evaluation of the objects movement's impact on the computation time.

Figure 5 :

 5 Figure 5: Tests performed from one thousand to one million of moving objects with three different algorithms. All objects are the same size. The computation time remains interactive (<25ms) in an environment composed of 100,000 objects with a Nvidia Quadro FX 3600M.

  The additional difficulty of comparison is that Lauterbach et al. have carried out their tests on a Nvidia GTX 285 while Pabst et al. have performed on a Nvidia GTX 295 and we have done ours on a Nvidia Quadro FX 3600M which makes it very tricky to compare.

Table 2 :

 2

	c The Eurographics Association 2012.

Table 3 :

 3 Percentage of errors of the Carmack's and Lomont's square root approximation when using two or three accuracy steps without correcting errors.

	Nb Tests Carmack	Lomont Carmack	Lomont
	1K	2.9%	2.9%	0.7%	0.7%
	10K	0.89%	0.89%	0.12%	0.12%
	100K	0.27%	0.27%	0.035%	0.035%
	1M	0.1800% 0.1803% 0.0110% 0.0111%
	10M	0.3681% 0.3685% 0.0040% 0.0057%
	100M	0.9847% 0.9855% 0.0134% 0.0142%

Table 4 :

 4 Table 4 presents the results obtained by comparing normal CUDA square root (sqrt()) with the Carmack's and Lomont's approximation techniques. GPU Computation time (in ms) of the Carmack's and Lomont's approximation compared to the CUDA square root function. The tests were made from 1 million to 1 billion sqrts to compute. The value in bold is the shortest computation time. The γ value correspond to the % of computation time reduction obtained with the square root approximation.

	Nb tests Carmack	Lomont	CUDA	γ
	1M	0.042ms 0.042ms 0.048ms 12.6%
	10M	0.409ms 0.409ms 0.478ms 14.7%
	100M	4.102ms 4.102ms 4.722ms 13.2%
	1B	41.01ms 41.01ms 48.01ms 14.6%

c The Eurographics Association 2012.

The

Table 6 :

 6 Computation timeα: AABB Updateβ: GPU Grid Updateδ: GPU Brute Force SAP

	Nb Objects α(ms)	β(ms)	δ(ms)	Total
	1K	0.03	0.13	0.09	0.25ms
	10K	0.31	1.45	0.59	2.35ms
	100K	3.38	11.12	13.41	27.91ms
	1M	14.26 134.91 142.72 291.89ms
	Nb Objects Without Grid With Grid
	100		0.046ms	1.020ms
	1K		0.872ms	1.840ms
	10K		31.450ms	3.257ms

Table 7 :

 7 Computation time with/without the grid use.

c The Eurographics Association 2012.

Acknowledgments

The authors want to thank Alison Day and Colin Moore (Wechsler-Reya Lab) for their help in the English review process. This research was supported by Bretagne Region under project GriRV N o 4295.