Measurements of $$\pi ^\pm $$ π ± , $$K^\pm $$ K ± and proton double 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 the π ±    π±  , K ±    K±  , and proton double differential yields emitted from the surface of the 90-cm   cm  -long carbon target (T2K replica) were performed for the incoming 31   31   GeV/c   GeV/c   protons with the NA61/SHINE spectrometer at the CERN SPS using data collected during 2010 run. The double differential π ±    π±   yields were measured with increased precision compared to the previously published NA61/SHINE results, while the K ±    K±   and proton yields were obtained for the first time. A strategy for dealing with the dependence of the results on the incoming proton beam profile is proposed. The purpose of these measurements is to reduce significantly the (anti)neutrino flux uncertainty in the T2K long-baseline neutrino experiment by constraining the production of (anti)neutrino ancestors coming from the T2K target.

section_content: NA61/SHINE (SPS Heavy Ion and Neutrino Experiment) [1] is a large hadron spectrometer at the CERN Super We dedicate this article to the memory of our colleague Professor Stepan Agaronovich Bunyatov. a e-mail: Boris. Popov@cern. ch Proton Synchrotron (SPS). The NA61/SHINE collaboration is pursuing several physics goals including hadron production measurements for T2K (Tokai to Kamioka) [2] an accelerator-based long-baseline neutrino experiment in Japan. The NA61/SHINE measurements are used to reduce the systematic uncertainties associated to the prediction of the (anti)neutrino fluxes in T2K. New measurements of π ±, K ± and p yields coming from the surface of the T2K replica target are presented here. These results aim at reducing the T2K (anti)neutrino flux uncertainties down to a 3-4% level. They can further be used to tune hadron interaction and transport models. 

The paper is structured as follows: a motivation for these hadron measurements is first presented. In Sect. 2, a brief overview of the NA61/SHINE setup is shown followed by the description of the analysis in Sect. 3. Section 4 gives a concise description of all systematic uncertainties. Results and comparison with Monte Carlo (MC) models and with the previous NA61/SHINE measurements are presented in Sect. 5. Finally, in Sect. 6, proper usage of these results in the T2K neutrino beam simulation is discussed, followed by a short conclusion.
section_name: Introduction

section_content: The T2K neutrino beam is produced at the Japan Proton Accelerator Research Complex (J-PARC) [3] by directing a 30 GeV (kinetic energy) proton beam towards a 90-cm-long graphite target. Produced mesons, mostly pions and kaons [4], are focused by a set of three magnetic horns [2] and decay to neutrinos in the decay volume. Focusing horns are aluminum conductors that create a toroidal magnetic field with respect to the beam direction. The polarity of the current in the horns can be changed, so positively or negatively charged pions (kaons) can be focused. In this way, T2K can produce either a neutrino-enhanced (ν μ ) or an antineutrino-enhanced (ν μ ) beam. Since direct hadron production measurements in situ are not possible, Monte Carlo models are used to predict the (anti)neutrino fluxes at the T2K near and far detectors. However, significant differences between hadron production models induce considerable systematic uncertainties on the (anti)neutrino fluxes. Without any experimental data to constrain the hadron production models, these systematic uncertainties would be larger than 25%. For this reason, T2K uses available hadron production measurements, but mostly relies on dedicated NA61/SHINE hadron production data [4].
section_name: The T2K neutrino beam

section_content: The NA61/SHINE experiment took data for T2K in 2007, 2009 and 2010. In 2007 and 2009, measurements were performed with a 30. 92 GeV/c proton beam and a 2-cm-thick carbon target (4% of a nuclear interaction length, λ i ) to measure hadron multiplicities (π ±, K ±, p, K 0 s, Λ) and the production cross section [5] [6] [7] [8]. Production events represent the fraction of inelastic events in which at least one new hadron in the final state is produced. By using these measurements, the uncertainty on the neutrino flux in T2K was reduced to about 10% [4]. However, the hadron production component of the uncertainty still dominated the total uncertainty due to the insufficient precision of the production cross section measurements and re-interactions inside the target and aluminum horns, which cannot be directly constrained by the measurements of the primary proton-carbon interactions. To further reduce the hadron production uncertainty of the neutrino flux, hadron production measurements with a replica of the T2K target were needed. Measurements with the T2K replica target can directly constrain up to 90% of the neutrino flux because hadron yields at the surface of the target are measured not just for primary interactions, but also for re-interactions inside the target [9]. These measurements were also done in 2007, 2009 and 2010. Measurements from the 2007 run [9] were used as a proof-ofprinciple, while measurements from the 2009 run [10, 11] are being incorporated in the T2K neutrino flux simulation. The expected (anti)neutrino flux uncertainty is around 5% [12], representing a significant improvement with respect to the previously published uncertainty of 10%. Measurements from the 2010 run are the main topic of this paper. With respect to the previous results, these measurements of π ± yields have smaller statistical and systematic uncertainties. Furthermore, K ± and p yields are measured for the first time. This is expected to further reduce the uncertainties on the (anti)neutrino fluxes in T2K from about 10% down to the level of 3-4%. Moreover, since (anti)neutrino fluxes in T2K in the energy region above ∼ 2 GeV are produced mainly by kaon decays, the corresponding uncertainties will be greatly reduced. 1. 3 The T2K requirements T2K imposes strict requirements on hadron production measurements with the T2K replica target [10]. Hadron yields must be measured as a function of the momentum, polar angle and longitudinal position of the hadron exit point on the target surface. In other words, vertices inside the target are not reconstructed. Tracks measured in a detector are only extrapolated to the target surface. Essentially, the target is treated as a black box and measured hadron yields are the sum of yields coming from primary interactions and re-interactions. In addition, T2K requires that the longitudinal position along the target surface is binned in at least six bins: five longitudinal bins 18 cm in size and the downstream target face as a sixth bin [13]. Since the T2K target is inserted in the first focusing horn, hadrons coming from different parts of the target will have different paths through the magnetic field and they will be focused differently by the horns. The previous NA61/SHINE measurements with 2007 and 2009 data have already proven that the NA61/SHINE setup is suitable for measuring hadron yields on the T2K replica target surface, although the target is placed outside of the tracking system.
section_name: The NA61/SHINE measurements for T2K

section_content: The NA61/SHINE setup during the T2K replica target datataking in 2010 consisted of five Time Projection Chambers (TPCs), three Time-of-Flight walls (ToF-F, ToF-L and ToF-R), three Beam Position Detectors (BPD-1, BPD-2, BPD-3), a 90-cm-long graphite target, five scintillator counters and two Cherenkov detectors. Two of the vertex TPCs (VTPC-1 and VTPC-2) are inside the magnetic field created by two superconducting magnets. For the study presented here the magnetic field of the dipole magnets was set to a bending power of 1. 14 T m (standard magnetic field), while a small subset of data was also taken with the full magnetic field of 9 T m. A schematic overview of the setup is presented in Fig. 1. More details can be found in Ref. [1]. The coordinate system is defined as follows: the z-axis is in the nominal direction of the beam, the x-axis in the horizontal plane is such that positively charged particles are bent in the positive x-direction, and the y-axis is perpendicular to the horizontal plane and points upward. The origin is located in the centre of the VTPC-2.
section_name: NA61/SHINE experimental setup

section_content: The NA61/SHINE spectrometer [1] is served by the so-called H2 beamline at the north area of the SPS. The beamline is designed to transport both primary and secondary hadrons and ions from the maximum SPS energy (400 GeV/c) down to 13 A GeV/c. Secondary hadron beams of various momenta
section_name: Beamline

section_content: Fig. 1 A top view (a) and a side view (b) of the NA61/SHINE experimental setup used in 2010 for the T2K replica target data-taking. The beam comes from the left. The orientation of the coordinate system is shown in the bottom left corner, while its origin is located at the centre of VTPC-2 are produced by impinging a primary 400 GeV/c proton beam on a beryllium target. Produced hadrons are selected by two spectrometers and a set of collimators according to their rigidity (momentum to charge ratio). The secondary beam is then transported towards NA61/SHINE. The beam is defined by three scintillator counters used in coincidence (S 1, S 2 and S 3 ) and two scintillator counters with holes used in anticoincidence (V 0 and V p 1 ). The S 3 counter has a radius equal to the target radius, and it was placed 0. 5 cm upstream of the target. The V 0 counter has a 1 cm diameter hole centered on the beam axis, which allows the selection of a narrower beam, if necessary. A Cherenkov Differential Counter with Achromatic Ring Focus [14] (CEDAR) and Threshold Cherenkov Counter (THC) were used in coincidence and anti-coincidence respectively to identify beam particles. By changing the gas pressure in these detectors, it is possible to estimate the beam composition. The 30. 92 GeV/c beam contains around 76. 3% π +, 1. 6% K +, and 12% protons [13]. The estimated purity of the selected proton sample is better than 99. 9%. Finally, the beam position is measured by a set of six multi-wire proportional chambers. Two chambers, one measuring position along the x-axis and another measuring the position along the y-axis, are part of one BPD. The precision of the position measurement in one BPD is around 200 µm [1]. 

The third BPD was mounted on the same support as the target, just upstream of the S 3 counter.
section_name: (b)

section_content: The T2K target [2] has a modular structure which includes a graphite core, a thin titanium case, cooling pipes filled with helium. The replica target design matches dimensions and material of the T2K graphite core. It was made of Toyo Tanso IG-43 graphite [15]. The target density was estimated by measuring its volume and mass, and it is equal to 1. 83±0. 03 g cm -3. This value is 1. 5% higher than the density of the T2K target, but well within the measurement uncertainty. The target shape was machined to form a 2 cm thick disk with an 88 cm long rod coming from its centre. The radii of the disk and rod are 3. 5 cm and 1. 3 cm, respectively. The disk was used to mount an aluminium flange which was, in turn, tightened to a target holder. The upstream side of the flange has a hole to minimize beam interactions with aluminium and provide space for the S 3 counter used in the trigger. A schematic overview of the target can be seen in Fig. 2. The position and tilt of the target were adjusted by the screws on the target holder. The downstream face of the target was placed around 67 cm from the upstream side of VTPC-1.
section_name: Target

section_content: The tracking system consists of two Vertex TPCs located in the magnetic field with an additional small TPC called the Gap TPC (GTPC) located between them, and two Main TPCs (MTPC-L and MTPC-R), located downstream symmetrically with respect to the beamline. All TPCs were inherited from the previous NA49 experiment [16]. The GTPC was introduced in order to measure forward-going particles with small polar angles, θ < 10 mrad, in the laboratory frame. The magnetic field was set so that the total bending power of the magnets was 1. 14 T m. This allowed a momentum resolution of σ p / p 2 = 5×10 -3 (GeV/c) -1, except for very forward tracks that only passed through the GTPC and MTPCs and had a momentum resolution of 22 × 10 -3 (GeV/c) -1 [13]. The angular resolution is around 3-4 mrad and it does not change much as a function of angle. The longitudinal position resolution changes with angle. For angles around 20 mrad it is around 15 cm. For larger angles, it goes down to around 2. 5 cm. Particles emitted from the downstream face of the target behave differently. Forward going particles are usually emitted from the center of the downstream target face. Therefore, extrapolation usually works quite well in this case. However, higher-angle particles coming from the downstream face of the target are usually emitted closer to the edge. In that case, migration effects increase. 

The measured energy loss (dE/dx) in the TPCs provided excellent particle identification capabilities; the achieved resolution is around 4%. In the region where the energy loss distributions for different particles cross, time-of-flight measurements (tof ) were used for particle identification.
section_name: Tracking system

section_content: Although the NA61/SHINE spectrometer has three different Time-of-Flight (ToF) detectors, for the T2K replica target measurements only the forward wall (ToF-F) was used for particle identification. The left and right walls (ToF-L and ToF-R) have superior granularity, but their coverage is insufficient. The ToF-F wall was constructed in 2007 and upgraded in 2009 for the T2K hadron production measurements. It consists of 80 scintillator bars with dimensions W×H×L = 10 × 120 × 2. 5 cm arranged in ten separate modules. The signal is read out by two PMTs placed at both ends of the bars. The estimated tof resolution was 115 ps [13].
section_name: Time-of-flight walls

section_content: The proton beam profile on target at NA61/SHINE is wider than that at J-PARC. Previous studies [10] have shown that longitudinal (z) distribution of the hadrons emitted from the replica target surface depends on the width and position of the incoming proton beam at the target upstream face. In the simplest case where the beam has some radial distribution and where beam divergence and possible re-interactions in the target are neglected, a hadron, produced at position z i along the target with polar angle θ and radial position Δr with respect to the target longitudinal axis, would escape the target at a minimal longitudinal position, z, given by:

Evidently, for narrower beams, the longitudinal (z) distribution of the hadrons exiting the target surface is pushed more downstream. With a sufficiently narrow beam, one can suppress the number of hadrons coming from the first 18 cm of the target by an order of magnitude or more. Although this effect seems to be purely geometrical, a Monte Carlo simulation suggest that different beam widths and positions can lead to slightly different numbers of low momentum pions produced in re-interactions. The difference in the number of pions is around 1%. The NA61/SHINE proton beam achieved in 2009 [10] is much wider than the Fig. 3 Radial distributions of the incoming proton beam with respect to the center of the target: T 3 beam profile (blue), T 2 beam profile (red) and T2K beam profiles for two different runs (green). Beam profiles in T2K change on a run-by-run basis and perfect agreement is not possible. All distributions are normalised to the corresponding total number of protons on target J-PARC beam. To avoid large differences between the NA61/SHINE beam and the J-PARC beam, four special triggers were used during the 2010 data-taking (see Fig. 1) :

The T 1 trigger is the only trigger that selects beam protons that did not necessarily hit the target or when a hit in the S 3 scintillator was not detected. The T 2 and T 3 triggers select only beam protons that hit the target. The only difference is that T 3 selects a narrower beam because the V 0 counter has a smaller hole. The fourth trigger selects all beam particles that hit the target. It is important to note that T 1, T 2 and T 4 were prescaled -only a fraction of them was recorded. The prescaling was applied to T 2 to achieve the desired radial distribution of the beam. If an event is a T 2 but not a T 3, this means that the triggered beam particle is in the tail of the radial distribution on the upstream target face. To reduce the beam tail, only one out of two such events was registered. 

The achieved radial distribution is shown in Fig. 3. Perfect agreement with the T2K beam profile is not possible since in T2K the beam profile changes on a run-by-run basis. Any differences must be studied before using these measurements in the T2K neutrino beam simulation. Such study is presented in Sect. 6. In total, 10. 2×10 6 triggers were recorded. Around 11. 7% of the triggers were recorded with the maximum magnetic field of 9 T m. This allowed the recording of beam protons that pass through the target without interacting and the measurements of their parameters in the TPCs. However, all sec-ondary hadrons below 6GeV/c were bent out of the TPCs and did not reach the ToF-F wall. The analysis of data recorded with the maximum magnetic field is a subject of an independent publication. In this paper, the data taken with the standard magnetic field configuration (89. 3% of triggers) were analysed, while the high magnetic field subset was used for alignment between the BPDs and TPCs (see Sect. 3. 1).
section_name: Triggers

section_content: The reconstruction and simulation procedures used for the replica target analysis are described in earlier publications [9, 10]. 

The analysis of the data taken with the standard magnetic field configuration was performed with the so-called tof -dE/dx method, which makes use of the TPC measurements of the specific energy loss (dE/dx) and the time-of-flight measurements [9, 10]. The specific energy loss of a track is calculated as a truncated mean [16] of the charges of the clusters (points) on the track traversing the TPCs. After event and track selection, selected tracks were binned in momentum ( p), polar angle (θ ), both measured in the laboratory frame, and longitudinal position of the exit point on the target surface (z), as required by T2K. The different hadron yields were estimated by fitting a two-dimensional dE/dx-m 2 tof distribution for every bin. These raw yields were then corrected for all inefficiencies by applying Monte Carlo and data-based correction factors. Before any event or track selection is applied, it was necessary to determine the target position with respect to the BPDs and TPCs in order to refine the original surveyor measurements.
section_name: Analysis

section_content: The target upstream face x and y positions with respect to the BPDs can be determined just by plotting the T 2 beam profile at the target upstream z position. The position of the target centre was obtained by maximising the number of entries in the circle with radius equal to 1. 3 cm. This was only possible because S 3 counter, which was part of the trigger, had the same radius as the target and it was placed just upstream of the target. The target upstream z position was measured by surveyors before the data-taking. Once the relative alignment of the target and BPDs was known, the alignment of the TPCs with respect to the BPDs could also be determined. To do so, the data taken with the high magnetic field configuration were used to select events with only one high momentum track corresponding to the beam track. Those tracks were extrapolated to the surveyed z position of the target using an analytical method (see Ref. [17] ) with error propagation that accounts for multiple scattering (see Ref. [18] ). The positions of the beam tracks extrapolated forward from the BPDs and 5 The distribution of z positions of closest approach for TPC tracks extrapolated to the beam track. Only tracks whose distance to the beam track is within extrapolation uncertainty were selected. The red line shows fitted rising edge and the position of the half maximum was taken to be the target upstream z position backward from the TPCs were compared and the mean value of the differences was taken to be the misalignment of the TPCs with respect to the BPDs (see Fig. 4 ). 

Because of the non-negligible probability of re-interactions and long extrapolation distance, it is difficult to reconstruct interaction vertices in the target. However, by extrapolating TPC tracks and selecting only those for which the point of closest approach to the beam track is within the extrapolation uncertainty, it was possible to check the surveyors' measurement of the target z position. The z distribution of the points of closest approach was plotted and a rising edge could be seen as shown in Fig. 5. The z position at the halfmax of this rising edge was taken to be the upstream target face. Full agreement with the surveyor position measurement was obtained. The final step was to determine the tilt of the target in x -z and yz planes. First, the target was assumed to be parallel to the z-axis. Then, only events with beam tracks passing through the whole length of the target were selected. The TPC tracks were again extrapolated backwards until the minimum distance from the beam track was reached. By plotting the x and y distributions for 18 slices in z and checking the width and position of the distributions, it is straightforward to extract any possible tilt of the target. For example, if the target is tilted in positive x direction, edges of the x distributions would be shifted towards the centre because the target is moved from the beam. The shift would increase with z. No appreciable tilt was found. In Table 1 the measured position of the target is summarized with uncertainties for all parameters.
section_name: Target position and alignment

section_content: Three different criteria were imposed on the recorded events. Beam particles were required to hit the target, in other words, an event must correspond to T 2 or T 3 trigger. Also, the position of each beam particle must have been measured by all three BPDs. And finally, the reconstructed path of the beam particle must have passed through the whole length of the target. 

Only 75% of the events passed this selection. Events were mostly discarded at the second requirement due to inefficiencies in the BPD measurements. The event selection should Fig. 6 The θ -φ distribution of selected tracks: regions of flat acceptance are selected. The selected regions in azimuthal angle φ depend on the polar angle θ interval Table 2 Fractions of (anti)neutrino fluxes at Super-Kamiokande produced by π +, π -, K +, K -, and p emitted from the surface of the T2K target (estimated using Fluka2011. 2c. 5 simulations) and covered by the measurements presented in this paper. The fractions are presented for the forward and reverse horn current configurations. The last column is a sum of the coverage for five hadrons presented in this paper and does not include other particles contributing to the neutrino flux (K 0 s, Λ, μ) also ensure that events containing off-time TPC tracks were not selected as these tracks can potentially bias the measured hadron yields. The beam intensity during the data-taking was around 8. 3 kHz which means that the mean time difference between beam particles was around 120 μs. However, the ToF-F wall has a time acquisition window of 100 ns, and because of the requirement that all tracks must have a tof hit for particle identification, most of the off-time TPC tracks were automatically discarded. The fraction of accepted offtime tracks was negligible.
section_name: Event selection

section_content: Track selection was rather simple and can be divided into two parts. First, a track was required to have fitted momentum, energy loss, and time-of-flight measurements. Afterwards, cuts were applied to increase the quality of the selected track sample. Tracks can have segments in different TPCs and therefore, can be divided into different topologies. For each topology, a different requirement on the number of clusters was applied. For example, if a track passes through the VTPCs, at least 25 points in the VTPCs were required. A cut on the number of clusters in MTPCs was not applied to such tracks. However, tracks with momentum measured only in the GTPC have between four and seven GTPC points. At least six points in the GTPC were required as well as additional 30 points in the MTPCs. This ensured that the track parameters were measured with sufficient precision and that the track did not pass very close to the MTPC walls where possible field distortions are the largest. Afterwards, the azimuthal angle distributions were plotted for polar angle intervals of 20 mrad or 40 mrad. Two different types of azimuthal regions were removed: the regions with rapidly changing acceptance and the regions in which the Monte Carlo acceptance does not correspond to the data. This may be caused by slightly different sizes of the TPCs or distortions of the magnetic field which are not present in the magnetic field maps used in the Monte Carlo simulation. The selected regions in (φ, θ) space are presented in Fig. 6. This reduces the possibility of bias once the Monte Carlo corrections were applied to the data. Finally, the selected TPC tracks were extrapolated towards the target surface. The track parameters and covariance matrices were saved for the positions where tracks hit the target surface or the minimum distance from the target surface was less than 3σ R, where σ R is the extrapolated radial uncertainty. This cut reduced the number of selected tracks that were not coming from the target, but were created in decays or interactions outside of the target.
section_name: Track selection

section_content: Selected tracks were binned in ( p, θ, z) space. As required by T2K [10], six longitudinal bins in total were used: the 90-cm-long graphite rod was divided into five bins of 18 cm each and the target downstream face was considered as an additional sixth bin. Different ( p, θ) binning was applied for extracting π ±, K ± and p yields. The reason for this is due to statistics: several times fewer kaons than pions are expected, so coarser binning must be used for kaons. Unequal bin sizes [19] [20] [21]. Also, the starting value for the momentum binning was carefully adjusted. The necessity for this comes from the fact that the ToF-F response was not simulated in the Monte Carlo. Instead, in the reconstruction chain applied to simulated events the tof measurements were just assigned to the tracks that hit the ToF-F wall, and later on the inefficiency of the ToF-F wall was corrected with a data-based correction. Actually, very low momentum particles, depending on their mass, cannot reach the ToF-F wall within the ToF-F acquisition window. Therefore, Monte Carlo corrections in this Fig. 16 Efficiency of the scintillator bars in the ToF-F wall as a function of the x-coordinate. It is clear that the efficiency drops if the bar is closer to the beamline (x = cm). This is due to an increasing track density with decreasing distance to the beam region would be heavily biased. For example, momentum bins for protons must start at 0. 5 GeV/c. The overall phase space binning covers more than 96% of (anti)neutrinos crossing the T2K far detector (Super-Kamiokande) which are produced by the charged hadrons emitted from the target. More than 97% of (anti)neutrinos produced by pions are covered, while the coverage drops for kaons and protons. This is summarised in Table 2. The phase space of the charged hadrons coming from the T2K target surface that contribute to the (anti)neutrino fluxes in Super-Kamiokande are plotted in Figs. 7, 8, 9, 10 and 11. The top row in each figure represents the coverage for the forward horn current (FHC) corresponding to the neutrino mode, while the bottom row represents the coverage for the reverse horn current (RHC) corresponding to the antineutrino mode. Each column represents a different longitudinal bin. The phase space coverage of the measurements presented in this paper is overlaid on top of the figures as a black line. Since NA61/SHINE added new forward TPCs in 2017 [22], possible future measurements could improve coverage for the forward-going high momentum K ± and p [23].
section_name: Phase space

section_content: Since large fractions of the phase space are covered, a robust particle identification method is needed. For the momentum range between 1 GeV/c and 3 GeV/c (see Fig. 12 ), the energy loss distributions cross, hence particle identification based on the energy loss alone is not possible. However, in these regions, tof measurements can be used to distinguish between particles. Kaons can be separated easily up to 3 GeV/c, while protons can be separated up to 8 GeV/c, as can be seen in Fig. 13. For higher momentum, the tof resolution becomes too poor, while energy loss measurements allow for better identification. It is clear that both approaches are complementary and therefore they can be combined to cover all the bins. The tof measurement was used to calcu- Fraction Fig. 19 Uncertainties on K + yields shown as a function of p and θ: statistical uncertainties (top row), total systematic uncertainties (middle row) and total uncertainties (bottom row). Each column corresponds to a different longitudinal bin late the particle mass squared (m 2 tof ) and it was combined with the energy loss. Particles were represented by islands in the m 2 tof -dE/dx space. Therefore, raw hadron yields could be obtained by fitting an appropriate function to the distribution for each phase space bin. Both, m 2 tof and dE/dx were assumed to be normally distributed. This assumption for dE/dx must be closely examined. In general, dE/dx is not normally distributed, as its distribution is similar to the Landau distribution. The measured dE/dx in NA61/SHINE is calculated as a truncated mean (see Ref. [16] ) of all energy depositions for all clusters. Therefore, we assume that truncated mean is normally distributed. This is only valid if all tracks have the However, the dE/dx resolution in NA61/SHINE saturates at around 3. 5% for tracks with more than 70 clusters. Around 98. 5% of the selected tracks have more than 70 clusters. This is a good justification for using only a single Gaussian for describing energy loss in a single phase space bin for one particle species. Our previous measurements prove this assumption (see Refs. [8, 10] ). The total fitting function was constructed from four two-dimensional Gaussians, one for each particle species (e ±, π ±, K ±, p (p)). The fitting was done in the RooFit framework [24] by using extended loglikelihood minimization which treats the number of observed Fig. 22 Systematic uncertainties for the measured yields as a function of momentum. Each column corresponds to a different particle type, from left to right: π +, π -, K +, K -, and p. Each row corresponds to a different longitudinal bin. Specific theta ranges are given on top of each plot. Presented ( p, θ, z) bins are selected to show typical sizes of the systematic uncertainties particles as a Poisson random variable. An example of a fit is shown in Fig. 14. In total, there were 20 parameters in the fit: eight mean values, eight standard deviations, and four particle multiplicities. The obtained raw multiplicities must be corrected for various inefficiencies with Monte Carlo and data-based corrections.
section_name: Particle identification

section_content: Interactions in the target were simulated by Fluka2011. 2c. 5 [19] [20] [21]. The kinematics of the particles at their exit position on the surface of the target was then fed to Geant3 /GCALOR [25] for propagation through the NA61/SHINE detector. The incoming proton beam information from the data was used to generate the Monte Carlo beam profile. The beam track parameters for each simulated proton were thrown randomly according to the beam divergence and position distributions shown in Fig. 15. The x-and y-positions and divergences were considered to be independent. In total, nearly 40 × 10 6 protons on target (POT) were simulated. 

A simple Monte Carlo correction was applied to the raw yields:

where N i jk sim is the number of simulated hadrons in the bin ( p, θ, z) = (i, j, k) and N i jk rec,sim is the number of simulated hadrons reconstructed in the same bin. This correction can be separated into several factors, which allow for a better determination of the associated systematic errors:

where ε φ i jk is the geometrical efficiency (because of the φ cut), ε loss i jk is the hadron loss efficiency due to decays and interactions in the detector (hadrons which are produced in re-interactions and by chance extrapolated to the target surface are included in this correction), ε rec i jk is the reconstruction efficiency, ε sel i jk is the selection efficiency, ε mig i jk is the bin migration efficiency and ε feed-down i jk is the feed-down efficiency. The geometrical efficiency gives the percentage of 4π solid angle covered by the NA61/SHINE detector after the selection. Some of the produced hadrons inside the detector coverage can decay or re-interact in the detector and therefore, they do not hit the TOF-F wall. The percentage of surviving hadrons is given by the hadron loss efficiency. The reconstruction efficiency is a convolution of the efficiency of software algorithms used for the reconstruction and the efficiency of the detector. The selection efficiency corrects for the hadrons lost in the quality selection. The values of momentum ( p), polar angle (θ ) and longitudinal position along the target surface (z) for each hadron are reconstructed with finite precision. Consequently, some of the hadrons will be placed in a wrong ( p, θ, z) bin. The ratio of the number of hadrons which are reconstructed in a given bin to the number of hadrons which are simulated in 123 the same bin represents the migration efficiency. The feeddown efficiency corrects for weak decays outside of the target whose decay products are extrapolated to the target surface. 

Although only the total Monte Carlo correction was applied to raw yields, it was useful to study each contribution separately in order to estimate possible biases. The systematic uncertainties for each correction are presented in the next section. As previously described, the efficiency of the ToF-F detector was not included in the Monte Carlo correction.
section_name: Monte Carlo correction factors

section_content: A different track selection was used to estimate the ToF-F efficiency from the data. Only tracks that reached the end of the MTPCs were selected. Since the ToF-F wall was placed just downstream of the MTPCs, it was necessary to make sure that tracks did not re-interact or decay before hitting the ToF-F wall. The efficiency of each scintillator bar was simply calculated as the number of tracks with hits in the ToF-F wall divided by the total number of selected tracks. Inefficiency was caused mainly by two effects. First, if two tracks hit a scintillator bar in the same event, the hit was discarded during reconstruction since it was not possible to distinguish between these tracks (close to the beamline there is a 4% probability another track will be matched to the same ToF-F hit, while this probability is below 1% further away from the beamline). Second, a hit is also discarded if the difference in time of flight measured by the top and bottom PMTs in a bar was larger than 2000 ps. During the data-taking period, four scintillator bars had only one PMT working or there was a DAQ problem for these channels. These scintillators did not suffer from the inefficiency caused by the quality cut during reconstruction. However, these bars had a worse time resolution by a factor of √ 2. The efficiency of each bar is shown in Fig. 16. Efficiency uncertainties were calculated as binomial errors. This information was then used to estimate the ToF-F efficiency for each phase space bin. Tracks in a single bin usually contain hits in several scintillators (usually around three or four). The error on the efficiency for each bin was taken as a ToF-F systematic uncertainty.
section_name: ToF-F efficiency correction factors

section_content: A 50% drop in the number of ToF-F hits in the data was observed for a couple of upstream longitudinal bins and high polar angle (θ ≈ 300 mrad). However, such a drop was not seen in the MC simulation and could not be related to the ToF-F inefficiency. Under further inspection, it was discovered that the corresponding tracks actually did not reach the MTPCs but were bent out or absorbed while pass-ing very close to the edge of the magnetic field map used in the simulation. It is thus possible that the field map in this region is biased. It appears to be biased asymmetrically as negatively-charged tracks were more affected than positively-charged ones, pointing out to asymmetric E × B corrections in the TPCs. Nevertheless, the affected bins represent only around 1% of the total π ± bins and other hadrons were not affected, since they were not measured in the affected area because of the low statistics. To account for this inefficiency, an additional ad hoc correction factor was applied:

where

is the number of selected tracks with tof hit in the Monte Carlo (data), n MC(data) sel

is the number of selected tracks in the Monte Carlo(data) without the tof requirement, and C tof i jk is the ToF-F efficiency for the given phase space bin. It is clear that this correction can be model-dependent and therefore this point must be addressed for a possible systematic bias.
section_name: Ad hoc correction factors

section_content: The systematic uncertainties of this analysis were carefully studied. A brief summary is given below, while more detailed information can be found in Ref. [13]. 

The following systematic uncertainties were considered:

(i) Hadron loss uncertainty. Apart from the ToF-F inefficiency, tracks can miss the ToF-F wall, re-interact in the detector, or decay before reaching the end of the MTPCs. Possible imperfections in the Monte Carlo description of the detector can create biases while correcting for the previously described effects. To check for potential biases, an attempt to impose more strict cuts and to select only tracks with a long segment in the MTPCs was made. After re-calculating and comparing the yields to the standard results, the difference was taken as a systematic uncertainty. In the majority of bins, this uncertainty is well below 1%. However, for tracks with a momentum of several GeV/c in the 20-40 mrad region in the most upstream longitudinal bin, the hadron loss uncertainty goes up to 20% (for a single bin it reached 50%). These tracks cross the beamline and pass very close to the walls of the MTPCs. (ii) Backward track extrapolation (bin migration) uncertainty. The uncertainty on the target position can induce systematic bin migration while extrapolating tracks backwards to the surface of the target. The position of the target was changed in the Monte Carlo and hadron yields were recalculated. Differences were taken as systematic uncertainties. For tracks with a very small angle originating from the first two longitudinal bins, this uncertainty goes up to 8%. Also, high angle tracks extrapolated toward the downstream face of the target are sensitive to the x and y shifts of the target, and for them the bin migration uncertainty is up to 5%. For all other regions of the phase space, the bin migration uncertainty drops below 1%. (iii) Reconstruction efficiency uncertainty. This uncertainty was estimated in the previous NA61/SHINE publications [8, 10]. Since the experimental setup and the reconstruction software did not change, the same value was used for the measurements presented in this paper. A conservative value of 2% is assigned for all bins. (iv) Particle identification (PID) uncertainty. The width of the energy loss distribution for a single bin depends on the number of clusters in tracks. The energy loss distribution should be a sum of Gaussians, where each Gaussian represents tracks with a distinct number of clusters. Selected tracks, in this case, are very long and the width of the energy loss distribution becomes nearly constant. For this reason, a single Gaussian was used to describe the energy loss distributions. However, at higher momenta and, thus, for larger momentum bins, the mean values will also shift and possibly create an asymmetry in the distribution. To check for potential biases, a sum of two Gaussians was used and the obtained yields were compared. The mean values of Gaussians were kept independent, while the width of the second Gaussian was forced to be larger than the width of the first Gaussian. For low momentum, no difference was found, however, for higher momentum and larger bins, a difference of up to 1% was observed for π ± and p. A constant value of 1% was assigned to all of these bins. Positively charged kaons are located under the proton peak in the energy loss distribution and the pion peak in the m 2 tof distribution for momentum larger than 3GeV/c. For this reason, the uncertainty for K + can reach 15% for the higher momentum bins. In case of K -, the obtained uncertainties are similar to the ones for K +. (v) Feed-down uncertainty. Produced K 0 s and Λ in the target can decay before reaching the VTPC-1 and some of the daughter particles reconstructed in the TPCs can be extrapolated to the target surface. Since the correction for this effect is calculated with the Monte Carlo simulation, the number of produced K 0 s and Λ is model dependent. A systematic uncertainty of 30% on the feed-down correction factor was assigned, following the standard approach in Refs. [8, 10]. The systematic uncertainty is up to 2% for π ± and protons. Mostly low momentum bins in the first and last longitudinal bins are affected. (vi) The ToF-F efficiency uncertainty. The uncertainty of the ToF-F efficiency was taken as a systematic uncertainty. It goes up to 2%, but for most of the bins it is well below 1%. (vii) Ad-hoc uncertainty. A constant conservative uncertainty of 25% was assigned to the bins that are corrected by the ad-hoc correction factor. The value was based on half the size of the typical ad hoc correction (50%). 

Total systematic uncertainties for individual analysis bins and for each particle type are shown in Figs. 17, 18, 19, 20 and 21 (middle row). 

Respective contributions to the total systematic uncertainties for different particle types are shown in Fig. 22. Different regions of the phase space have different dominant contributions. For example, the dominant contribution in the third longitudinal bin (z 3 ) for π ± yields comes from the reconstruction efficiency uncertainty, while the dominant contribution for K ± comes from the PID uncertainty.
section_name: Systematic uncertainties

section_content: Results are presented in the form of double differential yields normalised by the total number of protons on target:

where (i, j, k) are the (p, θ, z) bin numbers, N P OT is the number of protons on target, α = π +, π -, K +, K -, p is the hadron species, N α i jk is the number of extracted hadrons in a given bin, C i jk is the correction factor, Δp i jk is the momentum bin size, and Δθ i j is the polar angle bin size. Uncertainties on the hadron yields are shown in Figs. 17, 18, 19, 20 and 21. The top rows show statistical uncertainties, the middle rows show systematic uncertainty and the bottom rows show total uncertainties. For π +, π -and p statistical uncertainties dominate in the low (< 40 mrad) and high (> 200 mrad) polar angle region. While for 40 < θ < 200 mrad, the systematic uncertainties are comparable to or larger than the statistical uncertainties. For kaons, the total uncertainty is mostly dominated by the statistical uncertainty. Tables with numerical results are presented in Ref. [26]. QGSP_BERT (blue) physics lists from Geant4. 10. 03 [27, 28] are overlaid on top of the data
section_name: Results and comparison with hadron production models

section_content: The results were compared with NuBeam and QGSP_BERT physics lists from Geant4. 10. 03 [27, 28]. Detailed comparisons are presented in Figs. 23, 24, 25, 26, 27 better predictions of π +, π -, K +, and K -yields. However, both physics lists fail to accurately predict p yields. Similar behaviour was noticed in the previous NA61/SHINE measurements of primary proton-carbon interactions at 31 GeV/c [8], where it was shown that most of the models poorly reproduce proton yields. The earlier observation [10] that a better agreement with Fluka2011 predictions can be obtained by lowering the production cross section in Fluka2011 is confirmed.
section_name: Comparisons with models

section_content: As previously mentioned, NA61/SHINE provided measurements of the π ± yields emitted from the T2K replica target using data collected in 2009 [10]. Since the beam profiles for the data taken in 2009 and 2010 are different, a direct comparison between the 2009 measurements and the measurements presented in this paper is not possible. The beam profile for the 2009 measurements is wider, and the 2009 pion yields are expected to be larger for the low angle upstream longitudinal bins. The differences decrease for more downstream longitudinal bins and higher angles. But even for the downstream face of the target and high angles (> 100 mrad) perfect agreement is not expected since Monte Carlo simulations show small differences in re-interactions for different beam profiles. Detailed Monte Carlo studies can be found in Ref. [13]. However, under the assumption that major dif- QGSP_BERT (dashed black line) physics lists from Geant4. 10. 03 [27, 28] are overlaid on top of the data applied in the phase space binning as defined in Ref. [10]. Also, Monte Carlo and tof correction factors were calculated, but the ad hoc correction was not applied. It is important to note that systematic uncertainties were not reevaluated for the new phase space bins. Qualitative comparisons of π + and π -yields can be seen in Figs. 53, 54, respectively. Only two polar angle bins are selected and presented for all longitudinal bins. A quantitative agreement between the two sets of measurements is not expected due to the simplified assumptions described above, as well as differences in the beam divergence which were not accounted for. Indeed, some differences in the π + yields can be observed for small angles. However, in general, these two sets of measurements are consistent within about 20%. Note that the systematic uncertainties were not reevaluated for these special 2010 measurements and they are not shown on the plots. Additionally, only differences in the radial beam distribution were taken into account, while differences in the beam divergence were ignored. 

In order to demonstrate the reliability of the computed π -yields, an independent full chain of data analysis was performed. By means of the software tools and procedures QGSP_BERT (dashed black line) physics lists from Geant4. 10. 03 [27, 28] are overlaid on top of the data approaches provide yields in acceptable agreement. Some minor differences were observed only for the first longitudinal target bin, for small momentum (below 6 GeV/c) and polar angles (θ below 40 mrad), and in the θ interval [300; 400] mrad for the first three longitudinal target bins.
section_name: Comparison with previous NA61/SHINE measurements

section_content: The results presented in Sect. 5 will be used in the T2K neutrino beam simulation for re-weighting of the hadron yields on the target surface. Details of the currently used re-weighting procedure based on NA61/SHINE thin-target measurements are presented in Ref. [4]. 

The weights are simply calculated as the ratio of the data yields and the Fluka2011 yields for each ( p, θ, z) bin as follows:

These weights are applied to all neutrinos that have a hadron ancestor exiting the target surface in the phase space covered by NA61/SHINE. As previously mentioned, hadron QGSP_BERT (dashed black line) physics lists from Geant4. 10. 03 [27, 28] are overlaid on top of the data yields coming from the surface of the long target depend on the width and position of the incoming proton beam hitting the upstream face of the target. Since it was shown that this effect is mostly geometrical, there is the possibility that the reweighting factors are invariant under the beam profile change. To test this hypothesis, Fluka2011. 2c. 5 was used as input for data and the NuBeam physics list from Geant4. 10. 03 as input for the Monte Carlo. Around 40 × 10 6 protons on target were simulated for NA61/SHINE T 2 and T 3 beam profiles using both models. Re-weighting factors were calculated for both beam profiles and compared. Examples of the ratios for π + and p are shown in Fig. 55. Some differences can be noticed for the upstream z bins, while for other bins any differences are mostly below 2%. It is important to note that the differences between the T 2 and T 3 beam profiles are very large in comparison to differences between the T 2 and T2K beam profiles. and QGSP_BERT (dashed black line) physics lists from Geant4. 10. 03 [27, 28] are overlaid on top of the data
section_name: Dependence of the T2K re-weighting factors on the proton beam profile

section_content: Measurements of the double differential yields of π +, π -, K +, K -and p emitted from the surface of a 90-cm-long carbon target (T2K replica) with incoming 31 GeV/c protons were performed by the NA61/SHINE experiment at the CERN SPS. Yields of π + and π -were measured with improved precision compared to the previously published NA61/SHINE results, while K +, K -and p yields were obtained for the first time. These measurements are crucial for reducing the hadron production component of the T2K (anti)neutrino flux error, which is the dominant component in the flux uncertainty. Any reduction of the flux uncertainties will directly improve measurements of the (anti)neutrinonucleus cross sections as well as of the (anti)neutrino oscillation parameters in T2K. A simple method of using these results in T2K which avoids problems with the dependence on the incident proton beam profile is proposed.
section_name: Summary and conclusions

section_content: Acknowledgements We would like to thank the CERN EP, BE and EN Departments for the strong support of NA61/SHINE. This work was supported by the Hungarian Scientific Research Fund (Grants NKFIH 123842-123959 ), the János Bolyai Research Scholarship of the Hungarian Academy of Sciences, the Polish Ministry of Science and Higher Education (Grants 667/N-CERN/2010/ 0, NN 202 48 4339 and NN 202 23 1837 ), the Polish National Center for Science (grants 2011/03/N/ST2/03691, 2013/11/N/ST2/03879, 2014/13/N/ST2/02565, 2014/14/E/ST2/00018, 2014/15/B/ST2/02537 and 2015/18/M/ST2/00125, 2015/19/N/ST2 /01689, 2016/23/B/ST2/ 00692, UMO-2014/15/B/ST2/03752 ), the Russian Science Foundation, grant 16-12-10176, the Russian Academy of Science and the Russian Foundation for Basic Research (Grants 08-02-00018, 09-02-00664 and 12-02-91503-CERN ), the Ministry of Science and Education of the Russian Federation, Grant No. 3. 3380. 2017/4. 6, the National Research Nuclear University MEPhI in the framework of the Russian Academic Excellence Project (contract No. 02. a03. 21. 0005, 27. 08. 2013 ), the Ministry of Education, Culture, Sports, Science and Technology, Japan, Grant-in-Aid for Scientific Research (grants 18071005, 19034011, 19740162, 20740160 and 20039012), the German Research Foundation (Grant GA 1480/2-2 ), the Bulgarian Nuclear Regulatory Agency and the Joint Institute for Nuclear Research, Dubna (bilateral contract No. 4418-1-15/17 ), Bulgarian National Science Fund (Grant DN08/11 ), Ministry of Education and Science of the Republic of Serbia (grant OI171002 ), Swiss Nationalfonds Foundation (Grant 200020117913/1 ), ETH Research Grant TH-01 07-3 and the U. S. Department of Energy.

section_name: Data Availability Statement

section_content: This manuscript has no associated data or the data will not be deposited. [Authors' comment: All tables with numerical results are available from the link presented in Ref. [26]. ]
section_name: Data Availability Statement

section_content: Open Access This article is distributed under the terms of the Creative Commons Attribution 4. 0 International License (http://creativecomm ons. org/licenses/by/4. 0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3.