of , and proton
the surface of the T2K replica target for incoming 31 GeV/ protons with the NA61/SHINE spectrometer at the CERN SPS \PreprintIdNumberCERN-EP-2018-222 \ShineJournalEur. Phys. J. C \ShineAbstract Measurements of the , , and proton double differential yields emitted from the surface of the -cm-long carbon target (T2K replica) were performed for the incoming protons with the \NASixtyOnespectrometer at the CERN SPS using data collected during 2010 run. The double differential yields were measured with increased precision compared to the previously published \NASixtyOneresults, while the 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.
N. Abgrall, A. Aduszkiewicz, E.V. Andronov, T. Antićić, B. Baatar, M. Baszczyk, S. Bhosale, A. Blondel, M. Bogomilov, A. Brandin, A. Bravar, W. Bryliński, J. Brzychczyk, S.A. Bunyatov, O. Busygina, A. Bzdak, H. Cherif, M. Ćirković, T. Czopowicz, A. Damyanova, N. Davis, M. Deveaux, W. Dominik, P. Dorosz, J. Dumarchez, R. Engel, A. Ereditato, G.A. Feofilov, L. Fields, Z. Fodor, A. Garibov, M. Gaździcki, O. Golosov, M. Golubeva, K. Grebieszkow, F. Guber, A. Haesler, T. Hasegawa, A.E. Hervé, S.N. Igolkin, S. Ilieva, A. Ivashkin, S.R. Johnson, K. Kadija, E. Kaptur, N. Kargin, E. Kashirin, M. Kiełbowicz, V.A. Kireyeu, V. Klochkov, T. Kobayashi, V.I. Kolesnikov, D. Kolev, A. Korzenev, V.N. Kovalenko, K. Kowalik, S. Kowalski, M. Koziel, A. Krasnoperov, W. Kucewicz, M. Kuich, A. Kurepin, D. Larsen, A. László, T.V. Lazareva, M. Lewicki, K. Łojek, B. Łysakowski, V.V. Lyubushkin, M. Maćkowiak-Pawłowska, Z. Majka, B. Maksiak, A.I. Malakhov, D. Manić, A. Marchionni, A. Marcinek, A.D. Marino, K. Marton, H.-J. Mathes, T. Matulewicz, V. Matveev, G.L. Melkumov, A.O. Merzlaya, B. Messerly, Ł. Mik, G.B. Mills, S. Morozov, S. Mrówczyński, Y. Nagai, T. Nakadaira, M. Naskręt, K. Nishikawa, V. Ozvenchuk, V. Paolone, M. Pavin, O. Petukhov, C. Pistillo, R. Płaneta, P. Podlaski, B.A. Popov, M. Posiadała, S. Puławski, J. Puzović, W. Rauch, M. Ravonel, R. Renfordt, E. Richter-Wąs, D. Röhrich, E. Rondio, M. Roth, B.T. Rumberger, A. Rustamov, M. Rybczynski, A. Rybicki, A. Sadovsky, K. Sakashita, K. Schmidt, T. Sekiguchi, I. Selyuzhenkov, A.Yu. Seryakov, P. Seyboth, M. Shibata, M. Słodkowski, A. Snoch, P. Staszel, G. Stefanek, J. Stepaniak, M. Strikhanov, H. Ströbele, T. Šuša, M. Tada, A. Taranenko, A. Tefelska, D. Tefelski, V. Tereshchenko, A. Toia, R. Tsenov, L. Turko, R. Ulrich, M. Unger, F.F. Valiev, D. Veberič, V.V. Vechernin, M. Walewski, A. Wickremasinghe, Z. Włodarczyk, A. Wojtaszek-Szwarc, O. Wyszyński, L. Zambelli, E.D. Zimmerman, R. Zwaska,
The T2K Beam Group
L. Berns, G.A. Fiorentini, M. Friend, M. Hartz, T. Vladisavljevic, M. Yu
National Nuclear Research Center, Baku, Azerbaijan
Faculty of Physics, University of Sofia, Sofia, Bulgaria
Ruđer Bošković Institute, Zagreb, Croatia
LPNHE, University of Paris VI and VII, Paris, France
Karlsruhe Institute of Technology, Karlsruhe, Germany
Fachhochschule Frankfurt, Frankfurt, Germany
University of Frankfurt, Frankfurt, Germany
Wigner Research Centre for Physics of the Hungarian Academy of Sciences, Budapest, Hungary
University of Bergen, Bergen, Norway
Jan Kochanowski University in Kielce, Poland
H. Niewodniczański Institute of Nuclear Physics of the Polish Academy of Sciences, Kraków, Poland
National Centre for Nuclear Research, Warsaw, Poland
Jagiellonian University, Cracow, Poland
AGH - University of Science and Technology, Cracow, Poland
University of Silesia, Katowice, Poland
University of Warsaw, Warsaw, Poland
University of Wrocław, Wrocław, Poland
Warsaw University of Technology, Warsaw, Poland
Institute for Nuclear Research, Moscow, Russia
Joint Institute for Nuclear Research, Dubna, Russia
National Research Nuclear University (Moscow Engineering Physics Institute), Moscow, Russia
St. Petersburg State University, St. Petersburg, Russia
University of Belgrade, Belgrade, Serbia
University of Geneva, Geneva, Switzerland
Fermilab, Batavia, USA
University of Colorado, Boulder, USA
University of Pittsburgh, Pittsburgh, USA
Los Alamos National Laboratory, Los Alamos, USA
Institute for Particle and Nuclear Studies, KEK, Tsukuba, Japan
Tokyo Institute of Technology, Department of Physics, Tokyo, Japan
Department of Physics and Astronomy, York University, Toronto, ON, Canada
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, University of Tokyo, Kashiwa, Japan
TRIUMF, Vancouver, BC, Canada
Oxford University, Department of Physics, Oxford, United Kingdom
University of Bern, Bern, Switzerland
(SPS Heavy Ion and Neutrino Experiment)  is a large hadron spectrometer at the CERN Super Proton Synchrotron (SPS). The \NASixtyOnecollaboration is pursuing several physics goals including hadron production measurements for T2K (Tokai to Kamioka)  - an accelerator-based long-baseline neutrino experiment in Japan. The \NASixtyOnemeasurements are used to reduce the systematic uncertainties associated to the prediction of the (anti)neutrino fluxes in T2K. New measurements of , 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 Section 2, a brief overview of the \NASixtyOnesetup is shown followed by the description of the analysis in Section 3. Section 4 gives a concise description of all systematic uncertainties. Results and comparison with Monte Carlo (MC) models and with the previous \NASixtyOnemeasurements are presented in Section 5. Finally, in Section 6, proper usage of these results in the T2K neutrino beam simulation is discussed, followed by a short conclusion.
1.1 The T2K neutrino beam
The T2K neutrino beam is produced at the Japan Proton Accelerator Research Complex (J-PARC)  by directing a GeV (kinetic energy) proton beam towards a -cm-long graphite target. Produced mesons, mostly pions and kaons , are focused by a set of three magnetic horns  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 %. For this reason, T2K uses available hadron production measurements, but mostly relies on dedicated \NASixtyOnehadron production data .
1.2 The \NASixtyOnemeasurements for T2K
The \NASixtyOneexperiment took data for T2K in 2007, 2009 and 2010. In 2007 and 2009, measurements were done with a proton beam and a -cm-thick carbon target (4% of a nuclear interaction length, ) to measure hadron multiplicities (, , p, , ) 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 . 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 % 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 . These measurements were also done in 2007, 2009 and 2010. Measurements from the 2007 run  were used as a proof-of-principle, 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 % , representing a significant improvement with respect to the previously published uncertainty of . 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, 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 . 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 . 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 cm in size and the downstream target face as a sixth bin . 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 \NASixtyOnemeasurements with 2007 and 2009 data have already proven that the \NASixtyOnesetup is suitable for measuring hadron yields on the T2K replica target surface, although the target is placed outside of the tracking system.
2 NA61/SHINE experimental setup
The \NASixtyOnesetup during the T2K replica target data-taking 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 -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 Figure 1. More details can be found in Ref. . The coordinate system is defined as follows: the -axis is in the nominal direction of the beam, the -axis in the horizontal plane is such that positively charged particles are bent in the positive -direction, and the -axis is perpendicular to the horizontal plane and points upward. The origin is located in the centre of the VTPC-2.
The \NASixtyOnespectrometer  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 () down to . Secondary hadron beams of various momenta are produced by impinging a primary 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 \NASixtyOne. The beam is defined by three scintillator counters used in coincidence (, and ) and two scintillator counters with holes used in anticoincidence ( and ). The counter has a radius equal to the target radius, and it was placed cm upstream of the target. The counter has a 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  (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 beam contains around % , % , and % protons . The estimated purity of the selected proton sample is better than %. Finally, the beam position is measured by a set of six multi-wire proportional chambers. Two chambers, one measuring position along the -axis and another measuring the position along the -axis, are part of one BPD. The precision of the position measurement in one BPD is around .
The third BPD was mounted on the same support as the target, just upstream of the counter.
The T2K replica target was made of Toyo Tanso IG-43 graphite  which was also used for the original T2K target. The target density was estimated by measuring its volume and mass, and it is equal to g cm. This value is % higher than the density of the T2K target, but well within the measurement uncertainty. The target shape was machined to form a cm thick disk with an cm long rod coming from its centre. The radii of the disk and rod are and 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 counter used in the trigger. A schematic overview of the target can be seen in Figure 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 cm from the upstream side of VTPC-1.
2.3 Tracking system
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 . The GTPC was introduced in order to measure forward-going particles with small polar angles, mrad, in the laboratory frame. The magnetic field was set so that the total bending power of the magnets was T m. This allowed a momentum resolution of (), except for very forward tracks that only passed through the GTPC and MTPCs and had a momentum resolution of () . The measured energy loss () provided excellent particle identification capabilities; the achieved resolution is around %. In the region where the energy loss distributions for different particles cross, time-of-flight measurements () were used for particle identification.
2.4 Time-of-flight walls
Although the \NASixtyOnespectrometer 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 scintillator bars with dimensions WHL = cm arranged in ten separate modules. The signal is read out by two PMTs placed at both ends of the bars. The estimated resolution was ps .
The proton beam profile on target at NA61/SHINE is wider than that at J-PARC. Previous studies  have shown that longitudinal () 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 along the target with polar angle and radial position with respect to the target longitudinal axis, would escape the target at a minimal longitudinal position, , given by:
Evidently, for narrower beams, the longitudinal () 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 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 %. The \NASixtyOneproton beam achieved in 2009  is much wider than the J-PARC beam. To avoid large differences between the \NASixtyOnebeam and the J-PARC beam, four special triggers were used during the 2010 data-taking (see Figure 1):
The trigger is the only trigger that selects beam protons that did not necessarily hit the target or when a hit in the scintillator was not detected. The and triggers select only beam protons that hit the target. The only difference is that selects a narrower beam because the counter has a smaller hole. The fourth trigger selects all beam particles that hit the target. It is important to note that , and were prescaled - only a fraction of them was recorded. The prescaling was applied to to achieve the desired radial distribution of the beam. If an event is a but not a , 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 Figure 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 Section 6.
In total, triggers were recorded. Around % 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 secondary hadrons below 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 (% of triggers) were analysed, while the high magnetic field subset was used for alignment between the BPDs and TPCs (see Subsection 3.1).
The analysis of the data taken with the standard magnetic field configuration was performed with the so-called - method, which makes use of the TPC measurements of the specific energy loss () and the time-of-flight measurements [9, 10]. The specific energy loss of a track is calculated as a truncated mean  of the charges of the clusters (points) on the track traversing the TPCs. After event and track selection, selected tracks were binned in momentum (), polar angle (), both measured in the laboratory frame, and longitudinal position of the exit point on the target surface (), as required by T2K. The different hadron yields were estimated by fitting a two-dimensional - 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.
3.1 Target position and alignment
The target upstream face and positions with respect to the BPDs can be determined just by plotting the beam profile at the target upstream position. The position of the target centre was obtained by maximising the number of entries in the circle with radius equal to cm. This was only possible because 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 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. ) with error propagation that accounts for multiple scattering (see Ref. ). The positions of the beam tracks extrapolated forward from the BPDs and 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.
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 position. The distribution of the points of closest approach was plotted and a rising edge could be seen as shown in Fig. 5. The position at the half-max 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 and planes. First, the target was assumed to be parallel to the -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 and distributions for slices in 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 direction, edges of the distributions would be shifted towards the centre because the target is moved from the beam. The shift would increase with . No appreciable tilt was found. In Table 1 the measured position of the target is summarized with uncertainties for all parameters.
|[cm ]||[cm ]||[cm ]||[mrad ]||[mrad ]|
3.2 Event selection
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 or 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 % 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 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 kHz which means that the mean time difference between beam particles was around . However, the ToF-F wall has a time acquisition window of ns, and because of the requirement that all tracks must have a hit for particle identification, most of the off-time TPC tracks were automatically discarded. The fraction of accepted off-time tracks was negligible.
3.3 Track selection
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 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 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 mrad or 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 , where 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.
3.4 Phase space
Selected tracks were binned in () space. As required by T2K , 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 () binning was applied for extracting , 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 were used to reduce the large variability in the statistical uncertainty for low and high momenta. The appropriate size of the bins was estimated from Fluka2011.c. simulations [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 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 region would be heavily biased. For example, momentum bins for protons must start at .
|[%]||[%]||[%]||[%]||p [%]||Total [%]|
The overall phase space binning covers more than of (anti)neutrinos crossing the T2K far detector (Super-Kamiokande) which are produced by the charged hadrons emitted from the target. More than 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 - 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 anti-neutrino 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 \NASixtyOneadded new forward TPCs in 2017 , possible future measurements could improve coverage for the forward-going high momentum and p .
3.5 Particle identification
Since large fractions of the phase space are covered, a robust particle identification method is needed. For the momentum range between to (see Fig. 12), the energy loss distributions cross, hence particle identification based on the energy loss alone is not possible. However, in these regions, measurements can be used to distinguish between particles. Kaons can be separated easily up to , while protons can be separated up to , as can be seen in Fig. 13. For higher momentum, the 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 measurement was used to calculate the particle mass squared () and it was combined with the energy loss. Particles were represented by islands in the - space. Therefore, raw hadron yields could be obtained by fitting an appropriate function to the distribution for each phase space bin. Both, and were assumed to be normally distributed. This assumption for must be closely examined. In general, is not normally distributed, as its distribution is similar to the Landau distribution. The measured in NA61/SHINE is calculated as a truncated mean (see Ref. ) 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 same number of TPC clusters. The selected tracks can have a different number of clusters which can significantly affect the resolution if the number of TPC clusters is small. However, the resolution in \NASixtyOnesaturates at around % for tracks with more than clusters. Around % of the selected tracks have more than 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 (, , , p ()). The fitting was done in the RooFit framework  by using extended log-likelihood minimization which treats the number of observed particles as a Poisson random variable. An example of a fit is shown in Fig. 14. In total, there were 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.
3.6 Monte Carlo correction factors
Interactions in the target were simulated by Fluka2011.c. [19, 20, 21]. The kinematics of the particles at their exit position on the surface of the target was then fed to Geant3 GCALOR  for propagation through the \NASixtyOnedetector. 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 - and - positions and divergences were considered to be independent. In total, nearly protons on target (POT) were simulated.
A simple Monte Carlo correction was applied to the raw yields:
where is the number of simulated hadrons in the bin ()=() and 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 is the geometrical efficiency (because of the cut), is the hadron loss efficiency due to decays and interactions in the detector, is the reconstruction efficiency, is the selection efficiency, is the bin migration efficiency and is the feed-down efficiency. The geometrical efficiency gives the percentage of solid angle covered by the \NASixtyOnedetector 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 (), polar angle () and longitudinal position along the target surface () for each hadron are reconstructed with finite precision. Consequently, some of the hadrons will be placed in a wrong ) bin. The ratio of the number of hadrons which are reconstructed in a given bin to the number of hadrons which are simulated in the same bin represents the migration efficiency. The feed-down 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.
3.7 ToF-F efficiency correction factors
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. 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 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 . 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.
3.8 Ad hoc correction factors
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 (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 passing 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 corrections in the TPCs. Nevertheless, the affected bins represent only around 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 hit in the Monte Carlo (data), is the number of selected tracks in the Monte Carlo(data) without the requirement, and 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.
4 Systematic Uncertainties
The systematic uncertainties of this analysis were carefully studied. A brief summary is given below, while more detailed information can be found in Ref. .
The following systematic uncertainties were considered:
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 . However, for tracks with a momentum of several in the mrad region in the most upstream longitudinal bin, the hadron loss uncertainty goes up to (for a single bin it reached ). These tracks cross the beamline and pass very close to the walls of the MTPCs.
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 . Also, high angle tracks extrapolated toward the downstream face of the target are sensitive to the and shifts of the target, and for them the bin migration uncertainty is up to . For all other regions of the phase space, the bin migration uncertainty drops below .
Reconstruction efficiency uncertainty. This uncertainty was estimated in the previous \NASixtyOnepublications [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 is assigned for all bins.
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 was observed for and p. A constant value of 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 distribution for momentum larger than . For this reason, the uncertainty for can reach for the higher momentum bins. In case of , the obtained uncertainties are similar to the ones for .
Feed-down uncertainty. Produced 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 and is model dependent. A systematic uncertainty of on the feed-down correction factor was assigned, following the standard approach in Refs. [8, 10]. The systematic uncertainty is up to for and protons. Mostly low momentum bins in the first and last longitudinal bins are affected.
The ToF-F efficiency uncertainty. The uncertainty of the ToF-F efficiency was taken as a systematic uncertainty. It goes up to , but for most of the bins it is well below .
Ad-hoc uncertainty. A constant conservative uncertainty of 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 ().
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 () for yields comes from the reconstruction efficiency uncertainty, while the dominant contribution for comes from the PID uncertainty.
5 Results and comparison with hadron production models
Results are presented in the form of double differential yields normalised by the total number of protons on target:
where () are the () bin numbers, is the number of protons on target, is the hadron species, is the number of extracted hadrons in a given bin, is the correction factor, is the momentum bin size, and is the polar angle bin size. Uncertainties on the hadron yields are shown in Figs. 17 - 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 (mrad) and high (mrad) polar angle region. While for 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. .
5.1 Comparisons with models
The results were compared with NuBeam and QGSP_BERT physics lists from Geant4.10.03 [27, 28]. Detailed comparisons are presented in Figs. 23 - 28 for yields, in Figs. 29 - 34 for yields, in Figs. 35 - 40 for yields, in Figs. 41 - 46 for yields, and in Figs. 47 - 52 for p yields. The NuBeam physics list provides better predictions of , , , and yields. However, both physics lists fail to accurately predict p yields. Similar behaviour was noticed in the previous \NASixtyOnemeasurements of primary proton-carbon interactions at , where it was shown that most of the models poorly reproduce proton yields. The earlier observation  that a better agreement with Fluka2011 predictions can be obtained by lowering the production cross section in Fluka2011 is confirmed.
5.2 Comparison with previous \NASixtyOnemeasurements
As previously mentioned, \NASixtyOneprovided measurements of the yields emitted from the T2K replica target using data collected in 2009 . 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 (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. . However, under the assumption that major differences are only due to geometry - i.e. low angle secondaries cannot exit the target from the first longitudinal bin if the beam is very narrow, one can apply scaling to the radial beam distribution and achieve a similar radial distribution as in the previous measurement. The radial beam distribution was divided into bins, mm in size. A simple weight was applied to all selected tracks:
where is the number of beam protons hitting the target in the radial bin in data from and is the total number of beam protons hitting the target. After scaling, the particle identification procedure was applied in the phase space binning as defined in Ref. . Also, Monte Carlo and 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 and 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 . 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 used in the analysis of data collected in 2009  the differential yields of negatively charged hadrons () were extracted. The analysis was done with the standard selection and reconstruction cuts. The quantities computed with the measured data were corrected with a corresponding Monte Carlo correction. The yields were computed with the most appropriate and trigger definitions. The resulting yields were compared with the yields obtained within the analysis framework described above. Such a comparison is reasonable since at the beam momentum of 31 the negatively-charged hadrons are strongly dominated by negatively-charged pions. It was concluded that both approaches provide yields in acceptable agreement. Some minor differences were observed only for the first longitudinal target bin, for small momentum (below 6 GeV/) and polar angles ( below 40 mrad), and in the interval [300; 400] mrad for the first three longitudinal target bins.
6 Dependence of the T2K re-weighting factors on the proton beam profile
The results presented in Section 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 \NASixtyOnethin-target measurements are presented in Ref. .
The weights are simply calculated as the ratio of the data yields and the Fluka2011 yields for each 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 \NASixtyOne. As previously mentioned, hadron 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 re-weighting 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 protons on target were simulated for \NASixtyOne and 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 bins, while for other bins any differences are mostly below . It is important to note that the differences between the and beam profiles are very large in comparison to differences between the and T2K beam profiles.
Thus, one can use the measured \NASixtyOnedouble differential hadron yields in the T2K flux simulation in the following way:
Run the simulation of interactions inside the target for the \NASixtyOne beam profile.
Calculate the re-weighting factors.
Run the simulation with the T2K beam profile.
Re-weight the new simulation using the previously calculated re-weighting factors.
7 Summary and conclusions
Measurements of the double differential yields of , , , and p emitted from the surface of a -cm-long carbon target (T2K replica) with incoming protons were performed by the \NASixtyOneexperiment at the CERN SPS. Yields of and were measured with improved precision compared to the previously published \NASixtyOneresults, while , 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)neutrino-nucleus 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.
The results were also compared with the NuBeam and QGSP_BERT physics lists from Geant4.10.03 and show that none of the models predicts accurately the measured proton yields. This information can be used by model builders to improve hadronic Monte Carlo generators.
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), 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.
-  N. Abgrall et al., [NA61 collaboration Collab.] JINST 9 (2014) P06005, arXiv:1401.4699 [physics.ins-det].
-  K. Abe et al., [T2K Collaboration Collab.] Nucl. Instrum. Meth. A659 (2011) 106–135, arXiv:1106.1238 [physics.ins-det].
-  “J-parc tdr, kek-report 2002-13 and jaeri-tech 2003-044,” 2003. http://hadron.kek.jp/ accelerator/TDA/tdr2003/index2.html.
-  K. Abe et al., [T2K Collaboration Collab.] Phys. Rev. D87 (2013) 012001, arXiv:1211.0469 [hep-ex].
-  N. Abgrall et al., [NA61/SHINE Collaboration Collab.] Phys. Rev. C89 (2014) 025205, arXiv:1309.1997 [physics.acc-ph].
-  N. Abgrall et al., [NA61/SHINE Collaboration Collab.] Phys. Rev. C84 (2011) 034604, arXiv:1102.0983 [hep-ex].
-  N. Abgrall et al., [NA61/SHINE Collaboration Collab.] Phys. Rev. C85 (2012) 035210, arXiv:1112.0150 [hep-ex].
-  N. Abgrall et al., [NA61/SHINE Collab.] Eur. Phys. J. C76 no. 2, (2016) 84, arXiv:1510.02703 [hep-ex].
-  N. Abgrall et al., [NA61/SHINE Collaboration Collab.] Nucl. Instrum. Meth. A701 (2013) 99–114, arXiv:1207.2114 [hep-ex].
-  N. Abgrall et al., [NA61/SHINE Collab.] Eur. Phys. J. C76 no. 11, (2016) 617, arXiv:1603.06774 [hep-ex].
-  A. Hasler Ph.D. Thesis, University of Geneva (2015) . https://cds.cern.ch/record/2039148. Presented 22 Jun 2015.
-  L. Zambelli, A. Fiorentini, and T. Vladisavljevic, [T2K, NA61/SHINE Collab.] J. Phys. Conf. Ser. 888 no. 1, (2017) 012067.
-  M. Pavin, “Measurements of hadron yields from the T2K replica target in the NA61/SHINE experiment for neutrino flux prediction in T2K,” 2017. https://cds.cern.ch/record/2292904. Ph.D. Thesis, University of Paris VI.
-  C. Bovet et al. CERN-YELLOW-82-13.
-  “Toyotanso graphite.” http://www.toyotanso.co.jp/Products/Special_graphite/data_en.html.
-  S. Afanasiev et al., [NA49 Collaboration Collab.] Nucl.Instrum.Meth. A430 (1999) 210–244.
-  S. Gorbunov and I. Kisel Nucl. Instrum. Meth. A559 (2006) 148–152.
-  E. J. Wolin and L. L. Ho Nucl. Instrum. Meth. A329 (1993) 493–500.
-  G. Battistoni et al. AIP Conf. Proc. 896 (2007) 31–49.
-  A. Ferrari, P. R. Sala, A. Fasso, and J. Ranft CERN-2005-010, SLAC-R-773, INFN-TC-05-11 .
-  T. Böhlen et al. Nuclear Data Sheets 120, 211-214 (2014) .
-  A. Aduszkiewicz et al., [NA61/SHINE Collaboration Collab.] CERN-SPSC-2017-038 ; SPSC-SR-221 (2017) .
-  A. Aduszkiewicz et al., [NA61/SHINE Collaboration Collab.] CERN-SPSC-2018-008 ; SPSC-P-330-ADD-10 (2018) .
-  W. Verkerke and D. P. Kirkby eConf C0303241 (2003) MOLT007, arXiv:physics/0306116 [physics]. [,186(2003)].
-  R. Brun, F. Bruyant, M. Maire, A. McPherson, and P. Zanarini CERN-DD-EE-84-1 (1987) .
-  “Tables with numerical results for the t2k replica target 2010 paper.” https://edms.cern.ch/document/1828979/1.
-  S. Agostinelli et al., [GEANT4 Collaboration Collab.] Nucl. Instrum. Meth. A506 (2003) 250–303.
-  J. Allison et al., [GEANT4 Collaboration Collab.] IEEE Transactions on Nuclear Science 53 no. 1, (2006) 270–278.