The mass hierarchy with atmospheric neutrinos at INO

We study the neutrino mass hierarchy at the magnetized Iron CALorimeter (ICAL) detector at India-based Neutrino Observatory with atmospheric neutrino events generated by the Monte Carlo event generator Nuance. We judicially choose the observables so that the possible systematic uncertainties can be reduced. The resolution as a function of both energy and zenith angle simultaneously is obtained for neutrinos and anti-neutrinos separately from thousand years un-oscillated atmospheric neutrino events at ICAL to migrate number of events from neutrino energy and zenith angle bins to muon energy and zenith angle bins. The resonance ranges in terms of directly measurable quantities like muon energy and zenith angle are found using this resolution function at different input values of $\theta_{13}$. Then, the marginalized $\chi^2$s are studied for different input values of $\theta_{13}$ with its resonance ranges taking input data in muon energy and zenith angle bins. Finally, we find that the mass hierarchy can be explored up to a lower value of $\theta_{13}\approx 5^\circ$ with confidence level $>$ 95% in this set up.


I. INTRODUCTION
The recent experiments reveal that neutrinos have mass and they oscillate [1] from one flavor to another flavor as they travel. In the standard framework of oscillation scenario with three active neutrinos, the mass eigenstates |ν i and the flavor eigen states |ν α are related by a mixing matrix U [2]: The flavor mixing matrix in vacuum U (PDG representation [1]) is In this oscillation picture, there are six parameters: two mass squared differences ∆m 2 i j (∆m 2 12 , ∆m 2 32 ), three mixing angles (θ 12 , θ 23 , θ 13 ), and a CP-violating phase δ. At present, information is available on two neutrino mass-squared differences and two mixing angles: from atmospheric neutrinos one gets the best-fit values with 3σ error |∆m 2 32 | ≃ 2.50 +.52 −0.60 × 10 −3 eV 2 , θ 23 ≃ 45.0 •+10.55 • −9.44 • ; while solar neutrinos tell us ∆m 2 21 ≃ 7.9 × 10 −5 eV 2 , θ 12 ≃ 33.21 •+6.02 • −3.88 • [3]. Here ∆m 2 i j = m 2 i − m 2 j . The present experiments are sensitive to ∆m 2 32 only and the both sign of ∆m 2 32 fit the data equally well. The neutrino mass hierarchy, whether normal (m 2 3 > m 2 2 ), or inverted (m 2 2 > m 2 3 ), is of great theoretical interest. For example, the grand unified theory favors normal hierarchy. It can be qualitatively understood from the fact that it relates leptons to quarks and quark mass hierarchies are normal. The inverted mass hierarchy which is unquarklike, would probably indicate a global leptonic symmetry [4].
For nonzero value of θ 13 , the hierarchy can be determined from the matter effect on neutrino oscillation. The contributions of this effect to the effective Lagrangian during the propagation through matter are opposite for neutrino and antineutrino, which depends mainly on the value of θ 13 and the sign of ∆m 2 32 . So the total number of observed events is expected to differ for normal hierarchy (NH) and inverted hierarchy (IH). The hierarchy can be more prominently distinguished if the experiment is able to count neutrinos and antineutrinos separately. This needs charge identification of the produced leptons.
In case of θ 13 = 0, the mass hierarchy can also be observable in principle even in case of vacuum oscillation due to nonzero value of ∆m 2 21 [5]. Presently over the world, there are many ongoing and planned experiments: MINOS [6], T2K [7], ICARUS [8], NOVA [9], D-CHOOZ [10], UNO [11], SKIII [12], OPERA [14], Hyper-K [13], and many others. It is notable that all these experiments are planned in the north hemisphere of the Earth. Among them MINOS is the only magnetized one and has a good charge identification capability. A large magnetized ICAL detector is under strong consideration for the proposed India-based Neutrino Observatory (INO) [15] near the equator. Since it has high charge identification capability ( > ∼ 95%) [16], it is expected to have reasonably better capability to determine the neutrino mass hierarchy.
It is also notable here that the CERN-INO baseline happens to be close to the so-called 'magic' baseline [17,18] for which the oscillation probabilities are relatively insensitive to the yet unconstrained CP phase. This permits such an experiment to make precise measurements of the masses and their mixing avoiding the degeneracy issues [19] which are associated with other baselines.
The mass hierarchy with atmospheric neutrinos at the magnetized ICAL has been studied in [20,21,22]. In [20,21], the number of events are calculated in energy and zenith angle bins. Then the chi-square is calculated including the effect of possible uncertainties over a large range of zenith angle and energy. It was found in [20] that a precision of ≃ 5% in neutrino energy and 5 • in neutrino direction reconstruction are required to distinguish the hierarchy at 2σ level with 200 events for sin 2 2θ 13 = 0.1 and roughly an order of magnitude larger event numbers are required in case of resolutions of 15% for the neutrino energy and 15 • for the neutrino direction. However, in [21], the authors are able to relax the resolution width to 10 − 15% for energy and 10 • for direction reconstruction in distinguishing the type of hierarchy at 95% CL or better for sin 2 2θ 13 ≥ 0.05 with 1 Mton-yr exposure of ICAL.
In [22], an asymmetry parameter is considered as where, u, (ū), d, (d) are the number of up going and down going events for neutrino (anti-neutrino), respectively for a baseline L and an energy E. For down going events, the 'mirror' L is considered, which is exactly equal to that value when neutrino comes from the opposite direction. Here the asymmetry is calculated integrating over E ≥ 4 GeV. Here a statistically significant region with maximum number of events corresponds to the range 500 < ∼ L/E < ∼ 3000. However, the results depends crucially on the L/E resolutions.
In [23], the neutrino and the anti-neutrino events for direct as well as inverted hierarchy are considered for the E range 5 − 10 GeV and the L range 6000 − 9500 km.
In our work, we find the resonance ranges in terms of directly measurable parameters θ zenith µ and E µ for a given set of oscillation parameters, where the matter effect contributes significantly in distinguishing the hierarchy. We have also studied the possible observables and find the suitable ones for discrimination of the hierarchy. Finally, we make a marginalized χ 2 study over a range of oscillation parameters and find how low θ 13 can be explored in discriminating the hierarchy at INO.

II. THE INO DETECTOR
The detector is a rectangular parallelepiped magnetized Iron CALorimeter detector [15]. Here the simulation has been carried out for 50 kTon mass with size 48 m × 16 m × 12 m. It consists of a stack of 140 horizontal layers of 6 cm thick iron slabs interleaved with 2.5 cm gap for the active detector elements. For the sake of illustration, we define a rectangular coordinate frame with origin at the center of the detector, x(y)-axis along the longer (shorter) lateral direction, and zaxis along the vertical direction. A magnetic field of strength ≈ 1 Tesla will be applied along +y-direction. The resistive plate chamber (RPC) appears to be the best option for the active part of the detector. It is a gaseous detector consisting of two parallel electrodes made of 2 mm thick 2 m × 2 m glass/Bakelite plates with graphite paint on the out sides and separated by a gap of 2 mm. When a charge particle passes through this active part, it gives a transient and a very localized electric discharge in the gases. The readout of the RPCs are the 2 cm width Cu strips placed on the external sides of the plates. This type of detector has a very good time (∼ 1 ns) as well as spatial resolution.

III. ATMOSPHERIC NEUTRINO FLUX AND EVENTS
The atmospheric neutrinos are produced from the interactions of the cosmic rays with the Earth's atmosphere. The knowledge of primary spectrum of the cosmic rays has been improved from the observations by BESS [24] and AMS [25]. However large regions of parameter space have been unexplored and they are interpolated or extrapolated from the measured flux. The difficulties and uncertainties in calculation of the neutrino flux depends on the neutrino energy. The low energy fluxes have been known quite well. The cosmic ray fluxes (< 10 GeV) are modulated by solar activity and geomagnetic field through a rigidity (momentum/charge) cutoff. At higher neutrino energy (> 100 GeV), solar activity and rigidity cutoff are irrelevant [26]. There is 10% agreement among the calculations for neutrino energy below 10 GeV because different hadronic interaction models are used in the calculations and because the uncertainty in cosmic ray flux measurement is 5% for cosmic ray energy below 100 GeV [26]. In our simulation, we have used a typical Honda flux calculated in 3-dimensional scheme [26].
The interactions of neutrinos with the detector material are simulated using the Monte Carlo model Nuance (version-3) [27]. Here the charged and neutral current interactions are considered for (quasi-)elastic reactions, resonance processes, coherent and diffractive, and deep inelastic scattering processes.
In fig. 1, we show the variation of total Charged Current (CC) cross section with neutrino energy for both neutrino and anti-neutrino. In the third plot of this figure, we also show how rapidly the neutrino flux changes with energy.

IV. NEUTRINO OSCILLATION THROUGH THE EARTH MATTER
The effective Hamiltonian which describes the time evolution of neutrinos in matter can be expressed in flavor basis as where, A = G √ 2N e 2E is the effective potential of ν e with electrons, U is the flavor mixing matrix in vacuum (eq. 2), G is the Fermi constant, N e is the electron density of the medium, and E is the neutrino energy.
For a time t, the evolution of neutrino states is given by ν(t) = S(t)ν(0) with S(t) = e −iHt for constant matter density. The corresponding effective Hamiltonian for anti-neutrinos is obtained by U → U * and A → −A. To understand the analytical solution one may adopt the so called "one mass scale dominance" (OMSD) frame work: |∆m 2 21 | << |m 2 3 −m 2 1,2 | [28,29]. With this OMSD approximation, the survival probability of ν µ can expressed as The mass squared difference (∆m 2 31 ) m and mixing angle sin 2 2θ m 13 in matter are related to their vacuum values by From eq. 5 and 6 it is seen that a resonance in P m µµ will occur for neutrinos (anti-neutrinos) with NH (IH) when sin 2 2θ m 13 → 1 or, A = ∆m 2 31 cos 2θ 13 .
Then resonance energy can be expressed as Fig. 2 shows how the resonance energy and the average density of the Earth changes with zenith angle.

V. ROLE OF HIERARCHY IN OSCILLATION PROBABILITY AT DIFFERENT
The significant difference in survival probabilities between NH and IH arises due to matter resonance. For a neutrino energy (E ν ), the corresponding baseline (L ν ) or zenith angle (θ zenith ) where the resonance occurs, can be understood from fig. 2 for both cases of ν µ andν µ . There are two zones of energy: 1. low energy zone (E ≈ 4 GeV) for the events through the core of the Earth, and 2. high energy zone (E ≈ 6 − 8 GeV) for the events through the mantle of the Earth.
Here we have studied in detail the difference in survival probabilities at different E ν and cosθ zenith ν using full three flavor oscillation formula with the Earth density profile of PREM model [30]. The dependence on the values of the oscillation parameters are also studied in detail. The observations are the following. The survival probabilities for both NH and IH without matter effect produce a sinusoidal behavior in L for a fixed E. In case of non-zero θ 13 with NH, it is seen that when ν µ with energy ≈ 2 − 6 GeV passes through the core of the Earth (density ≈ 12gm/cc), a large depletion in survival probability P(ν µ ↔ ν µ ) with respect to IH arises due to the increase of P(ν µ ↔ ν e ) (see first row of fig. 3). In case ofν µ , it happens for IH (see second row of fig. 3). The important point is that the atmospheric neutrino flux is sufficiently high at this energy range.
Again, there is also a large depletion in survival probability of neutrino (anti-neutrino) for NH (IH) with respect to IH (NH) at the high energy (≈ 5 − 10 GeV) for −0.75 > ∼ cos θ zenith ν > ∼ −0.40 as seen in fig. 4, and 5.
This pattern remains almost same over the present allowed range of oscillation parameters.
Though all the above effects diminish rapidly as θ 13 goes to zero, there remains a finite difference in the survival probabilities for NH and IH at θ 13 = 0 due to nonzero value of ∆m 2 21 .

VI. POSSIBLE SYSTEMATIC UNCERTAINTIES
The systematic uncertainties may enter into the analysis of data in three steps: i) flux estimation, ii) neutrino interaction, and iii) event reconstruction.
We can divide each of them into two categories: I) overall uncertainties (which are flat with respect to energy and zenith angle), and II) tilt uncertainties (which are function of energy and/or zenith angle).
The uncertainty in flux of category II may arise due to the tilt in its shape with energy and zenith angle. This can be expressed in the following fashion for the case of energy: This arises due to the uncertainty in spectral indices. The uncertainty δ E =5% and E 0 = 2 GeV is considered in analogy with [33,34]. Similarly the flux uncertainty as a function of zenith angle can be expressed as However, this uncertainty is less (≈ 2%) [33, 34] than the energy dependent uncertainty. The experimental systematic uncertainties may come through reconstruction of events. Till now it is not estimated by the simulation of ICAL at INO. However the fact from simulation is that the wrong charge identification possibility is almost zero for the considered energy and zenith angle ranges. Since the detector is symmetric to up and down going events, it is expected that the up/down ratio cancels most of the experimental uncertainties. The quantitative details of the uncertainties of category I (except those arising from the event reconstruction) can be found in [35].

VII. MIGRATION FROM NEUTRINO TO MUON ENERGY AND ZENITH ANGLE
In the CC processes, a lepton of same leptonic family of neutrino is produced with addition of some hadrons. Among the leptons, only muon produces a very clean track in ICAL detector and the hadrons produce a few hits around the vertex of the interaction. Our analysis considers only muons, and no hadrons for simplicity. For a given neutrino energy, the scattering angle and the energy of the muon are related to each other. Again, the neutrino direction is specified by two angles: the polar angle and the azimuthal angle. Therefore, the actual resolution function should be a function of energy and these two angles.
We are only interested in zenith angle which is a complicated function of above two angles. We reduce the above complexity of two angles by constructing a two dimensional resolution function R(x E , x zenith ) for a given neutrino energy with To find the resolution function R(x E , x zenith ), we generate one thousand years un-oscillated atmospheric neutrino data for ICAL. Here, the two angles automatically come in proper way in zenith angle resolution and it should work well in analysis of the atmospheric data. We divide the whole data set into 17 neutrino energy bins (for E =0.8-18 GeV) in logarithmic scale. It is also notable that the zenith angle resolution slightly depends on the value of neutrino zenith angle due to the spherical shape of the Earth. To account for this dependence, we divide the whole zenith angle range (cos θ z = −1 to + 1) into 10 bins for every energy bin. In fig. 6 we show typical resolu-tion functions at low (≈ 1 GeV) and high (≈ 15 GeV) energy for both neutrinos and anti-neutrinos. Here, we see, the zenith angle resolution improves and energy resolution worsens with increase of energy. The zenith angle resolution is significantly better forν µ s than ν µ s. From the distribution of resolution function ( fig. 6), it is clear that the probability is high near zero of x zenith and it rapidly falls as we go far from it. Therefore we use a varying bin size in our analysis for both variables, small near zero and gradually bigger at far values.
How we obtained the number of events in muon energy and zenith angle bins from atmospheric neutrino flux using the above resolution function is described in the following steps.
1. For a given neutrino energy and zenith angle bin, we obtain a distribution of events in a plane with axes from a large number of events generated by Nuance as shown in fig. 6. From this distribution we find the probability of getting a muon at a given energy and zenith angle bin dividing the number of events of this bin by the total number of events generated for the particular neutrino energy and zenith angle. 3. To obtain how many events at a particular muon energy and zenith angle bin may come from the number of events of a given neutrino energy and zenith angle bin, we multiply this previously calculated number of events for the given neutrino energy with the probability at this muon energy and angle bin obtained in step 1.
4. The total contribution at a particular muon energy and zenith angle bin is obtained by adding the contribution from all neutrino energy and zenith angle bins. Finally we obtained the total muon events at different energy and zenith angle bins in the considered ranges.

VIII. CHOICES OF OBSERVABLES AND ITS VARIABLES
The effect of matter on the oscillation probability is a complicated function of E, L and the density of medium (ρ). From the study of survival probability P(ν µ ↔ ν µ ) and P(ν µ ↔ν µ ) in section V, it is clear that the difference between NH and IH is more significant for particular ranges of E and L rather than the range of the combination L/E. In case of the combination L/E, there remains some L(E) values for some values of E(L) where the difference between NH and IH is insignificant. For these values of L and E, it increases the magnitude of statistical errors and eventually kills the significance which comes from the interesting regions of L and E. For this reason, we need to consider only important ranges of the baseline L i for each range of energy E i . We may consider the following observables to determine the hierarchy: In case of observables with the combination of neutrinos and anti-neutrinos, we will encounter the uncertainty of the ratio ν µ /ν µ which is ≃ 1 at low energy and it increases with energy [26,31,32]. Since the matter effect and the cross sections for neutrinos and anti-neutrinos are dependent on E, the uncertainty in the ν µ /ν µ ratio may affect significantly in determination of hierarchy for the combined observables.
Firstly, we choose the observables as a ratio of u and d to cancel all the overall uncertainties (discussed in section VI). Secondly, since the magnetized ICAL detector has the capability to distinguish µ − and µ + , we consider the ratios separately for ν µ s andν µ s: In fig. 2, it clearly shows that there is a sudden fall in resonance energy due to a sudden rise in average density of the Earth. We choose two ranges of energy: one for core (vertical events) and another for mantle (slant events). The difference in survival probability between NH and IH can also be seen in figs. 3, and 4 at these two ranges of zenith angle for two typical resonance energy of 3.7 GeV and 7.5 GeV, respectively. The resonance range of energy for a typical value of the zenith angle is also shown in fig. 5. From the discussion of section III, V and VII, we observed the following : 1. The zenith angle resolution width is large with a peak around zero ( fig. 6).
2. From the distribution of the energy resolution, it is clear that a muon of energy E µ can come from any neutrino having energy E ν > ∼ E µ with almost equal probability for E ν > ∼ 3 GeV where the deep inelastic events dominate. Again, as the flux falls very rapidly with energy (see fig. 1) and E ν > ∼ E µ (this is obvious for E ν > 1 GeV for targets at rest), the lower value of the energy range will play a crucial role in the results.
The bonus point of these two ranges of zenith angle (vertical and slant) is that the vertical/horizontal uncertainty in flux is minimized and becomes negligible. Among the systematic uncertainties discussed in section VI, we remain with the energy dependent flux uncertainty. Here we will not consider any uncertainty coming from reconstruction of the events. This will be taken care of in near future in GEANT-based studies.

IX. χ 2 ANALYSIS
Here we are interested to find the hierarchy discrimination ability of ICAL@INO at low θ 13 using the resonance ranges.
It should be noted here that the resonance width falls very rapidly with θ 13 .
From the earlier works [20,21,22] it is found that as the energy and the angular resolution width increase the discrimination ability of hierarchy becomes worsened (see fig. 10 of [21]). In our actual simulation it is found that the width of these resolutions are large (see fig. 6).
We have tried to tackle the problem in the following way. We do not divide the selected ranges. We calculate the total up-going (u i ) and the total down-going (d i ) events and then their ratio (u i /d i ) for each selected range. Finally we perform the chi-square study with these u i /d i ratios of all ranges (see table I).
For a set of input values of the oscillation parameters, we generate the neutrino events for 1 MTon.year exposure of ICAL detector using the Nuance (see fig. 7). Next we calculate the total u i and the total d i for each considered range of E µ i − cos θ zenith µ i (see table I) and then their ratio u i /d i . We call it the "experimental data". The statistical error is estimated by the following formula We adopt the pull method for chi-square analysis. The theoretical data in chi-square is calculated with the following numerical method.
1. Considering a set of the oscillation parameters and a value of the pull variable for a considered systematic   uncertainty (using eq. 9), we obtain the neutrino flux in 200 bins of E ν and 300 bins of cosθ zenith ν . At low energy region the number of oscillation swings with E (L) for a fixed L (E) are very large (see fig. 3 and 5). Again the difference of the oscillation probabilities between NH and IH changes very rapidly for some ranges of L and E. To obtain the accurate oscillation pattern one needs finer binning of L and E. We have checked the oscillation plots varying the number of bins for E ν and cos θ zenith ν and finally fix 200 bins for E ν and 300 bins for cos θ zenith ν . We have used same number of bins to obtain both experimental and theoretical data. The number of events from neutrino to muon energy and zenith angle bins are obtained using the resolution func-tion and the cross section as discussed earlier.
2. Next we calculate the total u i and the total d i separately for each selected range of E µ i -cos θ zenith µ i (see table I). Then the ratio u i /d i is determined for each range. We call it the "theoretical data". It should be noted that the number of bins has been changed in chi-square from that of the flux. (The reason is described above.) 3. We calculate the standard chi-square from these "theoretical data", the above "experimental data" and the "statistical error" (eq. 13) for all ranges. Next the minimization of this chi-sqaure is done with respect to the pull variable to find its best-fit value keeping all other oscillation parameters fixed. Finally the "theoretical data" is again calculated using the best-fit value of the pull variable.
In the last the χ 2 is calculated from the final "theoretical data" and the above "experimental data". The whole process is done to find the χ 2 µ and χ 2 µ for ν µ s andν µ s, respectively. At the end, we find the total χ 2 = χ 2 µ + χ 2 µ for the considered set of oscillation parameters.
As discussed in section VI, we consider only neutrino flux uncertainty due to its tilt with energy[37].

X. RESONANCE RANGES IN TERMS OF E µ AND θ zenith µ
The optimization of the ranges by analytical calculation is almost impossible for the complexity in the behavior of survival probability, the smearing of resolution, and the rapid fall of flux with increase of energy.
We calculate the χ 2 varying the energy and zenith angle ranges for both NH and IH for a given input of θ 13 and type of hierarchy keeping all other oscillation parameters at their best-fit values. Then, we select the ranges where the χ 2 difference for NH and IH is the largest. This is done for the input values of θ 13 = 4 • , 5 • , 7 • , 10 • and tabulated in table I. Due to the effect of energy resolution function and the nature of flux, the higher end of the energy range can not be precisely determined. We see, the ranges squeeze rapidly with the decrease of θ 13 . It should be noted that the ranges determined here has a significant weight of the atmospheric neutrino flux and it is not independent of flux since it is not from the migration of neutrinos to muons only. However, there is no significant change in chi-square if the ranges are changed within few percent of their boundary values.

XI. MARGINALIZATION AND RESULTS
The neutrino oscillation parameters can not be measured at infinite accuracies. To find at what statistical significance the wrong hierarchy can be disfavored, we calculate the χ 2 for both true and false hierarchy varying the oscillation parameters in the following ranges (unless we specify) : 1. We set other oscillation parameters at their best-fit values and choose CP-violating phase δ = 0. For the data set of each input values of θ 13 , we used the corresponding resonance range of E µ − cosθ zenith µ . The marginalized χ 2 for both true and false hierarchy as a function of θ 13 is shown in fig. 8 for the input of IH and different values of θ 13 (5 • , 7 • and 10 • ). Here the marginalization is done over other two oscillation parameters |∆m 2 32 | and θ 23 . The similar plot is also shown for the input of NH in fig. 9, but with E µ − cos θ zenith ranges as optimized for IH. Finally we show the discrimination sensitivity ∆χ 2 = χ 2 (false)− χ 2 (true) marginalized over all three oscillation parameters |∆m 2 32 |, θ 23 and θ 13 as function of input values of θ 13 in fig. 10.
The changes reflected in χ 2 with the input of θ 13 and the type of hierarchy can be understood from the following facts: 1. Normally the number of events is larger for ν µ thanν µ due to larger ν µ cross section. So the ν µ contribution to the difference in chi-square due to hierarchy is larger thanν µ . Since ν µ s are suppressed in NH, the input with IH (larger statistics) will give a larger difference in chisquare and NH can be more easily disfavored than IH. (See fig. 8 and 9 for the input value of θ 13 = 10 • .) However, once any one type of the hierarchies is disfavored, there remains the other type only. So it is essential to disfavor any one of them.
2. The resolution is better for anti-neutrinos than neutrinos. Again as θ 13 decreases, the resonance ranges squeeze very rapidly. This leads to a competition between resolution and statistics resulting a comparable contribution of ν µ andν µ to chi-square difference for θ 13 < ∼ 7 • (compare fig. 8 and 9 for the input value of θ 13 = 7 • ).
3. Though there is an appreciable difference in survival probabilities between NH and IH below θ 13 value of 5 • , the tilt uncertainty with energy kills the difference in chi-squares.

XII. DISCUSSION
Here we will discuss some points for the estimation of detector efficiency and the errors that may occur in event reconstruction of the experimental data.
1. One can estimate the efficiency of the ICAL detector in the following procedure. The fiducial volume of the detector can be estimated so that an event with vertex inside this volume can be reconstructed. To reconstruct an event one may need the minimum number of hits > ∼ 9. Therefore, one can estimate the efficiency just neglecting 10 layers from top and bottom and 50 cms from Resonance ranges in terms of muon energy and zenith angle θ 13 Core Mantle cos θ zenith lateral sides. For this into consideration, the efficiency may be roughly ≈ 120 (layers) × 4700 cm × 1500 cm 140 (layers) × 4800 cm × 1600 cm ≃ 79%.
2. For the maximum difference between NH and IH, we find the lowest value of E µ ≈ 1 GeV at cos θ zenith µ = −1 and the highest value of E µ ≈ 15 GeV at cos θ zenith µ ≈ −0.3. For these cases, we will get number of hits > ∼ 12, and there will be a very high efficiency of filtration and reconstruction of these events in actual experiment.
3. For the muons below the energy of 15 GeV, the charge identification capability of ICAL@INO is > ∼ 95% with a magnetic field of 1 Tesla[15]. Here, we consider only the muons whose reconstruction efficiency is very high and the energy and angle resolutions are within 3-6% for the considered ranges. Moreover, the wrong charge identification possibility is very negligible [15]. So one may easily expect that the result estimated in this paper will not change appreciably after GEANT [36] simulation.

XIII. CONCLUSION
We have studied the neutrino mass hierarchy at the magnetized ICAL detector at INO with atmospheric neu-trino events generated by the Monte Carlo event generator Nuance(version-3). We have done the analysis in muon energy and zenith angle bins. We have adopted a numerical technique for the migration of number of neutrino events from neutrino energy and zenith angle bins to muon energy and zenith angle bins to find the "theoretical data" in chi-square analysis. The so called "experimental data" are obtained in muon energy and zenith angle bins from the Nuance simulation. The resonance ranges in terms of muon energy and zenith angle are found for different values of θ 13 . Then we find suitable variables and make a marginalized chi-square analysis using the pull method for both true and false hierarchy. The choice of variables and the pull method minimize the possible systematic uncertainties in the analysis. Finally we find the sensitivity as the difference ∆χ 2 = χ 2 (false) − χ 2 (true) as a function of θ 13 as shown in fig. 10. It is found that the neutrino mass hierarchy with atmospheric neutrinos can be probed up to a low value of θ 13 ≈ 5 • at ICAL@INO by judicious selection of events and observables. Here we have also shown that > 95% CL can be achieved for discrimination of hierarchy at θ 13 > ∼ 5 • in a 1 Mton-year exposure of ICAL detector.