text:
# DESY 01-053 April 2001 Multiplicity moments in deep inelastic scattering at HERA

## 1 Introduction

The hadronic final state in deep inelastic scattering (DIS) events is the result of a hard partonic scattering initiated at large momentum transfers, $`Q\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is a characteristic QCD scale of the order of a few hundred MeV. The subsequent parton cascade is followed by a soft fragmentation process. The latter occurs with small momentum transfers which may be considered to extend up to a value $`Q_0`$, which is a QCD cut-off above which perturbative methods can be applied. There are two main approaches for describing the hadronic final state. The first comprises analytic perturbative QCD calculations based on the hypothesis of Local Parton-Hadron Duality (LPHD) ; this hypothesis states that parton-level QCD predictions are applicable for sufficiently inclusive hadronic observables, without additional assumptions about hadronisation processes. Therefore, the hadronic spectra are proportional to those of partons if the cut-off $`Q_0`$ is decreased towards a small value of $`\mathrm{\Lambda }`$. The second approach is formulated in Monte Carlo (MC) programs, which generate a partonic final state according to the perturbative QCD picture and, below $`Q_01`$ GeV, hadronise the partons on the basis of non-perturbative models. The MC models are able to reproduce many detailed properties of the hadronic final state, but contain a large number of free parameters.

Multiplicity distributions and correlations between final-state particles are an important testing ground for analytic perturbative QCD in conjunction with the LPHD hypothesis, as well as for MC models describing the hadronic final state . In this paper, a recent analysis of two-particle angular correlations  is extended to a study of multiplicity distributions measured in restricted phase-space regions of neutral current DIS events. The particle multiplicities are studied in terms of the normalised factorial moments111It may be noted that, for $`q=1`$, $`F_q=1`$ by definition

$$F_q(\mathrm{\Omega })=n(n1)\mathrm{}(nq+1)/n^q,q=2,3,\mathrm{},$$

for a specified phase-space region of size $`\mathrm{\Omega }`$. The number, $`n`$, of particles is measured inside $`\mathrm{\Omega }`$ and angled brackets $`\mathrm{}`$ denote the average over all events. The factorial moments, along with cumulants  and bunching parameters , are convenient tools to characterise the multiplicity distributions when $`\mathrm{\Omega }`$ becomes small. For uncorrelated particle production within $`\mathrm{\Omega }`$, Poisson statistics holds and $`F_q=1`$ for all $`q`$. Correlations between particles lead to a broadening of the multiplicity distribution and to dynamical fluctuations. In this case, the normalised factorial moments increase with decreasing $`\mathrm{\Omega }`$. If the rise follows a power law, this effect is frequently called “intermittency” .

The LPHD hypothesis, originally suggested for single-particle spectra, presently has experimental support in $`e^+e^{}`$ annihilations . However, the success of LPHD is less evident for the moments of inclusive single-particle densities at HERA energies  and two-particle angular correlations . In contrast to the single-particle densities and global event-shape variables, the factorial moments are most directly affected by inclusive many-particle densities. This can be seen from the relation $`n(n1)\mathrm{}(nq+1)=_\mathrm{\Omega }\rho ^{(q)}(p_1,\mathrm{}p_q)𝑑p_1\mathrm{}𝑑p_q`$ between the factorial moments in a region $`\mathrm{\Omega }`$ and the inclusive $`q`$-particle density, $`\rho ^{(q)}(p_1,\mathrm{}p_q)=(1/\sigma _{tot})d^q\sigma /dp_1\mathrm{}dp_q`$. The latter represents the probability to detect $`q`$ particles with momenta $`p_1,\mathrm{}p_q`$, irrespective of the presence of any other particles. Thus the factorial moments are an important tool to study such densities, for which the applicability of LPHD remains questionable even for higher $`e^+e^{}`$ energies at LEP . Deviations of the measured moments from theoretical predictions would imply the necessity for further refinements of perturbative QCD calculations and/or the failure of the LPHD hypothesis. In the latter case, the magnitude of such deviations can shed light on details of the hadronisation processes.

## 2 Analytic QCD results

The dominant source of particle production inside a jet is gluon splitting in the QCD cascade, where the presence of a gluon enhances the probability for emission of another gluon nearby in momentum space. This leads to inter-parton correlations and non-Poissonian statistics for the multiplicity distributions in restricted phase-space intervals where partons are counted. Thus the factorial moments, $`F_q`$, are expected to deviate from unity.

There exist a number of analytic QCD predictions for the moments of the parton multiplicity distributions obtained in the Double Leading Log Approximation (DLLA), which are discussed below.

### 2.1 Multiplicity-cut moments

The normalised factorial moments of the multiplicity distributions of gluons which are restricted in either transverse momentum $`p_t<p_t^{\mathrm{cut}}`$ or absolute momentum $`p𝐩<p^{\mathrm{cut}}`$ are expected, respectively, to have the following behaviour :

$`F_q(p_t^{\mathrm{cut}})`$ $``$ $`1+{\displaystyle \frac{q(q1)}{6}}{\displaystyle \frac{\mathrm{ln}(p_t^{\mathrm{cut}}/Q_0)}{\mathrm{ln}(E/Q_0)}},`$
$`F_q(p^{\mathrm{cut}})`$ $``$ $`C(q)>1,`$
where $`E`$ is the initial energy of the outgoing quark that radiates the gluons and $`C(q)`$ are constants depending on $`q`$. The maximum transverse momentum, $`p_t^{\mathrm{cut}}`$, and the $`p_t`$ values are defined with respect to the direction of this quark (see Fig. 1(a) for a schematic representation in the case of DIS).

Equations (1) indicate that there are positive correlations between partons, because the factorial moments are larger than unity. For small $`p_t^{\mathrm{cut}}`$ values, however, the correlations vanish due to the presence of an angular ordering of the partons in the jet. The existence of the angular ordering leads to the absence of branching processes with secondary gluon emission at small $`p_t^{\mathrm{cut}}`$. This ultimately leads to the suppression of the correlations in this phase-space region and to independent parton emission222 Note that this effect is similar to multiple photon bremsstrahlung in QED. ($`F_q1`$ for $`p_t^{\mathrm{cut}}=Q_0`$). On the other hand, the distribution of soft gluons with limited absolute momenta, $`p<p^{\mathrm{cut}}`$, remains non-Poissonian for any small value of $`p^{\mathrm{cut}}`$.

Note that the theoretical formulae (1) can be regarded as asymptotic calculations, i.e. they are valid at small $`p_t^{\mathrm{cut}}`$ and $`p^{\mathrm{cut}}`$. Therefore, the calculations should be considered only as qualitative predictions when compared to the data using the LPHD hypothesis.

### 2.2 Angular moments

Gluon multiplicities in polar-angle rings around the outgoing-quark direction are expected to exhibit a distribution which is substantially wider than for Poisson statistics . The moments rise with increasing values of the phase-space variable $`z`$ as

$$\mathrm{ln}\frac{F_q(z)}{F_q(0)}=za(1D_q)(q1),$$
(2)
where $`z`$ and the constant $`a`$ are defined as:

$`z`$ $`=`$ $`a^1\mathrm{ln}{\displaystyle \frac{\mathrm{\Theta }_0}{\mathrm{\Theta }}},`$
$`a`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{E\mathrm{\Theta }_0}{\mathrm{\Lambda }}}.`$
The angle $`\mathrm{\Theta }_0`$ is the half opening-angle of a cone around the outgoing quark radiating the gluons, $`E`$ is its energy, and $`\mathrm{\Lambda }`$ is an effective QCD scale ($`\mathrm{\Lambda }\mathrm{\Lambda }_{\overline{MS}}`$). A schematic representation of the variables is shown in Fig. 1(b), in which a window in the angular half-width $`\mathrm{\Theta }`$ is subtended by the ring centered at $`\mathrm{\Theta }_0`$. The factors $`D_q`$ are the so-called Rényi dimensions  known from non-linear physics. A decrease of the angular window $`\mathrm{\Theta }`$ corresponds to an increase of the variable $`z`$. The maximum possible phase-space region ($`\mathrm{\Theta }=\mathrm{\Theta }_0`$), corresponding to $`z=0`$, is used to calculate the factorial moments, $`F_q(0)`$, which are necessary for normalisation in Eq. (2). Such normalisation removes the theoretical uncertainty in the absolute values of $`F_q(z)`$ and thus allows a comparison of the predicted $`z`$-dependence of the moments with experimental results.

For independent particle production in a restricted range of $`2\mathrm{\Theta }`$, $`D_q=1`$ and $`F_q(z)=F_q(0)`$ in Eq. (2). In contrast, the analytic QCD expectations for $`D_q`$ obtained in the DLLA are as follows :

* if the strong coupling constant, $`\alpha _s`$, is fixed at scale $`QE\mathrm{\Theta }_0`$, then for large angular bins (i.e. small $`z`$),

$$D_q=\gamma _0\frac{q+1}{q},$$
(4)
where $`\gamma _0=\sqrt{2C_\mathrm{A}\alpha _s/\pi }`$ is the anomalous QCD dimension, and $`C_\mathrm{A}=3`$ is the gluon colour factor. The $`\alpha _s`$ value is calculated at the one-loop level with $`n_f=3`$ (number of flavours). Such a behaviour of the moments corresponds to the expectations from intermittency, which can be interpreted in terms of the fractal structure of the gluon splittings .

For running-$`\alpha _s`$ values, the $`D_q`$ can be approximated by different expressions, depending on the treatment of the non-leading DLLA contributions to the moments:

* for Dokshitzer and Dremin :

$$D_q\gamma _0\frac{q+1}{q}\left(1+\frac{q^2+1}{4q^2}z\right);$$
(5)
* for Brax et al. :

$$D_q2\gamma _0\frac{q+1}{q}\left(\frac{1\sqrt{1z}}{z}\right);$$
(6)
* for Ochs and Wosiek :

$$D_q2\gamma _0\frac{qw(q,z)}{z(q1)},w(q,z)=q\sqrt{1z}\left(1\frac{\mathrm{ln}(1z)}{2q^2}\right);$$
(7)
* an estimate for $`D_q`$ which includes a correction from the Modified Leading Log Approximation (MLLA) has also been obtained by Dokshitzer and Dremin . In this case, Eq. (5) remains valid, but $`\gamma _0`$ is replaced by an effective $`\gamma _0^{\mathrm{eff}}(q)`$ which depends on $`q`$.

All the theoretical calculations quoted above are performed for partons. Therefore, for comparisons of the predictions with data, LPHD is assumed.

## 3 Experimental set-up

The data were taken in 1996 and 1997 using the ZEUS detector at HERA. During this period, the energy of the positron beam, $`E_e`$, was $`27.5`$ GeV and that of the proton beam was $`820`$ GeV. The integrated luminosity used for the present study was $`38.4\pm 0.6`$ pb-1.

ZEUS is a multi-purpose detector described in detail elsewhere . Of particular importance for the present study are the central tracking detector (CTD), positioned in a 1.43 T solenoidal magnetic field, and the calorimeter.

The CTD is a cylindrical drift chamber with nine superlayers covering the polar angle333 The ZEUS coordinate system is a right-handed Cartesian system, with the $`Z`$ axis pointing in the proton beam direction, and the $`X`$ axis pointing left towards the centre of HERA. The coordinate origin is at the nominal interaction point. The polar angle $`\theta `$ is defined with respect to the positive $`Z`$-direction. The pseudorapidity, $`\eta `$, is given by $`\mathrm{ln}(\mathrm{tan}\frac{\theta }{2})`$. region $`15^o<\theta 164^o`$ and the radial range 18.2-79.4 cm. Each superlayer consists of eight anodes. The transverse-momentum resolution for charged tracks traversing all CTD layers is $`\sigma (p_{})/p_{}=0.0058p_{}0.00650.0014/p_{}`$, with $`p_{}`$ being the track-transverse momentum (in GeV). The single hit efficiency of the CTD is greater than $`95\%`$.

The CTD is surrounded by the uranium-scintillator calorimeter (CAL), which is longitudinally segmented into electromagnetic and hadronic sections. The relative energy resolution of the calorimeter under test beam conditions is $`0.18/\sqrt{E}`$ for electrons and $`0.35/\sqrt{E}`$ for hadrons (with $`E`$ in GeV). The CAL cells provide time measurements with a time resolution below 1 ns for energy deposits larger than 4.5 GeV.

## 4 Data selection

### 4.1 Event kinematics

The basic event kinematics of DIS processes are determined by the negative squared four-momentum transfer at the positron vertex, $`Q^2=q^2=(kk^{^{}})^2`$ ($`k`$ and $`k^{^{}}`$ denote the four-momenta of the initial- and final-state positrons, respectively) and the Bjorken scaling variable $`x=Q^2/(2Pq)`$, where $`P`$ is the four-momentum of the proton. The fraction of the energy transfered to the proton in its rest frame, $`y`$, is related to these two variables by $`yQ^2/sx`$, where $`\sqrt{s}`$ is the positron-proton centre-of-mass energy.

The moments were studied in the Breit frame , which provides a natural system to separate the radiation from the outgoing struck quark from the proton remnants. In this frame, the exchanged virtual boson with virtuality $`Q^2`$ is completely space-like and has a momentum $`q=(q_0,q_{X_\mathrm{B}},q_{Y_\mathrm{B}},q_{Z_\mathrm{B}})=(0,0,0,Q)`$. In the quark-parton model, the incident quark has $`p_{Z_\mathrm{B}}=Q/2`$ and the outgoing struck quark carries $`p_{Z_\mathrm{B}}=Q/2`$. All particles with negative $`p_{Z_\mathrm{B}}`$ form the current region. These particles are produced by the fragmentation of the struck quark, so that the current region is analogous to a single hemisphere in $`e^+e^{}`$ annihilations. The variables $`p_t=\sqrt{p_{X_\mathrm{B}}^2+p_{Y_\mathrm{B}}^2}`$ and $`\mathrm{\Theta }_0`$ were determined with respect to the negative $`Z_\mathrm{B}`$-direction, which, for high $`Q^2`$ events, is a good estimate of the outgoing struck-quark direction. The moments in the current region of the Breit frame were measured in the kinematic region $`Q^2>1000`$ GeV2, which is sufficiently high to avoid major contributions from the boson-gluon fusion (BGF) processes (see Section 7).

### 4.2 Event and track selection

The kinematic variables $`x`$ and $`Q^2`$ were determined by a combination of methods: (i) the measurement of the energy and angle of the scattered positron (denoted by the subscript $`e`$); (ii) the double angle (DA) method  and (iii) the Jacquet-Blondel (JB) method , depending on which is most appropriate for a particular cut. In the DA method, the variables $`x`$, $`Q^2`$ and $`y`$ were reconstructed using the angles of the scattered positron and the hadronic energy flow. In the JB method, used only in the event selection, the kinematic variables were determined entirely from the hadronic system.

The boost vector from the laboratory to the Breit frame was determined using the DA method, which is less sensitive to systematic uncertainties in the energy measurement than the other methods. In the reconstruction of the Breit frame, all charged particles were assumed to have the pion mass.

The triggering and online event selections were identical to those used in a previous publication . To preselect neutral current DIS events, the following cuts were applied:

1. $`E_e^{^{}}10`$ GeV, $`E_e^{^{}}`$ being the energy of the scattered positron;
2. $`Q_{\mathrm{DA}}^21000`$ GeV2;
3. $`35\delta =E_i(1\mathrm{cos}\theta _i)\mathrm{\hspace{0.25em}60}`$ GeV, where $`E_i`$ is the energy of the $`i^{\mathrm{th}}`$ calorimeter cell and $`\theta _i`$ is its polar angle with respect to the beam axis;
4. $`y_e\mathrm{\hspace{0.25em}0.95}`$;
5. $`y_{\mathrm{JB}}\mathrm{\hspace{0.25em}0.04}`$;
6. $`Z`$ coordinate of the reconstructed vertex position, as determined from the tracks fitted to the vertex, in the range $`40<Z_{\mathrm{vertex}}<50`$ cm;
7. a timing cut requiring that the event time as measured by the CAL be consistent with an $`e^+p`$ interaction.

The resulting data set contained 7369 events.

The present study is based on the tracks reconstructed with the CTD and fitted to the event vertex. The scattered positron was removed from the track sample. In addition, the following cuts to select tracks were imposed:

1. $`p_t^{lab}>0.15`$ GeV, where $`p_t^{lab}`$ is the transverse momentum of the track;
2. $`\eta <1.75`$.

These cuts restrict the data to a well understood, high-acceptance region of the CTD.

## 5 Event simulation

Neutral current DIS events were generated with ARIADNE 4.10  with the high $`Q^2`$ modification developed during the HERA Monte Carlo workshop  using tuned parameters . This MC program is based on the colour-dipole model, in which gluons are emitted from the colour field between quark-antiquark pairs, supplemented with the BGF process. The hadronisation was simulated using the Lund string model as implemented in JETSET . Charged hadrons with a mean lifetime longer than $`10^8`$ seconds were considered as stable. The CTEQ4D  parameterisation of the proton parton distribution functions was used. The events were generated without any cuts on kinematic variables and track quantities. The MC events obtained by this procedure defined the generator-level sample.

To obtain a detector-level sample, events were generated with the DJANGO 6.24  program based on HERACLES 4.5.2 , in order to incorporate first-order electroweak corrections. The parton cascade and hadronisation were then simulated with ARIADNE. The events were then processed through a simulation of the detector using GEANT 3.13  to take into account particle interactions with the detector material, particle decays, event and track migrations, resolution, acceptance of the detector and event selection. The detector-level MC events were processed with the same reconstruction program as was used for the data.

In addition to ARIADNE, the LEPTO 6.5  and HERWIG 5.9  generators with tuned parameters  were also used to compare with the data. The parton cascade in LEPTO is based on a matrix-element calculation matched to parton showers according to the DGLAP equations. As in ARIADNE, the hadronisation in LEPTO is simulated with JETSET. HERWIG also has a parton shower based on DGLAP evolution, but parton emissions are ordered in angle, in contrast to LEPTO, which orders the parton emissions in invariant mass with an additional angular constraint to ensure coherence. The hadronisation in HERWIG is described by the cluster fragmentation model .

The Bose-Einstein (BE) interference between identical particles can produce additional correlations which distort the behaviour of the moments measured in small phase-space regions. By default, the BE interference is turned off in JETSET; it is absent in HERWIG. The effect of the BE correlations on the present results is considered in the next section.

## 6 Correction procedure

The bin-by-bin correction factors, $`𝒞_q=𝒜_q^{\mathrm{gen}}/𝒜_q^{\mathrm{det}}`$, for a given kinematic region in $`Q^2`$ were evaluated using the ARIADNE 4.10 MC event samples, separately for each observable, $`𝒜_q=F_q(p_t^{\mathrm{cut}}),F_q(p^{\mathrm{cut}})`$ and $`F_q(z)/F_q(0)`$, where $`𝒜_q^{\mathrm{gen}}`$ ($`𝒜_q^{\mathrm{det}}`$) was calculated at the generator (detector) level. The factors $`𝒞_q`$ correct the data for the detector effects and event-selection cuts as well as for the effects that are not included in the analytic calculations, namely (i) Dalitz decays of the $`\pi ^0`$; (ii) leptons in the final state and (iii) initial-state photon radiation.

The correction factors for the $`p_t^{\mathrm{cut}}`$ and $`p^{\mathrm{cut}}`$ moments are approximately flat and vary between about 1.0 and 1.3, except for the smallest momentum cuts and the highest-order $`p_t^{\mathrm{cut}}`$ and $`p^{\mathrm{cut}}`$ moments, where the correction factors are as large as $`1.7`$. The correction factors for the $`F_q(z)/F_q(0)`$ moments vary smoothly from 1.0 to 1.5. The contribution of Dalitz pairs to the correction factors is negligible, except for the smallest phase-space regions in $`p_t^{\mathrm{cut}}`$, $`p^{\mathrm{cut}}`$ and $`\mathrm{\Theta }`$, where it does not exceed $`20\%`$.

The analytic QCD calculations do not include BE interference. To take into account this phenomenon, the correction factor, BEq, was calculated for each data point

$$\mathrm{BE}_q=\frac{𝒜_q^{\mathrm{gen}}}{𝒜_q^{\mathrm{BE},\mathrm{gen}}},$$
(8)
where $`𝒜_q^{\mathrm{gen}}`$ ($`𝒜_q^{\mathrm{BE},\mathrm{gen}}`$) are the moments at the generator level without (with) the BE effect. For this correction, the BE effect according to a Gaussian parameterisation with JETSET default parameters was used. As a cross check, the BE correction was determined from the ZEUS data and also using the parameters obtained by the H1 Collaboration . The differences between the correction factors BEq for these three sets of parameters were small. The values of the BE correction factors are discussed in Section 8.

## 7 Systematic uncertainties

The overall systematic uncertainty for each moment was determined from the following uncertainties added in quadrature:

1. the choice of the MC models. HERWIG was used in place of ARIADNE to determine the correction factors. This uncertainty was typically $`45\%`$ of the total systematic uncertainty;
2. event reconstruction and selection. Systematic checks were performed by varying the cuts on $`\delta `$, $`y_e`$, $`y_{\mathrm{JB}}`$ and the vertex position requirement ($`40\delta \mathrm{\hspace{0.25em}55}`$ GeV, $`y_e\mathrm{\hspace{0.25em}0.85}`$, $`y_{\mathrm{JB}}\mathrm{\hspace{0.25em}0.05}`$ and $`35<Z_{\mathrm{vertex}}<45`$ cm). The overall uncertainties associated with the event selection were typically $`25\%`$ of the total systematic uncertainty;
3. track reconstruction and selection. Tracks were required to have transverse momenta larger than $`0.2`$ GeV and $`\eta <1.5`$. Tracks that reach at least the third CTD superlayer were used. These uncertainties were typically $`30\%`$ of the total systematic uncertainty.

As an example, the total systematic uncertainty of $`F_2(p_t^{\mathrm{cut}})`$ for the smallest $`p_t^{\mathrm{cut}}`$ value used in this analysis is approximately $`6\%`$. For the large $`p_t^{\mathrm{cut}}`$ regions, a typical value of the systematic uncertainty is $`1\%`$.

To ensure that the two particles in a pair were resolved, the angle between two tracks in the current region of the Breit frame should not be smaller than two degrees . This limit determines the minimum size of phase-space regions used in this analysis.

The kinematic variables $`p_t`$, $`\mathrm{\Theta }`$ and $`\mathrm{\Theta }_0`$ have to be measured with respect to the outgoing-quark direction, which coincides with the $`Z_\mathrm{B}`$-direction of the Breit frame only for the zeroth-order $`\gamma ^{}q`$ scattering. The first-order QCD effects, BGF ($`\gamma ^{}gq\overline{q}`$) and QCD Compton ($`\gamma ^{}qqg`$), can lead to two hard partons that are not collinear to the $`Z_\mathrm{B}`$-direction. For these events, the determination of the kinematic variables with respect to the $`Z_\mathrm{B}`$ axis of the Breit frame is inappropriate, especially when these two hard partons escape the current region. Using the LEPTO MC model, it has been determined that, at $`Q^2>1000`$ GeV2, the number of BGF and QCD Compton events without hard partons in the current region, relative to the total number of events, is $``$ 15%. As a systematic check, the factorial moments have been calculated in kinematic variables with respect to the thrust axis determined from the charged particles in the current region of the Breit frame. The results agree with those for which the negative $`Z_\mathrm{B}`$-direction of the Breit frame is used.

## 8 Results

The charged-particle multiplicities measured in the entire current region without additional momentum cuts on the charged particles were compared to the previous ZEUS measurements  and to $`e^+e^{}`$ data . For the present measurement, the corrected average value of $`Q^2`$ is $`Q^22070`$ GeV2, which corresponds to an energy $`\sqrt{Q^2}46`$ GeV. The corrected average charged-particle multiplicity in the current region, $`5.76\pm 0.29`$, where the error represents statistical and systematic uncertainties added in quadrature, agrees with the previous ZEUS measurement, $`6.01\pm 0.46`$, for a similar value of $`Q^2`$. This average multiplicity is lower than that found in a single hemisphere of $`e^+e^{}`$ at $`\sqrt{s}=44`$ GeV, which is $`7.58\pm 0.37`$ . The second-order factorial moment $`F_2`$ in a single hemisphere of $`e^+e^{}`$ annihilations is444 The second-order moment in a single hemisphere of $`e^+e^{}`$ annihilation can be estimated from the dispersion $`D=\sqrt{n^2n^2}`$ and the average multiplicity $`n`$  using the relation $`F_2=1+D^2n^2n^1`$. $`F_2=1.036\pm 0.026`$, which is consistent with the value $`1.051\pm 0.018`$ obtained in this study of the entire current region. This difference in the average multiplicity can be attributed to a migration of particles between the current and target regions due to BGF and QCD Compton processes, an effect which has been theoretically quantified  and directly measured by ZEUS . Such a migration leads to current-target anti-correlations which do not vanish even at the high $`Q^2`$ values studied here. This is unlike the $`e^+e^{}`$ annihilation processes, where weak and positive forward-backward correlations have been observed at a similar energy .

### 8.1 Multiplicity-cut moments

The multiplicity moments of order $`q=2,\mathrm{},5`$ as a function of $`p_t^{\mathrm{cut}}`$ are shown in Fig. 2. The $`p_t^{\mathrm{cut}}`$ moments were calculated with respect to the negative $`Z_\mathrm{B}`$-direction of the Breit frame. As $`p_t^{\mathrm{cut}}`$ decreases below 1 GeV, the moments rise. This disagrees with the DLLA theoretical calculations, which predict that the moments approach unity (see Eq. (1)). On the other hand, Monte Carlo models show the same trend as the data, although some differences in slopes are apparent. LEPTO exhibits smaller values of the moments than for the data at small $`p_t^{\mathrm{cut}}`$ values, while HERWIG overestimates the moments. ARIADNE gives the best overall description of the correlations.

Figure 2 also shows the factors BEq from Eq. (8) used to correct the data for the BE interference. The largest correction was found in the case of moments of high orders ($`q=4,5`$) for small $`p_t^{\mathrm{cut}}`$ cuts ($`<0.5`$ GeV).

To investigate contributions of the hadronisation processes to the moments, further samples of the ARIADNE MC were generated: (1) considering neutral hadrons in addition to charged particles; (2) using only hadrons coming from the Lund string fragmentation without resonance decays. For these two cases, the moments have an even steeper rise at low $`p_t^{\mathrm{cut}}`$. Replacing the Lund string-fragmentation scheme by independent fragmentation and removing resonance decays led to an increase of the moments for $`p_t^{\mathrm{cut}}0.5`$ GeV, but below this cut the moments approach constant values which are larger than unity. An MC sample generated without the colour-dipole parton emission and resonance decays, but retaining the Lund string fragmentation, leads to an increase of the moments with decreasing $`p_t^{\mathrm{cut}}`$.

Figure 3 shows the factorial moments for restricted $`p`$ values. As was the case for the $`p_t^{\mathrm{cut}}`$ moments, there is no indication that the values approach unity at small $`p^{\mathrm{cut}}`$ values. All MC models at the hadron level agree with the data. The BE interference has a large influence on this measurement, giving a $`20\%`$ correction.

The DLLA QCD calculations neglect energy-momentum conservation in gluon splittings, considering idealised jets at asymptotically large energy. This effect can be taken into account in the Modified Leading Log Approximation; however, such calculations have not yet been performed. Therefore, it is important to compare the experimental results with the parton shower of MC models that explicitly include energy-momentum conservation and the complete parton-splitting functions. In addition, such comparisons allow the investigation of the contributions from hadronisation.

The moments restricted in $`p_t`$ and $`p`$ variables were determined from the ARIADNE parton-level, the physics implementation of which strongly resembles the analytic calculations . The LPHD requires $`Q_0`$ to approach $`\mathrm{\Lambda }`$. The ARIADNE default value of $`\mathrm{\Lambda }`$ is 0.22 GeV, while the transverse momentum cut-off $`Q_0`$ is $`0.6`$ GeV. Therefore, to satisfy the LPHD requirement, the parton predictions were obtained after switching off the string fragmentation and lowering the $`Q_0`$ value in ARIADNE to $`Q_0=0.27`$ GeV. This value was chosen in order to insure that the average multiplicity of partons in the current region of the Breit frame equalled that of hadrons, which is important for normalisation purposes. This method is somewhat different to that used in a similar study of the ARIADNE parton cascade for $`e^+e^{}`$ annihilation events . In contrast to that approach, the present modification of the model retains the original default value of $`\mathrm{\Lambda }`$, to preserve the kinematics of first-order QCD processes in the Breit frame.

The parton-level ARIADNE for $`p_t^{\mathrm{cut}}`$ and $`p^{\mathrm{cut}}`$ moments is shown as thin solid lines in Figs. 2 and 3. Agreement between the parton-level predictions and the data for the lowest $`p_t^{\mathrm{cut}}`$ values can be excluded at greater than $`99\%`$ confidence level for $`q=2,3`$ and at greater than $`95\%`$ confidence level for $`q=4,5`$. The data for all $`p^{\mathrm{cut}}`$ moments disagrees with the parton-level predictions by many standard deviations. The MC predictions for partons have the trends expected from the analytic DLLA results, showing a decrease of the moments for decreasing $`p_t^{\mathrm{cut}}`$ and $`p^{\mathrm{cut}}`$. This is contrary both to the measurements and to the hadron-level MC predictions, which show a rise of the moments. Analogous differences between the hadron and parton levels of ARIADNE have also been observed in $`e^+e^{}`$ annihilation .

The parton predictions in the full current region (large $`p_t^{\mathrm{cut}}`$ and $`p^{\mathrm{cut}}`$) show larger values of the moments than the data, despite the tuning of the parton to the hadron multiplicities. In the case of $`e^+e^{}`$, the MC-tuned parton level describes the hadron-level predictions for large momentum cuts . It is likely that the discrepancy between partons and hadrons in DIS is due to additional non-perturbative effects related to the proton remnant.

### 8.2 Angular moments

Figure 4 shows the angular moments $`F_q(z)/F_q(0)`$. This represents the first measurement of these quantities in DIS. The measurements are compared to DLLA QCD calculations in Eq. (2), after substitution of different forms of the Rényi dimension (see Eqs. (4)-(7) ) and using the effective QCD parameter $`\mathrm{\Lambda }=0.15`$ GeV as in LEP studies . For DIS in the Breit frame, the energy $`E=Q/2`$ was chosen to define the $`\alpha _s`$ value, corresponding to the outgoing quark momentum in the quark-parton model. The factors used to correct the data for the BE interference are close to unity for all angular moments (not shown). A significant disagreement with the data is observed, which becomes smaller for higher orders of the moments. Although the error bars are clearly correlated, even for the high-order moments, agreement between the predictions and the data can be excluded for $`z=`$0.05-0.2 at greater than $`95\%`$ confidence level. For running-$`\alpha _s`$ calculations, the best description is observed for high-order moments in the case of the DLLA result with a correction from the MLLA. This result is similar to that obtained for $`e^+e^{}`$ annihilation . It should be noted that the present analysis, however, uses a larger angle $`\mathrm{\Theta }_0=45^o`$, in order to increase the statistics, and that the average energy ($`E23`$ GeV) of the outgoing quark giving rise to the parton shower is smaller than at LEP by a factor of two.

Figure 5 shows the comparisons of the angular moments with Monte Carlo models at the hadron level. All models reproduce the trends of the data but somewhat overestimate the strength of the correlations; however, LEPTO fails to give a reasonable description of the angular moments.

The parton-level ARIADNE predictions for the angular moments are shown in Fig. 5. The ratio $`F_q(z)/F_q(0)`$ predicted by ARIADNE at the parton level (thin lines) overestimates the data, but is closer than the analytic DLLA calculations for the second-order moment. This comparison with the MC parton level and similar studies of $`e^+e^{}`$ annihilations at LEP  imply that MLLA analytic calculations are necessary to give a better description of the measured angular moments.

## 9 Conclusions

Multiplicity moments have been studied in deep inelastic scattering at $`Q^2>1000`$ GeV2 with the ZEUS detector at HERA. The moments measured in the current region of the Breit frame show a strong rise as the size of the restricted intervals in $`p_t`$, $`p`$ and $`\mathrm{\Theta }`$ decreases. Monte Carlo models, which include hadronisation effects, reproduce the correlation pattern of the hadronic final state, although there are some discrepancies with the data.

The analytic QCD calculations for the angular moments agree qualitatively with the data, but there are discrepancies at the quantitative level. An inadequate treatment of energy-momentum conservation in gluon splittings, which is important at low energies, is likely to be responsible for the discrepancies between the analytic QCD expectations and the data. The parton level of ARIADNE, which does not suffer from the limitations of the DLLA calculations, indicates that the discrepancies can be reduced. MLLA calculations are required to provide an accurate description of the angular moments. Similar studies of the angular moments measured at LEP also show that the analytic perturbative calculations need to be improved.

In contrast to the angular fluctuations, the moments calculated in restricted $`p_t`$ and $`p`$ regions show significant discrepancies, both for the DLLA calculations and for the parton-level of ARIADNE. This is most clearly seen for the $`p_t^{\mathrm{cut}}`$ moments; while the analytic calculations and ARIADNE predictions for partons decrease and approach unity for small $`p_t`$ regions, the data show a steep rise. This is the first indication that perturbative QCD, in conjunction with the LPHD hypothesis, fails on a qualitative level to describe the hadronic multiplicities. Monte Carlo studies of the hadronisation effects indicate that this is due to violation of the LPHD hypothesis for many-particle densities measured in the regions restricted in the transverse momenta with respect to the direction of the outgoing quark.

## Acknowledgements

The strong support and encouragement of the DESY Directorate have been invaluable, and we are much indebted to the HERA machine group for their inventiveness and diligent efforts. The design, construction and installation of the ZEUS detector have been made possible by the ingenuity and dedicated efforts of many people from inside DESY and from the home institutes who are not listed as authors. Their contributions are acknowledged with great appreciation. We thank W. Ochs and J. Wosiek for helpful discussions.


format: text

text:
# 1 Introduction

## 1 Introduction

The Earth’s atmosphere is being bombarded continuously by a flux of particles (cosmic rays), coming from all directions. Their energies range from a few MeV to more than 1020 eV. Their spectrum follows a power law with a negative exponent which is almost constant over thirteen orders of magnitude in energy. The origin of the cosmic rays is still an open question. Those rays with energy below $``$ 1 GeV are likely to have a solar origin, but for higher energies their acceleration mechanism remains in mystery. It is believed that up to $``$ 4.1015 eV they can be accelerated by diffuse shock processes produced in supernova explosions. In this energy region, (usually called “the knee”), the exponent of the power law describing the cosmic ray flux per units of area, time, solid angle, and energy, suddenly steepens from $``$ -2.7 to $``$ -3.2, and this change is believed to be related with the maximum energy that can be transferred by a supernova shock to a single particle. If the kinetic energy of the cosmic ray is high enough, then secondary particles are produced as a consequence of hadronic or electromagnetic interactions with the upper atmosphere atomic nuclei. Those secondary particles will, in turn, produce more particles, yielding a cascade which is known as an Extensive Air Shower (EAS). Depending on the primary energy and zenithal angle, this cascade can be stopped in the atmosphere, or even reach the ground level.

A method which has been used for the observation of EAS is the detection of the light emitted by the Čerenkov effect in air, while fast charged particles, (mainly electrons), are crossing the atmosphere (WHIPPLE, CANGAROO). Alternatively, it is possible to observe the UV light emitted by decay processes occurring in the atmospheric molecules after excitation by the EAS’s secondary particles (Fly’s Eye, HiRes, Pierre Auger Project, Telescope Array, OWL Project). The amount of UV and Čerenkov light emitted by an EAS is extremely faint, and because of this it is possible to observe these processes only during moonless dark nights, and by using relatively large telescope mirrors as light concentrators and sensitive photomultiplier tubes.

Another (and perhaps more common) approach (Haverah Park, AGASA, Volcano Ranch, SUGAR) is the direct detection of the shower secondary particles reaching the ground level. The size of the footprint at ground level is several thousand square meters for showers produced by primary cosmic rays of energies near the “knee” or higher. Because of this these experiments are designed so as to observe only samples of the particle showers using an array of ground-based detector stations, where gas-filled chambers, plastic scintillators or Čerenkov-effect detectors are typical components.

The detector stations of these ground arrays are usually capable to measure particle densities. In the case of the array described in the present work, where water Čerenkov detectors (WCD) are used, this measurement is performed through a sample of the amount of light emitted when the shower particles traverse through the water radiator. Also, the precise relative times of the signals produced by each station are recorded, together with the Čerenkov light intensity information. By using the relative hitting times at each station and the known geometry of the array it is possible to determine the direction of arrival of the primary cosmic ray, assuming that the general development of the shower follows a rather flat front profile. Rigorously the shower front is a curved surface whose radius of curvature could in principle be determined if the number (and quality) of the sampling detectors is high enough.

The determination of the primary energy from EAS measurements using ground-based detectors is closely tied to shower reconstructions, based on Monte Carlo simulations. These simulations correlate the primary energy to the particle densities at a fixed distance from the shower “core” position, that is, the center of gravity of the air shower at ground. In a simplified model, the primary energy is simply estimated as a magnitude proportional to the total number of particles in a shower. Hence, the particle density measurements performed by each station is used to estimate the total number of particles of the shower.

In the following sections the design and construction of this new air shower experiment, which has been optimized for the “knee” region of the energy spectrum, are described. The necessary simulations, which were required to set the numerous design parameters of the array, are presented in detail.

## 2 The array

The TANGO (TANdar Ground Observatory) Array has been constructed in Buenos Aires, Argentina, at ($``$ 15 m a.s.l), 35 34’ 21” S and 58 30’ 50” W, in the Campus of the Constituyentes Atomic Center, belonging to the Argentinean Atomic Energy Commission (CNEA). The data acquisition (DAQ) room was set inside the TANDAR Accelerator Building. Three detectors are placed on the vertices of an almost isosceles triangle, and a fourth detector was installed on top of the building, in a convenient position close to the center of the triangle, as shown in Figure 1. The final positions of the detector stations were conditioned by the free space available between the existing buildings, and an effort was made to come up to an overall shape as close as possible to an equilateral triangle, which maximizes the effective collection area. The distances between surface stations were measured using a GPS and their error has been estimated in $`\pm `$ 1 m (the measurement was performed after release of the high precision GPS service). The final configuration encloses a geometrical area of 31286 m2. The array has a yearly average overburden of $``$ 1000 gr/cm2.

During an EAS event the DAQ system measures both, the intensities of the Čerenkov photons emitted by the water when crossed by the secondary particles of a high energy cosmic ray’s EAS, and also, the arrival time of the these particles to each station. The threshold energy of the array, resulting from the geometry and from the particular detector conditions, (present noise, trigger levels, etc) is close to 1014 eV for vertical showers. The detector stations are connected by low attenuation (RG-213) coaxial cables to the DAQ room, where the signals are recorded using a 4-channel digital oscilloscope connected to a computer. Depending on the selected trigger conditions, which are generated by standard NIM electronic modules, the number of accepted events ranges from $``$ 250 to $``$ 1500 per day.

### 2.1 The detector stations

This array is a project which grew up from the first 1:1 scale prototype of a WCD (See Figure 2) built in 1995 by members of the local Pierre Auger Project Collaboration . This first detector (labelled A in Figure 1) was construced in a tank, cylindrical in shape, made of 0.68 mm stainless steel plate, with a footprint area of 10 m2. The effective water depth is 120 cm. Three, 8-inch photomultipliers (Hamamatsu R1408), symmetrically placed at 120 cm from the tank axis were installed looking down on the top of the detector, having only the photocathodes immersed in the water working as Čerenkov radiator. Thus, a sample of the Čerenkov photons emitted when a charged particle crosses the tank are collected by the three PMTs. This detector, being a prototype, was designed as a flexible system, allowing the introduction of modifications in the photomultiplier positions, the effective water height, or the inner lining material. During the measurements as a component of the TANGO array the configuration of this detector was that of the Pierre Auger Project baseline design, with the dimensions mentioned above. In order to improve the optical properties of the inner surfaces all detectors were fully lined with Tyvek which is a highly UV-diffusive and reflective material.

The two other detectors sitting in the vertices of the triangle (B and C in Figure 1) have the same general dimensions quoted previously. They are made of 1 mm thick stainless steel, and the external walls are shaped as a dodecagon (see Figure 3). In these detectors we used Hamamatsu R5912, 8-inch diameter PMTs, arranged with the same geometry used in the first detector tank. The fourth detector (D) is smaller, it was made using a fiberglass-reinforced polyestyrene tank, with a footprint of 0.5 m2 and an effective water depth of 80 cm. The tank was also internally lined with Tyvek and only one, 3-inch PMT, was installed centered on the top of the tank.

The larger outer detectors are more adequate to measure lower particle densities generated by showers falling relatively far away from them, either close to the center of the array or very away the whole array. The relatively large ratio between the volume of the external detectors to that of the smaller central detector helps to improve the accuracy in the determination of the particle densities in those cases where the shower core falls close to the center of the array. This is so because of the larger dynamic range of the central detector, which admits higher particle densities without going into saturation (See Section 2.3, and the higher sensitivity of the larger detectors placed on the vertices of the triangle.

All PMTs used in the WCDs were mounted in water-tight enclosures that protect from moisture their voltage dividers and only the photocathode areas of the glass bulbs are immersed in the water radiator (See Figure 2. The glass bulbs of the PMTs were glued to the PVC housings using an elastic silicone compound to reduce mechanical stresses that could break the glass, as happened in the Milagrito experiment. Local high voltage power suplies, fed from the AC mains, were installed near each station. The bias configuration of all PMTs was adopted as grounded cathode, to prevent eventual noise produced by electrical leaks or discharges through the glass.

The water used to fill the tanks was treated in a reverse-osmosis plant, producing an average final water resistivity of about 1 M$`\mathrm{\Omega }`$-cm. Before filling the tanks they were carefully degreased, brushed with water and mild detergent and rinsed abundantly with the same water used as the detector material. These precautions, together with the darkness and the fact that the water used as detector material has a very low level of bacteria nutrients, virtually blocked any extensive biological activity . After more than 1 year since the filling of the detectors, no significative decrease in the signal strength has been observed.

### 2.2 Characterization of the photomultiplier tubes

The gains of the three (R1408) PMTs used in the first prototype were measured previously, and we have built a dark box, adequate for measuring the gain and dark current of the new Hamamatsu R5912 tubes, purchased for the new detector stations. We also used this darkbox to characterize the photocatode sensitivity profiles and the influence of the Earth magnetic field direction of the PMTs. The dimensions of the box are 50 x 50 x 100 cm and accepts one PMT, which is mounted in an axially rotatable holder.

Using electrons thermally emitted from the photocathodes at room temperature, we measured the gains of the PMTs by means of the single electron technique. In these measurements we tested different voltage divider configurations in a range of high voltage from 1100 to 1800 V. In all cases the tubes were kept in total darkness for at least 2 hours before collecting the single electron spectra, to reduce the rate of multi-electron emissions due to fast decaying fluorescence in the photocatodes. The dark current pulse rate (threshold = 1/3 p.e.) after one hour of storage in total darkness was about 900 Hz at 1500 Volts. A typical spectrum and a plot of the gain values are shown in Figure 4.

A study of the influence on the PMT gain of the Earth magnetic field direction relative to the dynode geometry, was also performed. The gains were carefully measured at a fixed HV setting (1.5 kV) for 8 different azimuthal orientations of the PMTs, covering 360 . The tubes were measured keeping always their axis in vertical position. No significative shifts (less than $`\pm `$ 1%) were observed in the peak positions of the single electron spectra obtained in this way. From these results we concluded that the local influence of the Earth magnetic field direction on the gain can be neglected, and thus no special care was taken to install the tubes in the detector stations at any special dynode orientation.

### 2.3 Gain settings and calibration of the detector electronics

After testing several possible configurations we adopted a pasive, grounded cathode design for the voltage dividers, as simple and reliable as possible, assuring both wide dynamic range, and linearity. The final configuration chosen is similar to that recommended by Hamamatsu for the R5912 tubes. Metal film resistors were used, and decoupling capacitors stabilize the 3 last dynodes. The finished printed cards were protected with water resistant varnish and installed with bags of dissicant material in the water-tight enclosures mentioned above.

It is important that the gains of the 3 PMTs installed in each detector station are matched to each other to avoid unbalances in charge collection. Unmatching could impair the homogeneity in the response of the detector, and even reduce the dynamic range. Because of availability of equipment and space, we used a common high voltage supply in each detector station. Because of this, and provided that the observed differences in gains were small (less than 15 %), we compensated the differences in gain by using passive, constant-impedance variable attenuators to reduce as necessary the output pulse amplitude of the two tubes having larger gains in each station.

In order to determine the relative gains we adopted a routine procedure based on the measurement of the signal from each PMT produced by background muons. The trigger for this measurement is taken from the signals from other PMT belonging to the same station, as described in. Once the average relative gains are obtained for the three tubes, the attenuators are set in the two PMTs with higher gains, matching the peak position produced by the PMT with lowest gain. This procedure for gain matching and calibration has been performed on a monthly basis during the complete period of measurements. The system proved to be very stable and very few adjustments were required along this time.

In addition to this periodic gain matching monitoring procedure, a daily routine for monitoring the overall gain of the 4 detector stations has been performed. It also uses the natural background muons falling in the detectors as the source of signals for calibration. It has been found that the spectra of the summed signals of the three PMTs within each peripheral detector station, and also the response of the 3” PMT installed in the small central detector show clearly a peak when they are triggered by themselves. Although it is somewhat broad, the position of this “background” muon peak is very closely the same as the position of the similar peak obtained when a pair of external plastic scintillators are used to select vertical and central muons for triggering. This is valid for both, voltage and charge spectra. This experimental result, which might be due to the remarkable uniformity in the light distribution produced by the Tyvek liners, provides a simple and reliable procedure for remote monitoring and calibration of the station gain . This peak value has been called VEM (Vertical Equivalent Muon), and is defined as the charge (or voltage) peak produced by singly charged, energetic particles, crossing vertically the detector along its axis. This VEM-value is a characteristic parameter of each detector, and depends on its components, geometry, construction, and also on its operation conditions (transparency of the water radiator, bias voltage, etc). The VEM-value provides a practical way to normalize the signals from different detectors and moreover, to express the total signal produced in each station by an EAS (i.e. muons, electrons, gamma rays, etc., hitting the station) in terms of an “equivalent reference particle”. Muons have been selected in this case as they are present everywhere and proved to be very convenient for calibration.

In previous studies performed with the prototype detector, very good homogeneity in charge collection was obtained by using the sum of the signals of the three PMTS. This behavior, which might be again attributed to the excellent light spread produced by the Tyvek liners, is kept almost independently of the entrance points and directions of the muons.

According to these results, our design included fast active adder circuits installed in the peripheral detectors. The operational amplifier employed (CLC 452) works also as the driver for the relatively long RG-213 cable, carrying the signals from each detector to the DAQ room. Although their response in speed is excellent (we require a 130 MHz bandwidth), these circuits introduce a limitation in the dynamic range as the maximum span voltage is less than 2 V. In order to reduce the reflection of the pick up noise signal the impedance of the cable was matched at the sending end, but this reduces further the available amplitude to only 1.4 V. On the other hand, an acceptable signal to noise ratio in the DAQ room asks for minimum signal amplitudes for single muons of $``$ 100 mV. These figures, together with the measured RF pick up and the signal attenuation produced along the cables, limit the final available dynamic range from 1 to 15 muons. Even when this dynamic range is limited, it has been found to be acceptable because only $``$ 30% of the events had to be rejected in the off-line data analysis due to electronic saturation in any station.

Because there is only one PMT in the central detector, and the cable length to the DAQ room is relatively short ($``$ 45 m), its anode signal was directly sent without summing circuit nor attenuator, and hence without the limitation in the dynamic range present in the outer detectors.

### 2.4 Trigger System and Data Acquisition

The signals from the four stations arriving to the electronics front panel are split using linear fan-in/fan-out (FIFO) modules (see Figure 5). Then, they are fed directly to the input connectors of a four-channel digital oscilloscope (Tektronix TDS 3034 set at 500 Ms/s). The oscilloscope is the core of our DAQ system and works as the digitizing stage for detector signals under control of a STROBE pulse. The detector signals are at this point unsynchronized after travelling different lengths of cables (206, 196, 310 and 44 m, for detectors A,B,C,D, respectively). Thus, in order to generate valid trigger conditions, it is essential to compensate these different transit times. With this purpose, we use the second signals from the FIFOs to generate logical pulses in analog discriminators (discrimination level $``$ 1 VEM), then these pulses are delayed accurately to compensate for these differences in time and then they are fed to a majority logic coincidence unit to select the desired trigger condition. The time window in this module is set to 1.1 $`\mu `$s, covering safely the maximum time used by the EAS front to go across the array, even for the case of almost horizontal directions.

The digital oscilloscope available does not feature external trigger. For this reason one of the analog channels had to be used for triggering purposes, in addition to its signal digitizing function. With this purpose the STROBE signal generated by the coincidence unit, indicating the production of an event of interest (in practice 3 or 4-fold coincidences) is delayed about 8 $`\mu `$s after arrival of the last detector signal and then summed to one of the detector channels (channel 4 in Figure 5 and 7). Because of the relatively low singles counting rates and with the introduced delay of 8 $`\mu `$s no overlaps are produced in practice. Provided that the STROBE pulse is summed with opposite polarity respect to the detector signals the Advanced Trigger feature of the oscilloscope could be safely used for triggering.

When the STROBE pulse is detected by the oscilloscope, the SAVE procedure is initiated, i.e., the traces stored in the four channel memories corresponding to the last 16384 ns, are frozen and transferred to the PC disk. This time slice allows us to obtain a good measurement of both, the desired detector signals and the unavoidable radio noise pick-up in the long cables carring the signals.

The system dead time (digitalization, data transfer and PC storage) is 22 seconds for event. This relatively long dead time is primarily produced by the transfers, through the RS-232 serial port working, at 19200 bps. This dead time is considered acceptable in comparison with the average time between events, which is of the order of 6 minutes.

The first time region of 8192 ns (up to the first cursor in Figure 7) is used to compute the bias level at the time of presentation of the detector signals. The typical pick up noise appears as a dominant oscillation with a period of the order of 1 $`\mu `$s, corresponding mainly to the AM broadcasting stations. The following 8192 ns region, between the cursors, is the time region where the detector signals are stored. The last region which contains the STROBE signal is not saved to disk. The internal 150 MHz bandwidth low-pass filter built in the oscilloscope is active in order to reduce the amplitude of higher frequency signals. A fast Fourier analysis of the detector signals indicated that their main harmonic components extend up to $``$ 100 MHz, thus little distortion in the detector signals is introduced by the filter.

A special program was written to drive the data acquisition system in a completely automatic way. Normal collection of shower events is performed when the program runs in “Survey Mode”. The detection of a STROBE pulse causes that the oscilloscope traces recorded in the 4 channels are saved to disk, together with information on the year, day, and local civil time, which allows to reconstruct the equatorial or galactic coordinates of the shower arrival direction.

In addition, every day the program switches at a predetermined time to the “Calibration Mode”. In this mode, the collection of data from each detector station is self-triggered, in order to record background events to calibrate the stations, i.e. to determine the daily VEM value for each station. The four detectors are measured sequentially in this mode, under program control. The data are stored to disk and analyzed off-line. Roughly, one hour and a half is required to acquire 3000 background events (found to be adecuate to obtain the VEM values with an error of $``$ 5%) for each station and to save the calibration data from the four detectors. The starting time for the calibration procedure, and the total amount of background events for each detector, are set in an ASCII file. Once the calibration is completed, the program automatically switches to the “Survey mode” described above. This mode of operation is kept until the “Calibration Mode” is called up again, at the programmed time next day.

The singles counting rates of the four detector stations are permanently recorded using a CAMAC scaler with a refreshing time of 1 s, and are also saved to disk. This information is valuable for monitoring the status of each station. It helped discarding particular data when the operating condition of a particular station became unstable due, for instance, to a high level of pick up noise, or gave an alert signal for the need of maintenance of a station, in occasional cases of light leaks. The recording of the counting rates is also program-controlled and does not require operator action to run, once it is launched.

## 3 Simulated performance of the array

In order to characterize the behavior of the array, detailed simulations were performed to estimate its efficiency for shower detection and its angular and energy resolutions. A special routine, simulating the detector response to the different shower particles, has also been written to provide an input for the reconstruction routines.

### 3.1 Shower database

The AIRES program  using the SYBILL hadronic package was used in the first step of the simulation pipeline: the construction of an adequate shower database containing detailed information about the secondary particles at ground produced by primary cosmic rays of the energies of interest.

The shower simulation starts with the injection of a primary particle in the high atmosphere ($``$ 100 km above sea level) and tracks down the different generations of secondary particles in the subsequent cascade. The technique known as thinning was used to reduce the CPU time and the disk storage requirement. This procedure consists in tracking explicitly all particles above certain energy threshold (or thinning energy), and those particles having an energy level below the threshold are computed using statistical weight. In our simulations we set a relative thinning energy level of 5.10-5 with respect to the primary particle energy.

To construct the shower database, twenty primary energies ranging from 1014 eV to 1018 eV were selected, and two nuclear species (protons and iron nuclei) were considered as primary particles. They were injected at zenithal angles from 0 to 60, in 15 steps. To reduce the artificial fluctuations due in part to the thinning method , and also to obtain representative values of the relevant parameters, batches of 100 showers were simulated under the same initial conditions (as described above), and their average and RMS values were used. All these simulations were performed considering a ground level of 15 m.a.s.l.

The AIRES program produces a set of tables written in ASCII code, in which the information referent to the secondary particles reaching ground level (after deconvolution of the thinning algorithm) can be expressed simply as particle densities as a function of the core distance. This function is called lateral distribution function (LDF). In addition, for those particles “reaching” the ground level, these tables provide the landing time as function of the shower core distance and also their energy distribution. These tables include the mean, RMS and extreme values for each computed variable bin.

Figure 8 shows a typical example of AIRES results for the LDF and for the energy distributions produced by a proton primary of 2.1015 eV impinging at a zenithal angle of 30.

Thus, by running protons and iron nuclei as primary particles, the AIRES database contains an amount of 20000 simulated showers covering the energy and zenithal range of interest for the TANGO Array. The shower database tables contains only particle densities, energies and arrival times (with respect to the core position particles arrival times) for muons, electrons, and gamma-rays.

### 3.2 Array simulation procedure

In order to predict the response of the array, the information on showers contained in the AIRES database tables was used to simulate events , i.e. the effect of individual showers falling relatively close to the array.

A simulated shower is a set of information describing in detail the calculated number and properties of secondary particles reaching the ground level. A simulated event is the set of information about the calculated effect of the shower on the array, taking into account the simulation of the detector, DAQ hardware, electronics, etc. With this purpose the showers are read from the AIRES database tables, establishing the number, energy and type of particles hitting each detector station, and their relative arrival time. The result of an event is a set of electronic signals from the detectors which are stored in computer memory. The simulation of each particular event allows to determine if it triggers or not the DAQ system. A total number of 360000 “events” , including all primary species and energies have been simulated from the information contained in the AIRES tables and they constitute the simulated events database. The procedure to simulate one event is described as follows:

* For each primary energy, 9000 shower core positions were selected landing at random in a area larger than the geometrical area of the array. In this way we can estimate the effect of showers falling outside the boundaries and obtain an estimation of their triggering efficiency. The size of this landing area was scaled logarithmically with the primary energy with the purpose to take into account the increase of the shower size at ground with energy.
* For each core landing position the zenithal and azimuthal angles of the event were chosen as follows: the azimuthal angle was uniformly distributed, and the zenithal angle distribution can be described by a cos3($`\theta `$) function, with a cut-off at 45 . This cut-off was selected accordingly with the atmospheric depth at Buenos Aires, where most EASs arrive within a cone of $``$ 40 . The exponent of the distribution was chosen so as to produce a distribution flatter than the flattest one reported up to date . This was done with the purpose of including in the database a number statistically significant of simulated events at higher zenithal angles.
* Once the event core position and the angles were established for each particular event, the distances from each detector station to the core were calculated. Then, from the AIRES tables the apropiate mean values and dispersions were extracted and interpolated to reproduce the simulated event.
* To include the shower-to-shower fluctuations, uniform random number generators were profiled (using the accept-reject technique)  to reproduce the mean value and dispersion of the AIRES particle density tables contained in the database, according to the particular secondary particle considered. With the AIRES tables interpolated to the particular conditions of each simulated event, and the modified random number generators, the densities (particles/m2) of muons (both charges), electrons (both charges) and gamma-rays hitting each detector neighbourhood, were obtained. Finally, these density values were scaled according to the geometrical area of each detector to obtain the number of particles falling over each detector in each shower.
* The energy and arrival time of every individual particle hitting each detector station were obtained using the same procedure (the accept-reject technique). This was made taking by into account the particle species and its distance to the core.
* Once the number of particles, energies and arrival times of all particle species falling on each station for the event, were obtained, the detector signal was obtained as is described in detail in 3.2.1.
* The next step in this calculation was the simulation of the response of the data acquisition electronics by performing a check to determine whether each particular simulated shower produces or not a valid trigger. With this purpose, the simulated traces for each detector station were scanned, searching for the threshold crossing times in each channel (there could exist multiple crossings in a single event). Then, these threshold crossing times, determined in each channel, were compared to those corresponding to the other channels to establish the presence of temporal coincidences between the traces (an EAS). If multiple crossing times were present in one or more channels, each one of them was searched for time coincidence with the other channels. The threshold levels in all channels were set as equivalent to the signal amplitudes produced by 1 VEM from each particular detector, and the time window for the coincidences was set to 1.1$`\mu `$s, in correspondence to the real situation during measurements. If a coincidence condition was found, then the event was classified accordingly to the number of stations involved in the coincidence.
* Finally, the behavior of the A/D converter stage was also simulated, featuring a FADC working at 500 Ms/s (like that used in the data acquisition system). An appropriate noise generator has been included. From noise spectrum measurements we concluded that the local AM radio stations are the main noise sources, contributing with $``$ 15 to 30 mV to the signal (the typical signal amplitude corresponding to one single particle is $``$ 100 mV). The noise spectrum can be described as a continuous distribution with superposed, strongly varying peaks corresponding to the well-known local AM broadcasting frequencies, ranging from $``$ 550 to $``$ 1650 kHz. The FM band is also seen in the noise spectrum. However, its amplitude is much lower and can be safely ignored. On this basis, in order to obtain a realistic simulation, we added to the simulated signals a noise spectrum which follows the description given above.

All computer programs required for the simulation pipeline (except AIRES) were especially developed in the present work.

#### 3.2.1 Detector simulation

A simple and very fast simulation program was written to emulate the detector response. In this program, instead of simulating in detail the production and transmission of the Čerenkov photons emitted during the passage of charged particles through the water, we used the detailed knowledge of the detector behavior achieved during the previous years of operation of the first prototype. On this basis, the detector response was reproduced accordingly to a large set of measured parameters. In the following, prior to describing the detector simulation program, we present a summary of the experimental data, obtained previously.

In previous experiments  the response of the WCD to vertical and tilted muons has been observed in detail. In these experiments, the entrance and exit points of the muons on the detector surface have been carefully selected to cover as much as possible all possible situations. A total of 38 different particle track lengths have been measured, corresponding to 162 different situations derived from the symmetry properties of the detector configuration. A particle track is considered to be “fully contained” when the entrance point of the muon is anywhere on the lid and the exit point is in the bottom, or when the entrance and exit points are near diametrally placed on the lateral cylindrical wall of the detector. Either an entrance or exit point on the side wall, and the other in the bottom or in the lid, are considered to produce a “clipping corner” track. As a result of these measurements (which are summarized in the Figure 9) we have found that the sum of the charges collected in the three PMTs of our WCD is, very approximately, directly proportional to the track length of the particle in the water radiator, and this is valid regardless of the entrance point position or the zenithal angle of the track.

For all measured tracks the digitized pulse shapes were recorded. The rise and fall times remain almost constant for the whole range of track lengths, which might be understood from the fact that these parameters are primarily determined by the highly diffusive properties of the Tyvek liner. These measurements have also shown that the fluctuations of the measured parameters (rise and fall times, voltage amplitude, and charge) are not larger than about $`\pm `$ 10 % of their mean value. These results were supported by GEANT  simulations, performed previously.

In addition to the response to fast muons, the response of the WCD to fast electrons and gamma-rays was obtained. Both, electrons and gamma-rays, produce also an amount of light proportional to their track lengths. It should be taken into account that gamma-rays are detected through their interaction with the water going essentially through pair-creation processes. This is the most probable case, given the relative cross sections at the typical gamma-ray energies present in the EAS. Therefore, the signals produced by gamma-rays are roughly the same as those produced by fast electrons, provided their energy distributions are similar (see Figure 8).

It should be taken into account that the signal produced by a muon with energy higher than $``$ 400 MeV becomes indistinguishable from the signal produced by an electron with energy higher than $``$ 250 MeV . Hence, these values were used to normalize the signals from electrons to those corresponding to muons.

In order to include in the simulations the effect of the signal distortions in the cables, we have recorded in a previous work the average pulse shape for vertical muons transmitted through 200 m of RG-213 cable.

By taking into account all this information, the simulation of the surface detector signal was carried out as described below:

* Muons: For each muon hitting a detector station, a zenithal angle is selected using a gaussian-shaped random number generator, with its mean value centered in the zenithal angle of the primary particle of the EAS, and a sigma value of 4. Hence the particles are restricted to an angular range of about $`\pm `$ 25 . Once the zenithal angle is established, the range of the particle in water is obtained according to its energy, and a peak amplitude is found as a function of its range. If the range of the muon exceeds the track length inside the surface detector, then the amplitude is made proportional to the track length. Finally, rise and fall times are selected with a gaussian shaped random number generator and the pulse shape is written to memory, considering the respective time delay from the AIRES results (again conveniently spread using a gaussian random number generator). The signal peak amplitudes, as well as the rise and fall times, are established also from gaussian-shaped random number generators, with their relative sigma-values obtained from measurements.
* Electrons: The general procedure is similar to that described for muons. The main difference occurs in the calculation of the range, which, in the case of the electrons, is assumed to be completely contained within the WCD, i.e. no backscattered electrons are simulated. The values of the peak amplitudes are obtained from electron simulations performed previously using the program GEANT.
* Gamma Rays: The energy of the $`\gamma `$-rays originated in an EAS range from $``$ 10 MeV to $``$ 100 MeV, and the main interaction channel in water goes through the pair creation process. In this energy regime, the mean interaction length of gamma-rays in water is about 80 cm. The track length for a specific gamma-ray (which depends on the zenithal angle selected as described above), determines the probability of creation of an electron-positron pair. In this case the electron simulation routine is called with two electrons, having a total energy balancing the gamma-ray energy. The energy of the recoiling nucleus is neglected.

The resulting final program is very fast; once the AIRES tables are locally available in a 233 MHz PC running under Linux, it simulates an average of 100 events/minute, and produces realistic pulse shapes as shown in Figure 10, where a real shower is compared with a simulated one.

As a summary, we have simulated the output of the digitizing electronic stage, reproducing the detector signal on the basis of previously known parameters relating the underlying physics with the detector response. The effect of the cables on the signal shape and the influence of the pick-up noise (only for the AM band), have been considered.

## 4 Shower Reconstruction

The reconstruction of the showers aims to find the direction of the shower axis, and to make an estimation of the energy of the primary cosmic ray. This reconstruction could be made in principle by performing a careful evaluation of a number of parameters which are measured from the oscilloscope traces, and from a comparison with the results of the simulations.

The reconstruction procedure is initiated by the obtention of the direction of the shower axis by fitting the arrival times to each detector, asuming a flat shower front. Once the direction is determined, the core position is found through minimization of the lateral distribution function using the particle density falling over each station. Then, using Monte Carlo simulations, it is possible to correlate the shower primary energy with the particle density measured by the detector stations.

### 4.1 Reconstruction of the shower direction

The reconstruction of the direction is based on the arrival times of the shower front particles to each detector station in the array. In order to determine the “trigger time” in each station, the voltage signal is time-integrated, and the crossing times of charge amplitude values equal to 10%, 50% and 90% of the maximum collected charge are determined, with the condition that no dynamic range saturation occurs. These times are called $`t_{10}`$, $`t_{50}`$ and $`t_{90}`$, respectively. These parameters behave like constant fraction discriminators crossing times, and they are valuable for the comparison of the overall time structure of the station’s signals when different particle densities are measured.

The $`t_{10}`$ are good indicators of the arrival time of the shower front to the detectors and they are used to obtain the shower axis direction, which is coincident with the primary cosmic ray arrival direction. On the other hand, the $`t_{50}`$ and $`t_{90}`$ are more closely related with the time structure and temporal width of the shower than with the shower direction, and can be used to estimate the primary mass composition and also the core distance .

If three ground stations detect a shower, its axis can be determined by triangulation with reference to the arrival times and to the positions of the stations. This is made by searching for a unique, downward-going shower front, which we assume to be a plane, and moves at the speed of light. When all four detectors are hit, then a least squares method is used to find the best fit to this plane shower front. More elaborate and detailed algorithms can be used to obtain the shower direction including, for instance, the radius of curvature of the shower front. However, if we use the available data obtained from our four stations and try to use more complicated algorithms they, very often, fail to converge.

#### 4.1.1 Angular resolution

As described above, the $`t_{10}`$ values were obtained from the simulated events database and used to obtain the arrival direction of each event. From this reconstruction the $`\theta `$ and $`\varphi `$, (zenithal and azimuthal angles) were obtained. These angles are the spherical angular components of a vector, normal to the (assumed) plane shower front. The accuracy in the reconstruction is determined by comparing these angles with the “true” angular direction of the particular simulated event, which is read from the events database.

As can be seen in Figure 11, the angular resolution ($`\sigma `$) of the array improves progressively with energy in the decade of 1014 eV, then remains almost constant in the decade 1015 eV and slowly decreases beyond $``$ 1016 eV.

The reconstructed plane is, actually, a plane parallel to the plane tangent to the shower curved front surface crossing the array at its center point. Beyond $``$ 1016 eV the shower front disk is much larger than the geometrical size of the array and the probability of having the shower core falling away from the array, and being still able to produce a trigger, is higher than the probability for the core to fall closer. Because of the finite radius of curvature of the shower front, the vector normal to the tangent plane is more tilted, respect to the shower axis, at points lying far away from the core than for points closer to the core.

Typical shower front curvature radius are in the order of 10 km, hence at a distance of 300 meters from the geometrical center of the array, the normal to the tangent plane is tilted $``$ 2 respect to the shower axis. This angular difference between the shower axis and the reconstructed direction have to be added to the intrinsic angular resolution of the array, which is of the same order of magnitude. This effect might explain the decrease in angular resolution at higher energies. Also, this may be interpreted as an energy limit for the validity of the flat shower front assumption, given the size of our array.

### 4.2 Reconstruction of shower energy

The axial symmetry assumption for an EAS is relevant for the energy analysis. It means that in a plane perpendicular to the shower axis, the particle density only depends on the radial distance from the axis. On the ground plane this symmetry is lost (unless the shower is vertical). However, for moderate zenithal angles ($``$ 40) the assumption of a symmetric distribution is a valid approximation. For instance, if the EAS has a zenithal angle of 40 and the shower front has a diameter of 300 meters (typical size of the TANGO array) the forward component of the shower front travels about 250 meters more than the backward component to reach the ground. Considering the measured attenuation length reported in the Haverah Park experiment of (780 $`\pm `$ 35) gr/cm2, the forward component of the shower front would be attenuated only $``$ 5% with respect to the backward component. This calculation shows the validity of the approximation. On this basis, we assume in the following that axial symmetry is a valid assumption.

A key parameter required to estimate the energy of an EAS is the LDF, i.e. the particle density as a function of the distance to the core position. From the results of previous experiments  it is possible to propose the functional dependence:

$$\rho =\frac{A}{r^{\eta +r/r_0}}$$
(1)
where $`\rho `$ is the particle density (\[VEM/m2\]), $`r`$ is the distance to the core (\[m\]), $`A`$ is a normalization constant (proportional to the primary particle energy) and $`\eta `$ and $`r_0`$ control the shape of the LDF. The last two parameters were obtained by fitting the previous expresion to simulated particle density distributions, which have included the detector response to different shower particle species ($`\mu ^\pm `$, $`e^\pm `$ and $`\gamma `$-rays) as described in 3.2.1.

Figure 12 shows the fits using Equation 1 to the simulated particle densities “measured” with the simulated WCD for several primary energies, where all zenithal angles included in the simulation were averaged within each core-distance bin. The small “plateau” observed in the leftmost part of the 2.1016 eV curve is produced by the (simulated) electronic dynamic range saturation.

It should be noted that the $`\eta `$ parameter is slightly sensitive to the primary particle mass  as was found by fitting the previous expresion to the simulated data events. The reproduction of such dependence from simulations is an encouraging result. Although different approaches were attempted to obtain, at least, a primary mass indicator from the reconstructed events in the simulated database, none of them were satisfactory, probably due to the simulated shower-to-shower fluctuations that might mask the small differences in the $`\eta `$ parameter for different primary species. Because of this we used in the following an average value for the $`\eta `$ parameter between the values corresponding to proton and iron primaries (see Table 1).

It is known from extensive Monte Carlo simulations that there exists a certain distance from the shower core for which the particle density of an EAS correlates with its primary energy and also their fluctuations are minimized. In the present experiment performed with only 4 detectors, we have used a simplified model where the normalization constant $`A`$ of LDF was correlated with the primary energy instead of the particle density at a fixed distance of the core position. The LDF was obtained from particle density measurements in each detector station, far away from the core.

The normalization constant of the LDF is found through minimization of the following equation

$$\chi ^2=\underset{i=1}{\overset{n}{}}\left(\rho _i\frac{A}{r_i^{\eta +r_i/r_0}}\right)^2$$
(2)
where $`\rho _i`$ and $`r_i`$ are the particle density and the distance between the core impact position and the i-th station, respectively, and $`\eta `$ and $`r_0`$ were obtained from simulations as mentioned before. Finally, the particle density measured by each station is obtained by the ratio of the time-integrated oscilloscope trace (bias subtracted) and the VEM value corresponding to that particular detector station. This ratio yields the number of equivalent particles falling in the station. Then, the equivalent particle densities (VEM/m2) is obtained by simple normalization to the respective detector area.

The miniminization of Equation 2 was performed through a grid search on the simulated data of the events database, yielding the $`x`$ and $`y`$ coordinates of the core position, as well as the normalization constant $`A`$. Table 2 shows the accuracy of the core position reconstruction obtained by this method for some energies. The accuracy is degraded at higher energies, probably because the shower front size becomes comparable with the array size.

#### 4.2.1 Correction by atmospheric attenuation

The axial symmetry asumption proved to be a valid approximation regarding the forward and backward components of the shower front for non-vertical showers. However, the effect of atmospheric attenuation in a tilted EAS development cannot be ignored. In order to get an evaluation of the magnitude of this effect we used the simulated events database to estimate the atmospheric attenuation on the shower propagation through the atmosphere.

By simple geometrical considerations it is possible to propose a functional dependence of the form

$$N=Ae^{\left[\beta (sec(\theta )1)\right]}$$
(3)
where $`N`$ is a normalization factor, proportional to the primary particle energy that includes the atmospheric attenuation correction factor and this constant includes the atmospheric attenuation correction factor.

For each primary energy, the simulated showers were divided in zenithal angle bins of 5 each, and a fit was performed to the data using the functional dependence shown in Equation 3, i.e., assuming only a dependence for the $`A`$ parameter on the zenithal angle. The average value obtained for $`\beta `$ by fitting the simulated data to Equation 3 is $`\beta =4.1\pm 0.1`$. For the zenithal range of interest ($`\theta `$ 30) this correction because of the atmospheric attenuation increases the estimated EAS’s energy up to $``$ 50%.

#### 4.2.2 Primary energy assignment

Finally, after minimization of Equation 2 and being performed the atmospheric attenuation correction (for which the directional reconstruction is required) it is possible to show the relationship between $`N`$ -a parameter obtained from the shower reconstruction routine- and the primary energy (obtained from the simulated events database). It should be noted that in this survey over the simulated events database we found that, beyond $``$ 2.1016 eV, $`N`$ fails to converge, and the linearity (in logarithmic scale) as a function of the primary particle energy is lost. Therefore, only data at lower energies are shown in Figure 14.

From these fits we obtain the following expressions, useful to correlate the parameter $`N`$ \[VEM/m2\] with the primary energy \[eV\]:

$$E_0=(4\pm 1)10^9N^{1.17\pm 0.03}$$
(4)
and

$$E_0=(5\pm 2)10^9N^{1.20\pm 0.03}$$
(5)
where Equations 4 and 5 correspond to proton and iron primaries, respectively.

From these equations it is possible to estimate the relative error in the energy reconstruction. Even when the expression has a dependence on the $`N`$ value, its dependence is logarithmic, and the $`\mathrm{\Delta }`$$`N/N`$ value was found to be $``$ 0.4 from the simulated database. This yields a relative error of 57% and 66% for protons and iron nuclei, respectively, in the energy range from $``$ 1014 eV to $``$ 1016 eV.

According to these results, a knowledge of the primary particle mass would be required to correctly correlate the N parameter with the primary particle energy by choosing the proper expresion. Strictly, this fact prevents us to make an unambiguous assignment of the primary energy. Furthermore, it should be recalled that both Equations 4 and 5, were obtained from surveys performed on the Monte Carlo simulations, which are dependent of the particular hadronic package utilized. On the other hand, however, the results obtained from both expresions are consistent within errors.

## 5 Summary

A new, Extended Air Shower Array has been constructed in Buenos Aires during 1999 and was commissioned in 2000. It consists of 4 Water Čerenkov Detectors, three of them are arranged in a triangular shape and the fourth is near the center of the triangle. The enclosed area is $``$ 30.000 m2. The detectors placed in the vertices of the triangle have a footprint area of 10 m2, the central detector has 0.5 m2.

Detailed Monte Carlo simulations of the showers were performed using the AIRES code with the SYBILL hadronic package. Various computer programs and routines were developed to simulate the array response including the surface detector, front end electronics, pick-up noise, and triggering. It should be noted that an effort was made to use experimental data whenever possible. The simulated events database contains a total of 360000 events.

A reconstruction routine has been developed from the simulated shower database. According to the simulations, the angular reconstruction resolution is better than 5 in the range 5.1014 eV to 1017eV. The accuracy expected in the energy resolution is roughly 60% in the range $``$ 1014 eV to $``$ 1016 eV. With respect to the primary mass determination it is concluded from the present simulations that no unambiguous assignement can be made, at present, from the showers measured with our array.

A fully automatic system for calibration, monitoring and data acquisition has been built using standard NIM and CAMAC modules and a 4-channel digital oscilloscope connected to standard PCs. Data have been continuously collected since September, 2000 and the shower reconstruction analysis will be published in an forthcoming paper.

## 6 Acknowledgements

We are very especially indebted to the late J.Vidallé for his unvaluable dedication and help in the early stages of the TANGO Array.

We are also deeply indebted to D.Simoncelli and E.Fisher for their outstanding work at the Mechanical Workshop in TANDAR Laboratory. Also, we would like to express our gratitude to P. Stoliar, H. Di Paolo, C.Bolaños, J.Fernández Vásquez and O.Romanelli, for their help with different aspects of the electronic system. Thanks are given to M.Figueroa and M.Wagner for their work in the characterization of the PMTs. We would like to mention H.Grahmann, O.Ruiz, E.Altmann, A.Ferrero and A.Etchegoyen. They helped us in many different ways. We also thank Prof. Ma Yu Quian, from Beijing University, for the donation of the 3-inch PMT used in the central detector and to Plásticos Industriales S.A., specially E.Carricondo and P.Martelli for the donation of the fiberglass-reinforced tank. Finally, we would like to thank Fermilab for the loan of some electronic modules used in the experiment (Fermilab Loan C96082).

The work of P. Bauleo, C. Bonifazi and A. Reguera was supported by different CNEA fellowships. This work was partially supported by a CONICET grant (PIP 4446/96).