Measurements of $\pi^{\pm}$ differential yields from the surface of the T2K replica target for incoming 31 GeV/c protons with the NA61/SHINE spectrometer at the CERN SPS

Measurements of particle emission from a replica of the T2K 90 cm-long carbon target were performed in the NA61/SHINE experiment at CERN SPS, using data collected during a high-statistics run in 2009. An efficient use of the long-target measurements for neutrino flux predictions in T2K requires dedicated reconstruction and analysis techniques. Fully-corrected differential yields of $\pi^\pm$-mesons from the surface of the T2K replica target for incoming 31 GeV/c protons are presented. A possible strategy to implement these results into the T2K neutrino beam predictions is discussed and the propagation of the uncertainties of these results to the final neutrino flux is performed.


Introduction
The NA61/SHINE (SPS Heavy Ion and Neutrino Experiment) experiment [1] at the CERN Super Proton Synchrotron (SPS) is pursuing a rich physics program in various fields. Precise hadron production measurements are performed for the T2K long-baseline neutrino experiment [2][3][4][5] and for more reliable simulations of cosmic-ray air showers for the Pierre Auger and KASCADE experiments [6,7]. The properties of the onset of deconfinement are studied with measurements of p+p [8], p+Pb and nucleus+nucleus collisions at the SPS energies [9,10].
In T2K, a high-intensity neutrino beam is produced at J-PARC by a 30 GeV (kinetic energy) proton beam impinging on a 90 cm long graphite target (1.9 nuclear interaction length, λ I ). Positively or negatively charged hadrons exiting the target (mainly π and K mesons) are focused by a set of three magnetic horns and decay along a 96 m long decay tunnel. The flavour content and energy spectrum of the neutrino interactions are measured at the near detector complex located 280 m away from the target station, and by the Super-Kamiokande (SK) detector at a distance of 295 km.
Accelerator-based neutrino beams provide well defined and controlled sources of neutrinos. However, intrinsic uncertainties on the fluxes predicted with Monte Carlo (MC) simulations arise from the models employed to simulate hadron emission from long nuclear targets used in accelerator-based experiments. In these types of experiments, a non-negligible fraction of the neutrino flux actually originates from particles which are produced in hadronic reinteractions in the long target. Up to now, neutrino flux predictions a present address: Department of Physics, COMSATS Institute of Information Technology, Islamabad 44000 Pakistan have been constrained by using either parametrizations based on existing hadron production data available in the literature, or dedicated hadron production measurements performed on thin nuclear targets.
This approach was also followed by the T2K experiment for (anti)neutrino flux predictions [11] with the NA61/SHINE measurements [12][13][14][15] performed using a thin (0.04 λ I ) graphite target and a 31 GeV/c proton beam. For the first time, the kinematical phase space of pions and kaons exiting the target and producing neutrinos in the direction of the near and far detectors was fully covered by a single hadron production experiment.
Such hadron production measurements performed with thin targets provide constraints on the production of secondary particles in the primary interaction of beam protons in the target. However, the lack of direct measurements of the production of tertiary particles in re-interactions, and hence the use of sparse data sets to cover these contributions, limits the achievable precision of the flux prediction. The main motivation for measurements of hadron emission from a replica of the T2K target is therefore to reduce the systematic uncertainties on the prediction of the initial neutrino flux originating from products of interactions in the target.
Thus, in addition to the published charged pion and kaon measurements in p+C interactions at 31 GeV/c on a thin target [12][13][14], already used for the T2K neutrino flux predictions [11], the NA61/ SHINE collaboration also performed studies of hadron emission in interactions of a 31 GeV/c proton beam with a full-size replica of the T2K target using data taken in 2007 [16]. New dedicated reconstruction and analysis techniques were developed and described. Uncorrected differential yields of positively charged pions at the surface of the replica target and their ratios with respect to the predictions of the model used to simulate hadronic interactions in the T2K target were published [16].
In this article we present new measurements of fully-corrected differential yields of π ± -mesons from the surface of the T2K replica target for incoming 31 GeV/c protons, denoted as p+(T2K RT)@31 GeV/c The results were obtained using data collected during a high-statistics run performed in 2009. The paper is organized as follows: Section 2 describes the T2K neutrino beam and discusses different sources of neutrinos. Section 3 is devoted to the NA61/SHINE experimental setup. Section 4 presents the analysis techniques used for the T2K replica target measurements, followed by Section 5 which describes the corresponding systematic uncertainties. Section 6 shows and discusses the results on the fully-corrected differential yields of π ± -mesons from the surface of the T2K replica target. A possible implementation of these results in the T2K neutrino flux predictions and the expected constraints on these predictions with the T2K replica target measurements are explained in Section 7.

T2K neutrino beams
The T2K neutrino beam is generated at the J-PARC complex by 30 GeV protons impinging on a target which is a 90 cm long graphite rod. The primary proton beam is monitored by a set of detectors which allows to precisely measure the characteristics of the beam. The produced hadrons are focused by magnetic horns.  [11] for the positive focussing at 250 kA horn current ('p250ka' configuration).
By choosing the polarity of the horn currents, it is possible to create either an enhanced neutrino beam or an enhanced antineutrino beam. In this article we concentrate on the case of the enhanced neutrino beam but the results of this paper can also be used for the prediction of the flux in the enhanced anti-neutrino configuration. A detailed description of the beam and its properties can be found in Ref. [11].
The neutrino beam predictions are based on a detailed Monte-Carlo simulation. The input parameters are given by the values describing the ellipsoid representing the primary proton beam impact points on the target upstream face as measured by the beam position detectors placed along the proton beam line. The FLUKA2011 [17][18][19] model is used to simulate the interactions of beam protons with the long graphite target. The propagation of the particles emerging from the surface of the target is modeled by a Z6 Z1 Z2 Z3 Z4 Z5 During the MC simulation, information on particle production and decay is stored, so the full history of neutrinos crossing either the near or far detector is available. This allows to study different components of the neutrino beam and the origin of the neutrino species. As shown in Fig. 1, the ν µ flux around the beam peak energy at the SK far detector arises mainly from pion decays, while it is mainly due to kaons at higher energies. This motivates the extraction of charged pion yields at the surface of the target, which is the subject of this paper.
It is important to note that not only the pion angular and momentum spectra are of interest, but also the longitudinal position where they exit the target. By dividing the 90 cm long graphite rod into 5 bins of 18 cm length each and considering the downstream face of the target as an additional sixth bin, as shown in Fig. 2, it is possible to study the contribution of each of these bins to the total neutrino flux. Figure 3 presents these different contributions as predicted at SK.

Requirements on the T2K neutrino flux prediction
The T2K experiment pursues three main physics goals [2] with an off-axis (essentially narrow band) neutrino or antineutrino beam peaked around the so-called atmospheric oscillation maximum (energy range from 0.2 to 1.2 GeV). These are: (i) the muon neutrino disappearance, (ii) the electron neutrino appearance (ν µ → ν e ), (iii) neutrino cross section measurements.
The muon-neutrino flux in the region of interest for the oscillation analysis is mainly generated by pion decays. For the oscillation measurements, the ratio of the flux of neutrinos at the near detector to the one at the far detector is the most important quantity and a desirable level of uncertainty is about 1-2%. Another quantity of interest for the electron neutrino appearance is the ratio between electron and muon neutrino cross sections, whose measurement in the near detector will require a knowledge of the electron to muon neutrino fluxes to better than about 2%. Failing to match this required precision might limit the precision of the results for the full expected T2K exposure. Existing data on (anti)neutrino cross sections in the energy range of interest are very limited, the precision of measurements ranging typically between 10% and 20%. A precision on the T2K neutrino flux with a 5% absolute normalization error would allow considerable improvement in the understanding of low energy neutrino interactions.
2.3 T2K flux predictions and the T2K replica target measurements in NA61/SHINE As already described in Ref. [16], the neutrino fluxes can be split into secondary and tertiary components. The secondary component originates from neutrino parents produced in the primary interaction of the beam protons in the target. The tertiary component refers to neutrino parents produced in interactions of secondary particles. The latter component is due to re-interactions in the target and re-interactions taking place in the elements of the beamline. Secondary and tertiary interactions occurring in the target are constrained by the measurements of identified hadron spectra from the surface of the T2K replica target.
The ν µ and ν e spectra around the most probable neutrino energy in T2K are predominantly produced by pions (see Fig. 1). Figure 4 shows the phase space (in the kinematic variables p and θ -the momentum and polar angle of particles in the laboratory frame) of the pions exiting the target surface and contributing to the ν µ flux at SK. The binning of the T2K replica target analysis is overlaid. The bins in polar angle and in z along the target were defined to ensure adequate sampling of the T2K beam focusing efficiency. The binning in momentum was then chosen to obtain roughly equipopulated bins. As can be seen, the T2K replica target analysis region covers most of the phase space of interest for T2K. Figure 5 presents the fractions of the ν µ and ν e fluxes at SK that can be constrained by the T2K replica target measurements presented in this article. The remaining flux originates from particles produced in interactions of primary protons or secondaries with the beam line elements, or by other particle species such as kaon decays, which are not included in the present analysis.

NA61/SHINE experimental setup
The NA61/SHINE apparatus is a wide acceptance spectrometer at the CERN SPS. Most of the detector components were inherited from the NA49 experiment and are described in Refs [1,22]. A more detailed analysis-oriented description of the NA61/SHINE setup can also be found in Ref. [12]-Sec. II. Only some features relevant for the 2009 running period are briefly mentioned here. The general layout of the detector is displayed in Fig. 6. The NA61/ SHINE right-handed coordinate system is displayed in the figure with the z axis along the beam line, the x axis in the horizontal plane and the y axis pointing upwards. The origin of the coordinate system is placed in the center of the VTPC-2.  Fig. 6: The NA61/SHINE experimental setup (horizontal cut). The beam is coming from the left, impinging on the T2K replica target shown in this figure. The chosen coordinate system is as follows: its origin lies in the middle of the VTPC-2, on the beam axis. The nominal beam direction is along the z axis. The magnetic field bends charged particle trajectories in the x-z (horizontal) plane. Positively charged particles are bent towards the top of the plot. The drift direction in the TPCs is along the y (vertical) axis.
The spectrometer is built around five Time Projection Chambers (TPCs): two Vertex TPCs (VTPC-1 and VTPC-2) placed in the magnetic field produced by two superconducting dipole magnets and two Main-TPCs (MTPC) located downstream symmetrically with respect to the beam line. A small additional TPC is placed between VTPC-1 and VTPC-2, covering the very-forward region, and is referred to as the GAP-TPC (GTPC).
The experimental setup is complemented by time-of-flight detectors. The ToF-F is located in the forward region, downstream of the MTPCs. It was used in the analysis presented in this paper. The detector consists of 80 scintillator bars read out at both ends by photo-multipliers. The time resolution of each scintillator is about 115 ps [23].
For the study presented here the magnetic field was set to a bending power of 1.14 Tm. This leads to a momentum resolution σ (p)/p 2 in the track reconstruction of about 5 × 10 −3 (GeV/c) −1 .
The T2K replica target is an isotropic graphite rod of density 1.83 g/cm 3 with a thickness along the beam axis of 90 cm, equivalent to about 1.9 λ I and a radius of 1.3 cm. Detailed descriptions of the minor differences between the target mounted in the T2K beam line and the T2K replica target are given in Ref. [16]. The downstream face of the target was placed 52 cm upstream of the VTPC-1. Aluminum flanges were mounted at the most upstream part of the target in order to fix it in the NA61/SHINE experimental set-up. The rod and the flanges can be seen in Fig. 2. Alignment screws, specially added for the 2009 data taking period and mounted on the flanges, allowed to precisely place and align the target with respect to the beamline.

Analysis
The analysis described in this paper is based on 2.8×10 6 reconstructed events collected during the 2009 data-taking period.

NA61/SHINE beam line
The NA61/SHINE beam line with different counters and positioning detectors used in the 2009 T2K replica target run is presented in Fig. 7. The beam line is instrumented with counters, S1-S3, and veto counters, V0 and V1 , which provide the beam definition. Furthermore S1 also sets the start time for all other detectors. The S3 scintillator counter has the same radius as the target and is positioned 5 mm upstream of the target upstream face. Hence, it ensures that the incident beam proton is impinging on the target. Two Cherenkov counters, CEDAR and THC, allow to trigger on protons as beam particles. The trigger used in the analysis includes the following combination of counters: (1) With the measurements of the three beam position detectors (BPDs), tracks of beam particles are reconstructed by fitting two straight lines in the x − z and y − z planes. Least square fits use the measured charge clusters in the BPDs.
The distribution of outgoing hadrons along the T2K replica target depends on the impact point of a primary proton at the target upstream face. This point is reconstructed using information from the BPD detectors. To ensure the quality of the beam track strict cuts on the BPD measurements are applied. All three BPDs must have properly reconstructed clusters in the x and y coordinates. Furthermore, a cut on the χ 2 of the fit of the reconstructed tracks allows to reach a resolution of 300 µm at the BPD-3 which is positioned 7 cm upstream of the target front face. Figure 8 shows the beam profile on the target upstream face as reconstructed with the BPD information.

S1
S2 VTPC1 V0 Fig. 7: A schematic layout of the NA61/SHINE beam line instrumentation [1] for the T2K replica target configuration. Horizontal cut in the beam plane is not to scale.

Reconstruction of particle parameters at the target surface
The final neutrino flux depends on the longitudinal distribution of the hadrons exiting the target surface. Hence, the exit position of the produced particles and their momentum vectors at the target surface were reconstructed.
The reconstruction procedure was similar to the one described in Ref. [15]. First, charge deposits by tracks traversing the TPCs measured by neighbouring readout pads are joined to form clusters. Local track segments are then reconstructed from the clusters in each TPC separately. The matching of the track segments from different TPCs allows to reconstruct global tracks. The track fitting through the magnetic field allows the determination of the track parameters at the first measured TPC cluster.
In order to finally determine the track parameters at the surface of the target, a backward extrapolation is performed from the first measured TPC cluster. A Runge-Kutta method is used to propagate track parameters and their uncertainties in the non-uniform magnetic field. An exit point of the track from the target surface is found when the backward extrapolated track intersects the target volume. If no intersection point can be found, the point of closest approach between the backward extrapolated track and the surface of the target is recorded and assigned as the exit position. The track is then considered as originating from the target if the distance between the point of closest approach and the target surface is within the one standard deviation uncertainty. Most particles point to the inside of the target, both in data and Monte-Carlo. The cut at 1 standard deviation keeps 96.6% of the selected tracks in the data and 96.5% in the Monte-Carlo. The resulting systematic uncertainty is considered negligible.

Determination of the target position and tilt
A precise knowledge of the position of the target and its alignment along the beam line is mandatory in order to reconstruct the exit point of the produced hadrons at the surface of the target. The alignment of the target with respect to the experimental setup is determined in three consecutive steps based on the data.
The first step consists of determining the angular alignment of the target. This is done by dividing the 90 cm long target into 17 slices. The reconstructed tracks are extrapolated backward from the TPCs to the center of each of the 17 slices. For each slice, a mean of the x (respectively y) position of the backward extrapolated tracks is determined and this mean is assigned as the central x (respectively y) position of the target. Fitting all the mean positions with a straight line allows to determine a tilt in the x − z (respectively y − z) planes. For the 2009 T2K replica target dataset, these tilts were found to be negligible as expected from the target positioning precision of better than 2.2 mrad.
The second step consists of determining the transverse position of the target with respect to the beam line. This position is determined with the help of the reconstructed beam tracks from the BPDs. Requiring a hit in the S3 scintillator counter ensures that the incident protons have reached the target upstream face and by drawing the x − y distribution of the beam particles under this requirement, the center of the target can be determined as the mean of the x − y distribution.
The third step consists of cross-checking the alignment between the beam line and the spectrometer and subsequently extracting the longitudinal position of the target. As this position is determined by using the beam tracks reconstructed from the BPD information, and as the position of the particles exiting the target surface is determined by using tracks from the TPCs, it is important to know the precision of the alignment between the BPDs and the TPCs. This is done by checking the consistency between the reconstructed vertices from the beam tracks and the reconstructed tracks from the TPCs for two independent track topologies. These topologies are defined as Then a potential misalignment in the transverse x − y plane can be reconstructed from a shift in the distribution of the longitudinal z coordinate. Figure 9 shows a schematic of this procedure. To reconstruct coordinates of the primary interaction point inside the long target, specific cuts have to be applied in order to reduce the number of tracks that originate from re-interactions inside the target and hence have trajectories not intersecting with the beam particle. The following cuts are applied to the events used for the construction of the distribution of the vertices: (i) the beam particle hits the target at least 0.5 cm inside the target surface (ii) the associated track in the TPCs has to be on the same side in the y − z plane as the beam particle, i.e. if the beam particle hits the target upstream face at the positive x, then the tracks have to exit the target surface at positive x as well.
In order to get sufficient precision only tracks with 100 < θ < 180 mrad were used. The procedure was tested and validated by applying it to a Monte-Carlo simulation for which the z position of the target is well known and the alignment of the beam is perfect with respect to the spectrometer.   There are several analysis techniques applied in NA61/SHINE for the extraction of raw yields of charged pions (see Refs. [12,15]). The analysis presented in this article, utilizes the TPC measurements of the specific energy loss (dE/dx) and the time-offlight measurements (to f ) using the ToF-F detector.
In this section, we will concentrate on specific details related to the extraction of charged pion yields at the surface of the T2K replica target.

Event and track selection
As mentioned in Section 4.1, the events of interest for this analysis are selected based on two requirements: (i) a hit in the S3 counter is required to make sure that the incoming proton is hitting the target upstream face. (ii) a strict cut on the χ 2 of the fitted beam track.
After applying these two cuts, a sample of 1.6×10 6 events remains available for the analysis.
The requirements on the reconstructed tracks depend on the analysis technique to be used (as explained later) but they are based on the same criteria: These criteria were also used for the target alignment procedure described in Section 4.3 as they select tracks with well fitted parameters.

Data binning
As mentioned in Section 2.1, different longitudinal sections of the target contribute differently to the final neutrino flux while the focusing of the horns will affect the particles depending on their momentum and polar angles. Hence, the analysis of the T2K replica target will be conducted in (p, θ , z) bins. The longitudinal z binning was determined by a study performed together with the T2K beam group. It was found that five longitudinal bins are sufficient to obtain a neutrino flux prediction that matches the nonbinned case, both in terms of shape and overall normalization, within a known and correctable bias of less than 2%. Hence the target surface is divided into five bins of 18 cm length and the downstream face of the target is taken as a sixth longitudinal bin, as shown in Fig. 2. The chosen (p, θ ) binning scheme is illustrated in Fig. 4.

The to f − dE/dx analysis for particle identification
Particle identification (PID) in NA61/SHINE relies on measurements of the energy loss dE/dx in the TPCs and the time-offlight that is used to compute the particle mass squared, m 2 . The method is illustrated in Fig. 11 (top panel) which depicts how the different particles (p, K, π and e) can be separated in the (m 2 , dE/dx) plane. A (m 2 , dE/dx) distribution, separately for positively and negatively charged tracks, is obtained for each bin (p, θ , z) determined at the surface of the replica target. The data distributions are then fitted to joint probability density functions (pdf) for the mass squared and the energy loss. Due to the independence of the dE/dx and m 2 variables, the joint pdf reduces to the product of the corresponding marginal distributions which are described by Gaussian distributions. The complete pdf is a sum of two-dimensional Gaussian distributions of four particle species, p, K, π and e. For the initialization of the fit, the resolution of the mass squared and the expected energy loss for each particle species is obtained from parametrizations of the data distributions shown in Figs. 12 and 13 as a function of the track momentum. The resolution of the energy loss measurement is a function of the number of reconstructed clusters on the track (1/ √ N). For the topology dependent cuts defined in this analysis the dE/dx resolution can be approximated by a constant value of 3% due to the sufficiently large number of clusters on each track. Independent normalization factors are introduced for each particle species. Since the individual pdfs are normalized to unity, particle yields are given by the normalization factors which are obtained from a two-dimensional log-likelihood minimization illustrated in Fig. 11. The projections on the m 2 and dE/dx variables better illustrate the quality of the fit results (see Fig. 11, bottom panels).

Simulation
Simulations were performed to generate events of 31 GeV/c protons interacting with the T2K replica target. In order to be consistent with the T2K neutrino beam simulation program [11], the simulation package FLUKA2011 [17][18][19] is used in NA61/SHINE to generate the interactions inside the graphite target. The GEANT3 [20] transport code was used to track the particles through the detector and GCALOR [21] handled the interactions in the spectrometer. More details can be found in Ref. [24]. The full simulation chain consists of three parts: (i) First, a stand-alone FLUKA simulation. The target geometry is described as a 90 cm long graphite rod with aluminum flanges and the S3 counter. The target was positioned at the location determined by the alignment procedure applied to the data as explained in Section 4.3. The incident proton beam profile was simulated following the shape of the distributions for the positions and divergences given by the data. The momentum of the beam was set precisely at 30.92 GeV/c to match the beam momentum measured in the data. Information on interactions happening inside the target was stored. Position, momentum as well as polar and azimuthal angles of the particles exiting the target were recorded at the surface of the target and saved as output of the FLUKA simulation. (ii) A GEANT3-based program used the kinematic parameters of particles produced by FLUKA at the surface of the target and propagated them through the NA61/SHINE experimental setup. The GCALOR model handled all hadronic interactions in the spectrometer. Moreover, a detailed simulation of various detector effects was included. (iii) The tracks were finally reconstructed following the same reconstruction procedure as the one applied to the data. All information from the FLUKA simulation and the simulated tracks until their reconstruction was stored in the final output files. This allows to get the full history of the simulation and to match the reconstructed to simulated tracks.
The simulation was used to calculate corrections for pions resulting from various sources: i) weak decay of heavier particles producing additional pions, ii) interactions in the detector material, iii) track reconstruction efficiency and resolution, iv) decay in flight. In total, 10 millions of protons on target were simulated.

Systematic Uncertainties
Six different sources of systematic uncertainties are considered. Their contributions to the systematic uncertainty are described in detail below. Five of them are similar to the thin target p+C@31 GeV/c analysis [15], only the last one (backward extrapolation) is specific to the T2K replica target analysis.

Particle Identification
The dE/dx distribution for the different particle species in each of the (p, θ , z) bins is approximated by a single Gaussian. In order to estimate the uncertainty related to this approximation two Gaussians with the same mean value but different widths were used to fit the dE/dx distributions. At low momenta, the particle identification is constrained by the ToF-F information and hence the magnitude of the uncertainty due to describing the energy loss by a single Gaussian is expected to be negligible. At higher momenta, when the resolution of the time-of-flight measurements does not allow to distinguish the different particle species, using two Gaussians instead of a single Gaussian in the fitting procedure leads to differences of up to 2% at momenta higher than 10 GeV/c.

Feed-down corrections
Pions not originating from the target but traversing the spectrometer and reconstructed as exiting the target surface represent the so-called feed-down contribution. The feed-down correction comes from particles of various origins: (i) interactions of particles outside the target, (ii) decays in flight of strange particles. The correction factor for the feed-down contribution is computed based on simulations produced with FLUKA as primary hadronic generator. This correction is model dependent and an uncertainty on this model prediction has to be assigned. As for the thin target analysis, 30% of the correction was assumed as the systematic uncertainty of the correction [12,15]. This uncertainty reaches values as large as 5% of the pion yield for momenta lower than 2 GeV/c. It decreases significantly at higher momenta.

Reconstruction efficiency
Following the thin target analysis, a constant 2% uncertainty on the efficiency of the reconstruction procedure is assigned [15].

ToF-F reconstruction efficiency
The correction for the ToF-F reconstruction efficiency is computed based on the procedure described in Ref. [25]. The efficiency was estimated on a sample of physics events with a strict cut on a time window around the triggered interaction. In the procedure a track traversing the ToF-F geometrical acceptance was required to make a hit in the corresponding slab that can be used later to compute a value of m 2 of the particle. The efficiency was parametrized as a function of slab position with respect to the beam and the track momentum. This dependence is small but not negligible. A constant 2% over the entire phase space is hence assigned as the systematic uncertainty of the ToF-F reconstruction efficiency.

π-loss
As mentioned above, the loss of pions can be regarded as tracks being measured in the TPCs and aiming towards the ToF-F acceptance but not having a recorded hit in the ToF-F due to decay and due to absorption or interactions of pions with the detector. The corrections related to the decay are computed via the precisely known pion decay which should be model independent. Hence, when varying the number of requested measured points in the MTPCs, one does not expect to see differences in the final spectra. Any variations would represent an uncertainty due to imperfections in the description of the spectrometer which can introduce a difference in the acceptance and reconstruction efficiency (merging track segments between VTPC-2 and MTPC-L/R) between simulated and real data. This uncertainty decreases fast with increasing particle momentum. Below 2 GeV/c this contribution can be larger than 5% but usually is not larger than 1% at higher momenta.

Backward track extrapolation
The uncertainty due to the backward extrapolation procedure is induced by the uncertainty on relative position of the target and TPCs, as presented in Section 4.3. The main goal of the backward extrapolation is to attribute to each track a specific longitudinal z bin as well as to determine the momentum and polar angle at the surface of the target. By shifting the target within the uncertainties on the different coordinates, the number of tracks exiting from each different (p, θ , z) bin will vary. This variation is used as binby-bin systematic uncertainty on the final spectra due to the backward extrapolation. This contribution is the most important one for the most upstream z bin at low polar angle, and for the most downstream bin at high polar angles. It can range up to 10% in these two specific phase-space regions.

Summary of systematic uncertainties
The systematic uncertainties are presented in Figs. 14 and 15 for positively charged pions and in Figs. 16 and 17 for negatively charged pions. They are displayed in z and θ bins as a function of momentum. The numerical values can be found in Ref. [25].
The momentum and angular resolutions are significantly smaller than the bin sizes, and bin-to-bin correlations are very small. The only significant bin-to-bin correlation concerns the first and second longitudinal bins along the target. These correlations were studied in the context of the systematic errors on alignment and backward extrapolation.
The dominant contribution for the most upstream and downstream z bins is due to the backward extrapolation. This can be understood by the fact that an uncertainty of the target position gets translated into a fake track migration between the longitudinal bins at the reconstruction level of the exit point of the tracks from the target surface. For the central part of the target (longitudinal bins z1 to z5) this effect gets averaged between the bins. At low momenta, the uncertainties of the feed-down corrections as well as of the pion loss are significant. In this region the systematic uncertainty due to particle identification is quite small. This can be explained by the fact that the to f − dE/dx approach separates well the particle species, as seen in Fig. 11.

Statistical uncertainties
Statistical uncertainties of the final corrected pion spectra receive contributions from the finite statistics of the data as well as the simulated events used to obtain the correction factors. The dominating contribution is the uncertainty of the normalisation factors returned by the maximum likelihood method applied for PID. The simulation statistics were higher than the data statistics by more than a factor of 6, hence the related uncertainties are by about a factor of 3 smaller. Added in quadrature to the statistical uncertainties of the data, these errors become negligible.

Cross-check with the h − analysis technique
The primary results obtained with the to f − dE/dx analysis technique were cross-checked using a complementary h − analysis method which was presented in previous publications [12,15].
The π − differential yields computed with the h − technique are, in general, in good agreement with those obtained with the to f − dE/dx approach. Some small deviations appear only in the upstream part of the target for the polar angle bin 0 < θ < 20 mrad and momentum below 7 GeV/c, where the π − differential yields from the h − analysis exceed those obtained with the to f − dE/dx approach by a few percent, actually remaining at the edge of the estimated uncertainties. In the other p − θ bins very good agreement between the two analysis techniques is observed. More details on the comparison of π − differential yields can be found in Ref. [10].
The agreement between the h − and to f − dE/dx analysis for the extracted π − yields gives confidence in the high precision of the results and in the reliability of the estimate of systematic uncertainties.
6 Results and comparison with the hadron production model used by T2K Differential multiplicities of positively and negatively charged pions emitted from the T2K replica target exposed to a 31 GeV/c proton beam are presented in Figs. 18 to 21 and in tabular form in Ref. [26]. The normalization is done to the number of protons on target.
To compare the replica target results with model predictions and with the results of thin target measurements [15], a procedure developed by the T2K beam group is used. A complete description of this procedure can be found in Ref. [11]. This procedure is based on the re-weighting of the model predictions with experimental measurements 1 . While running a Monte-Carlo simulation of the T2K neutrino flux predictions, the interactions of the incident protons with the 90 cm long graphite target and all subsequent processes leading to the creation of a neutrino are recorded. By considering only the particles exiting the target surface and reweighting only the interactions occurring inside the target, it is possible to directly use this procedure to compare model predictions with the T2K replica target measurements. Furthermore, the interactions taking place inside the target can be constrained by the available hadron production data. Those data are predominantly the NA61/SHINE measurements. Other external data are used as well as extrapolations to lower incident hadron energies [11] in order to constrain re-interactions inside the target. The constraint is done by re-weighting the produced Monte-Carlo particles based on two physics quantities: 1 Re-weighting is performed by assigning a multiplicative positive real number to each simulated particle, and if appropriate its interaction or decay products, in such a way that its entry in histograms, probability functions and flux calculations is multiplied by the product of all the weights it has received during the simulation. The weights are calculated in such a way that the Monte-Carlo simulation reproduces an input data set.
(i) the hadron differential multiplicities: at each interaction, the predicted hadron spectra normalized to mean multiplicities can be re-weighted with respect to the corresponding spectra measured by NA61/SHINE in p+C interactions at 31 GeV/c. (ii) the production interaction rate: at each interaction, the rate at which hadrons interact and produce new particles can be re-weighted to the production cross section extracted from the NA61/SHINE thin target measurements [15]; this production cross section is defined as σ prod = σ inel − σ qe , where σ inel is the inelastic cross section and σ qe is the quasi-elastic cross section which characterizes the break up of the carbon nucleus without production of new hadrons.
The weights are computed as the ratio between the measured values and the model predictions. The procedure developed by the T2K beam group allows to apply the two above-mentioned weighting schemes independently and hence to evaluate the effect of each of them.
In Figures 22 to 27 three comparisons are presented: (i) the T2K replica target results are compared with the nominal FLUKA [17][18][19] predictions. (ii) the T2K replica target results are compared with the FLUKA predictions re-weighted for the hadron multiplicities [15]. As can be seen, the weighting of the hadron multiplicities has a small effect which in general does not exceed a few percent. This is expected as the FLUKA predictions reproduce rather well the differential cross sections of pions measured by NA61/SHINE in p+C interactions at 31 GeV/c. (iii) the T2K replica target results are compared with the FLUKA predictions which are re-weighted for hadron multiplicities and production cross section [15]. In order to illustrate the sensitivity of the model predictions to the production cross section, a re-weighting procedure was applied such that it effectively decreases by 20 mb the production cross section in FLUKA 2 . As shown in Figs 22-27, this has a visible effect on the pion yields along the target. For the upstream part of the replica target and at lower pion momentum, the predictions are thus lowered by 5-10 %, bringing them closer to the replica target measurements. At higher momenta the predicted spectra stay unchanged. For the downstream part of the target, the agreement lies always within the uncertainty of the replica target data points.
The model whereby total cross sections are described by a sum of elastic scattering, inelastic scattering and production cross section, on which FLUKA is based, does not seem to reproduce well the longitudinal distribution of particles along the replica target. This point justifies further studies, first with the higher statistics data taken in 2010, and possibly with an upgraded experimental set-up in the future.   The implementation of the T2K replica target results within the T2K neutrino flux prediction can essentially follow the one for the thin target re-weighting procedure [11]. The major modification to be introduced is to save the relevant parameters of the particles exiting the surface of the target in the simulation process. This specific information is added to the already existing history chain for each produced neutrino. Recording these parameters allows, once the full neutrino beam simulation is completed, to come back to each pion exiting the surface of the target and being a parent of neutrinos. The neutrinos are then assigned a weight which is dependent on the (p, θ , z) parameters of the pions at the surface of the target and computed as where N meas NA61 are the hadron yields at the surface of the target measured by NA61/SHINE and N sim T 2K are the yields of emitted particles  at the surface of the target calculated within the T2K beam simulation program. In both cases the yields are normalized to the number of incoming protons.
Other particle species exiting the target surface are re-weighted following the already established thin target procedure. Re-interactions outside of the target (in the focusing horns and along the beam line) are also tuned with the nominal T2K neutrino flux calculation based on the thin target data.
In order to compare, in a consistent way, the neutrino flux prediction re-weighted with the thin target procedure and with the T2K replica target results, it is important to use the same settings for the entire flux computation. As the weights related to the T2K replica target dataset were computed with the NA61/SHINE beam profile, this same beam profile is used to run a full neutrino flux simulation within the T2K neutrino flux simulation code [11]. Using this code a comparison between the thin target and T2K replica target re-weighting procedures can then be performed.
In order to study different effects of the flux tuning with the thin target procedure, the neutrino flux at SK is first re-weighted for the multiplicity only and then for the multiplicity and production cross section for proton-Carbon interactions. For the latter, as in the case of the pion spectra, the production cross section is reduced by 20 mb as described in Section 6.
The ratio of ν µ flux re-weighted with the thin target procedure over the T2K replica target flux at SK is shown in Fig. 28. At low momentum the thin target procedure predicts a 10% larger value of the flux which is however consistent with the replica target prediction within one standard deviation. At high neutrino energy, the ratio tends to one. This is due to the fact that within this energy range the dominant part of the ν µ spectra comes from kaons and thus is not affected by the re-weighting procedure on pion spectra. At lower energy, where the pion parent particles dominate, the thin target tuning predicts a higher neutrino flux at SK. By lowering the proton-carbon production cross section, reduced spectra are   The beam profile of NA61/SHINE (space and angular distribution of incoming protons) were used for the neutrino flux calculations presented above. The profile of the T2K beam is about 50% narrower and better centered. Given the knowledge, on an eventby-event basis, of the incoming proton parameters, it is possible to re-weight the NA61/SHINE data to reproduce more precisely the situation of the T2K beam. Such an analysis was performed in Ref. [25]: it showed only a minor impact on the resulting neutrino fluxes that is small compared with the systematic errors.   propagation of uncertainties to the systematic uncertainties of the neutrino flux. For each neutrino energy bin i, one can write the contribution of the re-weighting factors, for each (p, θ , z) bin j, as a linear combination and the coefficients a i j are related to the contribution of each (p, θ , z) bin to each neutrino energy bin. Propagating the T2K replica target uncertainties means propagating the uncertainties on the weights ω j . The propagation is done via computation of covariance matrices for each source of uncertainties The off-diagonal elements σ x i ,x j , i = j of the matrix C X are the correlation coefficients for the T2K replica target results. The diagonal elements of the matrix C Y give the uncertainty on the final neutrino flux predictions with respect to the propagation of the replica target uncertainties. For instance, the first diagonal element of the matrix C Y can be written with the notation, Y 1 = f 1 (X), as: For each component of the systematic uncertainties related to the T2K replica target results, a full correlation between the anal-   ysis bins is considered, but no correlation is taken into account between the different contributions to the systematic uncertainties.
The statistical uncertainties are treated as uncorrelated between the different analysis bins as well as uncorrelated with the different components of the systematic uncertainties. Figure 29 shows the result of the propagation of the systematic and statistical uncertainties to the ν µ flux at SK. Each line for the systematic uncertainties corresponds to a component described in Section 5. It is important to note that these uncertainties correspond only to the component of the ν µ flux produced by pions exiting the target surface. Hence it does not represent the full uncertainties. As presented in Fig. 5, it covers about 87% of the flux at the most probable neutrino energy but only 10% at 4 GeV. Hence, with the T2K replica target results around 87% of the ν µ flux at SK can be predicted with 3.5% uncertainty at the most probable energy, while at 4 GeV only 10% of the flux can be predicted with a 4% uncertainty. The uncertainty on the remaining part of the flux will have to be estimated from the production of kaons off the surface of the target and/or the reinteractions along the beam line. These estimations are out of the scope of this paper.
The fraction of the SK neutrino flux that is coming from pions originating from the long target surface becomes very small at low    energies below 400 MeV, this explains the result that the relative flux uncertainty related to measurements made in NA61/SHINE is small. Probably the flux at these low energies originates from muon decays and tertiary pions.

Conclusions
The data taken in 2009 with 30 GeV protons impinging on a 90 cm long, 2.6 cm diameter carbon rod, were analysed and fullycorrected π + and π − spectra at the surface of the target were obtained. Using the beamline simulation program of T2K, these spectra were compared with the FLUKA2011 prediction re-weighted with cross section measurements obtained by NA61/SHINE with a thin target. A reasonable agreement was found, although an even better description was obtained when lowering the production cross section by 9%. A method allowing the direct implementation of the T2K replica target results, as well as the propagation of their uncertainties, in the T2K beam simulation was demonstrated. Fur-ther results will be obtained from a higher statistics dataset taken by NA61/SHINE in 2010.                  , FLUKA re-weighted for the multiplicities (green) and FLUKA re-weighted for multiplicities and production cross section σ prod (magenta) for the three upstream longitudinal bins, Z1-Z3, and in the polar angles between 220 and 340 mrad plotted as a function of momentum.    Z 5 dn/dp   Fig. 28: Ratio of the ν µ flux at SK re-weighted with the thin target procedure to the p+(T2K RT)@31 GeV/c flux. For the latter pion spectra presented in this paper have been used in the re-weighting. The dominant part of the ν µ spectra at high E ν comes from kaons and thus is not affected by the re-weighting of the pion spectra. The left plot shows the ratio calculated using the FLUKA production cross section, whereas the FLUKA cross section reduced by 20 mb was used to obtain the ratio presented in the right plot. Vertical error bars show the full uncertainties on the ratio which are dominated by systematical uncertainties.