Measurements of hadron production in π + + C and π + + Be interactions at 60 GeV = c

Precise knowledge of hadron production rates in the generation of neutrino beams is necessary for accelerator-based neutrino experiments to achieve their physics goals. NA61/SHINE, a large-acceptance hadron spectrometer, has recorded hadron þ nucleus interactions relevant to ongoing and future long-baseline neutrino experiments at Fermi National Accelerator Laboratory. This paper presents three analyses of interactions of 60 GeV =c π þ with thin, fixed carbon and beryllium targets. Integrated production and inelastic cross sections were measured for both of these reactions. In an analysis of strange, neutral hadron production, differential production multiplicities of K 0 S , Λ and ¯ Λ were measured. Lastly, in an analysis of charged hadron production, differential production multiplicities of π þ , π − , K þ , K − and protons were measured. These measurements will enable long-baseline neutrino experiments to better constrain predictions of their neutrino flux in order to achieve better precision on their neutrino cross section and oscillation measurements.


I. INTRODUCTION
The NA61 or SPS Heavy Ion and Neutrino Experiment (SHINE) [1] has a broad physics program that includes heavy ion physics, cosmic ray physics and neutrino physics. Accelerator-generated neutrino beams rely on beams of high energy protons which are directed towards a fixed target. The interactions of these protons result in secondary hadrons (especially pion, kaons, protons, neutrons and lambdas), some of which decay to produce the beam of neutrinos. As most neutrino beam lines use targets that are an interaction length or longer in length, many of the secondary hadrons can reinteract inside the target and other beam material (such as the decay pipe walls or material of the focusing horns). Thus, it is important to have accurate knowledge of not only the primary proton interactions in the target, but also of the reinteractions of secondary particles.
NA61/SHINE has previously measured hadron production in interactions of 31 GeV=c protons with a thin carbon target for the benefit of the T2K experiment [2][3][4][5]. The NA61/SHINE experiment is also well suited to making measurements of the beam line interactions that dominate the neutrino production in the Fermilab long-baseline accelerator neutrino program, including the existing NuMI beam [6], which is initiated by 120 GeV=c primary protons, and the proposed Long-Baseline Neutrino Facility (LBNF) beam line [7] that will supply neutrinos for the Deep Underground Neutrino Experiment (DUNE) [8], which will use 60-120 GeV=c primary protons. The current optimized beam line design for LBNF features a ∼2.2 m-long graphite target [9], but beryllium and hybrid targets have been considered as well.
In DUNE, near the oscillation peak at a neutrino energy of 3 GeV, roughly half of the neutrinos are produced from the decays of secondary particles generated in the interactions of primary protons (p → X → ν) [10]. The other half come from the decays of particles generated by the reinteractions of protons or hadrons (e.g., p → X → Y → ν). For the LBNF optimized beam, each neutrino in the near detector results from an average of 1.8 interactions in the beam line (including the interaction of the primary proton) [11]. After protons, the largest source of these interactions is pions with an average of 0.2 pion interactions contributing to each neutrino, and these pions typically have momenta in the range from roughly 10 to 70 GeV=c.
The current estimates of the flux uncertainties in DUNE [11] near the oscillation maximum are dominated by uncertainties on existing p þ C measurements such as those described in Ref. [12], proton and neutron interactions that are not covered by existing data and uncertainties on the reinteractions of pions and kaons. NA61/ SHINE seeks to improve on these uncertainties by making improved measurements of proton interactions with neutrino target materials (with more phase space coverage and larger statistics) and by making measurements of meson interactions with target and beam line materials. With the exception of the HARP measurements [13], there is little existing data on the particle production spectra from interactions of mesons in the incident momentum range of interest for long-baseline neutrino experiments. This paper presents new results on the yields of particles resulting from the interactions of 60 GeV=c π þ on carbon and beryllium targets recorded in 2016.
Three types of results are presented in this paper. Section IV presents measurements of the integrated production and inelastic cross sections for π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions, and describes the uncertainties on these measurements. Section V describes measurements of the differential multiplicity of neutral hadrons (K 0 S , Λ andΛ) produced in these interactions, in bins of the momentum and angle of the produced hadron. Section VI describes measurements of the differential multiplicity of the charged hadrons (π þ , π − , K þ , K − and p) in bins of the momentum and angle of the produced hadron. Section VII describes the systematic uncertainties on the results presented in Secs. V and VI.

II. DETECTOR SETUP
Located on a secondary beam line of CERN's Super Proton Synchrotron (SPS), NA61/SHINE probes the interactions of protons, pions, kaons and heavy ions with fixed targets. The 400 GeV=c primary protons from the SPS beam strike a target 535 m upstream of NA61/SHINE, generating the secondary beam. A system of magnets selects the desired beam momentum. Unwanted positrons and electrons are absorbed by a 4-mm lead absorber.
The NA61/SHINE detector [1] is shown in Fig. 1. In the 2016 operation configuration, the detector comprises four large Time Projection Chambers (TPCs) and a Time of Flight (ToF) system allowing NA61/SHINE to make spectral measurements of produced hadrons. Two of the TPCs, vertex TPC 1 (VTPC-1) and vertex TPC 2 (VTPC-2), are located inside superconducting magnets, capable of generating a combined maximum bending power of 9 T · m. Downstream of the VTPCs are the main TPC left (MTPC-L) and main TPC right (MTPC-R). Additionally, a smaller TPC, the gap TPC (GTPC), is positioned along the beam axis between the two VTPCs. Two side time-offlight walls, ToF-left and ToF-right, walls were present. Notably, the previously used ToF-forward wall was not installed during the 2016 operation. The Projectile Spectator Detector (PSD), a forward hadron calorimeter, sits downstream of the ToF system.
The NA61/SHINE trigger system uses two scintillator counters (S1 and S2) to trigger on beam particles. The S1 counter provides the start time for all counters. Two veto scintillation counters (V0 and V1), each with a hole aligned to the beam, are used to remove divergent beam particles upstream of the target. The S4 scintillator with a 1 cm radius (corresponding to a particle scattering off of the target at an angle of 2.7 mrad) sits downstream of the target and is used to determine whether or not an interaction has occurred. A Cherenkov differential counter with achromatic ring focus (CEDAR) [14,15] identifies beam particles of the desired species. The CEDAR focuses the Cherenkov ring from a beam particle onto a ring of eight photomultiplier tubes (PMTs). The pressure is set to a fixed value so that only particles of the desired species trigger the PMTs, and typically, a coincidence of at least six PMTs is required to tag a particle for the trigger.
The beam particles are selected by defining the beam trigger (T beam ) as the coincidence of S1 ∧ S2 ∧ V0 ∧ V1 ∧ CEDAR. The interaction trigger (T int ) is defined by the coincidence of T beam ∧ S4 to select beam particles which have interacted with the target. A correction factor is discussed in detail in Sec. IVA to correct for interactions that result in an S4 hit. Three beam position detectors (BPDs), which are proportional wire chambers, are located 30.39, 9.09 and 0.89 m upstream of the target and determine the location of the incident beam particle to an accuracy of ∼100 μm.
Interactions of π þ beams were measured on thin carbon and beryllium targets. The carbon target was composed of graphite of density ρ ¼ 1.80 g=cm 3 with dimensions of 25 mm (W) x 25 mm (H) x 14.8 mm (L), corresponding to roughly 3.1% of a proton-nuclear interaction length. The beryllium target had a density of ρ ¼ 1.85 g=cm 3 with dimensions of 25 mm (W) x 25 mm (H) x 14.9 mm (L), corresponding to roughly 3.5% of a proton-nuclear interaction length. The uncertainties in the densities of the targets were found to be 0.69% for the carbon target and 0.19% for the beryllium target.

III. EVENT SELECTION
Several cuts were applied to events to ensure the purity of the samples and to control the systematic effects caused by beam divergence. The same event cuts are used for the integrated cross section and differential cross section analyses in order to ensure that the normalization constants obtained from the integrated cross section analysis are valid for calculating multiplicities in the differential cross section analyses. First, the so-called wave form analyzer (WFA) cut was used to remove events in which multiple beam particles pass through the beam line in a small time frame. The WFA determines the timing of beam particles that pass through the S1 scintillator. If another beam particle passes through the beam line close in time to the triggered beam particle, it could cause a false trigger in the S4 scintillator and off-time tracks being reconstructed in the main interaction vertex. To mitigate these effects, a WFA cut of AE2 μs is used.
The measurements from the BPDs are important for estimating the effects of beam divergence on the integrated cross section measurements. To mitigate these effects, tracks are fitted to the reconstructed BPD clusters, and these tracks are extrapolated to the S4 plane. The so-called "good BPD" cut requires that each event includes a cluster in the most downstream BPD and that a track was successfully fit to the BPDs. Figure 2 shows the resulting BPD extrapolation to the S4 plane for the 60 GeV=c π þ beam. A radial cut was applied to the BPD tracks extrapolated to the S4, indicated by the red circles on Fig. 2, in order to ensure that noninteracting beam particles strike the S4 counter. This corresponds to a trajectory within 0.7 cm of the S4 center (compared to the S4 radius of 1 cm). It can be seen from these distributions that the beam, veto counters and S4 were well aligned during the data taking.
To begin the event selection, only unbiased T beam events are considered for the integrated cross section analysis. For the integrated cross section analysis of the π þ þ C at 60 GeV=c (π þ þBe at 60 GeV=c) data set, 191,099 (116,944) targetinserted and 86,022 (58,551) target-removed events were analyzed after the described selection. For the analysis of spectra, only T int events are considered. For the spectra analysis of the π þ þ C at 60 GeV=c (π þ þ Be at 60 GeV=c) data set, 1,496,524 (1,096,003) target-inserted and 86,764 (57,045) target-removed events were selected.

IV. INTEGRATED INELASTIC AND PRODUCTION CROSS SECTION ANALYSIS
The total integrated cross section of hadron þ nucleus interactions, σ tot , can be defined as the sum of the FIG. 2. Positions of BPD tracks extrapolated to the S4 plane in target-removed data runs from the π þ þ C at 60 GeV=c data set. The measured S4 position is shown as a black circle and the BPD radius cut is shown as a red circle in both figures. Left: Events taken by the beam trigger. Right: Events taken by the interaction trigger.
inelastic cross section, σ inel , and the coherent elastic cross section, σ el , Coherent elastic scattering leaves the nucleus intact. The sum of all other processes due to strong interactions makes up the inelastic cross section. The inelastic cross section can be divided into the production cross section, σ prod , and the quasielastic cross section, σ qe , In this paper, production interactions are defined as processes in which new hadrons are produced. Quasielastic interactions include processes other than coherent elastic interactions in which no new hadrons are produced, mainly fragmentation of the nucleus. In this paper, measurements of the production cross section, σ prod , and inelastic cross section, σ inel , are presented for π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions. These cross section measurements are important for accelerator-based neutrino experiments and are needed to normalize the differential cross section yields that are discussed in Secs. Vand VI. This analysis closely follows the method described in Ref. [16], but with some differences, which are discussed below.

A. Trigger cross section
For sufficiently thin targets, the probability P of a beam particle interacting is approximately proportional to the thickness, L, of the target, the number density of the target nuclei, n, and the interaction cross section, σ, The density of nuclei can be written in terms of Avogadro's number, N A , the material density, ρ, and the atomic mass, m a , The counts of beam (T beam ) and interaction triggers (T int ) that pass the event selection can be used to estimate the trigger probability with the target inserted (I) and with the target removed (R), Figure 3 shows an example of the trigger probabilities for each run for the π þ þ C at 60 GeV=c data set. The targetremoved runs were interspersed throughout the targetinserted data runs to ensure they represented comparable beam conditions. The trigger rates show consistency over the course of the runs, which were recorded over a period of about three days. Table I gives the trigger probabilities for both the target-inserted and target-removed samples of the π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c data sets. Taking into account the trigger probabilities with the target inserted and the target removed, P I T and P R T , the corrected trigger probability, P trig , can be obtained, FIG. 3. Trigger interaction probabilities for the π þ þ C at 60 GeV=c data set for target-inserted and target-removed runs.
Analogous to Eq. (3), the trigger cross section σ trig is defined as where the beam attenuation is taken into account by replacing L with L eff . The effective target length can be calculated using the absorption length, λ abs , where By combining Eqs. (7)-(9), σ trig can be rewritten as

B. S4 correction factors
The trigger cross section takes into account the interactions where the resulting particles miss the S4 scintillator. But even when there has been a production or quasielastic interaction in the target, there is a possibility that a forwardgoing particle will strike the S4 counter. Moreover, not all elastically scattered beam particles strike the S4. The trigger cross section must be corrected to account for these effects. Combining Eqs. (1) and (2), the trigger cross section can be related to the production cross section through Monte Carlo (MC) correction factors as follows: where f prod , f qe and f el are the fractions of production, quasielastic and elastic events that miss the S4 counter. The cross sections σ qe and σ el are also estimated from MC simulations. Equation (11) can be rewritten to obtain σ prod and σ inel as and A GEANT4 detector simulation [17][18][19] using GEANT4 version 10.4 with physics list FTFP_BERT was used to estimate the MC correction factors discussed above. The MC correction factors obtained for π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions are presented in Table II.

C. Beam composition
For the analyses of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions recorded in 2016, the beam composition could be constrained better than in the analysis of interactions recorded in 2015 by NA61/SHINE as discussed in [16]. Simulations of the H2 beam line show that the population of muons in the 60 GeV=c secondary hadron beam used to record these interactions is at the level of 1.5 AE 0.5% [20]. Nearly all of the muons come from decays of 60 GeV=c pions, so they have a minimum energy of 34 GeV=c. GEANT4 simulations were run to estimate the target-inserted and target-removed trigger rates due to muons, P I μ and P R μ . These simulations took the momentum distribution of muons into account. Additional H2 beam line simulations were run to more precisely estimate the level of positron contamination in the beam [21]. A conservative estimate of 0.5% AE 0.5% was attributed to this contamination. The trigger rates due to positrons, P I e and P R e , were also estimated with GEANT4 simulations. The effect of muon and positron contamination on the trigger cross section was estimated as follows: where f e ¼ 0.005, f μ ¼ 0.015 and f π ¼ 0.98. The resulting corrections applied to σ prod (σ inel ) were þ0.3% (þ0.3%) for π þ þ C at 60 GeV=c and þ1.1% (1.0%) for π þ þ Be at 60 GeV=c.

D. Systematic uncertainties
The integrated cross section results were evaluated for a number of possible systematic effects. The sources of TABLE I. This table presents the observed trigger interaction probabilities for both the target-inserted and target-removed samples of the π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c data sets.
Interaction TABLE II. Monte Carlo correction factors obtained for analyzing π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions.

Monte Carlo correction factors
Interaction uncertainty having a non-negligible effect on the results are the uncertainty in the density of the target, the uncertainty in the S4 size, the uncertainty on the beam composition and uncertainties on the S4 correction factors. The procedures used to evaluate these sources of systematic uncertainties were discussed in [16], so they are not discussed here.

Breakdowns of the integrated cross section uncertainties
The target density uncertainties, S4 size uncertainties, beam composition uncertainties and S4 correction factor uncertainties associated with the production and inelastic cross sections measurements for π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions are presented in Tables III and IV.

E. Integrated cross section results
Measurements of production cross sections for π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c are summarized in Table V along with statistical, systematic and physics model uncertainties. The production cross section of π þ þ C at 60 GeV=c interactions was found to be 166.7 mb, and the production cross section of π þ þ Be at 60 GeV=c interactions was found to be 140.6 mb. The result obtained for interactions of π þ þ C at 60 GeV=c with these 2016 data was lower compared to the result obtained with the 2015 data [16], but it is within the estimated uncertainty. Reasons for this difference could be due to the difference in the detector setup, the different target used and statistical fluctuations. These results, the results obtained by NA61/ SHINE from data recorded in 2015 and the measurements of Carroll et al. [22] are compared in Fig. 4.
The measurements of inelastic cross sections for π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c are summarized in Table VI along with statistical, systematic and physics model uncertainties. The inelastic cross section of π þ þ C at 60 GeV=c was found to be 182.7 mb, and the inelastic cross section of π þ þ Be at 60 GeV=c was found to be 154.4 mb. Again, the result obtained for interactions of π þ þ C at 60 GeV=c with these 2016 data was lower TABLE III. Breakdown of systematic uncertainties for the production cross section measurements of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions.
Systematic uncertainties for σ prod (mb) Interaction p (GeV=c) Density S4 size Beam purity MC statistical Total systematic uncertainties Model uncertainties TABLE IV. Breakdown of systematic uncertainties for the inelastic cross section measurements of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions.
Systematic uncertainties for σ inel (mb) Interaction p ( GeV=c) Density S4 size Beam purity MC statistical Total systematic uncertainties Model uncertainties Production cross section measurements of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions are presented.
The central values as well as the statistical (Δ stat ), systematic (Δ syst ) and model (Δ model ) uncertainties are shown. The total uncertainties (Δ total ) are the sum of the statistical, systematic and model uncertainties in quadrature.
Production cross section (mb) compared to the result obtained with the 2015 data [16], but it is within the estimated uncertainty. These results, the results obtained by NA61/SHINE from data recorded in 2015 and the measurements of Denisov et al. [23] are compared in Fig. 5.

V. ANALYSIS OF NEUTRAL HADRON SPECTRA
NA61/SHINE is able to identify a number of species of weakly decaying neutral hadrons by tracking their charged decay products. The simplest decay topology NA61/ SHINE can identify is the V 0 topology. This topology refers to track topologies in which an unobserved neutral particle decays into two child particles, one positively charged and one negatively charged, observed by the tracking system. This paper presents differential production cross section measurements of produced K 0 S , Λ andΛ in interactions of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c using a V 0 analysis.

A. Selection of V 0 candidates
To start with, every pair of one positively charged and one negatively charged track with a distance-of-closest approach less than 5 cm is considered as a V 0 candidate. Of course, many of these V 0 candidates are not true V 0 s. For example, a V 0 candidate might consist of two tracks that come from the main interaction point, the child tracks might come from two different vertices or the child tracks might come from a parent track, which is not a neutral particle. Additionally, photons converting to e þ e − pairs make up part of the V 0 sample.

Topological cuts
The topological cuts are designed to reduce the number of false V 0 s in the collection of V 0 candidates and to remove V 0 candidates that have poorly fitted track variables. Only V 0 candidates that have a reconstructed V 0 vertex downstream of the target are considered.
The second topological selection is the requirement that both child tracks have at least 20 reconstructed TPC clusters and that at least ten of those clusters belong to the VTPCs. This cut ensures that the reconstructed kinematics of the decay are reliable.
The third topological cut is the impact parameter cut, which removes many false V 0 candidates. This selection allows an impact parameter from between the extrapolated V 0 s track and the main interaction vertex of up to 4 cm in the x dimension and up to 2 cm in the y dimension.

Purity cuts
The purity cuts are designed to separate the desired neutral hadron species from other neutral species, as well as to remove additional false V 0 candidates. The first two purity cuts are applied in the same way to K 0 S , Λ andΛ. This first selection requires the reconstructed z position of the V 0 vertex to be at least 3.5 cm downstream of the target center. This cut removes many of the V 0 candidates coming from the main interaction vertex and neutral species that decay more quickly than K 0 S , Λ orΛ. Photons undergoing pair production (γ → e þ e − ) are present in the V 0 sample. Because the photon is massless, the transverse momentum of the decay is In order to remove most of these photons from the sample, the second purity cut requires a p T > 0.03 GeV=c.
3. Purity cuts for the selection of K 0 S At this point, it is necessary to assume a decay hypothesis. For K 0 S , the hypothesis is K 0 S → π þ π − . Therefore, it is assumed that the V 0 particle has a mass of m K 0 S ¼ 0.498 GeV=c 2 and the child particles have a mass of m π AE ¼ 0.140 GeV=c 2 [24].
To remove Λ andΛ from the K 0 S sample, cuts on the angles that the child particle tracks make with the V 0 track FIG. 5. Summary of inelastic cross section measurements. The results are compared to previous results from NA61/SHINE [16] and Denisov et al. [23]. VI. Inelastic cross section measurements of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions are presented.
The central values as well as the statistical (Δ stat ), systematic (Δ syst ) and model (Δ model ) uncertainties are shown. The total uncertainties (Δ total ) are the sum of the statistical, systematic and model uncertainties in quadrature. in the decay frame are applied to the sample. These angles are represented in Fig. 6. In order to remove Λ, cos θ þÃ < 0.8 is required and to removeΛ, cos θ −Ã < 0.8 is required. The next selection is an allowed range of the invariant mass. The invariant mass is calculated with the reconstructed momenta, assumed masses and energies of the child particles,

Inelastic cross section (mb)
The invariant mass range cut removes V 0 candidates with unreasonable values of M π þ π − , but is wide enough to allow a reliable fit to the background invariant mass distribution. For K 0 S , this range is chosen to be ½0.4; 0.65 GeV=c 2 . The final cut applied to the K 0 S selection is a cut on the proper decay length, cτ. The proper decay length can be calculated with the estimated momentum of the V 0 , p, the assumed mass, m, and the reconstructed length of the V 0 track, L, The purpose of this cut is to further reduce the number of false V 0 s and more quickly decaying neutral species. The chosen cut is cτ > 0.67 cm, which is a quarter of the proper decay length provided by the PDG [24], 2.68 cm.

Purity cuts for the selection of Λ andΛ
An invariant mass range cut and a proper decay length cut are used in the purity selection of Λ andΛ. The invariant mass hypothesis for the Λ decay is Λ → pπ − and the hypothesis for theΛ isΛ →pπ þ . An invariant mass range of ½1.09; 1.215 GeV=c 2 is used in both the Λ andΛ analyses.
A proper decay length cut is also applied to the Λ andΛ selection. The chosen cut is cτ > 1.97 cm, which is a quarter of the proper decay length given by the PDG [24], 7.89 cm.

Armenteros-Podolansky distributions
The effect of these selections on the V 0 candidates can be visualized with Armenteros-Podolansky distributions, which are distributions of α vs p T . The parameter α is the asymmetry in the longitudinal momenta of the child tracks with respect to the V 0 track, Figure 7 shows the V 0 candidates coming from π þ þ C at 60 GeV=c interactions before the V 0 selection cuts were applied and after the selection cuts were applied for the K 0 S , Λ andΛ analyses. It can be seen that the Λ andΛ candidates include part of the K 0 S spectra. These K 0 S are separated out from Λ andΛ during the fitting procedure discussed in the following section.

B. Fitting of invariant mass distributions
After applying the selection cuts for each particle species, the V 0 candidates are placed into the kinematic bins. For each of these kinematic bins, invariant mass distributions consist of both true K 0 S , Λ orΛ (signal) and the remaining background vertices. The objective of the fitting routine is to determine the number of true K 0 S , Λ andΛ in these invariant mass distributions. These fits are performed the same way on target-inserted and target-removed samples.

Signal model
In order to model the invariant mass distribution of K 0 S , Λ andΛ coming from the main interactions, template invariant mass distributions were derived from a GEANT4 MC production using the physics list FTFP_BERT. V 0 vertices are reconstructed, selected and binned in the same way as was done with the data. For each kinematic bin, MC templates are formed from the distributions of invariant mass from true K 0 S , Λ andΛ. These template distributions, g MC ðmÞ, are generated for both targetinserted and target-removed MC productions and were observed to peak at the known values of the K 0 S and Λ masses. In order to account for shifts in the invariant mass peaks and distortions of the signal shape due to misreconstruction of track variables and other possible effects, a mass shift, m 0 , and a smearing are applied to g MC ðmÞ. The smearing is applied by convolving g MC ðmÞ with a unit Gaussian distribution with width σ s . The parameters, m 0 and σ s , are allowed to vary for each kinematic bin and were observed to be small compared to the widths of the invariant mass distributions. The full signal distribution can be written as FIG. 6. This cartoon shows the relevant angles in V 0 decays in the rest frame of the V 0 . The child particles decay back to back in this frame. The angle at which the positively charged particle is emitted is θ þÃ , and the angle at which the negatively charged particle is emitted is θ −Ã .

Background model
It was observed that the shapes of the backgrounds in the invariant mass distributions vary among the K 0 S , Λ andΛ selection as well as among the kinematic bins. The background model was required to be flexible enough to account for the variation of background shapes in all of the kinematic bins for K 0 S , Λ andΛ. A second order polynomial was chosen to be used to fit the background distributions.

Fitting strategy
In order to fit for the signal and background contributions to the invariant mass distributions, a continuous loglikelihood function is constructed, to be part of the signal. The parameters, θ, include c s as well as the signal parameters, θ s , discussed in Sec. V B 1 and the background parameters, θ bg , which are the coefficients of the second degree polynomial. After obtaining c s from the fits, the raw yield of signal particles is calculated with y raw ¼ c s N V 0 Candidates . Figures 8 and 9 show example fits to K 0 S and Λ invariant mass distributions from the π þ þ C at 60 GeV=c data set. Averaging over the fit results for all kinematic bins, the observed K 0 S mass was 498.7 MeV=c 2 , which is slightly higher than the known value of 497.6 MeV=c 2 [24]. The average of the widths of the invariant mass distributions was observed to be 17 MeV=c 2 . The Λ andΛ masses were both observed to be 1; 117 MeV=c 2 , slightly higher than the known value of 1; 116 MeV=c 2 [24]. The widths of the Λ andΛ distributions were found to be 6 MeV=c 2 and 7 MeV=c 2 , respectively. These small discrepancies in the masses compared to the known values are likely due to small biases in the momentum reconstruction of tracks.

C. Corrections
The raw yields obtained from the fits discussed in the previous section must be corrected for systematic effects. These can roughly be categorized into several effects: branching ratio of the decay, detector acceptance, feeddown corrections, reconstruction efficiency and selection efficiency. The combined effect of these individual effects can be estimated as a single correction factor from Monte Carlo simulations. Using K 0 S as an example, the correction factor for kinematic bin i is given by The correction factors are calculated in the analogous way for Λ andΛ. The correction factors are obtained from the MC production using the FTFP_BERT physics list.

VI. ANALYSIS OF CHARGED HADRON SPECTRA
The analysis of produced charged hadrons is performed with a dE/dx analysis, which uses energy loss measured by the TPCs to separate particle species for both positively and negatively charged tracks. In particular, it was possible to measure spectra of produced π þ , π − , K þ , K − and protons with this method. Compared to past analyses of interactions of 31 GeV=c protons with a thin carbon target [2,3,5], in which the ToF-forward wall was used, in this analysis, proton and kaon spectra were not able to be distinguished for certain momentum ranges on the basis of dE/dx information alone.

A. Selection of tracks
The selection criteria are devised to remove off-time tracks and tracks coming from secondary interactions mistakenly reconstructed to the main interaction vertex. The selection cuts are also devised to filter out tracks with poorly determined track parameters, mainly p, θ and dE/dx. To start with, all tracks emanating from the main interaction vertex are considered for the dE/dx analysis.

Track topologies
There are a few ways tracks can be classified into different track topologies, including the initial direction of the tracks and which TPC chambers the tracks pass through. The most basic track topology classification used in NA61/SHINE analyses is the distinction between so-called right-side tracks (RSTs) and wrong-side tracks (WSTs) determined by the charge and direction emitted from the target. RSTs have a reconstructed p x that is in the same direction as the deflection by the vertex magnets. WSTs have a reconstructed p x opposite to the bending direction of the magnetic fields. This can be written more succinctly, For the same reconstructed momenta, RSTs and WSTs have very different detector acceptances, numbers of clusters and trajectories through different TPC sectors. Therefore, in this analysis, RSTs and WSTs undergo different selection criteria, are fit separately and have different corrections applied to them. This classification allows for a basic cross-check, since these two samples lead to two somewhat independent measurements. For the purposes of this analysis, the distinction between RSTs and WSTs is not made for the first angular bin ([0,10] mrad for pions and [0,20] mrad for kaons and protons), because it is difficult to accurately distinguish between RSTs and WSTs near θ ¼ 0 mrad.

Phi cuts
The azimuthal acceptance of the NA61/SHINE detector is highly dependent on the track topology and θ. In order to obtain samples of tracks with similar numbers of clusters, ϕ cuts were devised as a function of θ bin and track topology and applied to the selection.

Track quality cuts
The impact parameter of tracks (distance from the main interaction vertex and the extrapolation of the track to the plane of the target) is required to be less than 2 cm in order to remove off-time tracks and tracks produced in secondary interactions.
To ensure that the selected tracks have narrow enough dE/dx distributions to distinguish between particle species, at least 30 clusters are required in the VTPCs and MTPCs. In order to ensure tracks have good momentum estimations, there must be at least four clusters in the GTPC or ten clusters in the VTPCs. Additionally, to remove tracks resulting from secondary interactions that were falsely reconstructed to the main interaction vertex, a cut is applied to tracks with no reconstructed GTPC and VTPC-1 clusters. This cut requires there to be fewer than ten potential clusters in the VTPC-1 and fewer than seven potential clusters in the GTPC, where the potential clusters are calculated by extrapolating tracks through the tracking system.
Several dE/dx cuts were applied to remove tracks with nonsensical dE/dx values (MIP) and rare heavier mass or doubly charged particles, These cuts remove much less than 1% of tracks, so no correction is made to account for the dE/dx cuts. Figure 10 shows the dE/dx-momentum distribution of the selected positively charged and negatively charged tracks.

B. Fitting to dE/dx distributions
For each analysis bin, a fit is used to determine the yields of each particle species. Five particle species and their antiparticles are considered: e þ , π þ , K þ , protons and deuterons. Positively charged and negatively charged tracks are simultaneously fit to better constrain the parameters.

dE/dx model
The mean dE/dx, hϵi, of charged particles passing through NA61/SHINE's TPCs depends on the particles' values of β, which, for particles of the same momentum, depend on their masses. A Bethe-Bloch table provides initial guesses of hϵi for particle species within each bin.
The dE/dx distribution function describing the observed dE/dx of a charged particle passing through the TPCs depends on hϵi and the distance traveled through the TPCs. The distribution closely resembles an asymmetric Gaussian, where ϵ is the measured dE/dx of a track. The peak dE/dx of the distribution, μ, is related to hϵi through the relation where d is the asymmetry parameter, which controls the asymmetry of the distribution through the relation For a detector with uniform readout electronics, the width of the distribution for a single particle depends on the number of dE/dx clusters, N Cl , and on hϵi, where the parameter, α, controls how the width scales with hϵi and σ 0 is the base dE/dx width of a single cluster. However, in NA61/SHINE, nonuniform readout electronics leads to different base widths for clusters reconstructed in different areas of the detector. This effect is most apparent in three main areas of the NA61/SHINE TPC system: the MTPCs, the two most upstream sectors of the VTPCs, and the rest of the VTPCs. Different base widths characterize each of these regions: σ 0;M , σ 0;Up and σ 0;V . The dE/dx width of a single track can be parametrized more precisely by accounting for the numbers of clusters in each TPC region, N Cl;Up , N Cl;V and N Cl;M , At this point, some calibration and shape parameters need to be added in to account for imperfect dE/dx calibration, variation in pad response, variation in track angle and other effects that can cause hϵi and σ to deviate from the ideal model. Therefore, additional calibration parameters are added to allow the peaks and widths of the species distribution functions to vary slightly from the ideal model for each analysis bin.
The full form of the single species distribution function is then where σ i;j cal and μ i;j cal implicitly depend on the momentum p, the number of cluster variables and the calibration parameters.
With these single-species distribution functions the single-track distribution functions can be built for both charges, F þ and F − , where y i;j is the fractional contribution of species i to the sample of tracks with charge j. The yields for each charge are constrained such that they sum to 1.

Fitting strategy
To perform the minimization, a continuous log-likelihood function is constructed, logL ¼ analysis bin. In addition to the constraint that the yield fractions add up to 1 for each charge, soft constraints are applied to avoid the parameters converging to unreasonable values. For example, without constraints, it is easy for two species to swap the location of their dE/dx means. For fits to the target-removed data, all of the parameters are fixed to the fitted values from the target-inserted fits, except for the particle yields. Figure 11 shows a fit to the dE/dx distribution of an example bin. The estimated raw yield of a particle species in analysis bin k is obtained by multiplying the fractional yield obtained from the fit, y i;j k , by the number of positively or negatively charged tracks in that bin, N i k , For each of the π þ , π − , K þ , K − and proton analyses, a raw yield is obtained for each bin and for both the targetinserted and target-removed samples.
C. Corrections

Fit bias corrections
Simulated dE/dx distributions were generated in order to estimate the bias and the standard deviation of the particle yields obtained from the fitting procedure. Fifty simulated dE/dx distributions for each analysis bin were built from the dE/dx model discussed in the previous section. The kinematic variables of tracks from data and the resulting hadron yields were taken as inputs for the dE/dx simulation. The fit parameters are varied according to the spread of fit results observed in data.
The biases and standard deviations in the fitted yields are determined from the results of fits to these simulated dE/dx distributions. In general, the biases in the pion yields are small. The biases of the proton and kaon yields are larger in the high momentum regions and near the Bethe-Bloch crossing regions, where the particle distributions overlap significantly. The biases are used to correct the fit results with correction factors, c fit k , and the standard deviations are used to estimate the uncertainties related to the fitting procedure.

Monte Carlo corrections
The raw yields of particles obtained from the dE/dx fits must be corrected for a number of systematic effects. These can roughly be organized into detector acceptance, feeddown corrections, reconstruction efficiency, selection efficiency and in the case of pions, muon contamination. The combined effect of these individual effects can be estimated as an overall correction factor from Monte Carlo simulations, as was done in the V 0 analysis. A few of the forward kinematic bins contain particle trajectories that strike the S4. A further correction was applied to account for this effect, which reached about 7% for a few of the π þ bins, but did not exceed 2% for the other charged hadron species.
In the case of corrections for π þ and π − , because the dE/dx signal from muons is indistinguishable from pions, muon tracks that pass the selection criteria and are fitted to the main interaction vertex must also be accounted for, FIG. 11. An example fit to a dE/dx distribution is shown for the analysis of pions. On the top, the dE/dx distributions are shown for positively charged tracks (left) and negatively charged tracks (right) along with the fitted contributions due to the five particle species considered. On the bottom, the residuals of the fit with respect to the dE/dx distribution are shown.

Feed-down Reweighting
The feed-down correction, which can be as large as 20% for protons, is the main component of the MC correction factor that depends on the physics model. We cannot assume that the production of Λ,Λ and K 0 S is accurately predicted by the physics generators. This incurs an uncertainty on the MC corrections and subsequently, on the resulting multiplicity measurements.
We can constrain this uncertainty by reweighting our MC productions with the results of the V 0 analyses. When counting the number of reconstructed pions and protons passing the selection criteria, a weight is applied whenever that reconstructed track comes from a K 0 S , Λ orΛ, where m data β is the multiplicity measured in bin β of the V 0 analysis and m MC β is the multiplicity observed in the simulation in that bin.

VII. SYSTEMATIC UNCERTAINTIES ON SPECTRA MEASUREMENTS
A number of possible systematic effects on the multiplicity measurements have also been evaluated. These include biases and uncertainties incurred by the fitting procedures, uncertainties associated with the MC corrections, uncertainties incurred in the selection procedures and uncertainties associated with the reconstruction. On top of the uncertainties described in the following sections, an overall normalization uncertainty is attributed to all of the multiplicity measurements. It has been estimated to be AE 2 1 % by propagating the uncertainties on the normalization constants derived from the integrated cross section analysis through the multiplicity calculation, which is discussed in Sec. VIII.

A. Fit model uncertainty
In the V 0 analysis, it cannot be assumed that the fits to the invariant mass distributions perfectly separate the signal from the background. To check for biases in the fit results, the fitting procedure is performed on additional MC productions using GEANT4 physics lists QGSP_BERT, QBBC and FTF_BIC. With these samples, the numbers of true K 0 S , Λ andΛ are known, so the bias and the standard deviation of the fit result can be calculated. For K 0 S , Λ and Λ, the fitting bias, μ, on the signal fraction, c s , was found to be 3.3% AE 2.7%, 4.8% AE 4.2% and 11% AE 10%, respectively. The bias is not used as a correction for the fit results, but the values of μ AE σ are taken as upper and lower uncertainties on the signal fraction, which are propagated through the multiplicity calculation.
The fit model uncertainties on the charged spectra are obtained from the fits to simulated dE/dx distributions discussed in Sec. VI C 1. The standard deviations in the particle yields are propagated to the multiplicities and taken as the uncertainties associated with the fitting routine.

B. Physics uncertainties
Assuming different underlying physics can lead to different MC correction factors. For example, if the acceptance changes as a function of p and θ, different MC-predicted p and θ distributions can lead to different MC correction factors. This uncertainty is evaluated by applying correction factors obtained with additional MC productions using the physics lists: QGSP_BERT, QBBC and FTF_BIC. The upper and lower bounds on the uncertainties are taken as the maximum and minimum values of the multiplicity obtained using these additional MC correction factors for each analysis bin.

C. Feed-down uncertainties
The MC corrections account for a background of produced hadrons coming from heavier weakly decaying particles. However, it cannot be assumed that the physics generators correctly predict the production rates of these heavier weakly decaying hadrons. This uncertainty is evaluated by assuming a 50% uncertainty on the number of reconstructed feed-down particles when calculating the MC correction factors, unless the feed-down particle was a reweighted K 0 S , Λ orΛ. In this case, the upper and lower uncertainties on the associated neutral hadron spectra are assigned to the weight assigned to the feed-down particles. These uncertainties are then propagated to the multiplicities. This reweighting treatment results in a significant reduction of the uncertainties on the π þ , π − and proton spectra.

D. Selection uncertainties
Although the MC corrections account for the efficiency of the selection cuts, differences in data and MC could incur systematic biases in the result. It was found that tracks in data are typically composed of around 5% fewer clusters than tracks in MC for the same kinematics. To estimate the selection uncertainty, alternative sets of MC corrections were obtained by artificially decreasing the numbers of clusters in MC tracks by 5%. Higher multiplicities are obtained when applying these alternative correction factors, which are taken as the upper bounds of the selection uncertainty.

E. Reconstruction uncertainties
The MC corrections should account for inefficiencies in the reconstruction of tracks and V 0 s if the geometry and detector response are perfectly modeled by the simulation. Differences between the real detector and the simulated detector could lead to systematic effects on reconstruction efficiency component of the MC corrections. To estimate this uncertainty, the detectors were purposefully moved in the detector description model used by the reconstruction. Specifically, eight alternative productions were made after shifting the VTPC-1 and VTPC-2 by þ.2 mm and −2 mm in the x direction and þ5 mm and −.5 mm in the y direction. These shifts are considered to be rather large when compared to the alignment effects seen in the calibration of the data.
The numbers of selected charged tracks and V 0 candidates were calculated from these alternative productions. The maximum difference in the number of candidate tracks=V 0 s among the productions is calculated for the x shifts and the y shifts in each analysis bin. The effects of the x and y shifts are then added in quadrature to estimate the uncertainty for each bin. The resulting uncertainties are generally less than 1% and do not exceed 4%.

F. Momentum uncertainties
There is an uncertainty on the reconstruction of momentum due to uncertainties in converting the magnet currents to magnetic field strength. This uncertainty can be investigated by checking the invariant mass distributions fitted in the V0 analysis. The variation in the fitted means of the invariant mass distributions of K 0 S and Λ indicates an uncertainty in the reconstruction of momentum of up to 0.3%. Uncertainties on the measured multiplicities due to misreconstructed momenta was determined by varying the momenta of tracks by 0.3% and recalculating the numbers of selected tracks and V 0 candidates. This uncertainty was determined to be less than 1% for the majority of the analysis bins, but is on the level of the statistical uncertainty for some of the analysis bins at the edges of the phase space measured.

G. Breakdowns in uncertainties
The breakdowns in the uncertainties for π þ , K þ , proton, K 0 S and Λ spectra from π þ þ C at 60 GeV=c interactions are shown for representative angular bins in Fig. 12. These breakdowns include statistical uncertainties, fit uncertainties, physics uncertainties, feed-down uncertainties, selection uncertainties, momentum uncertainties and reconstruction uncertainties. The breakdowns of the uncertainties are largely similar for the measured hadron spectra from interactions of π þ þ Be at 60 GeV=c. Figures in Ref. [25] present breakdowns of the uncertainties for the complete set of spectra measurements for interactions of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c.
For the neutral spectra, the uncertainties are within 10% in the kinematic regions with good detector acceptance and high statistical power. In the low-momentum regions, uncertainties associated with the fitting routine tend to dominate the lower uncertainties and selection uncertainties tend to dominate the upper uncertainties. The physics model uncertainty is typically the largest component of the uncertainty in the high momenta regions.
For the charged spectra, the total uncertainties are generally around 5% or less except in the kinematic regions with poor acceptance or poor dE/dx separation. In spectra of π þ , the largest uncertainties tend to be reconstruction uncertainties at high momenta and dE/dx fit uncertainties at low momenta. In the case of π − , dE/dx fit uncertainties, physics model uncertainties and statistical uncertainties contribute the most to the total uncertainty. For kaons, dE/ dx fit uncertainties are dominant in the majority of the phase space measured. For protons, uncertainties related to the physics model and dE/dx fit uncertainties are dominant for the majority of the phase space measured.

VIII. DIFFERENTIAL PRODUCTION MULTIPLICITY MEASUREMENTS
The differential production multiplicity is the yield of particles produced per production interaction per unit momentum per radian in each kinematic bin k. The production multiplicity for neutral hadrons can be written as where ΔpΔθ is the size of bin k, and the yields, Y I;R k , are the total numbers of particles observed in bin k determined by the invariant mass fits for target-inserted and target-removed data. The constants σ trig , σ prod , f prod and ϵ are determined from the integrated cross section analysis and N I and N R are the numbers of selected events with the target inserted and target removed. The differential cross section is related to the multiplicity by a factor of σ prod , In order to calculate the multiplicity for produced charged hadrons (for each track topology, RST and WST), an additional correction factor is required for the fit bias corrections, c fit , For kinematic bins for which the detector acceptance and fit reliability is sufficient enough for multiplicity measurements in both RST and WST bins, the single-side FIG. 12. The breakdown of the fractional uncertainties on π þ , K þ , proton, K 0 S and Λ spectra from π þ þ C at 60 GeV=c interactions for select representative angular bins. The upper and lower uncertainties are shown on the positive and negative sides of the y axes.
FIG. 13. K 0 S multiplicity spectra from π þ þ C at 60 GeV=c interactions are shown for different regions of θ. The error bars represent total uncertainties except for the normalization uncertainty. The results are compared to the predictions of the GEANT4 physics lists QGSP_BERT and FTF_BIC as well as GiBUU2019 and FLUKA2011.
FIG. 14. Λ multiplicity spectra from π þ þ C at 60 GeV=c interactions are shown for different regions of θ. The error bars represent total uncertainties except for the normalization uncertainty. The results are compared to the predictions of the GEANT4 physics lists QGSP_BERT and FTF_BIC as well as GiBUU2019 and FLUKA2011. multiplicities, m R and m W , are merged by taking the weighted average, where the merged uncertainty, σ merged , is calculated with The uncertainties on the individual RST and WST multiplicities consider both the statistical uncertainties and the fit uncertainties, Measurements of spectra from π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c interactions are shown for produced π þ , K þ , proton, K 0 S and Λ for select representative angular bins. The error bars represent total uncertainties except for the normalization uncertainty.
In analysis bins for which the detector acceptance is only sufficient for either RSTs or WSTs, only the single-side multiplicity and uncertainty is taken as the result. Multiplicity spectra obtained for K 0 S , Λ andΛ in π þ þ C at 60 GeV=c interactions are presented in Figs. 13-15. The spectra are shown as one-dimensional momentum spectra for individual bins of θ. The error bars represent the total uncertainty except for the normalization uncertainty. The results are compared to the predictions of the GEANT4 physics lists: QGSP_BERT and FTF_BIC as well as GiBUU2019 [26] and FLUKA2011.2x.7 [27][28][29]. In general, the K 0 S spectra fall within the range of predictions of the models used. No single model describes the K 0 S spectra aptly for the full phase space. The models exhibit a large variability in their predictions of Λ and especiallyΛ spectra. QGSP_BERT seems to provide the best prediction of Λ spectra, while no single model seems to provide a satisfactory description ofΛ spectra. Tables in Ref. [25] present the numerical values of the multiplicity measurements of K 0 S , Λ andΛ along with statistical, systematic and total uncertainties for each kinematic bin analyzed. The normalization uncertainty of AE 2 1 % is not included in the values of the uncertainties shown in these tables but should be attributed to the multiplicity spectra of all hadron species analyzed.
Multiplicity spectra obtained for charged pions, charged kaons and protons in π þ þ C at 60 GeV=c interactions are shown in Figs. 16-20. The results are compared to the predictions of the GEANT4 physics lists: QGSP_BERT and FTF_BIC as well as GiBUU2019 and FLUKA2011. In general, no single model provides a good description of the charged hadron spectra for all particle species and for the full phase space, but FLUKA2011 seems to provide the best overall description. However, it should be noted that in the first few angular bins, an important region for neutrino beams, FLUKA2011 slightly overpredicts the production of π þ . Tables in Ref. [25] present the numerical values of the multiplicity measurements of charged pions, charged kaons and protons along with statistical, systematic and total uncertainties for each kinematic bin analyzed. The normalization uncertainty of AE 2 1 % is not included in the values of the uncertainties shown in these tables but should be attributed to the multiplicity spectra of all hadron species analyzed.
Measurements of spectra of produced π þ , K þ , proton, K 0 S and Λ from interactions of π þ þ Be at 60 GeV=c are shown in comparison to the results for interactions of π þ þ C at 60 GeV=c for representative angular bins in Fig. 21. The spectra are largely similar. The most notable difference in the spectra is that the multiplicities tend to be lower in the regions of low momentum and high production angle in interactions of π þ þ Be at 60 GeV=c. The full set of comparisons between the spectra results of π þ þ Be at 60 GeV=c and π þ þ C at 60 GeV=c is presented in Ref. [25].

IX. SUMMARY AND CONCLUSIONS
In summary, hadron production was studied in interactions of π þ þ C at 60 GeV=c and π þ þ Be at 60 GeV=c. For both of these reactions, the integrated production and inelastic cross sections were measured. Furthermore, differential cross sections were measured for produced π þ , π − , K þ , K − , protons, K 0 S , Λ andΛ. The inelastic cross section measurements are the first to be made at a beam momentum of 60 GeV=c. The production cross section of interactions of π þ þ Be at 60 GeV=c was measured for the first time as well. The differential cross sections were measured for the first time at this beam momentum scale, and compared to previous measurements at lower beam momenta, a larger kinematic phase space and more particle species were studied. These results will enable neutrino flux predictions to be constrained in neutrino experiments using the NuMI beam and future neutrino beam at LBNF. Specifically, these results can be used to reduce the uncertainties associated with secondary interactions of pions in the carbon targets and the beryllium elements in these beam lines.