{'fraction_non_alphanumeric': 0.032516784612593, 'fraction_numerical': 0.014044637996733805, 'mean_word_length': 4.740833333333334, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ALTAI: Computational code for the simulations of TeV air showers as observed with the ground-based imaging atmospheric Cherenkov telescopes

Ground-based atmospheric Cherenkov telescopes are proven to be effective instruments for observations of very high energy (VHE) gamma-radiation from celestial objects. For effective use of such technique one needs detailed Monte Carlo simulations of gamma-ray- and proton/nuclei-induced air showers in Earth atmosphere. Here we discuss in detail the algorithms used in the numerical code ALTAI developed particularly for the simulations of Cherenkov light emission from air showers of energy below 50 TeV. The specific scheme of sampling the charged particle transport in the atmosphere allows the performance of very fast and accurate simulations used for interpretation of the VHE gamma-ray observations.


Introduction
Recent exciting detections and observations of TeV γ-ray emission from a number of galactic and extragalactic objects (Ong 1998, Catanese & Weekes 1998) have shown the high performance of currently operating imaging airČerenkov telescopes (IACTs). Several projects for future detectors have been proposed lately. The major forthcoming instruments, such as CANGAROO IV, H.E.S.S., MAGIC and VERITAS, will have significantly better sensitivity to γ-ray fluxes in the energy range from 50 GeV up to 50 TeV.
For a lack of a collimated test beam of TeV γ-ray photons, the ground-basedČerenkov detectors heavily rely on the Monte Carlo simulations of theČerenkov light emission from air showers which are used to understand the performance of detector. Basic parameters of the instrument (detection area, angular and energy resolution, efficiency of cosmic ray rejection, etc) can be derived from the simulations. The crucial point is a precise measurement of the primary shower energy. For that one should include properly into the simulations all processes ofČerenkov light emission in the air shower as well as photon propagation on the way from the emitting shower particle to the telescope camera. Measurements of γ-ray fluxes, energy spectra, upper limits strongly depend on the absolute energy calibration of a telescope.
Here we give a description of ALTAI 2 computational code developed for detailed Monte Carlo simulations of theČerenkov light emission in TeV air showers. Among the other existing codes intended for such simulations, MOCCA (Hillas 1979), KASCADE (Kertzman & Sembroski 1994), CORSIKA (Heck et al. 1997), this code has a distinctive advantage -it's rather high speed of shower simulations due to a specific algorithm used for processing the multiple scattering of charged shower particles. This approach does not consume a lot of CPU time and allows to perform fast and accurate simulations. Together with the additional routine recently developed for the simulations of telescope response (Hemberger 1998) the ALTAI code was extensively used for production of a standard Monte Carlo database used in the HEGRA (High Energy Gamma Ray Astronomy) VHE γ-ray experiment (Konopelko The ALTAI code consists of two major programs which simulate the development of the electromagnetic (EMCCS) (hereafter we put in brackets the name of the corresponding subroutine of the code) and proton-nuclei (STEPAD, MULTIP, XPI) cascade in Earth atmosphere. We discuss in detail the procedure of simulating the electromagnetic and hadron-nuclei cascade in Sections 2 and 3, respectively. Section 2 also describes the algorithms ofČerenkov light emission by the shower charged particle. Section 4 deals with nucleus-nucleus interactions. In Section 5 we review the results of theČerenkov light simulations using ALTAI code compared with other relevant simulations and observational data.

Electromagnetic cascade
Here we overview the part of the code intended for the generation of an electromagnetic cascade in the atmosphere (EMCCS). Note that at present all cross-sections for electron and photon interactions are established and are well described in detail in the relevant energy range. An electromagnetic part of the code accounts for the following interaction processes: electronpositron pair production, Compton scattering for the primary photons and bremsstrahlung, ionization losses, Coulomb scattering for the primary electrons and positrons.
The cross-section of e − , e + pair production by primary photon was calculated according to the Bethe-Heitler formula taken from Motz et al. (1964) (MPRH). The emission angle of e − , e + components, defined with respect to the direction of initial photon, was sampled using the Bethe distribution (Motz et al. 1964). We used the Klein-Nishina formula (see e.g., Gaisser 1990) for the cross-section of Compton photon scattering (COMPH). The dependence of photon crosssections on energy and the atomic number of a medium was derived from data taken from Storm & Israel (1973) and Hubbel (1969) We simulated the bremsstrahlung interaction process (MRDH) according to the Bethe-Heitler formula (see Koch & Motz, 1959). The angle of emitted photon was sampled from the Schieff distribution (see Koch & Motz, 1959). The Bethe-Bloch formula, with the corrections for the density effect (Sternheimer 1952), was used for simulation of the mean ionization losses (MIONH) whereas for the differential cross-section of e − , e + ionization collisions we used the Möller formula (Möller 1932). Note that in the electromagnetic shower positrons were treated similarly to electrons.

Multiple scattering
The relativistic electron suffers an enormous number of interactions along its path length in the matter. Apparently, a direct simulation of all interactions is not possible. That is why all numerical codes for the shower simulation use a specific technique for grouping electron interactions (see for review Berger 1963, Akkerman 1991, Gaisser 1990. In this approach one can simulate only the so-called catastrophic e − , e + interactions. These interactions lead to emission of photon/electron of a relatively high energy, E > Q, which exceeds the intermediate energy threshold of Q. The same approach was used for the simulation of the e − , e + ionization process (Q = Q I ) and bremsstrahlung (Q = Q R ). The emitted electron must have a kinetic energy above a certain intermediate energy threshold, T > Q I .
A random path length between neighboring catastrophic collisions was divided into small seg-ments of ∆l ≤ 5 gr/cm 2 . We sampled the phase coordinates of a charged particle at the end of the segment ∆l according to the cumulative effect of all low energy collisions along the segment (SGMEH, SGMQH). Thus we simulated at the end of each segment the loss of the e − , e + kinetic energy, ∆T = T 1 − T 2 , the scattering angle, Θ = arccos( Ω 1 Ω 2 ), the azimuth angle of scattered particle, ϕ, the longitudinal Z and lateral ρ displacement of initial particle while passing over the segment (see Figure 1). Appropriate transport equations were used in order to derive the probability distributions for the phase coordinates at the end of a segment. We have obtained the analytical solutions for these distributions.

Energy losses
The energy losses of a charged particle, ∆E, in the segment, ∆l, were sampled according to the Landau-Vavilov formula (Landau 1944, Vavilov 1957). Plyasheshnikov & Kolchuzkin (1975) have tabulated this formula for a specific conditions of shower simulations. The distribution of the energy losses at the end of a multiple-scattering segment was described as follows: where is the mean energy loss per unit path length due to the inelastic collisions of a charged particle by production of low energy secondaries.
W I (Q), W R (Q) are the differential cross-sections for the ionization (I) and bremsstrahlung (R) interactions.
The standard inverse function technique was applied in order to simulate the energy loss ∆T according to the two-dimensional distribution f (µ, λ) (Eqn. (1)). For that one needs to compute the tables of the functionλ =λ(µ, α), which is inverse to the integral distribution These tables can be found, e.g., in Akimov et.al.
(1981) whereλ =λ(µ, α) was tabulated over a two-dimensional lattice with the steps of 0.025 and 0.01 over α and µ, respectively. First, it is necessary to generate the random number α uniformly distributed in the interval (0,1) and calculate parameter µ according to Eqn. (1). After that one should interpolate λ = λ(µ, α) using the above mentioned tables and finally determine ∆T as Formula (5) was derived using two asymptotics of initial distribution (1) for µ = 0 and µ = ∞, which were sewn up at µ = 0.1. The subroutine ELSH makes calculations according to these algorithms.

Angular deflection
We simulated the angle of multiple scattering using the Moliere theory (Moliere 1948, Bethe 1953). In addition we have improved the Moliere distribution by taking into account the energy losses at the multiple scattering segment. This distribution may be described as follows: Parameter B was determined by resolving the transcendental equation ). Functions f n (n = 0, 1, 2) were calculated as Quantities χ c and χ a are closely related to the differential cross-section for Coulomb scattering (see formula (10) for definition of the Coulomb cross-section). For a fixed segment length ∆l and kinetic energy of a particle T one can determine in a few iterations the parameter B, which defines the shape of the distribution P(Θ). Such calculations were done by use of subroutine ANGMH.

Space displacement
At the multiple scattering segments we simulated the longitudinal displacement Z of a charged particle. For that we used the Yang-Spenser distribution (Yang 1951, Spencer & Coune 1962). This distribution was adapted for the air shower simulations by Plyasheshnikov & Kolchuzhkin (1975). This distribution may be represented by the following expression where Here W c (T, θ) is a cross-section of the Coulomb scattering which determines parameters χ a and χ c .
To simulate the longitudinal displacement Z we use approach similar to that used in Section 2.1.1 for simulation of the energy loss ∆T . It is based on the interpolation of the two-dimensional function ξ = ξ(ν, α), which is the inverse function to the integral distribution One can find these tables, e.g., in Akimov et.al (1981). This approach includes (i) sampling of the random number α, (ii) calculation of ν on the basis of quantities T 1 and ∆l, (iii) determination of ξ on the basis of the two-dimensional interpolation using the above mentioned tables and, finally, calculation of Z according to the formula For the lateral displacement of a charged particles at the multiple scattering segment we defined the corresponding probability distribution using Fermi formula (see e.g., Kolchuzhkin & Plyasheshnikov 1975). This distribution was derived by solving the transport equations in the Focker-Planck approximation where the collision integral corresponding to the Coulomb scattering was replaced by a second order differential operator (see e.g., Kolchuzhkin & Uchaikin 1978). For the fixed angle of multiple-scattering Θ = (Θ x , Θ y ) = (Θ cos ϕ, Θ sin ϕ) this distribution was as follows where the parameter γ = η∆l defines the width of the distribution. As was shown by Kolchuzhkin & Plyasheshnikov (1975) more accurate numerical solution of the transport equations reduces by a factor of 1.5 the width of the final distribution. This difference may be corrected by introducing another definition for the parameter γ: whereB is determined fromB − lnB = ln((χ c /χ a ) 2 /1.80). The complete procedure for charged particle transport in multi-dimensional phase space of energy, angular and space coordinates was defined by a few parameters: two threshold energies Q I , Q R , and the length of a multiple-scattering segment ∆l. The extensive test calculations using the ALTAI code revealed the optimum values of parameters which allow to perform rather fast shower simulations without introducing systematic errors. Thus we used the following parameters for the calculational procedure: MeV is the threshold energy forČerenkov light emission by cascade electrons in air.
Apparently the length of a multiple scattering segment, ∆l, is one of the basic parameters of this method. By use of rather small segments one can reduce systematic error introduced by approximations of analytical solutions for the phase transformations at the multiple scattering segment. On the other hand this may slow down the procedure of shower simulation. We found that the optimum length of the multiple-scattering segment is within the range of ∼ 1 − 5 gr/cm 2 . Very accurate analytical solutions of a multiple-scattering process for a e − , e + transport allow us to use in calculations such segment length without introducing systematic errors in three-dimensional shower development.
In comparison to other computational codes we used more accurate distributions derived analytically from the multiple scattering theory. On the contrary the standard approach used for example in EGS-IV code does not include fluctuations of the energy losses at the multiple-scattering segment. In addition, in the EGS-IV code the lateral displacement of charged particle at the multiplescattering segment is taken as ρ = 0 and correspondingly the longitudinal displacement is of z = ∆l. All this necessitates a small segment length and makes significantly more time consuming the shower simulations. Thus for the simulations ofČerenkov light from the TeV air shower ALTAI code is faster by as much as a few times compared with the EGS-IV.

Emission ofČerenkov Light
A charged particle (e, µ) in an air shower can emitČerenkov light when its energy exceeds a certain threshold energy E > E th . The threshold where n is the refraction index in air, m o is the particle mass. Thus at sea level the energy threshold is E th ≃20 MeV for electron and ≃4 GeV for muons. To describe the altitude dependence of the refraction index we used the following expression (Beliaev et al. 1980) where η o = 2.9 · 10 −4 , ρ o = 1.22 · 10 −3 g/cm 3 , and ρ(H) is the air density at a height H above the sea. The model of a standard atmosphere (Elterman 1968) was used in the simulations. According to Frank & Tamm (1937) the mean number ofČerenkon photons emitted in a 1 cm pathlength by electron is described as where α is a fine structure constant α = 1/137 and θ c ≃ 1 o is the opening angle of theČerenkov light cone. The spectral region of emission is defined by wavelengths λ 1 and λ 2 . In simulations we included sampling of theČerenkov light attenuation in the atmosphere due to the Raleigh scattering, ozone and aerosol absorption. The crosssections for these processes were calculated using the data of Driscoll & Vaughan (1978).
In the shower simulation procedure (PARCHR, GSTCHR) (Konopelko 1990), we first define the total number ofČerenkov photons < n ph > emitted at the multiple scattering segment ∆l. As was discussed above we used rather small multiple scattering segments. Therefore, we can assume that allČerenkov photons are emitted from the center of the multiple scattering segment (see Figure 1). To a rather good approximation an intersection of theČerenkov light cone with the observation plane forms a circle of a radius r c which can be calculated as where Z and Z obs are the height ofČerenkov light emission and the height of the observation level above the sea, respectively. The number of photons hitting the telescope mirror can be calculated as where ∆L is length of the circle arc between points B and C (see Figure 2). p(Z, Z obs ) is a probability of theČerenkov light attenuation in the atmosphere along the way of a photon propagation. For specific response functions of the telescope camera one may calculate corresponding average number of photoelectrons as < n ph.e. >= ξ <ñ ph >, where ξ is a photonto-photoelectrons conversion efficiency. Finally one can use the Poisson distribution in order to simulate the random number of photons (photoelectrons) hitting the telescope. The photons reaching the telescope mirror were uniformly distributed over the circle arc BC (see Figure 2). In the latest version of the ALTAI code we save all, or a certain fraction, ofČerenkov photons hitting the telescope. Each photon has a set of 6 variables, υ(x, y, θ x , θ y , z o , t). z o is a height of the photon emission (in g/cm 2 ) and t is a time of photon arrival to the detector. Such database of the photon is used for detailed simulations of the telescope camera response (Hemberger 1998).

Hadron-nuclei cascade
For simulations of the hadron-nuclei cascade we used the phenomenological model of hadron interactions (Konopelko 1990) which, in the most part, is based on the available accelerator data. The energy spectra of secondary hadrons generated in a pA interactions were approximated by use of the radial scaling model (Hillas 1979). In this approach the energy spectra may be presented as where E o , E are the energies of the primary and secondary particles, respectively. Indeces p, q denote the type of primary and secondary particle (p, n, π + , π − , π o ). The basic functions are given below For the total cross-sections of inelastic hadron interactions we used the data given by Shabelski (1986,1987) in the following form GeV, E 2 = 10 4 GeV, y = ln(E/10 3 GeV ) (22) Table 1 Mean number of fragments created by different nuclei (Murzin 1988 The cross-sections from Eqn. (22) are measured in mBarn and correspond to the particle interaction with air. We assume that the mean value of atomic number for air nuclei is of 14.4. We used a simplified geometrical representation (Murzin 1988) for the total cross-section of the nucleus-nucleus interactions as follows where A 0 is the effective atomic mass of air (A o = 14.4) and b is the effective radius of nuclei overlapping zone (b = 1.17). σ N defines the cross-section of hadron interaction with air. The transverse momenta of secondary hadrons were calculated according to the distribution given by Ranft (1972) where the corresponding parameters (B,C,D) were defined for the interactions with air as described by Ranft et al. (1972). p ⊥ is measured in GeV/c. The ALTAI code was developed for the simulations in the energy range relevant for the very high energy γ-ray astronomy E ≤ 50 TeV. In this energy range the mean number of kaons (K ± , K o ) produced in hadron-nuclei interactions is very small compared with the number of emitted pions (π ± , π o ). Besides this, in this energy region kaons and pions exhibit similar properties of inelastic interactions (see, e.g., Grishin 1982). Thus in the shower simulations for this restricted energy range we can exclude the kaon production process by introducing appropriate corrections for the probabilities of pion production in the (p, n, π ± ) → X inelastic collisions.
For each inelastic hadron interaction (MUL-TIP) we sampled first the energy of the so-called leading particle (p, n, π ± ). This particle carries away the bulk of the primary energy. We assume that the energy of a leading particle E is uniformly distributed within the range (0, E 0 ) where E 0 is the energy of a primary interacting particle. Note that the type of the leading particle was selected randomly amongst proton and neutron for interaction (p, n) + A → (p, n), and amongst charged pion or neutral pion for interaction (π ± , π o ) + A → (π ± , π o ). Thus the total inelasticity coefficient was determined as κ = 1 − E/E 0 . At the next step we sampled multiple production of other secondary particles (π ± , π o ). The type of secondary particle was randomized assuming that π o pions forms, on average, one third of all pions generated in hadron interactions. The energy of secondary particle was simulated using the spectra given by Eqn. (20). We have completed the production of secondary particles when their total energy exceeds the energy of a primary particle. The energy of the last simulated particle was renormalized in order to conserve the total energy in each inelastic interaction.
The transverse momentum of the secondaries was simulated using Eqn. (24). By the renormalization of the transverse momentum of leading particles we allowed to fulfill the momentum conservation for each hadron interaction without distortion of stochastic properties of a multi-particle production process. The test calculations have shown that such an approach describes very well the initial inclusive spectra of the secondary particles if the number of secondaries is relatively high (≫ 1). Note that as almost identical algorithm was developed and tested by Barashenkov & Toneev (1972).
Note that we have included in the code the propagation of a single muon generated in the hadron-nuclei cascade due to the π ± decay (SGG,TRMUON). The emission of theČerenkov light from the muons was included in simulations according to the scheme described above.

Nucleus-nucleus interactions
We have implemented in the ALTAI code a model of independent nucleon interactions of colliding nuclei with nuclei fragmentation included (TAFRAC, TRAN, FRAG, OVERLAP). In this approach we assume that all nucleons of the projectile nucleus have the same energy determined as E 0 /A p . E 0 and A p are the energy and atomic number of projectile nucleus, respectively. All nucleons of the projectile nucleus which overlap with the target nucleus interacted independently with each other. The non-overlapping part of the projectile nucleus decayed into individual nucleons and heavier fragments. The energy of fragment with the atomic number A is defined as E f = (E 0 · A)/A p . We simulated a random number of a fragments according to the probabilities of different channels summarized in Table 1.

Comparison with other codes and data
An overview of Monte Carlo results on lateral, temporal, and angular characteristics of thě Cerenkov light in air showers of 10 GeV -1 TeV was recently given by Konopelko (1997). Most of theČerenkov light characteristics calculated with the ALTAI code reproduce rather well the results obtained with the MOCCA code (Hillas 1996). The stereoscopic observations of BL Lac object Mkn 501 in 1997 with the HEGRA system of imaging airČerenkov telescopes allowed the first measurements of the parameters of theČerenkov light images from the γ-ray-induced air showers. The HEGRA data are in excellent agreement with the results of Monte Carlo simulations obtained with ALTAI code . Using the Mkn 501 data sample comprising 38,000 γ-ray events we have tested in great detail the parameters ofČerenkov light image orientation and shape (Konopelko et al. 1999a) (see Figure 3) as Distribution of a shape parameter (mean scaled Width) ofČerenkov light images simulated with ALTAI code (histograms) and recorded with HEGRA system of IACTs (crosses).
predicted by the simulations. The HEGRA stereoscopic system of 5 IACTs was used for the measurements of the lateral distribution of theČerenkov light in the γ-rayinduced air showers. These measurements have been compared again with the Monte Carlo simulations using the ALTAI code (Aharonian et al. (1998)). The simulations using the ALTAI code reproduce very well the measured lateral distributions of theČerenkov light. Note that the measured shape of theČerenkov light lateral distributions is almost independent of the detector simulation procedure but strongly depends on the development in space of a multi-TeV γ-ray shower in the atmosphere.
We made several comparisons of the shape and size (total number of photoelectrons) of thě Cerenkov light images for the proton-and nucleiinduced air shower calculated with the ALTAI and CORSIKA codes (Heck et al. 1997 Figure 4. Distribution of image shape parameter Width for the ALTAI simulated events (histogram) and cosmic ray images recorded with a 10 m Whipple IACT (crosses) (Konopelko 1999b).
of the proton-nuclei interactions. This model is based on supercollider data and describes rather well the proton-nuclei air showers of energy below 10 TeV. In addition, the HDPM model of CORSIKA allows to perform relatively fast calculations with respect to other more modern shower generators, e.g. VENUS/QGSJET. Despite the different algorithms and schemes in the two codes the results appeared to be almost identical (Plyasheshnikov et al. 1997). The simulations with the ALTAI and CORSIKA codes well reproduce each other. Note that simulations with the ALTAI code are essentially less time consuming. The HEGRA data for the cosmic ray air showers were compared with the ALTAI simulations. Good agreement between simulations and data provided a precise measurement of the cosmic ray proton spectrum in the energy range ∼1-5 TeV (Aharonian et al. 1999b).
We show in Figure 3 the distribution of socalled mean scaled Width parameter for pri-mary γ-rays and cosmic rays extracted from the HEGRA data as well as from the Monte Carlo simulations. The simulations fit very well the data. Recently, we compared the relevant results of the ALTAI simulations with the cosmic ray data taken with a 10 m Whipple imaging airČerenkov telescope during the 1995/1996 Crab Nebula observations (Konopelko 1999b). For these simulations we have used the standard detector response functions offered by the Whipple collaboration. The simulations reproduced very well the shape of recordedČerenkov light images (see Figure 4).
As mentioned above the upper energy of the shower simulations is about 50 TeV. It is limited by the simplified phenomenological model of the hadronic cascade. Although in the case of pure electromagnetic showers one can extend simulations using the ALTAI code to much higher energies without breaking any model restrictions such extension may need the introduction of a number of changes into the code for the correct treatment of a larger number of particles.

Summary
Here we have presented a detailed description of the numeric code ALTAI, which was developed for simulations ofČerenkov light from extensive air showers. This code allows to make very fast and accurate calculations of the response of the ground-basedČerenkov detectors used in VHE γray astronomy. Although the code was designed for calculating the parameters and characteristics of the imaging airČerenkov telescopes it can be rearranged with minor changes in order to simulate the responses of showerfront sampling experiments like MILAGRO, or Tibet AS-γ.
We tested our computational code against the data taken with two currently operating detectors, the HEGRA system of imaging airČerenkov telescopes and the state-of-the-art single 10 m Whipple telescope.
These comparisons have proven the high precision of the simulations. The computational code ALTAI is an effective tool for producing the simulated data for VHE γ-ray astronomy. Note that the forthcoming imaginǧ Cerenkov detectors (CANGAROO IV, H.E.S.S., MAGIC, VERITAS) will need a large amount of simulated data. Therefore the performance of these instruments may benefit from the use of the ALTAI code.

Acknowledgments
The work for the ALTAI code was primarily initiated and substantially advanced at the Tomsk Technological University, Tomsk, Russia, where the major algorithms of the shower simulation scheme were developed. Authors thank Prof. We are grateful to an anonymous referee for detailed and helpful comments.

Availability
The ALTAI code is written in FORTRAN and may be easily installed at any computer platforms maintaining FORTRAN compiler. Regarding the availability of the code contact: alexander.konopelko@mpg.mpi-hd.de