A STATISTICAL METHOD FOR RECONSTRUCTING THE CORE LOCATION OF AN EXTENSIVE AIR SHOWER

A common method of finding the energy, type, and direction of a high-energy primary cosmic ray is to examine the extensive air shower (EAS) which it produces. Surface arrays, which can be used to study EASs, usually record the pulse height and arrival times of secondary charged particles. Pulse height can be a density measure of the passing secondary particles through any surface detector of a surface array. The beginning time of a pulse can be assigned to the first particle passing through a detector. Frequently, a minimum number of array detectors should be triggered to record more than a threshold number of particles. To find the arrival direction of a recorded EAS in a surface array, at first approximation, a plane wave front will be fitted to the forehead front of an EAS in the so-called plane wave front approximation (PWFA). Sometimes, in a second approximation, the curvature of the front of an EAS and its finite thickness (especially far from shower axis) are also taken into account. The errors for the arrival direction range from a few degrees for a small sparse array to less than a degree for a sophisticated designed array (Acharya et al. 1993). The next important parameter of an EAS is its size or total number of charged particles. In conventional methods, shower size estimation is intimately related to the core location determination and is usually performed together. In this procedure, a so-called lateral distribution function (LDF), which is defined to predict the density of the charged particles as a function of distance from shower axis, size, and age parameter, will be fitted to the measured number of particles in each triggered detector during every EAS event.

Traditionally, a variant of the NKG function (Greisen 1956; Greisen 1960 and Nishimura 1967) should be used as LDF. Normally, LDF depends linearly on the shower size, and therefore the shower size can be easily found from the derivation of the likelihood (chi-square) function:

$$\begin{eqnarray*}{\chi }^{2} & = & \displaystyle \sum _{i=1}^{N}\displaystyle \frac{{({N}_{e}f({{\boldsymbol{r}}}_{i})-\rho ({{\boldsymbol{r}}}_{i}))}^{2}}{{\sigma }_{i}^{2}},\\ \displaystyle \frac{\partial {\chi }^{2}}{\partial {N}_{e}} & = & 0,\end{eqnarray*}$$

$$\begin{eqnarray}\Rightarrow {N}_{e} & = & \displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}f({{\boldsymbol{r}}}_{i})\rho ({{\boldsymbol{r}}}_{i})/{\sigma }_{i}^{2}}{{\displaystyle \sum }_{i=1}^{N}f{({{\boldsymbol{r}}}_{i})}^{2}/{\sigma }_{i}^{2}},\\ \Rightarrow {\chi }^{2} & = & \displaystyle \sum _{i=1}^{N}(\displaystyle \frac{{\displaystyle \sum }_{j=1}^{N}f({{\boldsymbol{r}}}_{j})\rho ({{\boldsymbol{r}}}_{j})/{\sigma }_{j}^{2}}{{\displaystyle \sum }_{j=1}^{N}f{({{\boldsymbol{r}}}_{j})}^{2}/{\sigma }_{j}^{2}}f({{\boldsymbol{r}}}_{i})\\ & & -{\rho ({{\boldsymbol{r}}}_{i})\phantom{\displaystyle \frac{{\displaystyle \sum }_{j=1}^{N}f({{\boldsymbol{r}}}_{j})\rho ({{\boldsymbol{r}}}_{j})/{\sigma }_{j}^{2}}{{\displaystyle \sum }_{j=1}^{N}f{({{\boldsymbol{r}}}_{j})}^{2}/{\sigma }_{j}^{2}}}\phantom{\rule{-8.6em}{0ex}}{)}^{2}/}^{}\phantom{\rule{-6.5em}{0ex}}{\sigma }_{i}^{2},\end{eqnarray} \tag{ 1 }$$

where $\rho ({{\boldsymbol{r}}}_{i})$ is the measured density at location ${{\boldsymbol{r}}}_{i}$, ${\sigma }_{i}$ are the uncertainties related to the density measurements and particle fluctuations, and N_e is the shower size.

Consequently, if the arrival direction is assumed to be known, the likelihood is a function of three variables: distance to the axis (two components) and age parameter. Since LDF is usually highly nonlinear in distance to the axis and age parameters, we cannot find them analytically. Numerical techniques for finding the core location and age parameter require knowledge of their range and the initial estimation of their values. Frequently, the age parameter ranges from 0.4 to 1.6 and its initial trial value is taken as 1.0. So, the only remaining parameter whose initial trial value is important for estimating shower parameters is the core location. The initial trial core location (ITCL) is usually chosen at the location of the highest measured density (LHMD) or the density-weighted average of the detector coordinates (DWADC) that respond.

It may seem that finding the global optimum does not require any initial estimation. But finding the global optimum of a function (e.g., chi square) is far more difficult than local optima: analytical methods are frequently not applicable, and the use of numerical solution strategies often leads to very hard challenges. However, if we have a good estimate of the initial values of the parameters, we can be certain about reaching a global minimum by local search strategies (especially when our model function is a good description of the observed data). Therefore, the initial trial estimation of the shower parameters is essential for the reconstruction procedure.

Recently, a new statistical method named SIMEFIC has been established for finding a core location that uses arrival time information in addition to density information (Hedayati et al. 2009, 2011a). The SIMEFIC-refined version, called weighted center of gravity (WCG, from now on, we call it "previous WCG" to avoid ambiguity; when "previous" is not used, it means only the WCG, not a special method) allows us to find core location of an EAS in the internal region of a surface array better than other methods (Hedayati et al. 2011b). The only weakness of the previous WCG is its inaccuracy near the margin of the array.

In this paper, a refined form of the previous WCG is proposed that has the same precision as the previous WCG in the central region of the array and resolves its erroneous behavior near the array boundaries. The refined WCG (called SIMEFIC II) method is the best known method for finding core location (at least for a square array). Moreover, it can be safely used as the ITCL in all parts of a surface array.

The principles of SIMEFIC II are essentially the same as those of the SIMEFIC method and the previous WCG method, which have been explained in detail in Hedayati et al. (2011a, 2011b). Here, we bypass the physical principles and only describe the procedure. In the SIMEFIC II method, in addition to the density information of secondary particles in each responding detector of the array (during an EAS event), the so-called three-dimensional (3D) distances between them are also taken into account. In fact, the 3D distance between two detectors is the distance between the first crossing particles from each detector. In fact, surface detectors only record the arrival time of the first crossing particles.

For a vertical EAS, finding the 3D distance between two detectors is fairly straightforward: at first, detectors are indexed according to their responding times. If we know the horizontal distance between two responding detectors, i and j, as ${r}_{{ij}}=\sqrt{{({x}_{i}-{x}_{j})}^{2}+{({y}_{i}-{y}_{j})}^{2}}$ and the difference between the arrival times of the first crossing particles in each detector as ${\rm{\Delta }}{t}_{{ij}}={t}_{i}-{t}_{j}$, then the 3D distance between the first crossing particles in each detector will be

$$\begin{eqnarray}&&{d}_{3{ij}}=\sqrt{{r}_{{ij}}^{2}+c{\rm{\Delta }}{t}_{{ij}}^{2}},\end{eqnarray} \tag{ 2 }$$

where c is the speed of light, assuming that all of the secondary particles' directions are straight down.

For an inclined EAS, we only use the following transformations for finding 3D distances:

$$\begin{eqnarray}{x}_{i}^{\prime } & = & {x}_{i}\mathrm{cos}\theta \mathrm{cos}\varphi +{y}_{i}\mathrm{cos}\theta \mathrm{sin}\varphi ,\\ {y}_{i}^{\prime } & = & -{x}_{i}\mathrm{sin}\varphi +{y}_{i}\mathrm{cos}\varphi ,\\ {z}_{i}^{\prime } & = & -{x}_{i}\mathrm{sin}\theta \mathrm{cos}\varphi -{y}_{i}\mathrm{sin}\theta \mathrm{sin}\varphi ,\end{eqnarray} \tag{ 3 }$$

where θ and φ are the reconstructed shower zenith and azimuth angle, respectively. These transformations are the coordinates of detectors in a coordinate system with the reference plane perpendicular to the EAS axis. Then, we transform the particle arrival times as ${t}_{i}^{\prime }={t}_{i}+{z}_{i}^{\prime }/c$. Having performed these transformations, the 3D distances between two detectors will be

$$\begin{eqnarray}&&{d}_{3{ij}}=\sqrt{{({x}_{i}^{\prime }-{x}_{j}^{\prime })}^{2}+{({y}_{i}^{\prime }-{y}_{j}^{\prime })}^{2}+c{({t}_{i}^{\prime }-{t}_{j}^{\prime })}^{2}}.\end{eqnarray} \tag{ 4 }$$

If the number of detected particles in each detector is, respectively, n_i and n_j, a matrix element, ${m}_{{ij}}=({n}_{i}{n}_{j})/{d}_{3{ij}}$, will be attributed to the i and j pairs of detectors.

Next, we choose a pair of detectors with the largest matrix element, ${m}_{\gt }={m}_{{ij}}$. In this pair, if ${n}_{i}\gt ={n}_{j}$, then the ith detector will be selected as the nearest detector to the core location (note that, $i\lt j$). If ${n}_{j}\gt {n}_{i}$, then the jth detector will be chosen instead. Assuming the ith detector is selected, the ith detector is statistically the closest detector to the core location and we assign a weight ${w}_{i}={m}_{{ij}}$ to it (an indication of its closeness to the core position).

Then, the ith detector will be removed from the list of the responding detectors and the same procedure is repeated for the remaining responding detectors until all the detectors are exhausted. After this procedure, one of the detectors (the last one) will remain without weight. We set the weight of this last detector equal to that of the previous.

Assume that there are N responding detectors during an EAS event. We find the WCG of all the responding detectors:

$$\begin{eqnarray}{x}_{\mathrm{WCG}} & = & \displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{x}_{i}{w}_{i}}{{\displaystyle \sum }_{i=1}^{N}{w}_{i}},\\ {y}_{\mathrm{WCG}} & = & \displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{y}_{i}{w}_{i}}{{\displaystyle \sum }_{i=1}^{N}{w}_{i}}.\end{eqnarray} \tag{ 5 }$$

This location is an approximation of the true core location. Before going any further, we should compare the results of this reconstructed core with the results of the previous WCG and try to optimize the method. So, in the next section, we describe the simulation method and then return to the SIMEFIC II method.

In order to verify the SIMEFIC II method, 100,000 simulated showers were generated by CORSIKA version 7.4 (Heck et al. 1998). EAS simulations were performed for geographical longitude, latitude, and altitude of $51{\rm{E}}$, $35{\rm{N}}$, and $1200\;{\rm{m}}$. The Earth's magnetic field was assumed as ${B}_{z}=38.4\mu {\rm{T}}$ and ${B}_{x}=28.1\mu {\rm{T}}$ throughout the simulations. The low energy hadronic model was Fluka version 2011.2b (Ferrari et al. 2005), the high energy hadronic model was QGSJETII-04 (Ostapchenko 2011). For other parameters, we used CORSIKA default values. The primary particles for $90\%$ of the EASs were protons and the remaining primary particles were alphas. The assumed array is the same as in Hedayati et al. (2011b): ($21\times 21$) detectors distributed over a square array with a network constant of $10\;{\rm{m}}$. The total area of the array was $200\times 200\;{{\rm{m}}}^{2}$. Most of the comparisons among different methods were performed according to the following procedure.

The true core location of each EAS was assumed to be on the diagonal line $(i,i)$ (with i increments from 0 (array center) to $100\;{\rm{m}}$ at the corner of the array by steps of $1\;{\rm{m}}$). At every point, the distances from the reconstructed core to the axis for those showers which satisfied the threshold condition (see below) have been averaged.

In the calculations, the detection condition for an EAS (threshold condition) responded to at least 12% (53 detectors) of all detectors in the array (unless stated otherwise). To find the arrival direction of an EAS, PWFA was used. If the coordinates of the true core in the CORSIKA simulated showers were denoted as $({x}_{\mathrm{tc}},{y}_{\mathrm{tc}})$, then the distance between the true axis and the reconstructed core location, $({x}_{\mathrm{rc}},{y}_{\mathrm{rc}})$, would be $d=\sqrt{{({x}_{\mathrm{tc}}^{\prime }-{x}_{\mathrm{rc}}^{\prime })}^{2}+{({y}_{\mathrm{tc}}^{\prime }-{y}_{\mathrm{rc}}^{\prime })}^{2}}$. ${x}^{\prime }$ and ${y}^{\prime }$ are the new coordinates in a perpendicular plane to the shower axis:

$$\begin{eqnarray}{x}^{\prime } & = & x\mathrm{cos}{\theta }_{t}\mathrm{cos}{\varphi }_{t}+y\mathrm{cos}{\theta }_{t}\mathrm{sin}{\varphi }_{t},\\ {y}^{\prime } & = & {-x\mathrm{sin}\varphi }_{t}+y\mathrm{cos}{\varphi }_{t},\end{eqnarray} \tag{ 6 }$$

where ${\theta }_{t}$ and ${\varphi }_{t}$ are the true shower zenith and azimuth angle, respectively (provided by CORSIKA). Note that in the array plane, the true and reconstructed core locations do not have z coordinates. Distance d is a measure of the error which occurred during the reconstruction procedure.

When an EAS true core location crosses the peripheral part of the array, an asymmetry effect takes place. The reason behind this asymmetry is that the array registers more shower information on one side of the core than on the other side (beyond boundary of the array where there are no detectors). In the previous WCG method, two length scales were introduced to decrease this asymmetry effect. However, as can be seen in Figure 1, this correction cannot remove all of this asymmetry and, although the marginal EASs (because of they satisfy the threshold condition) have more average energy than the central EASs (Hedayati et al. 2011b; which can partly compensate for their asymmetry effect), the error in finding the core location in the previous WCG is worse for them.

Zoom In Zoom Out Reset image size

4.1. Directional Correction for the Finding Core Location

While the introduction of the length scales can remove some of the imbalance effect in the reconstruction of the core location of those EASs impacted in the outer region of the array, the error in the outer region is still notable. In fact, when the true core of an EAS crosses the outer region of the array, the triggered detectors far from the core, which are recorded only on one side of the core, displace the reconstructed core in their average directions and the length scales cannot completely remove this effect. A refinement to WCG of all the responding detectors, which can be used to improve the asymmetry, is to find the core location of an EAS with some of the responding detectors with the highest weights and then use the information from the lowest weighted detectors to correct the asymmetry effect.

Suppose the locations and weights of the responding detectors sorted from highest weight to lowest are $({x}_{1},{y}_{1},{w}_{1}),({x}_{2},{y}_{2},{w}_{2}),\ldots ,({x}_{N},{y}_{N},{w}_{N})$. Assume that we select the first M ($1\lt M\lt N$) of them to find an initial core location:

$$\begin{eqnarray}{x}_{\mathrm{rc}}^{M} & = & \displaystyle \frac{{\displaystyle \sum }_{i=1}^{M}{x}_{i}{w}_{i}}{{\displaystyle \sum }_{i=1}^{M}{w}_{i}},\\ {y}_{\mathrm{rc}}^{M} & = & \displaystyle \frac{{\displaystyle \sum }_{i=1}^{M}{y}_{i}{w}_{i}}{{\displaystyle \sum }_{i=1}^{M}{w}_{i}},\end{eqnarray} \tag{ 7 }$$

where $({x}_{\mathrm{rc}}^{M},{y}_{\mathrm{rc}}^{M})$ is the position of the reconstructed core location by the M highest weighted detectors. Then, we use the remaining $N-M$ to find the following correction terms:

$$\begin{eqnarray}{x}_{{cv}} & = & \displaystyle \frac{{\displaystyle \sum }_{i=M+1}^{N}({x}_{\mathrm{rc}}^{M}-{x}_{i}){w}_{i}}{{\displaystyle \sum }_{i=1}^{N}{w}_{i}}\\ & = & \displaystyle \frac{{\displaystyle \sum }_{i=M+1}^{N}({x}_{\mathrm{rc}}^{M}-{x}_{i}){w}_{i}}{{\displaystyle \sum }_{i=M+1}^{N}{w}_{i}}\displaystyle \frac{{\displaystyle \sum }_{i=M+1}^{N}{w}_{i}}{{\displaystyle \sum }_{i=1}^{N}{w}_{i}}\\ & = & ({x}_{\mathrm{rc}}^{M}-{x}_{\mathrm{rc}}^{N-M})\displaystyle \frac{{W}_{N-M}}{{W}_{N}},\\ {y}_{{cv}} & = & ({y}_{\mathrm{rc}}^{M}-{y}_{\mathrm{rc}}^{N-M})\displaystyle \frac{{W}_{N-M}}{{W}_{N}},\end{eqnarray} \tag{ 8 }$$

where $({x}_{\mathrm{rc}}^{N-M},{y}_{\mathrm{rc}}^{N-M})$ is the WCG of the remaining low weighted detectors. W_N is the sum of all the weights and ${W}_{N-M}$ is the sum of the weights of the low weighted detectors. Then, we find the reconstructed core location by

$$\begin{eqnarray}&&{x}_{\mathrm{rc}}={x}_{\mathrm{rc}}^{M}+{x}_{{cv}},\quad {y}_{\mathrm{rc}}={y}_{\mathrm{rc}}^{M}+{y}_{{cv}}.\end{eqnarray} \tag{ 9 }$$

Interpreting the above equations is simple. First, we find an initial approximation of core location by WCG of the highest weighted detectors and then use the WCG of the remaining low weighted detectors to displace this initial core location in the opposite direction of the WCG of low weighted detectors. The reason for the minus sign in this term is to compensate for the asymmetry effect, which mostly results from the distribution of the low weighted detectors. However, the effect of the low weighted detectors should be multiplied by the fraction of their weights to the total weight of all of the detectors, which is called asymmetry correction. Of course, this term is only important for those showers where the core is impacted far from the array center and has the reverse effect for those showers with the true core near the array center of symmetry.

Figure 2 shows the results of the asymmetry correction for $M=N/2$ on the reconstructed core location precision. Also, in this figure, the results of the previous WCG and the WCG of all responding detectors are shown. Although the precision of WCG with the correction term is slightly worse than the two others in the central region of the array, it has better results in the marginal part of the array which has a greater area.

Zoom In Zoom Out Reset image size

As can be seen in Figure 3, the correction term has a significant effect on the precision of the reconstructed core location. In the central part of the array, it has slightly worse results than WCG without the correction term; however, in the outer part of the array, it has significantly better results.

Zoom In Zoom Out Reset image size

4.2. Optimized SIMEFIC II

Based on what was mentioned, the following procedure can be recommended. At first, we find an initial core location based on some M values (e.g., $M=N/2$). Then, we find the distance of this initial core location from the array center (its precision is better than half of the array network constant). This distance can help decide which fraction of N is better to use for M. However, this procedure is not universal for every array and should be separately designed. For our assumed array, using $M=N/2$ in all parts of the array produces at most $4\;{\rm{m}}$ error, which is not too high. However, if we refine the procedure, we can reach a precision of better than $3.2\;{\rm{m}}$ (less than one-third of the array network constant) in all parts of the array for the threshold condition of 12% (Figure 4).

Zoom In Zoom Out Reset image size

Another optimization is possible if we set the coefficient of x_cv and y_cv to be different from 1. Namely, we use the following equations instead of Equation (9):

$$\begin{eqnarray}&&{x}_{\mathrm{rc}}={x}_{\mathrm{rc}}^{M}+{C}_{x}{x}_{{cv}},\quad {y}_{\mathrm{rc}}={y}_{\mathrm{rc}}^{M}+{C}_{y}{y}_{{cv}},\end{eqnarray} \tag{ 10 }$$

where C_x and C_y are two constants which have different values in different parts of the array. Figure 5 shows the effect of different C_x and C_y on the results of SIMEFIC II ($M=N/2$). In the central part of the array, it is better to use C_x and C_y of less than 1 and, for the peripheral region, values of more than 1 are better.

Zoom In Zoom Out Reset image size

Every method for the reconstruction of an EAS must be stable against timing errors of detectors. Unlike the conventional methods, time error has two effects for the SIMEFIC II method: decreasing the precision of the arrival direction estimation and decreasing the precision of 3D distances. In order to check the stability of our method against time error, we added a normally distributed random number with an average of 0 and a standard deviation of $30\mathrm{ns}$ to the particle arrival times. Figure 6 shows the effects of the time error on the results of SIMEFIC II (from now on, we set $M=N/2$, ${C}_{x}=1$, and ${C}_{y}=1$ in the SIMEFIC II method). As can be seen, this method is fairly stable against time error. Also, in this figure, the effect of eliminating 3D distances on the precision of the method is depicted (dashed–dotted (green line)). For this purpose, the value of ${w}_{i}={n}_{i}{n}_{j}$ is considered for the weights of the detectors (${d}_{3{ij}}$ is removed from the denominator of the weights).

Zoom In Zoom Out Reset image size

The effect of time error on the results in the central and intermediate areas of the array is greater than in the marginal area. The reason for this behavior, according to Figure 6, is that 3D distances between detectors are more effective in the central and intermediate areas of the array than in the marginal area. In the marginal area of the array, 3D distances have a negligible effect on the reconstructed core location precision. The explanation of this behavior is that, contrary to the density information, the 3D distances suffer from the asymmetry effect, especially in the marginal region.

A natural question may arise about the role of 3D distances. To what extent are they useful in determining the core position? In Figure 6, it is observed that they can increase the precision of the reconstructed core by more than $1\;{\rm{m}}$ (in some part of the assumed array). Figure 7 shows the results of using ${w}_{i}=1/{d}_{3{ij}}$ for the weights of detectors and removing densities from the numerator. In this figure, we do not use the correction term and only apply the WCG of the responding detectors (correction term is not useful for this situation). The error of this method is about $5\;{\rm{m}}$ (about half of the array network constant) in the central part of the array. Therefore, we can find the core location with timing detectors (which cannot measure density of the crossing particles) with a precision of about half the spacing between detectors in an array (at least in the central part of the array). However, this method should be optimized for different percentages of the chosen detectors.

Zoom In Zoom Out Reset image size

In Figure 8, the SIMEFIC II method with $M=N/2$ results is compared with the two common methods used for ITCL. SIMEFIC II (even in non-optimized version) has far better results than DWADC. The results confirm that DWADC has the worst results, even in the central region of the array, and should not be used as ITCL. Also, SIMEFIC II has better results than LHMD in most parts of the array.

Zoom In Zoom Out Reset image size

Our method must be independent of the shower simulation methods in order to be applicable for the real EASs. In order to confirm that the results of our method do not belong to a special kind of shower simulation, we use the proposed method for different kinds of showers. Also, more than 100,000 simulated EASs are used, which have the same parameters as the described showers, except that the low- and high-energy hadronic models are GHEISHA (Fesefeldt 1985) and QGSJET01c (Kalmykov et al. 1997), respectively. Figure 9 shows the results of comparing the prescribed showers and the new ones. Apart from a slight fluctuation, the results are the same.

Zoom In Zoom Out Reset image size

The distance to the axis of a shower is a positive quantity (bounded from one side, zero), and using its average as a measure of the error in finding the core location overestimates the error. Therefore, we prepared the distribution function of distance to the axis. To find this distribution function, the true core location of a simulated EAS is randomly distributed over the whole array area. Two sets of uniformly distributed random numbers between 0 and 100 for the x and y components of the true core are generated. Each EAS event is used 10 times (i.e., the true core of each CORSIKA simulated EAS is randomly chosen on the array surface 10 times). Then, we find the reconstructed core location of each accepted EAS (an EAS which triggers more than the threshold number of detectors).

Figure 10 shows the results for the core location estimation error in SIMEFIC II for two threshold conditions, corresponding to 12% and 24% of all detectors. Also, the results of DWADC and LHMD are demonstrated for the same threshold conditions. As can be seen, the most probable distance between the true core and the reconstructed core location is approximately half of the average distance between them. Also, unlike LHMD, by increasing the threshold number of detectors, the precision of SIMEFIC II will become considerably better (if we use the most probable value as a measure of the core location error).

Zoom In Zoom Out Reset image size

The results of this method should be checked against different arrays and different threshold conditions. In Figure 11, a different assumed array with $11\times 11$ (121) surface detectors distributed over an area of $150\times 150\;{{\rm{m}}}^{2}$ is used; so, the array network constant is $15\;{\rm{m}}$. The threshold condition of at least 12% and 24% (15 and 29 detectors, respectively) of all detectors responding during an EAS event is applied. For every threshold condition, we compare our method with other common methods. In this figure, the comparison of the SIMEFIC II (M = N/2) results with the LHMD and DWADC methods is shown. The SIMEFIC II method has significantly better results than LHMD for both threshold conditions for most of the array area. By increasing the threshold number of responding detectors, our method and LHMD have considerably better results. As before, the DWADC method is still not recommended for finding the core location.

Zoom In Zoom Out Reset image size

In this paper, a new model-independent method was revised for finding the core location. This method, called SIMEFIC II in its optimized form, can find the core location of an EAS with the precision of better than one-third of an array network constant (at least for a square array). Although, unlike the conventional methods, this method uses the arrival time information of secondary particles, it is stable against time error.

At first glance, the core location may not be important on its own. However, finding it by a model-independent method is important from two points.

The only parameter whose initial value is crucial for fitting an LDF is the core location. Other parameters can be found analytically (shower size) or their initial values are usually the same (age parameter). People often use one of these two methods for finding ITCL: LHMD or DWADC. SIMEFIC II has far better results than these methods.

The second important point of finding core location by a model-independent method is that we can estimate the accuracy of an LDF for the prediction of EAS parameters. It is unlikely that an LDF with statistically better results for the core location than another has worse results for the shower size or age parameter.

The last point that should be noted is that we think the SIMEFIC method is now in its primitive form and extensive work should be done on it. Examples include the following. Weights can be ${w}_{i}={n}_{i}^{a}{n}_{j}^{b}/{d}_{3{ij}}$ with different a and b whose values should be optimized for time order. Also, we considered the first and second moments of weights (${w}_{i}={n}_{i}$ and ${w}_{i}={n}_{i}{n}_{j}/{d}_{3{ij}}$), maybe considering other moments (e.g., third moment ${w}_{i}=({n}_{i}{n}_{j}{n}_{k})/({d}_{3{ij}}+{d}_{3{jk}})$ or other forms of it: e.g., ${w}_{i}=({n}_{i}{n}_{j}{n}_{k})/({d}_{3{ij}}\times {d}_{3{jk}})$) give better results. As an another example, using weights can be applied for arrival direction correction (Hedayati et al. 2011a) with different weights for different detectors in PWFA.

The authors would like to thank Dr. M. Khakian Ghomi for the CORSIKA-GHEISHA-QGSJET01c simulated showers which were used in this paper.