Measurement of the muon neutrino inclusive charged-current cross section in the energy range of 1-3 GeV with the T2K INGRID detector

We report a measurement of the $\nu_{\mu}$-nucleus inclusive charged current cross section (=$\sigma^{cc}$) on iron using data from exposed to the J-PARC neutrino beam. The detector consists of 14 modules in total, which are spread over a range of off-axis angles from 0$^\circ$ to 1.1$^\circ$. The variation in the neutrino energy spectrum as a function of the off-axis angle, combined with event topology information, is used to calculate this cross section as a function of neutrino energy. The cross section is measured to be $\sigma^{cc}(1.1\text{ GeV}) = 1.10 \pm 0.15$ $(10^{-38}\text{cm}^2/\text{nucleon})$, $\sigma^{cc}(2.0\text{ GeV}) = 2.07 \pm 0.27$ $(10^{-38}\text{cm}^2/\text{nucleon})$, and $\sigma^{cc}(3.3\text{ GeV}) = 2.29 \pm 0.45$ $(10^{-38}\text{cm}^2/\text{nucleon})$, at energies of 1.1, 2.0, and 3.3 GeV, respectively. These results are consistent with the cross section calculated by the neutrino interaction generators currently used by T2K. More importantly, the method described here opens up a new way to determine the energy dependence of neutrino-nucleus cross sections.


I. INTRODUCTION
Many recent long baseline neutrino oscillation experiments use muon-neutrino beams, with neutrino energies ranging from sub-GeV to a few GeV. The observed neutrino-nucleus charged-current (CC) interactions are then used to infer neutrino oscillation parameters. In this energy region, CC quasielastic and CC single-pion production reactions dominate the total cross section, and so understanding these channels is essential for precision measurements of the oscillation parameters. Measurements of exclusive cross sections, however, are complicated by reinteractions of the final-state hadrons as they exit the nucleus, known as finalstate interactions (FSIs). FSIs can absorb or produce particles, resulting in a different set of particles entering the detector than would be expected from the initial interaction. For example, the pion from a CC single-pion interaction might be absorbed in the nucleus, so the observable final state is similar to that from a CC quasielastic event. The CC inclusive channel is much less sensitive to these effects, since it requires only the detection of a charged lepton (muon) from the interaction. A precise measurement of this channel, combined with exclusive measurements, will help improve our understanding of neutrino interactions in this energy region.
So far, the MINOS and T2K experiments have measured the CC inclusive ν μ cross section on iron [1,2] using neutrino beams which cover the few-GeV region. The former experiment measured the CC inclusive cross section for neutrinos with energies above 3.5 GeV using the MINOS near detector. The latter used the T2K near detector, INGRID, to measure a flux-averaged CC inclusive cross section at a mean energy of 1.51 GeV, where the rms spread of the neutrino energy was 0.76 GeV (0.84 GeV) below (above) this mean energy. The MINERνA experiment also measured the CC inclusive cross section on iron, but only the cross-section ratio of iron to CH has been published [3]. The CC inclusive cross section on iron has not yet been measured in the 2-3 GeV energy range, and a measurement covering 1-3 GeV would provide a consistency check between the T2K and MINOS results.
The T2K experiment is a long baseline neutrino oscillation experiment in Japan [4]. T2K utilizes an almost pure ν μ beam produced as the decay product of π þ 's and K þ 's. The beam is first measured by the near detector, ND280, located 280 m downstream from the pion production target. After traveling 295 km, the neutrinos are then observed at the far detector, Super-Kamiokande. Oscillation parameters are determined by comparing the neutrino interactions observed at the near and far detectors.
T2K was built around the "off-axis beam method," where detectors are intentionally placed off the central axis of the neutrino beam (hereafter beam axis). The angle with respect to the beam axis is called the off-axis angle: θ OA . The direction of the neutrino parent particles is distributed around the beam axis, so θ OA is approximately equal to the angle between the parent particle and neutrino directions. In this case, the energy of a neutrino produced from the two-body decay (π → ν μ þ μ) can be expressed as follows: where m π and m μ are the masses of the pion and muon, respectively, while E π and p π are the energy and momentum, respectively, of the pion. The relation between E ν and a Also at J-PARC, Tokai, Japan. b Affiliated member at Kavli IPMU (WPI), the University of Tokyo, Japan.
p π for different θ OA 's is plotted in Fig. 1, showing the maximum neutrino energy reducing as θ OA increases. This indicates that the energy spectrum of the neutrino beam peaks at a lower energy and has a narrower width as θ OA increases. In the T2K experiment, θ OA ¼ 2.5°was chosen so that the neutrino flux peaks around 0.6 GeV, an energy which maximizes the oscillation probability of the muon neutrino at the far detector.
The T2K INGRID detector is installed on the beam axis at the near site. It consists of 14 identical modules, which are spread over a range of θ OA from 0°to 1.1°. Thus, the peak of the neutrino energy spectrum differs among the modules as in Eq. (1).
In this paper, we present a measurement of the ν μ inclusive CC cross section on iron in the energy range of 1-3 GeV with INGRID. This analysis uses data collected from 2010 to 2013, corresponding to 6.27 × 10 20 protons on target (POT). The neutrino interactions at different INGRID modules, which are distributed at different positions and thus observed different beam spectra, is used to extract the energy dependence of the cross section. The topology of each event, which is based on the kinematics of the outgoing muon, is also used to further improve the sensitivity of this measurement to the neutrino energy, since the two are directly related. The different energy spectra and event topologies are combined to construct a probability density function (PDF), which is used to measure the cross section using the least χ 2 method.
The paper is organized as follows. In Sec. II, we describe the T2K near detector INGRID. Section III introduces the Monte Carlo (MC) simulation used to predict neutrino event rates at the INGRID detector and describes the systematic uncertainties associated with this model. Next, the analysis method used to extract the energy dependence of the cross section is explained in Sec. IV with a discussion of the remaining systematics in Sec. V. Finally, the result of the analysis is presented in Sec. VI.

II. THE T2K NEAR DETECTOR: INGRID
Situated 280 m downstream from the pion production target, the INGRID detector monitors the neutrino beam direction and intensity. It consists of 14 identical modules, each of which is composed of nine iron target plates and 11 scintillator tracking planes. The iron plates and the tracking planes are stacked in alternating layers forming a sandwich structure, as shown in Fig. 2.
Each of the iron plates has dimensions of 124 × 124 cm 2 and a thickness of 6.5 cm, providing a total iron mass of 7.1 ton per module. The module is surrounded by scintillator veto planes, which detect charged particles coming from outside of the module. Each tracking plane has two layers of scintillator bars aligned orthogonally to one another, enabling particles to be tracked in all three dimensions as they pass through the plane. The veto planes are also formed from scintillator bars. The bars are coated in TiO 2 reflectors to help contain scintillation light, which is then captured by wavelength shifting (WLS) fibers which run through the center of the bars. This light is then read out by a multipixel photon counter (MPPC) [5,6], and the resultant signal is digitized by the Trip-t front-end board [7], to give the integrated charge and timing information. The integration cycle of the electronics is synchronized with the neutrino beam pulse structure, ensuring all data are captured.
The modules are installed in a cross shape centered on the beam axis. An overview of the INGRID detector is shown on the top in Fig. 3. An ID is assigned to each module as shown on the bottom in Fig 3. Further details of the detector and the basic performance of INGRID can be found in Ref. [8].

III. SIMULATING NEUTRINO INTERACTIONS IN THE INGRID DETECTOR
In this analysis, the cross section is measured by comparing data to a prediction, which is calculated using three sequential MC simulations: (1) prediction of the neutrino flux at the INGRID modules; (2) simulation of the neutrino-nucleus interactions in the iron target ( 56 Fe); (3) propagation of final-state particles through the detector and modeling of its response. After these processes, an event selection, which is detailed in Sec. IV B, is applied to the output in the same way as the data.
A. Neutrino flux

Flux prediction
A detailed description of the neutrino flux prediction can be found in Ref. [9]. In the simulation, protons are impinged upon the carbon target to produce hadrons, which decay into neutrinos. FLUKA2008 [10] and GEANT/GCALOR [11] are used to model hadron production in the target and surrounding material. Propagation of the resulting particles through the electromagnetic horns, which focus the charged hadrons along the beam axis, is simulated using dedicated, GEANT3-based [12] code, which also models the subsequent decay of the particles. For each hadron decay mode which produces neutrinos, the probability of the neutrinos to be emitted in the direction of the INGRID detector is calculated. The flux prediction is obtained by weighting the generated neutrinos with these probabilities. The flux is then tuned using hadron interaction data, primarily from the NA61/SHINE experiment [13]. Other hadron production data (Eichten et al. [14] and Allaby et al. [15]) are also used to tune the simulation in regions of the hadron production phase space that are not covered by the current NA61/SHINE measurement. In this analysis, the NA61/SHINE data taken in 2007 are used to correct the neutrino flux [16,17].
The neutrino flavor content across different energy regions at module 3 (one of the center modules) is summarized in Table I. This shows that muon neutrinos account for > ∼95% of the total flux for E ν < 3 GeV. The muon-neutrino flux fraction then falls to less than 90% for E ν > 3 GeV, whereν μ þ ν e account for > ∼10% of the total flux. Table II shows the ν e to ν μ flux ratio at modules 0 and 3, demonstrating that the difference in ν e contamination between the different modules is very small. Figure 4 shows the obtained muon-neutrino flux spectra at the INGRID modules. The neutrino energy spectrum changes with module position, with the spectrum at module 0 softer than that at module 3. This is because module 0 is located at θ OA ¼ 1.2°from the neutrino beam axis, and so the neutrino flux passing through it is shifted to lower energies due to the off-axis beam effect. This feature is, indeed, essential for this cross-section measurement.

Flux uncertainties
The systematic error on the neutrino flux prediction comes from uncertainties in hadron production and from errors in the measurement of the proton beam, the horn current, and the target alignment. The uncertainties in the  hadron production are mainly driven by uncertainties in the NA61/SHINE measurements and those in measurements by Eichten et al. and Allaby et al. The second category of flux errors is associated with inherent uncertainties and operational variations in the beam line conditions. They include uncertainties in the proton beam position, the beam direction, the absolute horn current, the horn angular alignment, the horn field asymmetry, the target alignment, and the proton beam intensity. The method used to estimate these flux uncertainties is described in Ref. [9]. Figure 5 shows the flux error at module 3, which includes all sources of uncertainty, and demonstrates that the systematic error on the neutrino flux is dominated by the uncertainties in the hadron production model. The propagation of the flux error in this analysis is described in Sec. VA.

B. Neutrino-nucleus interaction simulation
Neutrino-nucleus (Fe) interactions in the INGRID module are simulated by a neutrino event generator, which is a composite of different neutrino interaction models. NEUT (version 5.1.4.2) [18] is used as the primary event generator in this analysis. GENIE (version 2.8.0) [19], a different neutrino interaction simulation package, is also used for comparison. This section first describes the various interaction models simulated in the NEUT event generator and then explains the systematic uncertainties associated with each model. Details of the event generators used in T2K can be found in Ref. [20].

The NEUT neutrino event generator
Given a neutrino energy and a detector geometry, NEUT determines the interaction mode of an event and calculates the kinematics of the final-state particles. It also simulates FSIs for hadrons as they traverse the target nucleus. The following interaction modes are provided for both CC and neutral-current (NC) interactions by NEUT: (i) quasielastic scattering (CCQE or NCQE), (ii) resonant pion production (CC1π or NC1π), (iii) deep inelastic scattering (CCDIS or NCDIS), and (iv) coherent pion production.
Here N and N 0 denote nucleons, l is the lepton, and A is the target nucleus. Table III summarizes the parameters used for modeling neutrino interactions in NEUT. The systematic parameters listed in the table were evaluated in the previous analyses of neutrino oscillation from T2K [21,22] and fall into the following four categories:

Neutrino interaction model uncertainties
(i) M QE A , M RES A , and the nuclear model.-These parameters are used for modeling CCQE and CC1π interactions. A 16.5% error is assigned to the axial vector masses, which comes from a comparison of external measurements of these parameters. The uncertainty on the Fermi momentum and binding energy are estimated from electron scattering data [23]. The uncertainty in the nuclear model is evaluated by exchanging the Relative Fermi Gas (RFG) nuclear model with the spectral function model described in Ref. [24]. (ii) π-less Δ decay and W shape.-In the resonant pion production process, baryon resonances, mainly Δ's, can interact with other nucleons in the target nucleus and disappear without pion emission. The π-less Δ decay parameter is introduced to take into account the uncertainties in this process. The W shape parameter is introduced to modify the shape of the momentum distribution of pions from NC single pion production interactions so that it matches MiniBooNE data [25]. photon, kaon, or eta. According to the MINOS measurement [1], the relative uncertainty on the CC inclusive cross section on iron, which is dominated by CCDIS, is approximately 10% at 4 GeV. Using this as a reference point, the error on the CCDIS and CC resonant cross section is scaled using the following formula: Although the error goes to infinity as E ν approaches 0 GeV, this error does not have a significant contribution to the total cross-section error for the lower neutrino energy region because the interactions have a threshold energy of approximately 0.6 GeV and a small cross section in the 1 GeV energy region. Finally, the 1πE ν shape parameter is a weighting factor as a function of neutrino energy. This is introduced to cover the discrepancy between the MiniBooNE measurement of the CC1π cross section versus E ν and the NEUT prediction using the best fit parameters obtained from the fit to the MiniBooNE data [26]. This discrepancy is as large as 50% at 600 MeV, so an error of 50% is assigned to this weighting factor. In the nominal NEUT, this weighting is not applied ("off" as in Table III). (iv) Pion FSI.-There are three pion FSI processes of interest in the T2K energy range: absorption, charge exchange, and QE scattering. In addition to these interactions, the particle production process, defined as "inelastic scattering," was considered, since it is the dominant process at higher pion energies. Uncertainties on the FSI parameters are estimated using external data sets [27]. Propagation of these uncertainties is described in Sec. V B.

C. Detector simulation
The particle type and kinematic information provided by NEUT is passed to the detector simulation built within a GEANT4 framework [28]. All the detector components are modeled in the simulation. The energy deposited by each particle in the scintillator planes is converted into a number of photoelectrons, taking into account the nonlinear response of the scintillator, light collection efficiency and attenuation of the WLS fiber, and the nonlinearity of the MPPC response. The nonlinear response of the III. Neutrino interaction systematic parameters, nominal values, uncertainties (1σ), and interaction types (CC, NC, or CC þ NC). The first, second, and third groups represent the model parameters, the ad hoc parameters, and the pion FSI parameters, respectively. A 1 or 0 in the nominal value column means that the effect of the systematic parameter is implemented or not implemented by default [21,22].

Parameter
Nominal value Uncertainty (1σ) Interaction type ADCs on the front-end electronics is also modeled based on the results of charge injection tests. Particles generated in the wall upstream of the INGRID detectors are also propagated into the detector simulation and are treated as a background (BG) source.
The hadronic interaction of particles in the detector is simulated by GEANT4 using the FTFP_BERT physics list. In this physics list, the GEANT4 Bertini cascade model is used to simulate nuclear reactions by hadrons with kinetic energies below 5.5 GeV. For particles with kinetic energies above 5 GeV, the list uses the Fritiof model [29,30].

A. Overview
First, neutrino interactions are selected in each INGRID module. The different neutrino energy spectra sampled by the different modules provide a way of extracting the energy-dependent cross section. The selected events are further categorized according to the topology of their finalstate muon in order to improve the sensitivity of the samples to the energy of the incoming neutrino. A PDF is then constructed relating the different INGRID modules and event topologies to the neutrino energy. Finally, to extract the CC inclusive cross section, a χ 2 fit is performed between the selected events and the PDFs.

B. Neutrino event selection
A detailed description of the event selection for neutrino interactions at INGRID can be found in Refs. [2,31]. A brief explanation of each step of the selection is given here. A typical selected muon-neutrino interaction candidate is shown in Fig. 6.
(1) Preselection.-The integrated charge and timing of hits in each channel are recorded with a 2.5 photoelectron threshold. If there are more than three hit channels within a 100 ns time window, these hits are combined to form a hit cluster. The scintillator planes are then searched to find those with at least one hit in both their X-and Y-oriented layers. Such a plane is called an "active plane," and events with two or more active planes are selected.
(2) 2D track reconstruction.-A "cellular automaton" tracking algorithm [32] is applied to hits in the X and Y planes to obtain tracks in the XZ view and YZ view, respectively.
The difference in the most upstream layer hit between any two tracks in the XZ view and YZ view is used to determine if they originate from the same vertex. If this difference is greater than two planes, then the tracks are not matched.
The vertex is defined as the most upstream hit of the track. If there are two or more 3D tracks, a check is performed to see if they originate from the same vertex or not.
The T2K neutrino beam has eight bunches in each beam spill, and each bunch has a width of 58 ns. The selected events are required to lie within 100 ns of the expected time of each bunch.
The event is rejected if the reconstructed track has a hit in the upstream veto plane. The track is then extrapolated in the upstream direction until it intersects with the side veto plane. If there is a hit in the veto plane within 80 mm of this intersection point, the event is rejected.
The fiducial volume (FV) is defined as a cube with a ðAE50Þ × ðAE50Þ cm 2 transverse area, corresponding to the scintillator bars from the 3rd to 22nd channel in the X and Y direction, and from tracking planes 1-8 in the Z direction (see Fig. 7). Events with a vertex inside the FV are selected. After the event selection above, corrections are applied to the MC to account for differences in the individual iron target masses, the observed background rate, and the number of dead channels. We also apply a correction for event pileup, which depends on the beam intensity and results in a loss of efficiency [31]. Table IV summarizes the number of events predicted by the MC simulation and observed in the data.
The selected events contain the muon-neutrino signal as well as background events, such asν μ and ν e interactions. The other backgrounds come from muons, neutrons, and photons generated by neutrino interactions out of the FV or in the pit wall upstream of INGRID (hereafter called beamrelated BG). Since the contamination ofν e is negligible, it is not counted in the MC. The angular distribution of the lepton track for the selected events is shown in Fig. 8. Since the vertex is defined as the most upstream hit of the track, the angular acceptance is limited to between 0°and 90°.
The final selected event sample has > 70% efficiency for CC interactions from neutrinos with energies > 1 GeV, as shown in Fig. 9. Around 5% of the selection inefficiency in the higher energy region is due to events where the muon is produced at a large angle to the Z axis of the detector. This results in it escaping the module before passing through two active planes. Figure 10 shows the detection efficiency for CC events as a function of the true muon track angle, Track angle (degree) 20 40 60 80 Num. of events 6 10 ×    plotted for modules 0 and 3. Module 0 has a lower efficiency than module 3, especially at larger track angles. Module 0 is more off-axis than module 3 so samples a lower energy neutrino spectrum. This results in more low momentum, large angle leptons in module 0 than in module 3. This analysis requires tracks to have at least two active planes and so is less efficient at selecting lower momentum and larger angle tracks. The selection efficiency at module 0 is lower accordingly. The predicted energy spectrum of the reconstructed ν μ events at the INGRID modules is shown in Fig. 11.

C. Event topology
To improve the sensitivity of this analysis to the energy of the neutrino, the selected events are categorized into the following two topologies: (1) downstream (DS-) escaping and (2) non-downstream (non-DS-) escaping. If one of the tracks from the neutrino interaction penetrates the most downstream plane, as shown on the left in Fig. 12, that event is categorized as DS-escaping. All other events, i.e., both side escaping and fully contained events (see the right plot in Fig. 12), are categorized as non-DS-escaping.
Events are then further categorized according to the reconstructed vertex-Z position. Vertex-Z is defined as the most upstream active plane number and ranges from 1 (most upstream) to 8 (most downstream). Events whose vertex-Z is in the most downstream plane are greatly affected by uncertainties in the GEANT4 hadron production model, so only events with vertex-Z in the range 1-7 are used in this analysis.
In total, there are 14 event topologies: (i) DS-escaping: vertex-Z ¼ 1-7 and (ii) non-DS-escaping: vertex-Z ¼ 1-7. Figure 13 shows the energy spectra of "DS-escaping and vertex-Z ¼ 1" events and "non-DS-escaping and vertex-Z ¼ 7" events for module 0. The former have a more energetic μ track and are generally produced by higher energy neutrinos. The latter, on the other hand, tend to have muons produced at a larger angle to the neutrino beam or with a lower energy, and so the majority come from lower energy neutrinos.

D. Module grouping
A shift of the neutrino beam direction changes the peak of the neutrino energy spectra at the INGRID modules. In order to reduce this effect, for the horizontal and vertical directions separately, the two modules at beam-axis symmetric positions are grouped together. This results in seven module groups in total, defined in Table V. Note that the   off-axis angle for the vertical modules increases slightly, because those modules are located in front of the horizontal modules and are therefore closer to the target. The number of selected events for each module group and each topology is then defined as where the indices j and g denote the jth topology and the gth module group (g ¼ 1; 2; …; 7), respectively. The m and m 0 indices stand for the module numbers corresponding to each module group.

E. Detector response uncertainties
This section introduces two different kind of detector response uncertainty: those producing correlation among the event topologies and those that do not. They are summarized as follows: (1) uncorrelated error sources (a) mass of iron plate and (b) pileup correction; (2) correlated error sources (a) event selection and reconstruction and (b) MPPC noise rate. The treatment of these systematic uncertainties in this analysis is described in Sec. V.

Uncorrelated errors
Iron mass.-The error on the measurement of the mass of each iron plate and the machining tolerance for the plate area are taken into account in the systematic error on the iron mass. Since these errors are independent for each iron plate, an uncorrelated error of 0.09% is assigned to the number of selected events. The mass of the iron plates was measured using a crane scale. The typical measurement accuracy of the crane scale is better than 1 kg for a 1 ton load (0.1%). Assuming a fully correlated error resulting from this accuracy. the effect on this analysis is negligible.
Pileup correction.-The number of selected events is corrected to account for event pileup. The correction factor is estimated using data sets at different beam intensities. The uncertainty on the correction factor comes from the statistical error on the number of events, so the uncertainty is uncorrelated between the event topologies. An error of 0.5%-2.0%, depending on the event topology, is assigned in this analysis.

Correlated errors
MPPC noise.-MPPC noise hits sometimes result in misreconstruction of the event vertex or a miscounting of the number of active planes, which produces a variation in the neutrino event selection efficiency. The systematic error caused by the variation in the measured noise rate over time is evaluated by altering the noise rate in the MC. As a result, a 0.1%-1% error is assigned to the number of selected event in each topology.
Event selection.-In this analysis, uncertainties in the following event selection steps are taken into account: (i) 3D track matching, (ii) vertexing, (iii) veto cut, (iv) FV cut. The systematic error on the number of selected events is evaluated by varying the selection threshold for each step, picking the loosest or tightest threshold. The change from the nominal threshold to the loosest (tightest) is defined as the þ1σð−1σÞ change. The resultant fractional variation in the number of selected events (≡ΔN=N) due to the AE1σ change is computed for both the data and MC. Any difference in ΔN=N between the data and MC is then taken as a systematic error.
Uncertainties on the hit efficiency of the tracking planes, the contamination due to beam-related BGs, the hit inefficiency of the upstream veto plane, and the tracking efficiency were found to be negligible for this analysis and are not included in the final result.

F. Cross-section extraction
This analysis uses the least χ 2 method to fit to the observed number of events at each module group (gth bin: 1-7) and for each event topology (jth bin: 1-14): where N obs is the observed number of events, N cc , N nc , and N bg are the expected numbers of CC events, NC, and BG events, respectively, and Δf k and V k are the systematic parameter and the covariance for the kth error source, respectively. For the covariance term uncertainties in the neutrino flux, neutrino interaction model, and the detector response are taken into account and are described in Sec. V. The denominator in the χ 2 statistical term is composed of the statistical error on the observed number of events (N obs jg ), the MC statistical error (σ N mc jg ), and the error on the detector response, which is uncorrelated among event topologies (σ N det jg ): The expected number of CC events in the gth module group and for the jth event topology is expressed as where (i) Δf d is the systematic parameter for the detector response, (ii) Δf cc is the systematic parameter for the CC interaction model, (iii) Δf b is the systematic parameter for the ν μ flux, (iv) ϕ is the ν μ flux, (v) σ cc ccis the ν μ CC cross section, (vi) ϵ cc is the detection efficiency for the CC interaction, and (vii) T is the number of nucleons in the fiducial volume of the INGRID module. Here i goes over the neutrino energy bins, described in Sec. IV G, and Δf i is the parameter being fit, which is used to represent fractional deviations of the CC inclusive cross section. The systematic parameters Δf d j , Δf cc j , and Δf b ig are also fit to include the effect of these systematics into the cross-section result. The INGRID modules are formed from both iron and a plastic scintillator (CH). The effect on this result coming from the different target nucleus for CH interactions is found to be small, so the event rate per unit weight on CH is assumed to be equal to that on iron. Δf d j and Δf cc j are systematic parameters representing the uncertainty on the detector response and the CC cross-section model, respectively, for the jth topology bin. These uncertainties change the detection efficiency as a function of neutrino energy, resulting in a variation in the number of selected events. The difference in these uncertainties between the module groups is very small; therefore, the same parameters are applied to all module groups. Δf b ig parameterizes the flux uncertainty, changing the normalization of the neutrino flux, in the ith energy bin of the gth module. The Δf d j , Δf cc j , and Δf b ig parameters describe fractional deviations from the nominal MC and change the number of events in each event topology, module group, and energy bin.
Since the fraction of NC events in the selected sample is very small, the NC events are summed over the entire energy region and an averaged flux systematic parameter is applied to them: We express the number of the NC events as follows: The other BG events (the beam-related BG andν μ and ν e beam flux components) are summed for each module and for each topology (N bg jg ). The number of BG events in the sample and their associated errors are both small, so the errors on these BGs are neglected in this analysis.

G. Energy binning
The cross section is required to be continuous at the bin boundaries and is linearly interpolated between boundaries, so Eq. (6) is modified as follows: To interpolate, each ith energy bin is divided into fine bins (l ¼ 0; 1; …; L i ). The ith energy bin is defined as the "global bin" and the lth energy bin as the "local bin," respectively. The energy range of the global bins and the binning for the local bins are summarized in Table VI . Taking the average of neighboring parameters produces a measurement at the central energy between the bin boundaries, since, as in Eq. (9), a linear interpolation is applied between the neighboring Δf i parameters. As a result of this averaging, the cross section is measured at 1.1, 2.0, and 3.3 GeV. The final error on the cross section is smaller than those on the Δf i parameters due to the negative correlations between the Δf i 's. This feature is a result of the crosssection continuity requirement.
As seen in Fig. 9, the INGRID detection efficiency for CC interactions falls rapidly for neutrinos with energies less than 0.5 GeV. Since the event samples are not sensitive to the cross section in this region, Δf 0 is not used in the final result. For E ν > 4.0 GeV there is only a small  difference in the neutrino energy spectra between the INGRID modules. Therefore, sensitivity to the cross section for E ν > 4.0 GeV is expected to be worse compared to the lower energy regions. For these reasons it was decided, before fitting the data, to use Δf 1 -Δf 4 to measure the cross section at 1.1, 2.0, and 3.3 GeV. The PDF of CC events in the global energy binning is shown in Fig. 14. Here, the 'fraction" described by the z axis of the figure is obtained for each energy region by dividing the number of CC events in each bin by the total number of CC events in that energy regions. At lower neutrino energies most of the CC events are selected in the downstream vertex-Z bin for the DS-escaping topology, whereas at higher energies the DS-escaping events are distributed uniformly in vertex-Z. Non-DS-escaping CC events are selected in all vertex-Z bins for low energy neutrinos, but higher energy neutrinos tend to be located in upstream vertex-Z bins. In addition, more high-energy neutrino events are selected in modules closer to the beam axis.

V. PROPAGATION OF SYSTEMATIC UNCERTAINTIES
As described in Sec. IV F, the χ 2 has terms with covariance matrices for systematic parameters. In this section, we describe how the covariance matrices for the neutrino flux, the neutrino interaction model, and the detector response, which were introduced in Secs. III A, III B and IV E, respectively, are constructed.

A. Neutrino flux uncertainties
The covariance matrix for each source of error on the neutrino flux prediction, such as the horn current uncertainty, is calculated by taking the variation of the flux due to that error. The total covariance matrix is obtained by summing all these matrices, and Fig. 15 shows the correlation matrix derived from it. The energy binning in the covariance matrix is the same as that used to define the "global bin" in Table VI. One can see that it is largely positively correlated. This correlation comes mainly from the uncertainty associated with hadron production at the target, which varies the neutrino flux in the same way at all INGRID modules.

B. Neutrino interaction model uncertainties
Any systematic error in the CC interaction model would, by definition, alter the CC inclusive cross section itself. This is not the case for NC interaction model uncertainties, and so the systematic errors associated with these two processes are evaluated separately.

Systematics uncertainty on NC interactions
The uncertainty in the number of NC interactions in each bin is expressed by the normalization parameter: where N nc jg is the predicted number of NC events for the jth topology and the gth module group. N 0nc jg is the predicted number of events in the same bin but for the case where one of the NC systematic parameters has been changed by 1σ. The number of events is altered not only by the change in the cross section but also by changes in the event detection efficiency. The same normalization parameter f nc jg is used for all module groups, so the number of predicted events, N nc jg , becomes f nc j is estimated by combining NC events from all module groups. The fractional covariance for the topology bins is then calculated by varying each NC systematic parameter by AE1σ: where N nc iðjÞ is obtained by summing over all module groups: We found that the total NC error is fully correlated between all of the non-DS-escaping bins. Therefore, the seven non-DS-escaping bins are merged into a single bin. Figure 16 shows the correlation matrix for the NC interaction uncertainty, demonstrating that the event topology bins are almost fully correlated with each other. This correlation FIG. 15. Correlation matrix between module groups for the flux error. The energy binning used in this matrix is same as the global bin defined in Table VI. comes mainly from the NC normalization error (see Table III). The total normalization error on the number of NC events is 27%-30%, which is dominated by the uncertainty in the NC other normalization (see Table III), shown in Fig. 17. This gives a maximum error size of 5% on the total (CC þ NC) number of events.

Systematic uncertainty on CC interactions
Varying the CC interaction parameters results in a change in both the neutrino cross section and the detection efficiency. In this analysis, for the systematic uncertainty on CC interactions, only the latter change is taken into account, since the result of the analysis will be the cross section itself. A 1σ variation is applied to a CC interaction parameter, and the new selection efficiency ϵ 0cc ij is calculated. The change in the detection efficiency is then given by the ratio of the new efficiency to the nominal one: where the indices i and j denote the ith energy bin and the jth topology bin, respectively. The predicted number of CC events is then modified using w ij : A fractional covariance between topology bins is then calculated for each CC interaction parameter using the modified (N 0cc ) and nominal number of CC events. This is performed in the same way as for the NC interaction uncertainty, described earlier. The total covariance matrix is computed by summing up the individual matrix from each CC interaction parameter. The obtained correlation matrix and fractional error between the event topology bins are shown in Figs. 18 and 19.

Systematic uncertainty on FSI
For the pion FSI parameters, uncertainties on the absorption, charged exchange, quasielastic, and inelastic scattering cross sections of the pion are taken into account. These systematic errors are treated in a different way than the previous interaction uncertainties, because there are correlations between them. The INGRID data is fitted with N 0cc , obtained by changing each FSI parameter by 1σ, and the difference between the fitted and the nominal result is taken as the systematic error. The effect of each FSI parameter on the measured cross section was found to be negligible, except for the pion absorption uncertainty, which was then added in quadrature to the final result.

C. Detector response uncertainties
As described in Sec. IV E, each error source is categorized according to whether it produces correlation among topology bins or not. For the uncorrelated error sources, the iron mass and pileup correction, the individual systematic errors are summed quadratically and the total inserted into the denominator of the χ 2 statistical term [σ N det jg , in Eq. (5)]. For the correlated error sources, the size of the error does not vary between each module, so a covariance matrix is constructed from the topology bins using the average change over all modules. Namely, the errors are assumed to be fully correlated between module groups. For uncertainties from the event selection, the systematic error is evaluated by varying each selection threshold by 1σ, and the resultant fractional variation in the number of selected events (≡ΔN=N) computed for data and MC. For example, the systematic error on the track matching is estimated by changing the tolerance of ΔZ by AE1 plane (AE1σ change) from its nominal value (ΔZ ¼ 2), where ΔZ is the difference between the most upstream layer hit in the X − Z and Y − Z view. The change by AE1 plane covers the 1σ variation in the position of the upstream edge of the track. The difference in ΔN=N between data and MC is then taken as the systematic error, calculated as where the index j denotes the jth topology bin. If both a þ1σ and −1σ variation are applied to the event selection, then the following covariance is calculated: If only a þ1σ change can be applied, then the covariance becomes The statistical error of Δ i is also calculated and added to Eq. (17) [or Eq. (18)]. Finally, the total covariance is calculated by summing each individual covariance. Figures 20 and 21 show the correlation matrix obtained and the size of the fractional error for each event topology bin. Table VII summarizes the size of the detector systematic error for each error type.

D. Uncertainty in pion multiplicities and secondary interactions
Uncertainties associated with pion multiplicities and pion secondary interactions (SIs) are treated in a different way to the systematics described above. These uncertainties are evaluated by comparing the underlying model with external data. Any observed difference is used to correct the nominal MC sample, and then the χ 2 fit is performed using this "corrected" MC. The differences in the fitted values between the nominal and corrected MC are then taken as the systematic error on the final result.
Pion multiplicity.-In this analysis, the number of events is determined from the number of reconstructed vertices, which are sometimes missed due to the pileup of tracks from multiple neutrino interactions. Events with large numbers of tracks usually contain pions; therefore, the uncertainty associated with the pion multiplicity in these events needs to be considered. This uncertainty is estimated by following the method described in Ref. [33]. The probability of an event having a pion multiplicity of n is expressed as FIG. 20. Correlation matrix from the uncertainties in the detector response. The binning on the y axis is identical to that on the x axis.
where z ¼ n=hni; hni is the mean pion multiplicity and can be expressed by an approximate formula as in Eq. (20). The W is the hadronic invariant mass and is expressed as with E ν μ ðp ν μ Þ and E μ ðp μ Þ denoting the energy (3momentum) of the ν μ and μ, respectively;p N denoting the Fermi momentum of the nucleon and E N the nucleon energy. A, B, and c are derived by fitting two external data sets [34,35] with Eq. (19). These fitted parameters are compared with those used in NEUT, and the differences assigned as the systematic uncertainty. The parameters are used to produce a corrected MC sample which is input to the χ 2 fit. The differences in the fitted values coming from the corrected and the nominal MC are taken as the systematic error on the final result due to pion multiplicity uncertainties. Pion SI.-Hadrons produced in neutrino interactions can also interact while traveling through the detector, a process known as "secondary interaction." In the INGRID simulation the pion SI processes are controlled by GEANT4. In order to evaluate the uncertainty in the pion SI model, the following interaction modes are considered: (i) Quasielastic scattering (QEL) The final-state pion is the same type as the incoming pion.

(ii) Absorption (ABS)
The incident pion is absorbed by the nucleus, resulting in there being no pions in the final state. (iii) Single charge exchange (SCX) A charged pion interacts such that there is only one π 0 and no other pions in the final state. The existing experimental data used to evaluate this uncertainty are summarized in Table VIII. In the table, the reactive cross section is defined as σ total − σ elastic , where σ total is the total cross section and σ elastic is the elastic cross section. As seen in the table, Ashery et al. provide various cross sections across a range of pion momenta and target nuclei, including iron. The other data do not include measurements on iron. For these data, an A-dependent scaling is applied in order to extract the cross section on iron. To evaluate the systematic uncertainty coming from pion SI, we first tune the pion cross section in the momentum region covered by the data in Table VIII. Second, for the lower energy region not covered by data, the ABS π þ ðπ − Þ cross section < 20 (30) MeV is kept constant, motivated by the microscopic calculation from Ref. [36]. The QEL cross section is linearly extrapolated to 0 at 0 MeV. For the higher energy region, each of the cross sections above is tuned on the basis that the size of the reactive cross section is conserved, since the cross sections predicted by GEANT4 are in agreement with the experimental data in this energy region. This study gives four corrected MC samples in total, each of which is then used to fit the T2K data. Finally, the size of the systematic error on the f k parameters due to the uncertainty on pion FSI is calculated as follows: where the index i denotes the ith corrected MC sample and f k is the fitted normalization parameter for the kth energy bin.

VI. RESULT
In this section, we present the result of this ν μ inclusive CC cross-section measurement. Section VI A shows the data set used in this analysis, while Secs. VI B and VI C describe the output from the χ 2 fit and give a summary of this result, respectively.
A. Data set Figure 22 shows the observed and predicted topology distributions in all module groups for the data set used in this analysis, corresponding to 6.27 × 10 20 POT. The number of observed events for the non-DS-escaping topology is 3%-10% smaller than expected. TABLE VIII. Summary of pion-nucleus scattering data used to evaluate the pion SI uncertainty. The reactive cross section is defined as the sum of all the inelastic cross sections.

B. Cross-section fit
The topology distributions after the data fit are shown in Fig. 23. As seen in the figure, the predicted topology distributions with the fitted cross-section normalization parameters applied agree well with the observed distributions. Table IX shows the cross-section normalization param- The fitted values for the flux, detector, CC interaction, and NC interaction uncertainty parameters are shown in Fig. 24. A large deviation from 0 is seen for all the NC interaction uncertainty parameters. As described in Sec. VI A, the prediction overestimates the number of non-DS-escaping events by 3%-10%.
The fitter preferentially reduces the number of NC events to match the predicted topology distribution to the observed one. Since the NC interaction uncertainty parameters are almost fully correlated among topologies, as shown in Fig. 16, all the parameters move toward negative values. There are jumps seen in the CC interaction and detector uncertainty parameters. Both of the jumps appear at the boundary between parameters for non-DS-escaping and DS-escaping events.
In order to derive the normalization factor for the cross section, we take the average of the neighboring fitted f i parameters. The obtained cross-section normalizations are Table X summarizes the uncertainty on the fitted crosssection normalization parameters, broken down by error source. The errors on the combined normalization parameters are summarized in Table XI. The largest systematic error source is the flux uncertainty, which gives a 8%-9% error on the cross-section normalization. The cross-section normalization at 2.0 GeV is less affected by most of the systematic uncertainties than the other normalizations, as shown in Table XI. The reason for this is as follows. In Fig. 14, one can see that the PDF is well differentiated for both module group and event topology at  which results in weak correlation with the other crosssection normalizations. Thus, sensitivity to the crosssection normalization is good for E ν ¼ 2.0 GeV (¼ f 2 þf 3 2 ) but limited for E ν ¼ 3.

C. Summary
In the previous section, five individual fitting parameters were extracted with the least χ 2 method and used to calculate the following measured cross sections at energies of 1.1, 2.0, and 3.3 GeV:  Figure 27 shows these results compared to other measurements [1,2] and the neutrino event generators, NEUT (version 5.1.4.2) and GENIE (version 2.8.0). These measurements are consistent with the energy-dependent cross section measured by the MINOS near detector and the previous, fluxaveraged, cross section measured by INGRID, which used a subset of the data included in this analysis.
This analysis utilizes the different off-axis angle technique together with the final-state kinematics of the outgoing lepton to enhance sensitivity to the neutrino energy. Using final-state lepton kinematics in this analysis makes the result sensitive to uncertainties in the neutrino interaction model, which increases the uncertainty on the final measurement. In order to see how the result depends on the simulation, fits were performed to subsamples (DSescaping only or non-DS-escaping only) of the data. If the final-state lepton kinematics are very different between the data and the simulation, then the fit result should change depending on the choice of subsample. The results were found to be in agreement with each other, showing that the kinematic differences in the regions INGRID is sensitive to are small enough that the analysis is unaffected by them.
The errors achieved in this analysis are not small enough to distinguish between the neutrino models used in the different event generators. Nevertheless, the result seems to prefer NEUT for E ν ≤ 2 GeV and agrees with the MINOS data point and GENIE at E ν ¼ 3.3 GeV. Further reduction of the systematic error size would help in differentiating the neutrino models at the higher energy transition. This reduction could be made by utilizing neutrino beams covering a wider range of off-axis angles, which would provide a modelindependent way to infer the neutrino energy. This analysis demonstrates the feasibility of using different off-axis samples from the same neutrino beam to measure the energy dependence of neutrino interactions, which will provide useful information for future neutrino oscillation analyses.

VII. CONCLUSION
In this paper, we have reported the measurement of the energy-dependent inclusive ν μ charged-current cross section on iron using the T2K INGRID detector and the T2K neutrino beam. The unique variation in the neutrino flux across the INGRID modules, along with event kinematic information, was used to produce event samples FIG. 25. Error (left) and correlation (right) matrices for the five fitted parameters. In both of the matrices, the binning on the y axis is identical to that on the x axis.
FIG. 26. Error (left) and correlation (right) matrices for the cross-section normalization at 1.1, 2.0, and 3.3 GeV. In both of the matrices, the binning on the y axis is identical to that on the x axis. sensitive to neutrinos with energies from 1 to 3 GeV. These were used to extract the inclusive CC muon-neutrino interaction cross section on iron at 1.1, 2.0, and 3.3 GeV, using data corresponding to 6.27 × 10 20 POT. This result is consistent with the predictions of the NEUT and GENIE neutrino interaction generators. This is the first measurement to use the off-axis effect to measure neutrino cross sections as a function of energy and demonstrates the potential of this technique to provide useful information for future neutrino experiments.