First measurement of $\overline{\nu}_{\mu}$ and $\nu_{\mu}$ charged-current inclusive interactions on water using a nuclear emulsion detector

This paper reports the track multiplicity and kinematics of muons, charged pions, and protons from charged-current inclusive $\overline{\nu}_{\mu}$ and $\nu_{\mu}$ interactions on a water target, measured using a nuclear emulsion detector in the NINJA experiment. A 3-kg water target was exposed to the T2K antineutrino-enhanced beam with a mean energy of 1.3 GeV. Owing to the high-granularity of the nuclear emulsion, protons with momenta down to 200 MeV/$c$ from the neutrino-water interactions were detected. We find good agreement between the observed data and model predictions for all kinematic distributions other than the number of charged pions. These results demonstrate the capability of measurements with nuclear emulsion to improve neutrino interaction models.

This paper reports the track multiplicity and kinematics of muons, charged pions, and protons from charged-current inclusive νµ and νµ interactions on a water target, measured using a nuclear emulsion detector in the NINJA experiment. A 3 -kg water target was exposed to the T2K antineutrino-enhanced beam with a mean energy of 1.3 GeV. Owing to the high-granularity of the nuclear emulsion, protons with momenta down to 200 MeV/c from the neutrino-water interactions were detected. We find good agreement between the observed data and model predictions for all kinematic distributions other than the number of charged pions. These results demonstrate the capability of measurements with nuclear emulsion to improve neutrino interaction models.

I. INTRODUCTION
In accelerator-based long-baseline neutrino oscillation experiments, neutrino interactions with the nuclei are essential processes for measuring the neutrino oscillation parameters and searching for CP violation in the lepton sector [1][2][3][4][5][6]. However, an overall explanation of the neutrino interactions has not been provided thus far [7,8]. The charged-current quasi-elastic (CCQE) interactions, which excite one-particle-one-hole states, constitute the dominant interaction process in the energy region of the T2K neutrino oscillation experiment [9]. The CCQE interaction has one lepton and one nucleon in the final state. In addition, there are interactions with two-particle-two-hole (2p2h) excitations [10,11]. The charged-current 2p2h interaction has one lepton and two nucleons in the final state. The T2K far detector, Super-Kamiokande (SK) [12], is insensitive to most neutrons and protons. Events with a single lepton and no other visible particles are selected as the signals, and the incoming neutrino energies are reconstructed from only the outgoing leptons assuming the two-body kinematics of the CCQE interaction. Therefore, the 2p2h interactions involved in the selected events bias the reconstructed neu- * E-mail: hiramoto@scphys.kyoto-u.ac.jp † now at KEK trino energy. In T2K, neutrino interactions are measured and studied using the near detectors [13][14][15][16][17][18][19][20]. However, at present, the measurement of the 2p2h interaction is poor because the momentum threshold for protons is not sufficiently low to detect all the protons from the neutrino interactions. In addition to the proton measurements, precise measurements of interactions including low-momentum charged pions in the final state are important, as they also contaminate the signals at SK when the pions fall short of the Cherenkov threshold in water. Measurements of protons and pions from neutrino interactions with low momentum thresholds play an important role in constructing reliable models of the neutrino-nucleus interactions and reducing the systematic uncertainties in T2K.
Low-momentum hadrons produced by neutrino interactions, especially protons with momenta down to 200 MeV/c, have been measured using bubble chambers containing hydrogen or deuterium [21][22][23] as well as liquid argon time projection chambers [24]. By contrast, recent long-baseline experiments use carbon and oxygen as their targets. The proton momentum thresholds achieved for these nuclei are down to only around 400 MeV/c [17,25]. Hence, a new experiment using a nuclear emulsion detector was proposed to measure protons from neutrino-water interactions with a momentum threshold as low as 200 MeV/c. A nuclear emulsion detector is a high-granularity three-dimensional tracking de-arXiv:2008.03895v1 [hep-ex] 10 Aug 2020 vice. Emulsion detectors have contributed to advances in fundamental particle physics such as the discovery of the charm particles in cosmic rays [26], the direct observation of ν τ [27], and the discovery of ν τ appearance in neutrino oscillation [28]. The detection of extremely short tracks was key to these observations. The high granularity allows clear observation of short-range tracks from neutrino interaction vertices. The charged track multiplicity is determined by preparing an alternating structure of emulsion films and thin water-target layers.
A series of pilot experiments has been carried out by the NINJA collaboration beginning in 2014 [29,30]. This paper reports the results of a pilot run with a small-mass water target (J-PARC T68). A 3 -kg water target was exposed to the T2K antineutrino mode beam from 2017 to 2018. The signals are the charged-current (CC) inclusive ν µ and ν µ interactions on water, and muons, charged pions, and protons are detected as the outgoing particles. We measure the distributions of multiplicity, angle, and momentum of the outgoing particles. In particular, we focused on the measurement of protons in the 200-400 MeV/c range from neutrino-water interactions.
The remainder of this paper is organized as follows. Section II describes the experimental apparatus. Section III discusses the Monte Carlo (MC) simulation. Section IV describes the event reconstruction. Section V addresses the momentum reconstruction and particle identification (PID). Section VI describes the selection of neutrino events. Section VII discusses the estimation of the systematic uncertainties. Section VIII presents the results. Finally, Section IX concludes the paper.

II. DETECTOR CONFIGURATION AND DATA SAMPLES
Three detectors were installed in the T2K neutrino near-detector hall. Figure 1 shows a schematic view of the detectors. The main detector that records all charged particles from neutrino interactions is the watertarget emulsion cloud chamber (ECC). The ECC was installed upstream of one of the modules of INGRID, which is a T2K near detector [31]. In this measurement, INGRID is used to detect muons from the neutrino interactions in the ECC. The emulsion accumulates all the tracks after production without timing information, whereas INGRID records the tracks with timing information. The angular and position resolutions of IN-GRID are not sufficient to identify corresponding tracks between the ECC and INGRID. Therefore, a scintillating fiber tracker (SFT) was newly developed and installed between them. The ECC and SFT were placed in a cooling shelter to maintain the temperature at around 10 • C and the humidity below 60%. This is done to prevent the emulsion tracks from fading under high temperature and humidity as well as to prevent the films from warping due to fluctuations in the ambient temperature.

A. J-PARC neutrino beamline
The J-PARC accelerator provides a high-intensity 30 -GeV proton beam. The proton beam spill is delivered to a graphite target every 2.48 s. The spill has an eight-bunch structure, and each bunch has a full width of around 58 ns and separation of about 580 ns. Hadrons produced by the impinging protons are focused into a decay volume by three electromagnetic horns, where they decay mainly into muons and neutrinos. By changing the polarity of the horns, the charge of the focused hadrons and thus the production of either a neutrino or an antineutrino beam are selected. This measurement is performed with the antineutrino mode beam created by the decay of negatively charged hadrons, predominantly π − . For further details of the neutrino flux prediction, see Ref. [32] B. INGRID INGRID is a T2K on-axis near detector located 280 m downstream of the graphite target. It has 14 modules placed along the vertical and horizontal axes to measure the neutrino event rate and beam profile. We use one horizontal module next to the central module (Fig. 2 top) as a muon range detector. An INGRID module has a sandwich structure consisting of 9 iron plates and 11 scintillator tracking planes (Fig. 2 bottom). The thickness of each iron plate is 6.5 cm. INGRID measures the muon momentum up to around 1 GeV/c with a resolution of around 10%. Each scintillator tracking plane consists of 24 plastic scintillator bars aligned horizontally and 24 vertically. Each scintillator bar has dimensions of 120 cm × 5 cm × 1 cm, and photons are collected by a wavelength shifting (WLS) fiber inserted in a hole made along the longitudinal direction of the scintillator. A silicon photomultiplier (SiPM) is attached to one end of the WLS fiber with an optical connector. The angular and position resolutions of the reconstructed tracks are around 2.7 cm and 3.8 • , respectively [31].

C. ECC
The ECC is an emulsion-based detector composed of alternating layers of emulsion films and target materials. The target materials and their thickness can be selected flexibly. In addition, the alternating structure of emulsion films and thin target layers enables us to achieve a low momentum threshold. Figure 3 shows the structure of the ECC. The components of the ECC are placed in a desiccator. The desiccator is constructed with 2 -cm-thick walls, and it has inner dimensions of 21 cm × 21 cm × 21 cm. The structure formed by two emulsion films and a 500 -µm-thick iron plate vacuumpacked in a 115 -µm-thick aluminized packing film is referred to as a tracking layer. The iron plate is sandwiched between the two emulsion films, each of which consists of a plastic base film that has been coated with an emulsion gel on both sides. These iron plates are employed as supporting structures for the emulsion films and also used for the momentum measurement described in Sec- tion V. The tracking layers are placed at 2 -mm intervals using acrylic frames with a hollow square shape. The desiccator is filled with water, and 2 -mm water layers are formed inside the acrylic frames. Thus, charged particles from neutrino interactions occurring in the water layers make tracks on the upstream or downstream emulsion films. As the tracks are required to pass through at least one iron plate and two emulsion films, the momentum threshold for proton tracks is around 200 MeV/c, while that for pion tracks is around 50 MeV/c. An iron ECC is placed downstream of the water region to measure the momentum of the charged particles using multiple Coulomb scattering (MCS) at the iron plates. In addition, two special sheets (SSs) and one changeable sheet (CS) are installed in the most downstream region. SS1 is placed outside the desiccator, while SS2 is placed inside it. Each SS has four emulsion films with a 2 -mm-thick acrylic plate inserted between the emulsion films. This structure enables us to achieve a good angular resolution. The CS contains two emulsion films. They are replaced every month to separate the tracks into several time periods.

D. SFT
Although the emulsion detector has excellent angular and position resolutions, it does not provide any time information. For track matching between the ECC and IN-GRID, another device with time and position resolutions is required because the angular and position resolutions of INGRID are not sufficient to select a track candidate in the ECC to be connected to an INGRID track. In some cases, an emulsion shifter [30,33] is used to apply timestamps to tracks in the ECC. However, in this pi- lot run, the SFT is employed as a timestamper because it can provide more precise time information than the emulsion shifter. By arranging square fibers in a slanting lattice pattern as shown in Fig. 4, the ratio of the light yields at neighboring fibers can be used to obtain a precise track position. As the light yield at each fiber is proportional to the path length of a charged particle, the ratio of the light yields changes with the position of the particle. The track position d is calculated as where R is the fiber interval and N 1 and N 2 are the light yields from each fiber. The expected position resolution is proportional to 1/ √ N 1 + N 2 when the ratio of the light yields is used. Thus, with the same number of fibers a position resolution better than the typical A/ √ 12, where A is the fiber cross section, can be obtained. Although the position resolution is degraded as the injection angle of the particle increases, this effect is not significant for most muons from the neutrino interactions in the ECC. In this pilot run, 1 -mm square fibers (Kuraray, SCSF-78) are aligned at 0.725 -mm intervals to cover an area of 37 cm × 37 cm. A horizontal layer and a vertical layer are constructed, and each layer consists of 512 fibers.
Hamamatsu S13361-3050AE-04 16-channel SiPM arrays are used for the readout of the scintillation light, and NIM EASIROC modules [34] are used as the readout electronics. The total light yield in a layer is around 60 photoelectrons (p.e.). In this pilot run, the SFT recorded only one event per spill without timing information inside the spill. Therefore, only the first hit is recorded even if there are several hits in a spill.
To reduce the number of readout channels, one SiPM channel reads out four fiber signals, and the signals are read out from both ends of the fibers. As the combinations of the four fibers at the two ends are different, the hit fibers can be identified. Therefore, the total number of readout channels is 512 and the total number of fibers is 1024.

E. Data samples
There are two periods of beam exposure, which correspond to T2K Run 9. The first period (Run-a) is from October to December 2017 and the second period (Run-b) is from March to May 2018. Both Run-a and Run-b are separated into three periods using different CS films, and each period corresponds to roughly one month. Considering periods in which both the SFT and INGRID are collecting data, this analysis is performed with 7.1 × 10 20 protons on target (POT) of the antineutrino mode beam.

III. MONTE CARLO SIMULATION
The expected signals and backgrounds are generated by MC simulations. Three software packages are used: JNUBEAM [32] for the neutrino flux simulation, NEUT [35] for the neutrino-nucleus interactions, and a GEANT4 [36]-based framework for the detector response simulation. In this analysis, ν µ and ν µ interactions on H 2 O and Fe in the antineutrino mode beam are generated by JNUBEAM and NEUT. As the ν e and ν e components of the flux are less than 1%, ν e and ν e interactions in the ECC are not simulated. The MC predictions are normalized by POT and corrected by the detector efficiencies estimated using the data and the MC simulations.

A. Neutrino flux
JNUBEAM is a GEANT3 [37]-based neutrino flux simulator developed by T2K. Interactions of the primary protons from the accelerator and the graphite target are simulated by FLUKA [38,39]. Secondary particles produced are transferred to JNUBEAM, which simulates the propagation, interaction, and decay of these secondary particles. Neutrinos are generated by the decay of the hadrons. The hadron interactions are tuned by external measurements of hadron production such as CERN NA61/SHINE [40][41][42]. Figure 5 shows the predicted flux of the antineutrino mode beam at the location of the NINJA detector. The mean energy of the ν µ components is 1.3 GeV and that of the ν µ components is 2.0 GeV.

B. Neutrino interaction
Using the neutrino flux calculated by JNUBEAM, ν µ and ν µ interactions on H 2 O and Fe targets are gen-  Semi-classical intra-nuclear cascade model [35] erated by NEUT. In addition, neutrino interactions in the upstream wall and INGRID are generated as background sources. Table I summarizes the neutrino interaction models used in this analysis. The nominal MC predictions are generated using NEUT version 5.4.0, which uses the 1p1h model by Nieves et al. with correction by random-phase approximation (RPA) [43,44], and the axial mass M QE A is set to 1.05 GeV/c 2 for the CCQE interactions. The local Fermi gas model (LFG) is used as the nuclear model, while the Spectral Function (SF) [45,46] is prepared as an alternative model. For the 2p2h interactions, the model of Nieves et al. [10] is used. The single pion production is modeled by the Rein-Sehgal model [47], and the axial mass M RES A is set to 0.95 GeV/c 2 . The Berger-Sehgal model [48] is used for the coherent pion production, and the deep inelastic scattering (DIS) is described by the parton distribution function GRV98 and the cross-section model modified by Bodek and Yang [49]. The final state interactions (FSI) in the nuclear medium are simulated using a semi-classical intra-nuclear cascade model [35]. Samples with other parameters are studied for comparison and estimation of the systematic uncertainties as discussed in Section VII.

C. Detector response
The behavior of the particles from neutrino interactions is simulated by a GEANT4-based detector MC framework. The detectors and the wall of the detector hall are modeled. QGSP BERT is used as the default physics list, and muons, charged pions, and protons from the neutrino events generated by NEUT and their secondary particles are simulated. In addition to the neutrino interactions in the ECC, interactions in the INGRID modules and the upstream wall of the detector hall are generated for the background study. The background from cosmic rays is evaluated using the off-beam timing track data instead of the MC simulation.

IV. RECONSTRUCTION
This section describes the track reconstructions in the ECC and INGRID, the hit position reconstruction at the SFT, and the track matching between all the detectors.

A. Track reconstruction in ECC
After the beam exposure, all the emulsion films are developed and several steps of film treatment are performed. Then, they are scanned using Hyper Track Selector (HTS) [50] and the tracks are reconstructed automatically for each film [51]. The current scanning angle is limited to |tanθ| 1.5, where θ is the angle of a track with respect to the direction perpendicular to the emulsion films. Tracks satisfying |tanθ| < 1.3 are used in the analysis. The track density is O(10 3 ) per cm 2 and the detection efficiency of a single emulsion film is 98%-99%. The main components of the tracks in the emulsion films are cosmic rays and environmental radiation. Following the track reconstruction in each film, connections between the films are established by the autoreconstruction process [51]. The track connection process is applied not only to adjacent films but also those separated by one or two other films. The angular and position tolerances are defined as functions of the track angle, and they are determined on the basis of the scattering angle of the minimum ionizing particles (MIPs). The connection efficiency between two films for the MIPs is more than 99.8%. Therefore, by connecting the tracks between both adjacent films and films separated by one or two other films, the connection efficiency becomes greater than 99.99%.

B. Track reconstruction in INGRID
Channels with more than 2.5 p.e. are counted as hits. At least three continuous planes are required to have hits on both horizontal and vertical layers. The tracks are reconstructed using a cellular automaton algorithm [52], which is the same as that used in the event rate and the profile measurements of the T2K neutrino beam. In our analysis, the tracks are required to start at the most upstream plane of INGRID.

C. Hit reconstruction of SFT
As described in Section II D, the SFT fiber hits are identified on the basis of the combination of channels at both ends of the fibers. The hit threshold of the SFT is set at 2.5 p.e. and at least one hit is required in each layer. The hit position is reconstructed from the ratio of the light yields of neighboring fibers. If there is only one hit fiber, the particle is considered to have passed through exactly the center of the fiber because it is likely that the particle penetrated areas that are insensitive due to fiber cladding. To evaluate the track reconstruction efficiency, the effects of accidental noise hits as well as those of missing the second or later hits in a spill by the SFT are calculated using sand muon events which are from neutrino interactions on the upstream wall of the near-detector hall.

D. Track matching
After reconstructing the tracks and hit positions at each detector, a track matching process is performed. To connect tracks between the ECC and INGRID, matching between INGRID and the SFT is carried out first. Then, using the SFT hit position and INGRID angle, matching between the SFT hits and the ECC tracks is performed.
The track matching between the SFT and INGRID is performed using the position and timing information recorded at each detector. The INGRID tracks are extrapolated to the SFT position. If the extrapolated position is within ±10 cm from the reconstructed SFT hit position in the same spill, they are regarded as belonging to the same track. If there are several INGRID track candidates for one SFT hit, the INGRID track in the earliest bunch is selected. This is because the SFT records only the first hit in a spill owing to the limitation of the data acquisition system. By contrast, when one INGRID track has several SFT hit candidates, all of them are put forward to the neutrino event selection. The matching efficiency depends on the track angle, but it is higher than 95% in most regions. The efficiency of the SFT hit reconstruction is included in this matching efficiency.
After matching the INGRID tracks and SFT hits, track matching between the SFT and the ECC is carried out. The tracks recorded on the SS emulsion films are extrapolated to the SFT position. Hits are required to be recorded on at least one of the two films on both sides of SS1. In addition, hits are also required to exist on both CS films. To extrapolate tracks from the SS films, the angle is reconstructed not by a track angle on one film, but by two recorded tracks on the films over the 2 -mm-thick acrylic plate as they give a better angular resolution of around 1 mrad. The angular resolution of a track reconstructed in one emulsion film is typically 2 mrad. If the difference between the SFT hit position and the position of the extrapolated SS track is less than 600 µm, and the difference of their angles is less than 0.2 in terms of tan θ in the horizontal and vertical directions, the track is regarded as a matched track. If there are several candidates, all candidates remain until the neutrino event selection. The angular and position resolutions after the matching are 330 µm and 0.05 in terms of tan θ, respectively. Figure 6 shows the total muon detection efficiencies in Run-a and Run-b. The efficiency in Run-b is lower than that in Run-a because the CS films were slightly bent in Run-b. Thus, the distance between the SS and CS films varied depending on the position in the films, and the accuracy of the matching between these films was degraded.

V. MOMENTUM RECONSTRUCTION AND PARTICLE IDENTIFICATION
The INGRID matched tracks are considered to be muons, while the other tracks recorded in the ECC are considered to be protons or charged pions. This section describes the momentum reconstruction of each particle and separation of protons and pions. In emulsion detectors, the momentum of a charged particle is measured using MCS as P β, where P is the momentum and β is the velocity of the particle. As the angular and position resolutions of the emulsion detectors are sufficient to measure the MCS of the particles, the momentum of the particles can be measured without using a magnetic field. There are two methods for measuring P β: coordinate method [53] and angular method [54]. The coordinate method uses the positional displacement of a track on three films, while the angular method uses the scattering measured by the angular difference of a track between films. In this analysis, the coordinate method is used for the reconstruction of muon momenta, while the momenta of protons and pions are obtained by the angular method, as described later. Moreover, the track range is used to measure the momenta of protons fully-contained in the ECC.
A. Momentum reconstruction of muon tracks As described in Section II B, INGRID measures muon momentum up to 1 GeV/c. However, most of muons have higher momenta and they penetrate INGRID. Therefore, we adopt the MCS method using the ECC to reconstruct higher momenta. To reconstruct the muon momentum using the coordinate method, three films are used to calculate a positional displacement. The first and second films are used to reconstruct the track angle. Using the reconstructed angle, the track on the second film is extrapolated to the third film. Then, the positional displacement at the third film is measured. The upper limit of the measurable momentum is determined by the measurement error because the scattering angle of a high-momentum particle is small. The positional displacement is proportional to x 3/2 due to the nature of multiple scattering while its measurement error is proportional to x, where x is the thickness of the material between the second and third films. Hence, two films placed further apart distance can measure higher momentum compared to adjacent films because the measurement error becomes sufficiently small compared to the scattering angle. In this analysis, the second and third films are separated by five iron plates, which corresponding to around 1.5 cm and it enables us to measure momentum up to around 5 GeV/c. Films separated by a single water gap are used as the first and second films to reconstruct the track which is extrapolated to a film placed over five iron plates away. This is applied for all available combination of three films. Then, the positional displacement from the predicted position at each combination, y i (i=1, 2, 3...), is measured. The quadrature sum of y i is taken as y 2 meas , which includes both the positional displacement by MCS (y 0 ) and the measurement error (y err ): Therefore, the measurement error needs to be subtracted. The y err is mainly caused by the alignment error and estimated to be less than 10 µm depending on the distance between the segments of the track on the second and third films. Finally, P β is calculated from the following rela- tion [55]: where z is the distance between the second and third films, x is the thickness of the material, X 0 is the radiation length, and t is a correction factor for the effect of passing through several materials. It is assumed that only the iron plates affect the scattering of a particle when we assign values to x and X 0 . If the ECC has a simple structure of a single target material and emulsion films, and the mass of the emulsion films is sufficiently small in comparison to that of the target material, the scattering in the emulsion films is negligible. However, the water ECC contains several types of layers of different materials such as iron, water, emulsion film, and vacuum packing film. Thus, scattering in each material is considered and t is estimated to be 1.4 using the MC simulation. Figure 7 shows the relation between the true and reconstructed momenta of muons from the neutrino interactions in the MC simulation. With this method, our detector can reconstruct the muon momentum with a resolution of 30%-40%.

B. Momentum reconstruction of non-muon tracks
Another way to reconstruct momentum is the angular method. The angular difference between two segments of a track on different films is measured instead of the positional displacement. This method enables us to increase the statistic of the combination of films and to reconstruct the momentum of short tracks. However, the angular method cannot be used for the muon momentum measurement because the measurable momentum is limited by the angular resolution of the films, which is typically 2 mrad for the forward angle tracks. In this analy-sis, the maximum P β measured by the angular method is around 1.5 GeV/c, while P β up to 5 GeV/c can be measured by the coordinate method.
The root mean square of the scattering angle is denoted as θ 0 and it is related to P β as follows [55]: As discussed in the coordinate method, the ECC has a complex structure of several materials. The total scattering angle is considered as the quadrature sum of the scattering angle in each material. The measurement error is also considered and subtracted from the measured scattering angles. The angular method reconstructs the momentum of protons and pions with a resolution of 30%-40%.
In addition, the momenta of protons stopping in the ECC are measured by the track range in the ECC. This method is used only for the tracks identified as protons. In this method, the momentum is reconstructed with a resolution of 5%. Most of the low-momentum protons (typically below 400 MeV/c) are measured by the range.

C. PID
Muon-like tracks are identified by the track matching with INGRID. This section describes the PID of the other tracks. After the P β estimation, all the tracks are separated into nine angle and nine P β bins, and separation of the proton-like and pion-like particles in each bin is performed on the basis of the energy deposit in the emulsion films. The emulsion can measure the energy deposit by the blackness of the track, referred to as the volume pulse height (VPH) [56,57]. Tracks with sufficiently large energy deposits, such as the proton tracks, have large VPH, while MIPs have small VPH. Therefore, the distribution of the VPH shows two peaks. The smaller peak is called the MIP peak and the larger one is called the black peak. Each peak is fitted by a Gaussian function to obtain the mean (µ MIP , µ black ) and the deviation (σ MIP , σ black ). Then, the proton-like likelihood L proton and pion-like likelihood L pion are defined as follows: where v is the VPH of the track. The pion-like likelihood ratio R is defined as Particles with R greater than 0.5 are regarded as pions and those with R less than 0.5 are regarded as protons.
The proton selection efficiency is evaluated as 76.0% with 98.5% purity, while the pion selection efficiency is evaluated as 98.7% with 78.8% purity for the tracks from the neutrino interactions. The VPH decreases over time due to the fading of emulsion. The degree of fading in each film is measured using the sand muon samples, and the correction is applied to the beam timing events in the data. (1) INGRID matching: Track matching between the ECC and INGRID is performed using the SFT. After the selection, a total of 14495 events remain as CC interaction candidates.
(2) Fiducial volume cut: Most of the tracks identified as muons are sand muons. To select neutrino interactions occurring in the ECC, a fiducial volume (FV) is defined as an area of 16 cm × 17 cm from the water-gap next to the most upstream one to the most downstream one. The starting points of the muon candidates are required to be in the FV. After the FV cut, 350 events remain as candidates for the interactions in the ECC.
(3) Viewer check: Track segments might fail to be connected or wrong track segments might be connected due to the inefficiency of track detection on the emulsion films or the failure of the auto-reconstruction process. To find the misconnections and properly determine the muon starting position, all the event candidates are checked by the viewer from the event display. Track segments near the muon starting point are checked whether they are connected to the starting point of the muon track candidate.
There is also a possibility of misidentifying largekink sand muons as neutrino events having a backward-going pion. Additional selections are applied to such kink event candidates found in the viewer check based on the angle, momentum, and VPH of tracks upstream and downstream of the kink position. Backgrounds from kink events are evaluated using sand muons in the MC simulation.

(4) Manual check with a microscope:
After the viewer check, the interaction position is confirmed using a microscope manually. The vertex position of the multiple track events can be determined by extrapolating the track data in an emulsion film obtained by HTS. By contrast, the starting positions of single-track events are not determined by the data. Moreover, the selected events are contaminated by the interactions on the emulsion films and packing films because the scanned data do not contain track segments starting in the middle of the emulsion films. To exclude interactions on the emulsion films, the upstream emulsion film of the vertex position is manually checked. If a track starts in the emulsion film, it is considered as an interaction in the emulsion film and excluded. However, the interactions on the packing films cannot be excluded by this check. Therefore, the background from the interactions on the packing films is evaluated by the MC simulation.
After the viewer and manual checks, 97 events remain as interactions on water, 182 events remain as interactions on iron, and 71 events are excluded as misconnected tracks or interactions in the emulsion films.

(5) Momentum consistency check:
Cosmic rays coming from the downstream region may stop in the ECC, and could be connected to the INGRID tracks induced by the neutrino interactions by chance. The protons and pions from the neutrino interactions also contaminate the muon candidates. To exclude such tracks, the consistencies of the muon momentum measured by MCS in the ECC and that measured by the INGRID range are compared event by event. As cosmic rays stopping in the ECC have low momenta, the momentum measured by the INGRID range becomes larger than that measured by MCS. By contrast, in the case of protons and pions, the INGRID range becomes shorter than that expected by MCS. If the momentum measured by MCS is greater (smaller) than 175% (25%) of that measured by the INGRID range, these events are excluded. In the case of the INGRID-penetrating tracks, the maximum limit is not set because a momentum above 1 GeV/c cannot be measured. By this selection, 11 events are excluded from the neutrino-water event candidates.
The events selected above are considered as the candidates for muons from the neutrino interactions on the ECC water target. The number of selected events after each step is summarized in Table II. In this pilot run, a total of 86 candidate events of CC interactions on the water target are selected, while the MC prediction is 91.6 events. The observed number of events is consistent with the MC prediction within the statistical uncertainty. In the MC prediction, 58.7% of the events are ν µ interactions, while 18.0% are ν µ interactions. Although 23.4% are expected to be background events, cosmic rays are the dominant components and the amount is precisely predicted. The detection efficiency of the CC neutrino interactions within the acceptance of INGRID matching is 63.2%. It corresponds to a detection efficiency of 26.8% with respect to all the CC neutrino interactions on water target in the ECC FV. Following the selections above, the vertex position and multiplicity of each event are determined as below.

(6) Determination of the vertex:
After confirming the muon candidates, the vertex positions are determined. First, tracks that have a minimum distance of less than 600 µm are clustered. For each track, the midpoint of the closest point of that track and the muon candidate is calculated. The center of mass of those midpoints becomes a temporary vertex. Then, tracks that have a minimum distance less than 100 µm from the temporary vertex are clustered again, and their center of mass is regarded as the reconstructed vertex.

(7) Partner track determination:
Finally, partner tracks that make a vertex with the muon track are selected. To determine the track multiplicity, tracks with a minimum distance less than 50 µm from the vertex calculated in the previous step are selected. Track length selections are also applied to exclude very short tracks from nuclear spallations. It is required that the length of tracks with large VPH (black) are two or more layers, and that with small VPH (MIP) are nine or more layers.
After the determination of the multiplicity, the momentum reconstruction and the PID processes are applied. Figure 8 shows the distribution of the pion likelihood ratio. The data distributions are consistent with the MC predictions.
The selection efficiencies of protons and pions from the neutrino interactions are evaluated by the MC simulation. The efficiency is defined as the number of selected tracks divided by the number of tracks within the scanning angular acceptance. Figure 9 shows the result. Blank bins around 90 • are the region outside the acceptance. More than 50% of protons are expected to be detected in all angle regions, even in the 200-400 MeV/c regions.

VII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainty sources are classified into three categories: the neutrino flux, the detector response, and the background estimation. The uncertainty from each source is evaluated from the data and the MC simulation as follows.

Neutrino flux:
The neutrino flux uncertainty and correlations between each neutrino energy bin of both ν µ and ν µ components at the detector position are evaluated as a covariance matrix. The matrix is obtained from the uncertainty of the hadron interaction and the various configurations of the J-PARC neutrino beamline. Figure 10 shows the total flux uncertainty of ν µ and ν µ components in the antineutrino mode beam. Systematic uncertainties from the neutrino flux are calculated using a set of toy MC simulations. Weighting factors on flux bins are thrown according to the flux covariance matrix. Then, the change in the number of predicted neutrino interactions from the nominal value is estimated at each bin. This process is repeated 10 5 times, and the 68% range of the distribution is regarded as the size of the systematic uncertainty.

Detector response:
The uncertainties from the detector response are evaluated by the sand muons and the MC simulation. The reconstruction efficiency and matching efficiency between the detectors are evaluated using the sand muons, and its statistical error is taken as the uncertainty. The uncertainty of the muon momentum reconstruction is evaluated by varying the measurement error in the MC simulation. The uncertainty from the efficiency of the partner track search, momentum reconstruction, PID performance, and uncertainty of the GEANT4 physics list and detector modeling are checked by varying the selection criteria in the MC simulation on the basis of measurement error. The uncertainty of the target mass is also considered. The connection between the ECC and the SFT has around 3% uncertainty, and it is the dominant uncertainty source for muon detection. Around 7% uncertainty is assigned to the PID performance, which is the dominant uncertainty source for the kinematics measurements of the protons and pions.

Background estimation:
The beam-induced background, cosmic rays, and chance coincidence of those tracks at the vertex are considered as background sources. The uncertainty of the beam-induced background mainly originates from the normalization of the sand muons. The number of sand muons in the MC simulation is normalized with the data. There is a difference of around 30% between the MC prediction and the data; this is taken as the uncertainty. The cosmic background is induced by the misconnection of the track matching; however, the uncertainty is less than 1%, as the rate of chance coincidence is precisely estimated using mock data. The uncer- tainties of these backgrounds are sufficiently small compared to other uncertainties in most regions.
Besides the uncertainties above, uncertainties of the neutrino interaction modeling are evaluated as follows, and they are compared to the other uncertainties.

Neutrino interaction:
The neutrino interaction model and the FSI in NEUT have many uncertainties. Such uncertainties change the expected number of events in each measurement. Uncertainties from these sources are evaluated by changing parameters in the model based on the current understanding of the neutrino interaction model and the FSI [18]. Table III shows the nominal value and the 1σ error size of each parameter. In this analysis, the uncertainty of the nuclear binding energy is not evaluated. However, the uncertainty is covered by the comparison with an alternative nuclear model discussed in Section VIII. After evaluating the uncertainty induced by each parameter, the normalization of 2p2h interaction is found to be about 8%, making it the dominant uncertainty source in this analysis. The change of the detection efficiency by the change of the interaction model is separately estimated. The typical value is around 1-2% in each bin. It is added to the uncertainty of the detector response. Figures 11 and 12 show the systematic uncertainties of each measurement with a breakdown by category. The fractional uncertainty of the expected number of selected events in each bin is plotted. The quadrature sums of the uncertainties from the neutrino flux, the detector response, and the background estimation are smaller than the uncertainty of the neutrino interaction model in almost all bins. An uncertainty of only 5%-8% is derived from the flux uncertainty owing to the great improvement in the hadron interaction modeling by NA61/SHINE. The current detector uncertainty is slightly larger than the flux uncertainty, and desired to be improved in future analysis.

VIII. RESULTS
In this measurement, the statistics is not sufficient to precisely extract the cross section from the reconstructed distributions. Thus, distributions of the charged particles from the neutrino interactions on the water are compared with model predictions to get insights into the model validity and to demonstrate the feasibility of the NINJA detector for future measurements. It should be noted that our results include low-momentum charged particles, especially protons with momenta of 200-400 MeV/c, owing to the high granularity of the emulsion films.
First, raw data distributions are compared with sum of the neutrino event prediction estimated with the MC simulation and the cosmic-ray background prediction estimated with the off-beam timing data. Figure 13 shows the multiplicity of the charged particles and the number of pions and protons. The red boxes on the prediction correspond to the quadrature sum of the uncertainties of neutrino flux, the detector response, and the background estimation. Figure 14 shows distributions of the reconstructed kinematics of muons, pions, and protons. Though the angular resolution for all particles is sufficiently small compared to the bin width in the angle plots, the momentum resolution is typically larger than the momentum binning, especially for high momentum muons. In the proton momentum distribution, protons with momentum 200-400 MeV/c are successfully detected for the first time in measurements of neutrinowater interactions.
In order to extract the signal distributions from our data and to compare them with the prediction, backgrounds from neutral-current interactions, interactions on the packing films, cosmic rays, and chance coincidence of the off-beam timing tracks are subtracted from the data using the background prediction. Figures 15  and 16 show the results, in which the signal distributions are compared with the predictions. In Fig. 16, contaminations by misidentification between protons and pions are also subtracted. In these plots, the flux, detector response, and background estimation uncertainties are included in the error bars on the data points, while the hatched regions correspond to the uncertainty of the neutrino interaction model. The measurement's systematic error is smaller than the current model uncertainty. Hence measurements with the NINJA detector can be expected to constrain neutrino interaction models given more statistics.
Although the statistical uncertainty is large, the measurement result shows a slightly lower multiplicity of charged particles compared to the prediction. As shown in the bottom plots of Fig. 15, there is a tendency for the prediction to overestimate the number of pions. The number of detected pions is 4.9 ± 3.6 (stat.) ± 0.6 (syst.), while 14.5 are expected in the prediction. By contrast, 15.2 ± 4.2 (stat.) ± 0.8 (syst.) protons are detected, which is consistent with the prediction of 17.7 protons. The overestimation of charged pions may be induced by the inaccuracy of the modeling of either neutrino interactions or FSI. Besides this overestimation, the predictions explain the data well.
In addition to the one-dimensional kinematics distributions, Fig. 17 shows the two-dimensional distribution of the angle and the momentum for protons or pions. Although the statistics is limited, these plots show general agreement between the data and the predictions.
Finally, an alternative model of NEUT using SF, and another generator, GENIE [58,59], are studied for comparisons with the nominal model of NEUT using LFG. Figures 18 and 19 show the results. The interaction models used in the nominal MC simulation is summarized in Tab. I. Interaction models used in the alternative model of NEUT using SF is almost the same as those in Tab. I, but the nuclear model is changed to SF, and M QE A is set to 1.21 GeV/c 2 . GENIE v3.0.6 with G18 10b 02 11a tuning is used as an alternative generator. In GENIE, different axial mass values are used, and the Berger-Sehgal model [60] is used for the single pion production. For the FSI simulation, GENIE hN cascade model [61] is employed. With more statistics, our measurement can discriminate the models.

IX. CONCLUSION
The first results of the NINJA pilot run using a water-target emulsion detector are reported in this paper. Multiplicity, angle, and momentum distributions of the outgoing muons, charged pions, and protons from neutrino-water interactions are reported. Protons from neutrino-water interactions are measured with a 200 MeV/c threshold for the first time. Although the statistical uncertainty is large, we found that the current neutrino interaction models predict the kinematics distributions well within the measurement uncertainty, including for low momentum protons down to 200 MeV/c. Besides, in the measurements of pion kinematics, we found that there is a tendency to overestimate the number of charged-pions in the MC simulation. The related data shown in this paper can be found in [62].
The first physics run of the NINJA experiment concluded its beam exposure with the neutrino mode beam in early 2020. We expect 15 times more neutrino interactions, thus the statistical uncertainty will be as small as the current systematic uncertainty. In the current analysis, relatively large systematic uncertainty is applied to the measurements of pions and protons due to the uncertainty of the PID performance. This uncertainty can be reduced through further understanding of our detector response, and the size of the total uncertainty will be similar to that of the muon measurements. Then, we will measure the differential cross-section with about 10% uncertainty. We aim to characterize the nature of 2p2h interactions, which has the largest uncertainty in the current model. Moreover, differential cross-section measurements with respect to the number of protons with  FIG. 14. Distributions of muon, pion, and proton kinematics from neutrino-water interactions and backgrounds. The left column shows angle distributions, while the right column shows reconstructed momentum distributions. Though the angular resolution for all particles is sufficiently small compared to the bin width in the angle plots, the momentum resolution is typically larger than the momentum binning, especially for high momentum muons. The data points are shown by marker points with statistical error bars and the predictions are shown by histograms with systematic uncertainties as red boxes, which are the quadrature sum of the uncertainties of neutrino flux, the detector response, and the background estimation.   FIG. 16. Distributions of muon, pion, and proton kinematics. Backgrounds are subtracted from both the data and the prediction. The left column shows angle distributions, while the right column shows momentum distributions. Though the angular resolution for all particles is sufficiently small compared to the bin width in the angle plots, the momentum resolution is typically larger than the momentum binning, especially for high momentum muons. The flux, detector response, and background estimation uncertainties are included in the error bars on the data points, while the hatched regions correspond to uncertainty of the neutrino interaction model.  Backgrounds are subtracted from the data using the nominal prediction. The left column shows angle distributions, while the right column shows momentum distributions. Though the angular resolution for all particles is sufficiently small compared to the bin width in the angle plots, the momentum resolution is typically larger than the momentum binning, especially for high momentum muons. The flux, detector response, and background estimation uncertainties are included in the error bars on the data points.