Measurement of the single $\pi^0$ production rate in neutral current neutrino interactions on water

The single $\pi^0$ production rate in neutral current neutrino interactions on water in a neutrino beam with a peak neutrino energy of 0.6 GeV has been measured using the P{\O}D, one of the subdetectors of the T2K near detector. The production rate was measured for data taking periods when the P{\O}D contained water ($2.64\times{}10^{20}$ protons-on-target) and also periods without water ($3.49 \times 10^{20}$ protons-on-target). A measurement of the neutral current single $\pi^0$ production rate on water is made using appropriate subtraction of the production rate with water in from the rate with water out of the target region. The subtraction analysis yields 106 $\pm$ 41 (stat.) $\pm$ 69 (sys.) signal events, which is consistent with the prediction of 157 events from the nominal simulation. The measured to expected ratio is 0.68 $\pm$ 0.26 (stat.) $\pm$ 0.44 (sys.) $\pm$ 0.12 (flux). The nominal simulation uses a flux integrated cross section of $7.63\times{}10^{-39}$ cm${}^2$ per nucleon with an average neutrino interaction energy of 1.3 GeV.


(The T2K Collaboration)
The single π 0 production rate in neutral current neutrino interactions on water in a neutrino beam with a peak neutrino energy of 0.6 GeV has been measured using the PØD, one of the subdetectors of the T2K near detector. The production rate was measured for data taking periods when the PØD contained water (2.64 × 10 20 protons-on-target) and also periods without water (3.49 × 10 20 protonson-target). A measurement of the neutral current single π 0 production rate on water is made using appropriate subtraction of the production rate with water in from the rate with water out of the target region. The subtraction analysis yields 106 ± 41 (stat.) ± 69 (sys.) signal events, which is consistent with the prediction of 157 events from the nominal simulation. The measured to expected ratio is 0.68 ± 0.26 (stat.) ± 0.44 (sys.) ± 0.12 (flux). The nominal simulation uses a flux integrated cross section of 7.63×10 −39 cm 2 per nucleon with an average neutrino interaction energy of 1.3 GeV.

I. INTRODUCTION
The Tokai to Kamioka (T2K) long-baseline neutrino experiment is designed to make precision measurements of the neutrino oscillation parameters θ 23 and ∆m 2 32 via ν µ disappearance and to search for the mixing angle θ 13 via ν e appearance in a ν µ beam. An intense, almost pure beam of ν µ is produced by colliding 30 GeV protons with a graphite target at the J-PARC facility in Tokai-mura, Ibaraki [1]. The resultant neutrino beam is directed 2.5 • away from the axis between the target and the far detector, resulting in a narrow band beam with peak energy near 0.6 GeV. The direction, stability and flux of the beam are measured using a suite of near detectors (ND280) located 280 m downstream of the target. At this distance, the neutrino beam is not expected to have been affected by oscillations. The far detector, Super-Kamiokande, is located 295 km downstream of the target, a distance consistent with the oscillation maximum. Super-Kamiokande uses water as both a detection medium and target to measure the amount of ν e and ν µ present after oscillation has occurred. As neutral current π 0 events can cause an irreducible background to the ν e appearance signal, it is important to provide a constraint using measurements of the production rate on water using the near detector.
This paper details the first measurement of neutral current single π 0 production (NC1π 0 ) on water, using a neutrino beam with peak energy of 0.6 GeV [2]. The mean neutrino energy for the NC1π 0 interactions selected in this analysis is 1.3 GeV.
Two processes dominate neutral current single π 0 production by neutrinos: resonant production and coherent scattering. In resonant production, a neutrino interacts with a nucleon to produce a baryonic resonance, usually ∆(1232), which subsequently decays to a nucleon and a π 0 . Coherent scattering occurs when a neutrino interacts with the entire nucleus, exchanging little energy and leaving the nucleus in its ground state. The dominant decay mode for a π 0 is to two photons [3] and if one decay photon is not detected, an NC1π 0 event can be indistinguishable from a charged current ν e interaction, leading to an irreducible background in ν µ → ν e oscillation measurements. Whilst previous measurements performed using the T2K near detector have improved our knowledge of sub-GeV neutrino interactions, the rate of NC1π 0 production on water is still relatively unknown at the neutrino energies of the T2K beam. Measurements of NC1π 0 production on a variety of targets and different neutrino energy distributions have been made at other experiments [4][5][6][7][8].
In this analysis, the signal is defined by the final state particles, with an NC1π 0 interaction defined by a single π 0 particle exiting the nucleus along with any number of protons and neutrons but no charged leptons or other mesons. The rate of signal events on water is determined using event samples with a two photon signature from exposures with water in and out of the target region. Using the presence of a muon decay tag, two photon candidate events are divided into signal-enriched and background-enriched samples. The number of signal events, number of background events, energy scale, and shape of the background are then determined for each sample using a simultaneous maximum likelihood fit to the invariant mass distribution of the signal-enriched and background-enriched samples. Finally, the number of interactions on water is determined using a weighted subtraction of the rate of signal events determined during the exposure with water in the target region and the exposure with the water removed.
The major sections of this paper are as follows. Section II describes the T2K ND280 π 0 detector, as well as the simulation of the expected neutrino interactions and detector response. Section III describes the event selection efficiencies and reconstruction resolutions for signal-enriched and background-enriched event samples and the selected event samples are described in Section IV. The extraction of the number of signal events is described in Section V followed by a discussion of the systematic uncertainty in Section VI. Section VII describes the calculation of the event rate on water and compares it with the expectation.

II. T2K ND280 PØD DESCRIPTION AND SIMULATION
The T2K ND280 π 0 Detector (PØD) is a scintillator-based tracking calorimeter optimized to measure NC1π 0 production in the momentum range that contributes backgrounds to ν e appearance measurements [9]. The PØD is composed of layers of plastic scintillator alternating with water bags and brass or lead and is one of the first large scale detectors to use Multi-Pixel Photon Counters (MPPCs). Relative to the neutrino beam, it sits upstream of a tracking detector made up of two fine grain scintillator modules placed between three time projection chambers. Both the PØD and tracking detector are in a 0.2 T magnetic field and surrounded by electromagnetic calorimeters and muon range detectors [10][11][12][13].
The PØD comprises 40 scintillator modules, each 38 mm thick, formed from two layers of scintillating bars with the long axis oriented either horizontally, or vertically, and instrumented using wavelength shifting fibers with an MPPC on one end and mirrored on the other [9]. The triangular scintillating bars used to produce each of the two layers in each module have a height of 17 mm and a base of 32 mm and are interlocked to form a layer that is 17 mm thick. Two views are formed of an event, commonly labeled the X-Z, and the Y-Z view, where the z axis is horizontal and points downstream, the y axis points in the vertical direction, and the x axis is perpendicular to the Y-Z plane. A minimum ionizing particle will typically generate a charge in the MPPC equivalent to approximately 38 photoelectrons/cm, or an average of about 30 photoelectrons in a single bar. The scintillator modules are arranged in three regions. The most upstream and downstream regions are made of seven modules interleaved with 4.5 mm thick sheets of stainless steel-clad lead that function as 4.9 radiation length electromagnetic calorimeters to improve the containment of photons and electrons. The central region serves as a target containing water. It has 25 water target layers that are 28 mm thick interleaved with 26 scintillator modules and 1.3 mm brass sheets. When water is in the detector, the target fiducial region contains approximately 1900 kg of water and 3570 kg of other materials. Data collected with and without water in the PØD are analyzed separately. This analysis utilizes data collected with a predominantly ν µ beam generated between January 2010 and May 2013, see [2] for a detailed description. The neutrinos are generated using a fast extracted 30 GeV proton beam with a spill of 6-8 bunches that are separated by 582 ns. The proton beam strikes a graphite target producing pions and kaons which, after magnetically focusing the positive mesons, decay in flight to neutrinos. The magnetic focusing can be altered to focus negative mesons. The T2K runs, the configuration of the PØD, and the corresponding protons on target (POT) are summarized in Table I.
The simulated data set used in this analysis corresponds to 4.01 × 10 21 POT (water-out configuration), and 7.18 × 10 21 POT (water-in configuration). Neutrino interactions are simulated using the NEUT [14] event generator, version 5.1.4.2, with the interactions distributed within the full ND280 volume, as well as the surrounding hall. Interactions on all nuclear targets present in ND280 are simulated. Details of the neutrino interaction simulation process are described in [14], [15] and [16]. The T2K run periods are simulated using the nominal detector and beam configurations, and then combined using the appropriate POT normalization to form the final expectation. External, non-beam associated, backgrounds are not simulated, but are limited in the data sample by the duty cycle of the neutrino beam. Particles produced in neutrino interactions are simulated using GEANT 4.9.4 [17]. The standard GEANT physics list for electromagnetic interactions is used in the simulation.
Neutral current single π 0 production in the T2K neutrino beam is dominated by resonant ∆(1232) production, which is simulated using the Rein-Sehgal [18] model for neutrino-induced resonant pion production. The simulated NC1π 0 cross section on water integrated over the T2K neutrino beam flux is 7.63 × 10 −39 cm 2 or 4.24 × 10 −40 cm 2 nucleon −1 while the NC1π 0 cross section for the fiducial region in the water-out configuration is 4.20 × 10 −40 cm 2 amu −1 [2,14]. There is an additional 12% uncertainty in the neutrino flux integrated over the energy of neutrinos generating an NC1π 0 interaction [2]. This uncertainty is larger than presented in other T2K analyses because of the higher average neutrino energy, and, for this analysis, is unconstrained by other near detector measurements to allow direct comparison between the data and simulation.

III. EVENT RECONSTRUCTION AND SELECTION
Events are reconstructed in the PØD using scintillation light signals that occur in time windows containing the neutrino bunch arrival. A hit is constructed from the integrated charge during each time window, and the time relative to the start of the window at which the integrated light signal crosses a threshold equivalent to approximately 2.5 photoelectrons. Activity in different time windows is independently reconstructed as separate events. The reconstruction proceeds by selecting groups of hits consistent with a track-like signature that is classified as a light-ionizing track, such as a muon or charged pion, a heavy-ionizing track, such as a proton, or a non-track object such as a portion of an electromagnetic shower. Hits from non-track objects, as well as any hits not gathered into a track-like object are then used to form groups that are consistent with shower-like particles such as photons or electrons coming from a single vertex. In events without a track-like signature, the vertex is estimated by assuming that the particle signatures in the event emanate from a single point, with the particle directions going away from the vertex. While events with track-like objects will generally be rejected in the later analysis, if a track-like object is found, then the vertex is fixed at the upstream end of the longest track. After vertex reconstruction, all reconstructed non-track like objects, are classified as either EM-like, or shower-like. The shower-like objects primarily comprise interacting pions, interacting protons, or misidentified light-ionizing tracks. The result of the reconstruction is a single vertex with an associated collection of objects corresponding to light-ionizing tracks, heavy-ionizing tracks, EM-like, and other shower-like objects. A muon decay tag is associated with the reconstructed vertex when energy deposition consistent with a Michel electron is found.
A signal-enriched sample of exactly two photon candidates with invariant mass less than 500 MeV/c 2 is selected using eight selection criteria: event quality, vertex in the fiducial volume, energy containment in the PØD, lack of a muon decay signature, fraction of energy in the two most energetic photon candidates, particle identification, reconstructed direction and object separation. In comparison to the signal, a distinguishing characteristic of the background is that it contains either a µ, or a charged pion, both of which can generate a muon decay signature, so a separate background-enriched sample is selected by applying all criteria with the exception of the muon decay criterion which is reversed.
To be considered in this analysis, an event must occur during a neutrino beam spill and have a single reconstructed vertex as well as good data quality. The vertex must be in the fiducial volume defined as at least 25 cm from the edge of the active volume and inside the water target region of the PØD [19]. The containment criteria requires that all reconstructed objects are contained inside the PØD by requiring that no reconstructed objects have hits in the last layer of the PØD or in the outer two bars of any layer. This limits external background and improves the photon energy reconstruction.
The signature of interest is two reconstructed photons from the π 0 decay with no evidence of a muon-like object. To ensure that selected events have two reconstructed photon candidates containing most of the recorded energy deposition, a "charge-in-shower" requirement is placed on the fraction of energy in the two most energetic EMlike objects. The required fractions of 92% (water-in), and 80% (water-out) were chosen to optimize the statistical significance of the selected number of signal events using simulated samples. Due to the planar nature of the PØD, and the shape of the scintillator bars, the performance degrades for particles at an angle of more than approximately 75 • from the z axis. As such, the direction of the reconstructed total event momentum must be less than 60 • from the z axis, limiting the phase space covered by this measurement.
Two well-separated decay photon candidates are required to limit the background from particles with overlapping energy deposits. The object separation in each projection is calculated by finding the distance between the two closest hits of the reconstructed objects. Due to the planar nature of the PØD, it is possible for two objects to overlap in one projection, but not in the other and separation is only required in one of the two projections. The object separation is required to be greater than 9 cm (14 cm) in at least one projection for the water-in (water-out) configuration. Figure 1 shows the position of the reconstructed vertex relative to the true vertex position along the z axis for the water-out configuration of the PØD for simulated NC1π 0 events that have passed all selection criteria. For the water-in (water-out) configuration, the biases are −0.06 cm (0.08 cm) along the x axis, 0.06 cm (0.20 cm) along the y axis, and 1.67 cm (1.72 cm) along the z axis. The expected vertex residual distribution is asymmetric due it's dependence on the reconstructed photon position and direction, and is characterized by half the distance between the 16% and 84% quantiles. For the water-in (water-out) configuration, the resolutions are 5.5 cm (6.8 cm) along the x axis, 6.1 cm (8.0 cm) along the y axis, 8.6 cm (11.2 cm) along the z axis.
The momentum resolution of the π 0 is a combination of the energy and angular resolution for two reconstructed photons. The total energy for the reconstructed photons is determined calorimetrically, and the fraction of the total energy carried by each reconstruct photons is calculated using the projections where the photon objects are geometrically distinct. Figure 2 shows the fractional momentum residual (the difference of the reconstructed and true momenta divided by the true momentum) for the water-out configuration of the PØD for NC1π 0 events passing all selection criteria. A Gaussian is fit to the central region to determine the shift and width of the momentum distribution. The fractional momentum residual distribution has a mean of −3.2% with a width of 18.7% for the PØD water-in configuration. For the water-out configuration, the mean is −0.8% with a width of 21.1%. The reconstructed opening angle distribution for simulated NC1π 0 events passing all selection criteria in both the water-in and water-out configurations has a mean of −0.01 rad from the nominal value and an RMS of 0.06 rad. The reconstruction efficiency, , for an NC1π 0 event and signal-enriched sample purity are summarized in Table II. The efficiency is defined as the number of true NC1π 0 events reconstructed in the fiducial volume divided by the number of NC1π 0 interactions occurring within the same volume, while the purity is defined as the fraction of the selected events which result from a true NC1π 0 interaction. The average efficiency is 6.10% for the water-in configuration, and 4.79% for the water-out configuration. There is a small location dependence in the efficiency for the water-in configuration, so the efficiency is tabulated separately for interactions which occur in the water target, and for interactions which occur on another material. The average purity for the water-in (water-out) configuration is 48.7% (46.1%) for all events with a two photon invariant mass less than 500 MeV/c 2 corresponding to a rejection of more than 99.5% of the background events. Figure 3 shows the efficiency of the NC1π 0 selection as a function of the true π 0 momentum.

IV. SELECTED EVENT SAMPLES
Tables III and IV show the number of observed and expected events found in the signal-enriched and backgroundenriched samples. The expectation for each sample is broken down into the number of expected signal and background events, and the number of background events is further broken down by the presence of charged leptons with and without a π 0 in the final state of the neutrino interaction. Categories are also included for simulated events containing multiple neutrino interactions, and background entering from outside the PØD. Approximately 10% of the events in the background-enriched sample are due to signal interactions. In the data, 775 events were selected as an NC1π 0enriched sample for the water-in configuration and 555 events were selected for the water-out configuration of the PØD compared to an expectation of 893 (629) for the water-in (water-out) configuration. The distribution of the true neutrino energy for the selected sample of simulated events is shown in Figure 4 separated by event topology with the mean neutrino energy for the NC1π 0 signal being 1.3 GeV.
The event signature for this analysis is two reconstructed photons with an invariant mass, M , close to that of the π 0 . The reconstructed invariant mass is M = 2E 1 E 2 (1 − cos θ) where E 1 and E 2 are the reconstructed energies of each photon candidate and θ is the angle between the photon candidates. Figure 5 shows the distribution of the selected events in the signal-enriched and background-enriched samples where the expectation for each distribution has been normalized to the observed number of events. The reconstructed energy distribution of the signal-enriched samples is shown in Figure 6. The expected composition for each distribution is shown using the same breakdown as in Tables III and IV, however, the contributions from external and multiple interactions have been combined into a single category.

V. EXTRACTING THE SIGNAL EVENT RATE
The number of NC1π 0 events is found using a six parameter unbinned extended maximum likelihood fit to the invariant mass distribution of the signal-enriched and background-enriched samples. Four of the parameters, N SE Sig , N SE Bkg , N BE Sig , and N BE Bkg , are related to the number of signal and background events in the signal-enriched (SE) and background-enriched (BE) samples. The remaining two parameters control the energy scale of electromagnetic particles relative to minimum ionizing tracks, and the shape of the expected background.
The likelihood is extended by assuming that the probability of the observed numbers of signal-enriched (N SE γγ ) and background-enriched (N BE γγ ) events is given by the product of Poisson distributions and, in each case, the expected number of two photon events is a sum of the signal events (N Sig ) and the background events (N Bkg ). The expected invariant mass distribution for the signal (background) events in the signal-enriched (background-enriched) sample is generated using the simulation after event reconstruction. The distributions are normalized such that the sum over the signal-enriched bins is equal to the total number of events, N SE γγ , and, likewise, for the background-enriched sample, N BE γγ . The ratio of the number of signal and background events in the signal-enriched sample to the numbers in the background-enriched sample (N SE Sig /N BE Sig and N SE Bkg /N BE Bkg ) are determined by the efficiency of the muon decay tag and the probability of a muon decay tag false positive. Both relations are estimated using a sample of stopping muons from neutrino interactions occurring upstream of the PØD. This sample is selected by requiring a single track-like object entering the upstream face of the detector, and stopping in the water target region.
The number of background events in the signal-enriched sample is related to the background events in the  background-enriched sample by the muon decay reconstruction efficiency and is allowed to vary within the uncertainty on the muon decay tag efficiency. For the water-in (water-out) configuration, the expected efficiency is 45.6% (43.9%) and the observed efficiency is 44.1 ± 0.5% (46.2 ± 0.6%). The fractional difference between data and expectation is combined with its statistical error and used as a Gaussian constraint in the likelihood on the ratio between the number of background events in the signal-enriched and background-enriched samples. The constraint is 3.4% for the water-in configuration and 5.2% for the water-out configuration. The fitted number of signal events in the background-enriched sample relative to the number in the signal-enriched sample is allowed to vary within the uncertainty on the probability of a false positive muon decay tag. The uncertainty in modeling the false muon decay tag rate has been estimated by fitting the time distribution of muon decay tags occurring after a stopping muon, but within the same trigger window, to an exponential plus a constant. The fitted exponential lifetimes are consistent with the expectation for muon decay, and the constant term estimates the probability of incorrectly finding a muon decay tag. For the water-in (water-out) configuration there is a 1.1 ± 0.5% (1.3 ± 0.7%) difference in the constant term between the data and expectation which provides a ±1.6% (±2.0%) constraint on the false tag probability.
Since there is an uncertainty in the shape of the background underneath the π 0 invariant mass peak, an extra shape parameter has been added to the fit. The deviation from the expected background shape is assumed to have the same shape as the signal probability distribution, while the normalization is constrained by the background-enriched sample. The shape factor is allowed to be positive or negative meaning that the amount of background in the region of the π 0 invariant mass can be either increased or decreased. Two cases are considered in the fit. In the first instance, the number of signal and background events are determined by using the nominal shape for the background which is equivalent to fixing the shape parameter to a value of zero. In the second case, no prior constraint is placed on the shape parameter, and the uncertainty in the rate is estimated by constraining it with both the signal-enhanced and background-enhanced samples. See Section VI where this case is used to estimate the systematic uncertainty due to the unknown background shape.
The overall energy scale in the PØD is set using penetrating muon tracks, and must be translated to an energy scale for electromagnetic particles with uncertainty introduced due to the relative response of the detector to different particle types. The final electromagnetic energy scale is determined using the position of π 0 invariant mass peak. The mean difference between the reconstructed and true photon opening angle in the simulation is small (0.01 rad) and has negligible effect on the invariant mass distribution so the difference between the measured and simulated invariant mass scales is assigned to the energy scale uncertainty. No prior constraint is placed on the energy scale parameter, however, based on a survey of the detector material distribution and the uncertainties in the particle propagation model, the prior uncertainty is approximately 10% relative to the energy scale determined using penetrating muons. In the water-in (water-out) configuration, the fitted value for the electromagnetic energy scale parameter is 89.5±3.4% (96.7 ± 0.6%).
The best fit values for the number of signal and background events with the energy scale parameter unconstrained, while using the nominal shape for the background, are shown in Table V. Figure 7 compares the invariant mass expectation to the data, where the energy scale correction determined during the fit has been applied to the data. Because the criteria requiring events have a reconstructed invariant mass less than 500 MeV/c 2 is applied prior to determining the best fit energy scale, the mass bins above 440 MeV/c 2 are not fully populated with data. The goodness of fit is calculated as a binned χ 2 of the invariant mass distributions between 0-440 MeV/c 2 where the range has been limited to the region that the data populates. Considering only statistical uncertainty, the χ 2 value for the PØD water-in configuration is 40.4 for 39 degrees of freedom, leading to a p-value of 0.41. The χ 2 value for the PØD water-out configuration is 53.5 for 39 degrees of freedom, leading to a p-value of 0.06.

VI. SYSTEMATICS
The systematic uncertainties are summarized in Table VI, and are described below. The detector systematics are separately estimated for the water-in and water-out configurations. Since the detector performance is different, and run periods do not overlap in time, the systematic uncertainty related to detector performance is assumed to be uncorrelated between the water-in and water-out configurations.
Because the event reconstruction proceeds in two stages, first reconstructing track-like signatures in each event, and then reconstructing the remaining activity assuming showering signatures, reconstruction efficiencies primarily affect the result in two ways. First, an inefficiency is introduced when an electromagnetic object is reconstructed as track-like, because the object will then not be considered by the shower reconstruction. Other efficiencies are more closely related to the shower reconstruction, including efficiencies related to the particle identification of showering signatures, the reconstructed distance between showering objects, and the fraction of the visible energy assigned to each showering object.
The track particle identification efficiency uncertainty is estimated using the sample of stopping muons described in Section V. The uncertainty for each input parameter to the particle identification procedure is estimated and propagated through the particle identification likelihood to determine the effect on the identification efficiency. For simulated muons, there is a 5.40 ± 0.05% (5.06 ± 0.03%) uncertainty in the misidentification rates for the water-in (water-out) configuration. Combining the difference in quadrature with the statistical uncertainty leads to a total track object particle identification systematic of 5.4% (5.1%) for the water-in (water-out) configuration.
The efficiencies of the charge-in-shower, object separation and shower particle identification criteria are related to the properties of a showering particle and are studied using control samples selected by reversing these cuts to create double "side-band" distributions. For example, to estimate the uncertainty in the efficiency of the chargein-shower criterion, events that fail the object separation and shower particle identification criteria, but which pass all other criteria, are selected to create a control sample with low signal purity. The estimated uncertainty for the efficiency of the charge-in-shower criterion is then the relative difference between the percentage of control sample events passing the criterion relative to the expectation. For the water-in configuration, 54.0% of the simulated, and 51.3 ± 2.2% of the data control sample events are selected, combining the difference and the statistical error in quadrature leads to a systematic uncertainty of 6.6% in the efficiency due to the charge-in-shower criterion. A similar calculation is done for the water-out configuration. The procedure is then repeated for the object separation and shower particle identification criteria. Because these three uncertainties are estimated using statistically limited data sets collected during independent water-in and water-out run periods they are assumed to be uncorrelated between the configurations, and the uncertainty will likely be reduced by the collection of additional data. The on-water uncertainty for these uncertainties is estimated by combining the water-in and water-out uncertainties in quadrature (summarized in Table VI). After the best fit values were found in Section V, the fitted value and uncertainty on the energy scale were used to estimate the effect of the energy scale on the estimated efficiency. This effect was modeled by scaling the NC1π 0 reconstruction efficiency shown in Figure 3 using many trials of the energy scale parameter distributed according to the best fit parameter and statistical uncertainty. The shifted efficiency curve represents a new expectation for the trial energy scale parameter and is used to estimate the expected number of saved signal events in the simulation. The fractional shift and RMS of the distribution of the expected number of signal events are then added in quadrature to estimate the uncertainty due to the energy scale for the estimated efficiency. The time variation of the energy scale was tracked using through-going minimum ionizing particles, and it introduces an efficiency uncertainty of 1.8% for both the water-in and water-out configurations.
Several systematic uncertainties were associated with the fiducial volume. Uncertainties on measurements of the detector mass were used to reweight the selected events to extract the uncertainty due to fiducial mass. The effect of alignment between detector elements on the efficiency was studied and found to be negligible (< 0.1%). Additionally, there are two fiducial volume uncertainties. One reflects how the result is affected by changing the fiducial volume definition, while the other quantifies the uncertainty due to a systematic shift between the simulated and true detector volumes. When combined, the fiducial volume uncertainty is 1.9% for water-in and 2.6% for water-out.
The systematic uncertainty due to the background shape is estimated by comparing the effect on the fitted signal rate with the shape parameter fixed to when it is unconstrained. Following the procedure outlined in Section V, the best values for the shape parameter with the water-in and water-out configurations are found to be statistically consistent with the χ 2 value for the invariant mass distribution being 47.5 for 38 degrees of freedom for the water-in configuration and 38.7 for 38 degrees of freedom for the water-out configuration. Because the backgrounds in both the water-in and water-out configuration arise from the same physical processes, it is assumed that this uncertainty is fully correlated between the water-in and water-out configurations and the uncertainty is directly applied to the on-water signal event rate. The change in the water-in and water-out event rates between the case where the shape parameter is fixed to zero and where the shape parameter is free leads to a 19.8% uncertainty in the on-water rate due to the background shape parameter.
While this analysis uses background-enriched control samples to minimize the uncertainty due to the model of the cross sections, changes in the generated event kinematics have an effect. This was studied using the neutrino flux shape and cross section uncertainties detailed in [2] and [16]. For the water-in configuration, the flux shape and cross section introduce a 2.5% uncertainty, and for the water-out configuration, a 3.3% uncertainty. Because rate of observed events is consistent with the expectation, the full 12% flux uncertainty is included and given as a separate uncertainty.

VII. NUMBER OF EVENTS ON WATER
The number of NC1π 0 events in the PØD measured using both the water-in and water-out configuration (Table V) can be used to determine the number of events occurring directly on water. The measured number of NC1π 0 events with water in the PØD is found to be 342 ± 33 (stat.) ± 88 (sys.) during an exposure of 3.49 × 10 20 POT, where the systematic uncertainty includes effects that are correlated between the water-in and water-out configurations. The ratio between the observed and expected rate is 0.79 ± 0.08 (stat.) ± 0.20 (sys.). Similarly, with water out of the PØD, the measured number is 246 ± 26 (stat.) ± 61 (sys.) for an exposure of 2.64 × 10 20 POT, and the ratio between the observed and expected rate is 0.85 ± 0.09 (stat.) ± 0.21 (sys.). To allow direct comparison to the expected NC1π 0 event rate, the quoted ratios include neither the 12% flux normalization uncertainty nor the NC1π 0 cross section uncertainty The total number of events on water is found using a statistical subtraction by relating the event rate during the water-in exposure to the event rate during the water-out exposure. The total number of signal events in the water-in (WI) configuration can be divided into two parts, N WI = N On-Water + N NW , where N NW is the number of signal events that occur on targets other than water (referred to as "not-water" events) in the water-in configuration. The number of not-water (NW) events, which is proportional to the number of water-out (WO) events, can be subtracted from the total number of on-water events by where the efficiencies, NW and WO , are given in Table II, and the POT is given in Table I. After the subtraction in Equation 1, 106 ± 41 (stat.) ± 69 (sys.) events were found, where the uncertainties that are correlated between the water-in and water-out configurations have been taken into account. The simulation predicts 157 true NC1π 0 events on water. The ratio of the number of measured to number of predicted on-water events, including the correlated and uncorrelated uncertainties described in Section VI is 0.68 ± 0.26 (stat.) ± 0.44 (sys.) ± 0.12 (flux), where the NC1π 0 cross section uncertainties are excluded.

VIII. CONCLUSION
An on-water NC1π 0 rate measurement has been performed by combining data from a 2.64 × 10 20 POT neutrino beam exposure of the T2K ND280 PØD using a water-in configuration with a 3.49 × 10 20 POT exposure using a water-out configuration. This is the first use of the subtraction method to measure neutral current event rates with the T2K near detector.
The signal event rates are found using an extended maximum likelihood fit to the reconstructed invariant mass for each sample in a range of 0-500 MeV/c 2 . As described in Section III, the phase space of the analysis has been limited to the region where the PØD has acceptance. The analysis finds 342 ± 33 (stat.) ± 88 (sys.) (246 ± 26 (stat.) ± 61 (sys.)) signal events in the PØD water-in (water-out) data compared to an expectation of 433 (290) events. Excluding the 12% normalization and NC1π 0 cross section uncertainties, the resulting observed to expected ratios are 0.79 ± 0.08 (stat.) ± 0.20 (sys.) for water-in and 0.85 ± 0.09 (stat.) ± 0.21 (sys.) for water-out configurations. Subtracting the water-in and water-out samples after correcting for the different POT and reconstruction efficiencies yields 106 ± 41 (stat.) ± 69 (sys.) signal events on water compared to an expectation of 157 events and an on-water NC1π 0 production rate of 0.68 ± 0.26 (stat.) ± 0.44 (sys.) ± 0.12 (flux) relative to the NEUT expectation. As noted in Section VI, the largest systematic errors, for example, the uncertainty in the background shape, are determined using statistically limited data sets and additional exposure is expected to reduce these uncertainties. The observed event rates are consistent with the expectation and indicate that the event rate from neutral current π 0 production is not underestimated. This provides confidence that the neutral current π 0 background to electron neutrino appearance in T2K is not underestimated.