Study of neutrino oscillation parameters at the INO-ICAL detector using event-by-event reconstruction

We present the reach of the proposed INO-ICAL in measuring the atmospheric-neutrino-oscillation parameters θ23\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{23}$$\end{document} and Δm322\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta m^2_{32}$$\end{document} using full event-by-event reconstruction for the first time. We also study the fluctuations arising from the low event statistics and their effect on the precision measurements and mass-hierarchy analysis for a 5-year exposure of the 50 kton ICAL detector. We find a mean resolution of Δχ2≈2.9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta \chi ^2 \approx 2.9$$\end{document}, which rules out the wrong mass hierarchy of the neutrinos with a significance of approximately 1.7σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.7\sigma $$\end{document}. These results on mass-hierarchy are similar to those presented in earlier studies that approximated the performance of the ICAL detector.


Introduction
In the Standard Model (SM), neutrinos are massless fermions which interact only via the weak interaction through the exchange of W ± or Z 0 bosons. A series of experiments dedicated to neutrinos [1][2][3][4][5][6][7][8][9] have proved the existence of neutrino flavour oscillations, which implies that neutrinos are massive. The neutrino flavour states produced along with charged leptons are linear superpositions of the mass eigenstates. Due to the difference in phase between the wave packets of each of the mass eigenstates, neutrino oscillations occur [10,11].
The mixing angle θ 13 is found to be non-zero from reactor [6,8,9] and accelerator [15,16] neutrino oscillation experiments. The relatively large value of θ 13 ∼ = 8.6 • has intensified the search for C P violation effects in neutrino oscillations, and also the determination of the sign of Δm 2 32 via matter effects [17,18]. Matter plays an important role in enhancing the effect of sin θ 13 via resonance, which is sensitive to the sign of Δm 2 32 and is different for neutrinos and antineutrinos [19,20]. Determination of the sign of Δm 2 32 would help us understand the correct mass hierarchy (MH) of the neutrinos, i.e., whether the MH is (m 1 < m 2 < m 3 ) normal (NH) or (m 3 < m 1 < m 2 ) inverted (IH) hierarchy. A series of experiments with complementary approaches have been proposed using accelerator, reactor and atmospheric neutrinos to determine the MH [21,22]. The intermediate and long baseline, off-axis accelerator neutrino experiments T2K [5] and NOνA [23,24], search for the appearance of ν e in an intense beam of ν μ , wherein the appearance probability depends on the MH of the neutrino states. Recent results from NOνA [25] have disfavoured the IH by more than 93% C.L. in the lower octant for all values of CP phase δ, and it is in good agreement with the T2K result [26] giving a near maximal value of θ 23 and a weak preference of NH. Liquid scintillator detectors proposed at RENO-50 [8] and JUNO [27] could unravel the MH using reactor neutrinos. Atmospheric neutrino experiments using water or ice Cherenkov detectors, such as Hyper-K [28,29], MEMPHYS [30], ORCA and PINGU [31,32], make use of different cross-sections and different ν andν fluxes to study the MH.
The proposed magnetized Iron Calorimeter (ICAL), to be built at the India-based Neutrino Observatory (INO) [33], will study interactions involving primarily atmospheric muon neutrinos and anti-neutrinos. It will consist of three identical modules, each with dimension 16 m × 16 m × 14.5 m (length × width × height) placed in a line and separated by a small gap of 20 cm. Each module will consist of 151 layers of 5.6 cm thick iron plates interleaved with 4 cm air gap containing the active detector elements, glass Resistive Plate Chambers (RPCs). This huge size of 48 m × 16 m × 14.5 m magnetized detector, with a mass of 50 kton, provides the target nuclei to achieve a statistically significant number of neutrino interactions within a reasonable time frame.
One of the main goals of INO is to study the MH via earth matter effects, and to determine the octant of θ 23 . ICAL is designed to have very good muon detection efficiency of greater than 85% for muons greater than 2 GeV (with zenith angle cos θ z ≥ 0.4), combined with excellent angular resolution. The ICAL will prominently observe the matter effects on the upward going muon neutrinos, where it uses the magnetic field to distinguish between ν μ andν μ events by discriminating the charge on the final state muon. The probability of oscillations are different for ν μ andν μ in the presence of matter, and depends on the MH (sign of Δm 2 32 ). The seen relative difference between normal and inverted hierarchy is 3.9% (−5%) for ν μ (ν μ ) for cos θ z > 0.2, assuming NH to be the true hierarchy (where cos θ z = +1 is the up-going direction). Hence the ICAL could study the MH by observing earth-matter effects independently on ν andν. This paper shows the precision reach of ICAL in the sin 2 θ 23 − |Δm 2 32 | plane for a 5 year run of ICAL. Event-byevent reconstruction and fluctuations arising from low event statistics, using an analysis technique that will be suitable to be employed on the actual data. In the previous studies [33], a lookup table was obtained for the reconstruction efficiencies and resolutions from the studies based on single muons of fixed energy and direction. Later the generated data is folded with the detector efficiencies and smeared by the resolution functions within the lookup table. Hence, the previous methods used parameterizations of the efficiency and resolution that do not reflect the tails of these distributions. As well as those methods used very large sample sizes to negate the effect of low statistics. We also apply a few event selection criteria, as presented in Ref. [34], but within the framework of low event statistics and present its effect on the outcome of this analysis. This paper also compares the results from simulated unfluctuated data, by simulating data sets corresponding to 5-years of data. This paper is organized as follows. In Sect. 2, we outline the methodology describing the event detection in the INO detector and the software framework used to simulate and reconstruct the events in the detector. In Sect. 3, we describe the event generation and discuss the fluctuations in the data. We describe how the events are reconstructed as well as the event selection criteria applied to obtain a sample of events.
We also describe the oscillation analysis including the Earthmatter effects and discuss how data collected by the ICAL are sensitive to the MH. We also describe the χ 2 analysis and the binning scheme used, and also discuss the types of systematic uncertainties used in this analysis.
In Sect. 4, we present the results of our simulated analysis, showing the reach of the ICAL for atmospheric oscillation parameters θ 23 and Δm 2 32 . Event selection reduces the event statistics, hence we also present the results with and without event selection to see its effect. We also discuss the effect of fluctuations on the precision measurements and show the different possible outcomes resulting from the low event statistics. We also discuss the results on the MH of the neutrinos and the effects of fluctuations in determining it. Finally, in Sect. 5 we present the summary of our results and conclusions.

Methodology
The NUANCE [35] neutrino event generator, along with the Honda neutrino flux [36] at the Kamioka site, is used to generate neutrino interactions within the ICAL detector. The proposed ICAL geometry, containing mainly iron and glass components of the detector, is given as input to NUANCE. It generates secondary particles from interactions with these materials, and calculates event rates integrated over the weighted flux and cross sections of all chargedcurrent (CC) and neutral-current (NC) interactions, at each neutrino energy and angle. The output from NUANCE contains vertex and timing information, as well as energy and momentum of all initial and final state particles in each event.
In the ICAL, atmospheric neutrinos will interact with an iron nucleus, undergoing NC and CC interactions. The main CC interactions taking place in the detector are quasi-elastic (QE) and resonance (RS) at low energies and deep-inelastic scattering (DIS) at higher energies. All neutrinos interacting via CC interaction produce an associated lepton. The DIS events produce a number of hadrons along with the lepton, while RS interactions produce at most one hadron. In the QE process, no hadrons are produced and the final-state lepton takes away most of the energy of the incident neutrino.
A C++ code developed by the collaboration using the GEANT4-based [37] simulation toolkit, containing the full ICAL detector geometry, magnetic field map and RPC characteristics, is used to propagate the secondary particles. The output from GEANT4 contains the information on energy loss and momentum of the particles all along its path, which are then digitized and reconstructed to obtain the event observables. The detailed information of which are discussed in Sect. 3.2. This paper is based only on the CC neutrino events, where we consider the muon information alone. The information on the energy and direction of the muons is used to study the sensitivity to atmospheric neutrino oscillation parameters θ 23 and Δm 2 32 at INO-ICAL and resolving the MH. Including hadron information is beyond the scope of this paper.

Analysis procedure
The first step in the procedure is event generation. The generated events are reconstructed in the GEANT4-simulated ICAL and oscillations applied event-by-event after event selection. The oscillated events are binned and used in the χ 2 analysis to determine the oscillation parameters. Each of these procedures are described in detail in the subsections below.

Event generation
In this analysis, NUANCE data for an exposure of 50 kton × 1000 years is generated, out of which sub-samples corresponding to 5 years of data are used as the experimentally simulated sample and the remaining 995 years of data are used to construct probability distribution functions (PDF) that are used in the χ 2 fit. Hence the data are uncorrelated with the PDFs that are used to fit the data. This paper is based only on the CC neutrino events with energies less than 50 GeV which corresponds to 98.6% of the sample. The idealized case, where the NUANCE data is folded with detector efficiencies and smeared by the resolution functions obtained from GEANT-based studies of single muons with fixed direction and energy, has been presented previously [38]. In the earlier analysis, although the data was analysed for an exposure of 5 or 10 years, it was scaled down from the 1000 year sample. Hence, the reconstructed central value was always practically the same as the input value. Here we examine in detail the more realistic case, where the data size and central value are both subject to fluctuations.

Event reconstruction
The generated NUANCE data is simulated in GEANT4based detector environment. The energy deposited due to the energy loss of the charged particle in the RPCs are converted to signals, where they are detected by the mutually orthogonal copper strips (along x and y directions parallel to the global detector coordinates described in Sect. 3.3) on the RPCs. Hence the measured data is digitized to form (x, z) or (y, z) and time t of the signal, referred to as hits. Here the z position is given by the layer number of the RPC. A recursive optimal state estimator -the Kalman Filter [39,40], uses the local geometry and magnetic field information to fit the muon hits, where the muons passing through a minimum of three layers are fit to form the track. The direction and the momentum of the muon is obtained from the best fit values of the track. Also the timing information from the RPCs, with a resolution of approximately 1 ns enables the distinction between upward and downward going particles. More details can be found in Ref. [40].
Hence for the first time we have done this analysis performing event-by-event reconstruction, where each event is simulated through the detector and reconstructed to obtain the observables. Therefore the tails of the resolution functions, which have been approximated by single Gaussians and Vavilov [41] functions in the previous studies [33], are also taken in to account in this analysis.
The μ ± leave one or two hits per layer on average, forming a well-defined track, whereas the hadrons leave several hits per layer forming a shower of hits. Rarely (less than 1% of the time) a pion may also leave a well-defined track in the ICAL and may be misidentified as a muon. In this case the longest track is identified as the muon. The iron plates will be magnetised to produce a field upto 1.5 T and this will be used in the ICAL to probe the charge and momentum of the muon. The direction and the curvature of the muon trajectory, as it propagates through the magnetized detector, gives its charge and momentum, respectively. Figure 1 shows the zenith angle (θ z ) distribution before (true) and after reconstruction. Note that in the current analysis cos θ z = +1 is the up-going direction.
The energy of hadrons is obtained by calibrating the number of hits not associated with the muon track, in the event [42]. The incident neutrino energy (E ν ) can be reconstructed from the energies of the muons and hadrons produced in the detector. The poor energy resolution of hadrons [42] affects the reconstruction of the incident neutrino. Hence for ICAL physics analysis hadron and muon energies are used separately, without losing the good energy and angular resolution of muons [33,38,43,44].

Event selection
The reconstruction of muons is badly affected by the nonuniform magnetic field and dead spaces such as coil slots and support structures. Also the horizontal events which pass through very few layers giving very few hits are reconstructed poorly. Partially contained events, where the μ ± leaves the detector volume, are typically harder to reconstruct, as most of the time they leave a short track within the detector. To remove these badly reconstructed events and obtain a better reconstructed sample of data, we investigated applying selection cuts as used in the previous Monte Carlo (MC) studies [34].
Events with χ 2 /ndf < 10 are used in the analysis, where the χ 2 is the chi-square estimate of the fit for the track obtained from the Kalman filter, and ndf is the number of degrees of freedom. Here ndf = 2N hits − 5, where the Kalman filter fits five parameters to form the track and N hits are the number of hits associated with the track, with each hit having two degrees of freedom as they are either (x, z) or (y, z) in coordinates. The badly reconstructed horizontal events are removed by applying a cut of | cos θ z | ≥ 0.35. Also, to keep a check on the events leaving from top and bottom of the detector, a cut on the z-position of the event vertex is applied. Events with vertices lying below z = 6 m and above z = −6 m are the ones selected from up-going and down-going events respectively.
The entire ICAL detector was divided into three regions, depending on the magnitude of the magnetic field. Considering three modules of size 16 m × 16 m × 14.4 m each and choosing an origin at the centre of the central module, the ICAL will have conventionally 24 m, 8 m and 7.2 m on either side of the origin along x, y and z directions respectively. The region |x| ≤ 20 m, |y| ≤ 4 m, with z unconstrained is defined to be the central region. Here the magnetic field is highest and uniform in magnitude (with ≈ 12% coefficient of variation) despite the fact that the direction of the magnetic field would flip along y in the regions |x| < 4 m, 4 m ≤ |x| < 12 m and 12 m ≤ |x| < 20 m. In contrast, the region |y| > 4 m, termed peripheral region, has maximally varying magnetic field in both magnitude (with ≈ 28% coefficient of variation) and direction. Finally the third region |x| > 20 m and |y| ≤ 4 m, termed the side region, has a magnetic field smaller by ≈ 11% and opposite in direction to the central region. The side region has better uniform magnetic field among the three regions (with less than 5% coefficient of variation).
All the events with the interaction vertices in the central region, and with N hits > 0, are selected as they are either contained within the detector or can form a reasonable length of track to identify the direction and momentum. The rest of the events in the peripheral and side regions are classified into partially (PC) and fully (FC) contained events according to The response of ICAL can be quantified in terms of reconstruction efficiency ( rec ), the relative charge ID (CID) efficiency ( cid ), the energy (σ E μ ) and angular resolutions (σ cos θ z ). The Fig. 2a shows the comparison of these quantities with (WS) and without (WOS) event selection, as a function of true muon energy E true μ . The reconstruction efficiency increases with increase in muon energy, as the energetic muon pass through large number of layers. At higher energies the reconstruction efficiencies become almost a constant. The resolution improves considerably at low and high energies with the selection of events. An overall improvement of 23% and 19% is noticed in energy and zenith angle resolution of muons (averaged over muon energy, zenith and azimuthal angles) after the event selection. The CID efficiency, i.e., the fraction of events identified with correct muon charge among the total reconstructed events, shows ∼ 6% to 10% improvement at all muon energies, after the event selection.
The previous studies without the event-by-event reconstruction obtained an angular resolution less than 1 • for near horizontal events (cos θ z = −0.35) with E μ > 5 GeV [38], whereas it deteriorates to less than 8 • after realistic reconstruction. Similarly the CID efficiency was calculated to be above 95% for all the events, whereas it deteriorates in between 60% to 80% after event-by-event reconstruction. The resolutions and efficiencies in the low energy (E μ < 5 GeV) are relatively worse, whereas the near vertical events show a good agreement with the previous results.
The reconstruction efficiency decreases with event selection as expected, where approximately 40% of reconstructed events are lost with event selection. Figure 2b, shows the reconstructed cos θ z distribution and compares the reconstructed events with (WS) and without selection (WOS) criterion. The dip at cos θ z = 0 results from the difficulty to reconstruct horizontal events. This paper also studies the effect a δ cp is assumed to be zero of the applied event selection on the sensitivity of oscillation parameters in ICAL. Hence, the parameter sensitivity is studied with and without applying any selection criterion.

Applying oscillations
The muon signal in the ICAL will have contributions from the component of the ν e flux (Φ ν e ) that has oscillated to ν μ and the component from ν μ flux (Φ ν μ ) that has survived. As the 5-year pseudo-data would have negligible contribution from Φ ν e compared to Φ ν μ , we have only used the Φ ν μ events in our analysis. Hence, neglecting oscillations from Φ ν e , the total number of events appearing in the detector for an exposure time T is obtained from where N D is the number of targets in the detector and P μμ is the ν μ survival probability. Oscillation probabilities are calculated by numerically evolving the neutrino flavor eigenstates [19] using the equation, where [ν α ] denotes the vector of flavor eigenstates, ν α , α = e, μ, τ , U is the PMNS mixing matrix, and M 2 is the mass squared matrix. Here, A is the diagonal matrix, diag(A, 0, 0), with matter term A given by where the sign is positive for ν and negative forν. Here, G F is the Fermi coupling constant, E is the neutrino energy in GeV, and n e is the electron number density which is related to the matter density ρ in gcm −3 . The density profile of the Earth's matter, given by Preliminary Reference Earth Model (PREM) [45], is used to calculate oscillation probabilities for ν andν. The difference in sign of A for ν andν leads to differing oscillation probabilities, which in turn are sensitive to the sign of Δm 2 32 . The ICAL has an advantage because it can differentiate between ν andν events and observe the matter effects separately.
Oscillations are applied on the 5-year data sample using the accept or reject method. First, the survival probability P μμ is calculated for each ν orν with a given energy and direction. To decide whether an un-oscillated ν μ survives  oscillations to be detected as ν μ , a uniform random number r is generated between 0 and 1. If P μμ > r , the event is accepted to have survived the oscillations. Otherwise it is considered to have oscillated into another flavor and is rejected. Here, we have used the true values of the oscillation parameters, from Ref. [46]; see Table 2. The zenith angle distribution of muons before and after applying oscillations via the accept or reject method is shown in Fig. 3, where the effect of change in sign of Δm 2 32 (MH) is shown forν μ (Fig. 3a, c) and ν μ (Fig. 3b, d) events. It also compares the zenith angle distributions with (WS) and without (WOS) event selection. Note that the oscillation signatures are different inν μ and ν μ events for normal and inverted hierarchies. This difference is solely due to the matter effects, where it depends on the sign of Δm 2 32 as we have assumed no C P violation. (It has been clearly established that CC μ events in the ICAL are insensitive to δ cp [33].) Hence ν μ events are separated fromν μ events while binning, to have maximum sensitivity to the MH.

Binning scheme
During reconstruction, the positively charged particles are tagged with positive momentum and the negatively charged particles are tagged with negative momentum by convention. The reconstructed muons of positive and negative charges are binned separately in Q μ E μ and cos θ z bins after applying oscillations, where Q μ = ±1 for μ + /μ − . The events with negative and positive Q μ E μ indicate those identified as ν μ andν μ events respectively. The atmospheric neutrino flux falls rapidly at higher energies. Hence, wider bins were chosen in those energy regions to ensure adequate statistics. Also, within the frame work of low event statistics, increasing the number of high energy bins is not feasible due to the limited statistics. Table 3 summarises the binning scheme used in the current analysis.
The effect of finer binning is studied previously [38] for energies less than 11 GeV, and is known to marginally improve the precision in both sin 2 θ 23 and |Δm 2 32 |. Also, increasing the range of energies beyond 11 GeV is known to improve the result [47]. Increasing the number of high energy bins can improve the precision. The optimization of bin widths at higher energies will be a part of the future work, and the current analysis will focus on the effects of fluctuations arising from the low event statistics.

The χ 2 analysis and systematics
The pull approach [48] is used in defining the χ 2 such that systematic uncertainties are incorporated. The pull approach is equivalent to the covariance approach, but is computationally much faster. After binning the oscillated events, the 5 year simulated data set is fit by defining the following χ 2 [38,49]: where, Here, N data i j and N pdf i j are the observed and the expected number of muon events respectively, in a given (cos θ i z , E j μ ) bin, while n cos θ z and n E μ are the total number of cos θ z and E μ bins respectively. N data i j is calculated for true values of oscillation parameters, summarised in Table 2, whereas N pdf i j is obtained by combining ν μ andν μ PDFs as in Eq. 5, where T ν i j and Tν i j are the ν andν PDFs respectively normalized from 995 year sample, with f as a free parameter describing the relative fraction ofν μ and ν μ in the sample, with R being the normalization factor in the fit which scales the PDF to 5 years. Theν μ and ν μ PDFs with (WS) and without (WOS) selection criterion are shown in Fig. 4a, b.
The systematic errors and the theoretical uncertainties are parametrized in terms of variables {ξ k } called pulls. The value ξ k = 0 corresponds to the expected value, and the variation, ξ k = ±1 corresponds to a one standard deviation for each source of systematics resulting in an uncertainty of π k i j for the k th source.
In this analysis we have considered two systematic uncertainties, a 5% uncertainty on the zenith angle dependence of the flux and another 5% on the energy dependent tilt error [49], parametrized by ξ flux zenith and ξ flux tilt respectively. There is no systematic uncertainty related to the flux normalization as R and f are the fit parameters which fixes the overall and relative flux normalizations. To calculate the energy tilt error i.e., the possible deviation of the energy dependence of the atmospheric fluxes from the power law, we use the standard procedure as given in, for example Ref. [49], and define: Neglecting the effect of oscillations, the expected number of events Φ 0 (E) is calculated for each (i j) th bin. Then we compute Φ δ (E), where δ is the 1σ tilt error taken to be 5% and E 0 = 2 Gev, and find the relative change in flux to obtain the coupling π tilt i j . The coupling π zenith i j in each bin is calculated in proportion to the zenith angle value of that particular bin. The parameters in the fit are marginalized as given in Table 4, where f and R are always marginalized over the given ranges. The parameters sin 2 θ 13 , sin 2 θ 12 and Δm 2 21 had minimal effect when marginalized hence were kept constant in the fit without any prior constraint.

Parameter determination
The fluctuated pseudo-data set is first fit to determine sin 2 θ 23 , marginalizing over |Δm 2 32 |, for an input value of sin 2 θ 23 = 0.5. The comparison of Δχ 2 with (WS) and without (WOS) event selection is shown as a function of sin 2 θ 23 in Fig. 5a, where the octant degeneracy in sin 2 θ 23 , which stems from the leading term sin 2 2θ 23 in the oscillation probability, is broken due to the relatively large value of θ 13 = 8.5 • . Hence, the asymmetrical curve in sin 2 θ 23 shows the effect of matter oscillations in breaking octant degeneracy. The significance of the fit, i.e., how far the observed value (best fit value) is away from the parameters true value (input value), is defined as where Δχ 2 input and Δχ 2 min are the Δχ 2 values at the true and observed values of the parameter respectively. The fit to sin 2 θ 23 without event selection converges to a value of 0.586 +0.060 −0.093 with a significance of 0.86, i.e., within 1σ of the input value, whereas the fit after event selection converges to 0.676 +0.063 −0.072 within 2σ of the input value. The fit to sin 2 θ 23 with event selection shows relatively larger uncertainty at 2 and 3σ range.
The comparison of Δχ 2 with and without event selection is shown as a function of |Δm 2 32 | in Fig. 5b, where the data is fit to determine |Δm 2 32 |, marginalizing over sin 2 θ 23 , for an input value of Δm 2 32 = 2.32 × 10 −3 eV 2 . The fit without event selection converges to a value of (2.38 +0. 11 −0.39 ) × 10 −3 eV 2 , within 1σ of the input value with a significance of 0.51, whereas the fit after event selection converges to (2.184 +0. 23 −0.37 ) × 10 −3 eV 2 also within 1σ of the input value. The fit to |Δm 2 32 | also shows relatively larger uncertainty at 2 and 3σ range after applying event selection. The multiple local minimas in Δχ 2 function is due to the statistical uncertainty on the PDF, and it is observed to reduce with fits to PDFs constructed from larger MC samples. The free parameter f converges to 0.26±0.01 and 0.27 ± 0.01 before and after applying the selection, where as the parameter R converges to 5901 ± 79 and 3628 ± 62 respectively within 2% of uncertainty. Further adding prior constraints on sin 2 θ 13 and Δm 2 12 , was observed not to make any difference in the fit results or in terms of the coverage in the sin 2 θ 23 − Δm 2 32 plane. The fit with event selection shows larger coverage at 99% CL. Note that the data sets are fluctuated, hence many different sets were also studied and all showed a larger coverage for the fits with event selection.
To decouple the effect of fluctuations from the effect of event selection, a 1000 year sample is scaled to required 5-year sample to minimize the fluctuations. The Fig. 7a, b shows the comparison of Δχ 2 obtained with and without event selection as a function of sin 2 θ 23 and |Δm 2 32 | respectively for the fit without fluctuations. The fits converge near the input value, and shows a relatively worse precision after the event selection as expected. The deterioration of precision in determining sin 2 θ 23 and Δm 2 32 after the event selection is due to larger statistical uncertainty, as the sample size was reduced by 40%. Hence, the rest of the analysis in this paper mainly focus on the fits and the effects of fluctuations without event selection.

Effect of fluctuations
Earlier analyses [33] scaled the 1000 year sample to a size corresponding to 5 years to generate the pseudo-data set. The process of scaling nullifies the effect of fluctuations so that the best-fit is always close to the input value. Then the parameter sensitivities are to be understood as the median value when averaged over a large number of randomly generated samples. In order to see this, we generate an unfluctuated 5 year sample by scaling the 1000 year set and a similar analysis is performed as in Sect. 3. The comparison of Δχ 2 with (WF) and without (WOF) fluctuations is shown as a function of sin 2 θ 23 and |Δm 2 32 | in Fig. 8a Fig. 9 Significance of convergence obtained from sixty different data sets. The solid (blue) and the broken (orange) line represents the significance for the fit to sin 2 θ 23 and Δm 2 32 respectively and the dotted (black) line shows the significance obtained for the simultaneous fit to sin 2 θ 23 and Δm 2 32 in Δm 2 32 . Note that three fits with fluctuations (WF :1, WF : 2 and WF : 3) from three independent fluctuated data sets are used in comparison, and each of them differ in the parameter sensitivities and the best fit values.
The fit to sin 2 θ 23 for the data without fluctuation gives smaller uncertainty in the lower octant, and the fluctuations in the data leads to the fluctuations in the octant sensitivities (see Fig. 8a). Also the uncertainties in the parameter determination changes along with the significance of convergence, with each independent fluctuated pseudo-data set. The correlated precision is also observed to change with each fluctuated data set, where they show different coverages in sin 2 θ 23 − Δm 2 32 plane. The analysis was repeated for sixty different fluctuated data sets, performing separate (one parameter) and simulta-neous (two parameter) fits to determine sin 2 θ 23 and |Δm 2 32 |. Figure 9 shows the significance of convergence in terms of standard deviation σ . Almost 68% of times, the fit to sin 2 θ 23 converges within 1σ of the input value sin 2 θ 23 = 0.5. A similar trend is observed in the fit to Δm 2 32 , where 59% of times it converges within 1σ of the input value Δm 2 32 = 2.32 × 10 −3 eV 2 . Also 95% of the time the fit converges within 2σ , which evidently shows the Gaussian nature of the fit. The simultaneous fit to sin 2 θ 23 and Δm 2 32 also shows a similar behaviour in significance. Figure 10 shows the average coverage area with 99% CL in sin 2 θ 23 − Δm 2 32 plane, obtained by averaging the coverages from simultaneous fit to fifty different pseudo-data sets. The orange band signifies the 1σ uncertainty in calculating the average, where the asymmetrical widths from the best fit point of each data set were used. The precision reach for the fit without fluctuations is within 1σ of the average coverage area calculated.
Previous studies [33,47] have obtained a better precision in sin 2 θ 23 and |Δm 2 32 |, and have quantified the precision on these parameters as where P max and P min are the maximum and minimum values of the concerned parameter determined at the given C.L. Figure 11 compares the Δχ 2 obtained from the previous [47] and the current analysis methods. It is to be noted that we have used the same binning scheme and the same set of NUANCE data in both the analysis methods for comparison, but the previous method incorporates the smearing of resolution functions whereas the current method incorporates event-by-event reconstruction. The precision in sin 2 θ 23 for the previous method is 19.4% at 1σ , whereas for the current  Fig 11(a)). The parameter |Δm 2 32 | also shows a similar behavior, where the precision deteriorates from 5.9% to 12.9% at 1σ for the current method (see Fig 11(b)). The drop in precision for the current method is more predominant in |Δm 2 32 | and is clearly seen with a 30% difference in precision at 3σ .
The noted difference in precision is due to the realistic approach of the event-by-event reconstruction, where the tails of the resolution functions, which were approximated in the previous studies, have been included. In the previous methods, the NUANCE data was folded with the detector efficiencies and were smeared by the resolution functions obtained from GEANT-based studies of single muons with fixed energy and direction [38].

Mass hierarchy determination
The 5 year pseudo-data set is oscillated via accept or reject method assuming NH (IH), which is then fit with true NH (IH) and false IH (NH) PDFs. The parameters in the fit are marginalized as given in Table 4, and the Δχ 2 resolution to differentiate the correct hierarchy from the wrong hierarchy is defined as: where χ 2 true and χ 2 false are the minimum values of χ 2 from the true and false fits respectively. Figure 12a shows the true and false hierarchical fit to sin 2 θ 23 for a particular pseudo-data set with fluctuations, wherein the resolution of Δχ 2 MH = 7.2 rules out the wrong hierarchy with a significance greater than 2σ for this set. The procedure was repeated for sixty independent 5-year fluctuated data sets to see the effect of fluctuations on the mass hierarchy significance. Figure 12b shows the distribution of Δχ 2 MH obtained from the fit to sixty different sets. The mean resolution of Δχ 2 MH = 2.9 rules out the wrong hierarchy with a significance of ≈ 1.7σ for a 5 year run of 50 kton ICAL detector. The large uncertainty in Δχ 2 MH is due to the fluctuations in the data and the negative values signifies the identification of the wrong mass hierarchy. Note that an earlier analysis that excluded the effect of fluctuations [47] gave a value Δχ 2 MH = 2.7 and our mean value is compatible with this, as expected, given minor differences in the analysis procedures. A 13 year exposure of ICAL will give a 3σ separation to differentiate the correct MH. The classical statistical analysis used here is the same as that in the earlier works by the INO collaboration hence allowing a direct comparison. However, using a Bayesian approach [50] could improve the hierarchy sensitivity.

Discussion and summary
One of the main aims of the proposed ICAL at INO is to measure the atmospheric neutrino oscillation parameters sin 2 θ 23 and |Δm 2 32 |, and also to measure the mass hierarchy (MH) of neutrinos. The moderately large value of θ 13 and the ability of the magnetised ICAL to distinguish a neutrino event from an antineutrino event, allows the observation of earth matter effects separately in ν andν and helps identify the MH of neutrinos. In this analysis we focus on the precision measurements and the mass hierarchy resolutions that ICAL could attain within a period of 5 years.
Incorporating a realistic analysis procedure, it is for the first time we have applied event-by-event reconstruction and have considered the tails of angular and energy resolution which were approximated by single Gaussians and Vavilov functions in previous such studies [33]. We show that incorporating non-Gaussian resolutions are likely to affect the parameter sensitivities, whereas reduction in precision by 4% and 7% is obtained for sin 2 θ 23 and Δm 2 32 respectively at 1σ . Also for the first time we study the effect of low event statistics on the precision and MH measurements, by introducing fluctuations in the data. It is also for the first time within the framework of low event statistics, we show the effect of event selection criterion on the parameter sensitivities, and show that we can include all reconstructed muons to get better sensitivity of parameters. Hence within the framework of low event statistics, we show that the fit without any selection criterion (WOS), where we include all the reconstructed events, is the baseline to obtain better constraints on the parameters.
We start by using 5 year fluctuated pseudo-data set for our analysis and apply oscillations via the accept or reject method. The oscillated data is binned in E μ and cos θ z for the χ 2 analysis, where we have used the energy and direc- The constraints on sin 2 θ 23 and Δm 2 32 are compared with and without event selection criterion. Statistically we lose 40% of events after selection, hence we find large uncertainty in parameter determination after applying event selection. We use an ensemble of independent fluctuated data sets to study the effect of low event statistics on the precision measurements of the oscillation parameters. The constraints on sin 2 θ 23 and Δm 2 32 are compared with and without fluctuations, and we find a reasonable agreement between the unfluctuated and the average fluctuated precision reach obtained in sin 2 θ 23 − Δm 2 32 plane. As far as the mass hierarchy of the neutrinos is concerned, we find a mean resolution of Δχ 2 M H = 2.9 from an ensemble of sixty experiments. This rules out the wrong hierarchy with a significance of ≈ 1.7σ , consistent with earlier analysis obtained without considering fluctuations. We also find a significant deviation in the mean value of Δχ 2 M H , and roughly a 15% probability of obtaining the wrong hierarchy due to the fluctuations in the data. In the near future, combining the CP independent measurements at ICAL with the measurements from NOνA, T2K and other reactor experiments, will help to determine the correct MH [22,51].
This paper presents an analysis procedure which can be used on the real ICAL data where the fluctuations are inbuilt as the PDFs are uncorrelated. However, in this analysis we have only used muon information from CC ν μ events and we have ignored the small contribution from ν e to ν μ oscillated events, that is known to slightly dilute the sensitivities [19]. The ICAL can also measure the hadron energy via proper calibration of hits, and including the hadron energy information in CC events is expected to improve the sensitivity of the detector towards the oscillation parameters and also improve the MH significance [44]. Note that there are CC ν e as well as neutral current (NC) events in the detector. Separation of ν μ CC events from the others is quite robust for E μ 1GeV and has been discussed elsewhere [33]. Separation of low energy CC ν μ events from CC ν e and NC events is an ongoing effort of the INO-ICAL collaboration. A combined analysis including all the CC events along with the hadron information will give us the maximum sensitivity the ICAL can attain, and is likely to improve the results presented in this paper.