Large-angle production of charged pions with incident pion beams on nuclear targets

Measurements of the double-differential pi+/- production cross-section in the range of momentum 100 MeV/c<= p<= 800 MeV/c and angle 0.35 rad<= theta<= 2.15 rad using pi+/- beams incident on beryllium, aluminium, carbon, copper, tin, tantalum and lead targets are presented. The data were taken with the large acceptance HARP detector in the T9 beam line of the CERN Proton Synchrotron. The secondary pions were produced by beams in a momentum range from 3 GeV/c to 12.9 GeV/c hitting a solid target with a thickness of 5% of a nuclear interaction length. The tracking and identification of the produced particles was performed using a small-radius cylindrical time projection chamber (TPC) placed inside a solenoidal magnet. Incident particles were identified by an elaborate system of beam detectors. Results are obtained for the double-differential cross-sections d2sigma/dpdtheta at six incident beam momenta. Data at 3 GeV/c, 5 GeV/c, 8 GeV/c, and 12 GeV/c are available for all targets while additional data at 8.9 GeV/c and 12.9 GeV/c were taken in positive particle beams on Be and Al targets, respectively. The measurements are compared with several generators of GEANT4 and the MARS Monte Carlo simulation.


I. INTRODUCTION
In many particle and astroparticle physics experiments external knowledge of hadron production is required to make optimal use of the recorded data and to help design the experimental facilities. The HARP experiment [1] is motivated by this need for precise hadron production measurements. It has taken data with beams of pions and protons with momenta from 1.5 GeV/c to 15 GeV/c hitting targets made of a large range of materials. To provide a large angular and momentum coverage of the produced charged particles the experiment comprises two spectrometers, a forward spectrometer built around a dipole magnet covering the angular range up to 250 mrad and a large-angle spectrometer constructed in a solenoidal magnet with an angular acceptance of 0.35 rad ≤ θ ≤ 2.15 rad.
The main objectives are to measure pion yields for a quantitative design of the proton driver of future superbeams (high-intensity conventional beams) and a neutrino factory [2], to provide measurements to improve calculations of the atmospheric neutrino flux [3][4][5][6] and to measure particle yields as input for the flux calculation of accelerator neutrino experiments [7], such as K2K [8,9], MiniBooNE [10] and SciBooNE [11]. In addition to these specific aims, the data provided by HARP are valuable for validating hadron production models used in simulation programs. These simulations are playing an important role in the interpretation and design of modern large particle-physics experiments. In particular, the simulation of calorimeter response and secondary interactions in tracking systems needs to be supported by experimental hadron production data.
This paper presents the final measurements of the double-differential production cross-section, d 2 σ π /dpdθ for π ± production at large angles by positively and negatively charged pions of 3 GeV/c, 5 GeV/c, 8 GeV/c, 8.9 GeV/c (π + -Be), 12 GeV/c and 12.9 GeV/c (π + -Al) momenta impinging on a thin beryllium, carbon, aluminium, copper, tin, tantalum or lead target of 5% nuclear interaction length. The data are taken with the large-angle spectrometer of the HARP detector. The results of a similar analysis of π ± production data taken with proton beams on the same set of nuclear targets can be found in Ref. [12].
The results presented in this paper, covering an extended range of solid targets in the same experiment, make it possible to perform systematic comparisons of hadron production models with measurements at different incoming beam momenta over a large range of target atomic number A. The performance of these models can be further tested by extending the comparisons with our results on pion production in p-A interactions obtained with the forward spectrometer. These results are the subject of other HARP publications [13][14][15][16][17]. The analysis of results taken with charged pion beams on the full range of targets presented in this paper, but measured with the forward spectrometer, can be found in Ref. [18]. Pion production data in this energy region are extremely scarce, especially in pion beams, and HARP is the first experiment to provide a large data set taken with many different targets, full particle identification and a large-acceptance detector [44].
Data were taken in the T9 beam of the CERN PS. The collected statistics, for the different nuclear targets, are reported in Tables I and II. The analysis proceeds by selecting tracks in the time projection chamber (TPC) in events with incident beams of charged pions of both polarities. Momentum and polar angle measurements and particle identification are based on the measurements of track position and energy deposition in the TPC. An unfolding method is used to correct for experimental resolution, efficiency and acceptance and to obtain the double-differential pion production cross-sections. The analysis follows the same methods as used for the determination of π ± production by protons. These analysis methods are documented in Ref. [21] with improvements described in Ref. [12] and will be only briefly outlined here.

II. EXPERIMENTAL APPARATUS AND DATA SELECTION.
The HARP detector is shown in Fig. 1 and is described in detail in Ref. [22]. The forward spectrometer, mainly used in the analysis for the conventional neutrino beams and atmospheric neutrino flux, consists of a dipole magnet, large planar drift chambers (NDCs) [23] , a time-of-flight wall (TOFW) [24], a threshold Cherenkov counter (CHE), and an electromagnetic calorimeter (ECAL).
In the large-angle region a cylindrical TPC with a radius of 408 mm is positioned inside a solenoidal magnet with a field of 0.7 T. The TPC detector was designed to measure and identify tracks in the angular region from 0.25 rad to 2.5 rad with respect to the beam axis. The target is placed inside the inner field cage (IFC) of the TPC such that, in addition to particles produced in the forward direction, backward-going tracks can be measured. The TPC is used for tracking, momentum determination and measurement of the energy deposition dE/dx for particle identification [25]. A set of resistive plate chambers (RPCs) form a barrel inside the solenoid around the TPC to measure the arrival time of the secondary particles [26]. Charged particle identification (PID) is achieved with the TPC by measuring the ionization per unit length in the gas (dE/dx) as a function of the total momentum of the particle. Additional PID can be performed through a time-of-flight measurement with the RPCs; this method is used to calculate the small components of kaons and deuterons and heavier ions. Its composition depends on the selected beam momentum. The proton fraction in the incoming positive-particle beam varies from 35% at 3 GeV/c to 92% at 12 GeV/c. This explains the relatively low statistics for some of the 12 GeV/c data sets taken with the π + beam. The negatively charged particle beam is mainly composed of pions with small background components of muons and electrons.
At the first stage of the analysis a favoured beam particle type is selected using the beam time-of-flight system and the two Cherenkov counters. A value of the pulse height consistent with the absence of a signal in both beam Cherenkov detectors distinguishes protons (and kaons) from electrons and pions. We also ask for time measurements to be present which are needed for calculating the arrival time of the beam proton at the target. The beam TOF system is used to reject ions, such as deuterons, and to separate protons from pions at low momenta. At 3 GeV/c, the TOF measurement allows the selection of pions from protons to be made at more than 5σ. In most beam settings the nitrogen pressure in the beam Cherenkov counters was too low for kaons to be above the threshold. Kaons are thus counted in the proton sample. However, the fraction of kaons has been measured in the 12.9 GeV/c beam configuration and are found to contribute less than 0.5%, and hence are negligible in the pion beam sample. Electrons radiate in the Cherenkov counters and would be counted as pions. In the 3 GeV/c beam, electrons are identified by both Cherenkov counters, since the pressure was such that pions remained below threshold. In the 5 GeV/c beam electrons could be tagged by one Cherenkov counter only, while the other Cherenkov counter was used to tag pions. The e/π fraction was measured to be 1% in the 3 GeV/c beam and < 10 −3 in the 5 GeV/c beam. By extrapolation from the lower-energy beam settings this electron contamination can be estimated to be negligible (< 10 −3 ) for the beams where it cannot be measured directly. More details on the beam particle selection can be found in Refs. [22] and [13].
In addition to the momentum-selected beam of protons and pions originating from the T9 production target one expects also the presence of muons from pion decay both downstream and upstream of the beam momentum selection. Therefore, precise absolute knowledge of the muon rate incident on the HARP targets is required when measurements of particle production with incident pions are performed. The particle identification detectors in the beam do not distinguish muons from pions. A separate measurement of the muon component has been performed using data-sets without target ("empty-target data-sets"). Since the empty-target data were taken with the same beam parameter settings as the data taken with targets, the beam composition can be measured in the empty-target runs using the forward spectrometer and then used as an overall correction for the counting of pions in the runs with targets. Muons are recognized by their longer range in the beam muon identifier (BMI). The BMI is a small instrumented stack of iron absorbers at the downstream end of the spectrometer. The punch-through background in the BMI is measured counting the protons (identified with the beam detectors) thus mis-identified as muons by the BMI. A comparison of the punch-through rate between simulated incoming pions and protons was used to determine a correction for the difference between pions and protons and to determine the systematic error. This difference is the dominant systematic error in the beam composition measurement. The aim was to determine the composition of the beam as it strikes the target, thus muons produced in pion decays after the HARP target should be considered as a background to the measurement of muons in the beam. The rate of these latter background muons, which depends mainly on the total inelastic cross-section and pion decay, was calculated by a Monte Carlo simulation using GEANT4 [31]. The muon fraction in the beam (at the target) is obtained taking into account the efficiency of the BMI selection criteria as well as the punch-through and decay backgrounds. The analyses for the various beam settings give results for R = µ/(µ + π) of (4.2±1)% and (5.2±1)% for the low-momentum beams (3 GeV/c and 5 GeV/c) and between (4.1±1)% and (2.8±1)% for the highest momenta (from 8 GeV/c to 12.9 GeV/c). The uncertainty in these fractions is dominated by the systematic uncertainty in the punch-through background. The fact that the background does not scale with the decay probability for pions is due to the limited acceptance of the beam-line to transport the decay muons. The muon contamination is taken into account in the normalization of the pion beam and adds a systematic error of 1% to the overall normalization.
A set of MWPCs is used to retain events with only one beam particle for which the trajectory extrapolates to the target. An identical beam particle selection was performed for events triggered with the incident-beam trigger in order to provide an absolute normalization of the incoming pions. This trigger selected every 64 th beam particle coincidence outside the dead-time of the data acquisition system.
The length of the accelerator spill is 400 ms with a typical intensity of 15 000 beam particles per spill in the positive beam and lower for the negative beam. The average number of events recorded by the data acquisition ranges from 300 to 350 per spill for the different beam momenta and beam polarities. The analysis proceeds by first selecting a beam pion hitting the target, not accompanied by other beam tracks. Then an event is required to give a large-angle interaction (LAI) trigger to be retained. After the event selection the sample of tracks to be used for analysis is defined. Tables I and II show the number of events and the number of π ± selected in the analysis for the π + and π − data, respectively. The large difference between the first and second set of rows ("Total DAQ events" and "Accepted pions with LAI") is due to the relatively large fraction of protons in the beam and to the larger number of triggers taken for the measurements with the forward dipole spectrometer ("forward triggers"). The entry "fraction of triggers I: Total number of events and tracks used in the various nuclear 5% λI target data sets taken with the π + beam and the number of pions on target as calculated from the pre-scaled incident pion triggers.
Data set (GeV/c) 3   used" shows the part of the data for which the dynamic distortions in the TPC could be calibrated reliably.

III. PARTICLE IDENTIFICATION
The particle identification in the large-angle region uses the dE/dx information provided by the TPC. The electron, pion and proton populations are well separated at most momentum values. Simple momentum-dependent cuts are used to separate the different populations. The pions are identified by removing electrons and protons. The kaon population is negligible. The cuts were optimized to maximize the purity of the pion sample, accepting a lower efficiency in the selection. More details are given in Ref. [21].
The measurement of the velocity β of secondary particles by the time-of-flight determination with the RPC detectors drop in efficiency toward higher momenta is caused by the need to make a hard cut to remove protons. The migration of pions and protons into the wrong sample is kept below the percent level in the momentum range of this analysis (p < 800 MeV/c). This is important for the measurement of the π + production rate since the proton production rate is significantly larger in some of the bins. The small differences in efficiency (up to ≈ 5%) which are visible between the data and the simulation are dealt with in the analysis by an ad hoc correction to the cross-sections. It has been checked that the angular dependences of the PID efficiency and purity are negligible. As already stated, electrons and positrons with a momentum above 125 MeV/c cannot be separated cleanly from pions. A simulation study shows that the dominant source of electron background is due to π 0 production. The approach chosen is not to attempt to identify these electrons but to consider these as background to be subtracted.
The assumption is made that the π 0 spectrum is similar to the spectrum of charged pions. With this assumption, the electron background is negligible above 250 MeV/c-300 MeV/c. The subtraction uses an iterative approach. Initial π − and π + spectra are obtained in an analysis without π 0 subtraction. The spectra of pions with opposite charge compared to the beam are then used in the Monte-Carlo (MC) simulation for the π 0 distributions. A full simulation of the production and decay into γ's with subsequent conversion in the detector materials is used to predict the background electron and positron tracks. Most of these tracks have a momentum below the threshold for this analysis or low enough to be recognized by dE/dx. The tracks with a PID below the expected value for pions can be rejected as background. In the region below 125 MeV/c, a large fraction of the electrons can be unambiguously identified. These tracks are used as relative normalization between data and MC. The remaining background is then estimated from the distributions of the simulated electron and positron tracks which are accepted as pion tracks with the same criteria as used to select the data. These normalized distributions are subtracted from the data before the unfolding procedure is applied. The initial pion spectra obtained in the first pass are not subtracted for this background. This is not a large problem, since this overestimation occurs in a momentum region where the majority of electrons from these π 0 's are too soft to disturb the analysis. Uncertainties in the assumption of the π 0 spectrum are taken into account by an alternative assumption that their spectrum follows the distribution of pions with the same charge as the beam. An additional systematic error of 10% is assigned to the normalization of the π 0 subtraction using the identified electrons and positrons. At low momenta and small angles the π 0 subtraction introduces the largest systematic uncertainty. It is in principle possible to reject more electrons and positrons by constructing a combined PID estimator based on dE/dx and TOF. Indeed, such an analysis was performed and gave consistent results. However, its systematic errors are more difficult to estimate.

IV. DATA ANALYSIS
Only a short outline of the data analysis procedure is presented here: for further details see Refs. [12,21]. The double-differential cross-section for the production of a particle of type α can be expressed in the laboratory system as: where d 2 σα dpidθj is expressed in bins of true momentum (p i ), angle (θ j ) and particle type (α). The "raw yield" N α ′ i ′ j ′ is the number of particles of observed type α ′ in bins of reconstructed momentum (p i ′ ) and angle (θ j ′ ). These particles must satisfy the event, track and PID selection criteria. Although, owing to the stringent PID selection, the background from misidentified protons in the pion sample is small, the pion and proton raw yields (N α ′ i ′ j ′ , for α ′ = π − , π + , p) have been measured simultaneously. It is thus possible to correct for the small remaining proton background in the pion data without prior assumptions concerning the proton production cross-section.
Various techniques are described in the literature to obtain the matrix M −1 ijαi ′ j ′ α ′ . In this analysis an unfolding technique is used. It performs a simultaneous unfolding of p, θ and PID, with a correction matrix M −1 computed mainly using the Monte Carlo simulation.
The matrix M −1 ijαi ′ j ′ α ′ corrects for the efficiency and the resolution of the detector. It unfolds the true variables i, j, α from the reconstructed variables i ′ , j ′ , α ′ with a Bayesian technique [32] and corrects the observed number of particles to take into account effects such as trigger efficiency, reconstruction efficiency, acceptance, absorption, pion decay, tertiary production, PID efficiency, PID misidentification and electron background. The central assumption of the method is that the probability density function in the ("true") physical parameters ("physical distribution") can be approximated by a histogram with bins of sufficiently small width. A population in the physical distribution of events in a given cell ijα generates a distribution in the measured variables, M ijαi ′ j ′ α ′ . Thus the observed distribution in the measurements can be represented by a linear superposition of such populations. The task of the unfolding procedure consists then of finding the number of events in the physical bins for which the predicted superposition in the measurement space gives the best description of the data. The method used to correct for the various effects is described in more detail in Ref. [21].
To predict the population of the migration matrix element M ijαi ′ j ′ α ′ , the resolution, efficiency and acceptance of the detector are obtained from the Monte Carlo. This is accurate provided the Monte Carlo simulation describes these quantities correctly. Where some deviations from the control samples measured from the data are found, the data are used to introduce (small) ad hoc corrections to the Monte Carlo. Using the unfolding approach, possible known biases in the measurements are taken into account automatically as long as they are described by the Monte Carlo. In the experiment simulation, which is based on the GEANT4 toolkit [31], the materials in the beam-line and the detector are accurately described as well as the relevant features of the detector response and the digitization process. The time-dependent properties of the TPC, such as pulse-height calibration per channel and the presence of dead channels, were reproduced for each individual data set by running a dedicated set of high-statistics simulations corresponding to each data set. In general, the Monte Carlo simulation compares well with the data, as shown in Ref. [21]. For all important issues, physical benchmarks have been used to validate the analysis. The absolute efficiency and the measurement of the angle and momentum was determined with elastic scattering. The momentum and angular resolution was determined exploiting the two halves of cosmic-ray tracks crossing the TPC volume. The efficiency of the particle identification was checked using two independent detector systems. Only the latter needs a small ad hoc correction compared to the simulation.
The factor A NAρt in Eq. 1 is the inverse of the number of target nuclei per unit area (A is the atomic mass, N A is the Avogadro number, ρ and t are the target density and thickness). We do not make a correction for the attenuation of the beam in the target, so that strictly speaking the cross-sections are valid for a λ I = 5% target. The result is normalized to the number of incident pions on the target N pot . The absolute normalization of the result is calculated in the first instance relative to the number of incident-beam particles accepted by the selection. After unfolding, the factor A NAρt is applied. The beam normalization using down-scaled incident pion triggers has uncertainties smaller than 2% for all beam momentum settings.
The background due to interactions of the primary pions outside the target (called the empty target background) is measured using data taken without the target mounted in the target holder. Owing to the selection criteria which only accept events from the target region and the good definition of the interaction point this background is negligible (< 10 −5 ). To subtract backgrounds generated by π 0 's produced in hadronic interactions of the incident-beam particle, the assumption is made that the π 0 spectrum is similar to the spectrum of charged pions. In an iterative procedure the production spectra of pions with opposite charge compared to the beam particles are used for the subtraction, while the difference between π + and π − production is used to estimate the systematic error as explained in Section III. The absorption and decay of particles is simulated by the Monte Carlo. The generated single particle can re-interact and produce background particles by hadronic or electromagnetic processes. These processes are simulated and additional particles reconstructed in the TPC in the same event are taken into account in the unfolding procedure as background.
The effects of the systematic uncertainties on the final results are estimated by repeating the analysis with the relevant input modified within the estimated uncertainty intervals. In many cases this procedure requires the construction of a set of different migration matrices. The correlations of the variations between the cross-section bins are evaluated and expressed in the covariance matrix. Each systematic error source is represented by its own covariance matrix. The sum of these matrices describes the total systematic error. The magnitude of the systematic errors will be shown in Section V.

V. EXPERIMENTAL RESULTS
The measured double-differential cross-sections for the production of π + and π − in the laboratory system as a function of the momentum and the polar angle for each incident beam momentum are shown in Figures 3 to 16 for targets from Be to Pb. The error bars shown are the square-roots of the diagonal elements in the covariance matrix, where statistical and systematic uncertainties are combined in quadrature. Correlations cannot be shown in the figures. The overall scale error (< 2%) is not shown. The correlation of the statistical errors (introduced by the unfolding procedure) are typically smaller than 20% for adjacent momentum bins and even smaller for adjacent angular bins. The correlations of the systematic errors are larger, typically 80% for adjacent bins. The results of this analysis are also tabulated in Appendix A.
One observes that only for the 3 GeV/c beam is the statistical error similar in magnitude to the systematic error, while the statistical error is negligible for the 8 GeV/c and 12 GeV/c beam settings. The statistical error is calculated by error propagation as part of the unfolding procedure. It takes into account that the unfolding matrix is obtained from the data themselves [45] and hence contributes also to the statistical error. This procedure almost doubles the    statistical error, but it avoids an important systematic error which would otherwise be introduced by assuming a cross-section model a priori to calculate the corrections. The largest systematic error corresponds to the uncertainty in the absolute momentum scale, which was estimated to be around 3% using elastic scattering [21]. At low momentum in the relatively small angle forward direction, the uncertainty in the subtraction of the electron and positron background due to π 0 production is dominant (∼ 6 − 10%). This uncertainty is split between the variation in the shape of the π 0 spectrum and the normalization using the recognized electrons. The target region definition and the uncertainty in the PID efficiency and background from tertiaries (particles produced in secondary interactions) are of similar size and are not negligible (∼ 2−3%). Relatively small errors are introduced by the uncertainties in the absorption correction, absolute knowledge of the angular and the momentum resolution. The correction for tertiaries is relatively large at low momenta and large angles (∼ 3−5%). As expected, this region is most affected by this component.
As already mentioned above, the overall normalization has an uncertainty of 2% and is not reported in the table. It is mainly due to the uncertainty in the efficiency that beam pions counted in the normalization actually hit the target and the muon contamination in the beams, with smaller components from the target density and beam particle counting procedure.
One observes the weak beam energy dependence of pion production by incoming pions. This is particularly striking for low A targets and enhanced for same-charge pion production. The energy dependence is larger for the more forward bins in the range covered in this analysis; backward production also shows little beam energy dependence over the whole range in A. The 3 GeV/c data tend to show a markedly different behaviour from all higher beam energy data already at moderate A (Cu) in the most forward angular bins, a trend which extends to larger angles for larger A targets.
The dependence of the averaged pion yields on the incident-beam momentum is shown in Fig. 17. The π + and π − yields are averaged over the region 0.350 rad ≤ θ < 1.550 rad and 100 MeV/c ≤ p < 700 MeV/c (pions produced in the forward direction only). Whereas the beam energy dependence of the yields in the Be and C data differs clearly from the dependence in the Ta and Pb data one can observe that the Al, Cu and Sn data display a smooth transition between them. The dependence in the Be and C data is much more flat with a saturation of the yield between 8 GeV/c and 12 GeV/c with the Al, Cu and Sn showing an intermediate behaviour. The momentum dependence is larger for opposite-charge pion production than for equal-charge production.
The integrated π − /π + ratio in the forward direction is displayed in Figs. 18 and 19 as a function of the secondary momentum for π + and π − beams, respectively.
For the π + beams, in the momentum range covered in most bins more π + 's are produced than π − 's. In the π + -Ta and π + -Pb data the ratio is closer to unity than for the π + -Be and π + -C data. The π − /π + ratio is larger for higher incoming beam momenta than for lower momenta and drops with increasing secondary momentum. Comparing the ratios in the data taken with the π + beam with the ratios previously published for proton beam data [12] one observes a large similarity. The ratios in the p-A data display a somewhat larger beam momentum dependence than the π + data. Especially the 3 GeV/c data are closer to the higher beam momenta in the π + data.
Turning now our attention to the π − /π + ratios in the π − beam data, one observes a much larger beam momentum dependence with a ratio which is larger than unity everywhere. In this case, the 3 GeV/c displays, not surprisingly, the highest π − /π + ratio. Again, the ratio falls with increasing secondary momentum, except for the lowest beam momentum for Pb, a trend which is already visible for Ta. Turning back to the double-differential spectra, one can conclude that this feature is not due to an increase in the production of π − , but rather to a suppression of π + production.
A striking feature is the large π − /π + ratio in the lowest bin of secondary momentum (100 MeV/c-150 MeV/c) for the heavy nuclear targets (Pb and Ta) in the beams with 8 GeV/c and 12 GeV/c momentum. This effect had already been observed in the p-Ta data [21] and is visible in the p, π + and π − beams. For the positive beams these are the only bins where more π − 's are produced than π + 's. The E910 Collaboration had made a similar observation for their lowest momentum bin (100 MeV/c-140 MeV/c) in p-Au collisions at 12.3 GeV/c and 17.5 GeV/c incoming beam momentum [33]. They offered as explanation Λ 0 production at rest which would enhance the π − yield at low secondary momenta. Perhaps a more plausible explanation has been put forward in Ref. [34], where it was shown that this effect was predicted in the model described in that paper. These authors invoke an asymmetry in the production of ∆ resonances due to the large neutron excess in these heavy nuclei. They tested this hypothesis by simulating data on hypothetical n-p symmetric heavy nuclei and found that in that case the effect was not predicted. Perhaps the same explanation holds for the marked difference between the π + -C data and the other targets, carbon being the only n-p symmetric target reported here. The π + -C data show a significantly smaller increase of the π − /π + ratio towards small momenta, a feature also visible in the p-C data of Ref. [12], while not present in the π − -C data.
The dependence of the averaged pion yields on the atomic number A is shown in Fig. 20. The π + yields averaged over the region 0.350 rad ≤ θ < 1.550 rad and 100 MeV/c ≤ p < 700 MeV/c are shown in the left panel and the π − data averaged over the same region in the right panel for four different beam momenta. One observes a smooth behaviour of the averaged yields. The A dependence is slightly different for π − and π + production. The production of opposite-sign pions displays a steeper A dependence than same-sign pion production.
The experimental uncertainties are summarized for π + in Tables III and IV, and for π − in Tables V and VI for all used targets. The relative sizes of the different systematic error sources are very similar for π − and π + and for the different beam energies. Going from lighter (Be, C) to heavier targets (Ta, Pb) the corrections for π 0 (conversion, concentrated at low secondary momentum) and absorption/tertiaries are bigger. Since the production cross-sections are not very different in the pion beams from the proton beam data, the discussion and figures shown in Ref. [12] give a reliable indication of the momentum and angular dependence of the systematic error components and need not be repeated here.
FIG. 3: Double-differential cross-sections for π + production (left) and π − production (right) in π + -Be interactions as a function of momentum displayed in different angular bins (shown in mrad in the panels). In the figure, the symbol legend 9 refers to 8.9 GeV/c nominal beam momentum. The error bars represent the combination of statistical and systematic uncertainties.
FIG. 5: Double-differential cross-sections for π + production (left) and π − production (right) in π + -Al interactions as a function of momentum displayed in different angular bins (shown in mrad in the panels). In the figure, the symbol legend 13 refers to 12.9 GeV/c nominal beam momentum. The error bars represent the combination of statistical and systematic uncertainties.
FIG. 8: Double-differential cross-sections for π + production (left) and π − production (right) in π + -Ta interactions as a function of momentum displayed in different angular bins (shown in mrad in the panels). The error bars represent the combination of statistical and systematic uncertainties.
FIG. 13: Double-differential cross-sections for π + production (left) and π − production (right) in π − -Cu interactions as a function of momentum displayed in different angular bins (shown in mrad in the panels). The error bars represent the combination of statistical and systematic uncertainties.
FIG. 18: From top-left panel to bottom-right panel, the ratio of the differential cross-sections for π − and π + production in π + -Be (top-left), π + -C (top-right), π + -Ta (bottom-left) and π + -Pb (bottom-right) interactions as a function of the secondary momentum integrated over the forward angular region (shown in mrad). In the figure, the symbol legends 13 and 9 refer to 12.9 and 8.9 GeV/c nominal beam momentum, respectively.
FIG. 19: From top-left panel to bottom-right panel, the ratio of the differential cross-sections for π − and π + production in π − -Be (top-left), π − -C (top-right), π − -Ta (bottom-left) and π − -Pb (bottom-right) interactions as a function of the secondary momentum integrated over the forward angular region (shown in mrad).

A. Comparisons with MC predictions
In Figures 21-36, comparisons with a selected set of Monte Carlo generators of GEANT4 [31] and MARS [35] codes are shown. We stress that no tuning to our data has been done by the GEANT4 or MARS teams. The comparisons are shown for the C and Ta targets as examples of a light and a heavy target.
At intermediate energies (up to 5 GeV/c-10 GeV/c), GEANT4 uses two types of intranuclear cascade models: the Bertini model [36,37] (valid up to ∼ 10 GeV) and the Binary model [38] (valid up to ∼ 3 GeV). Both models treat the target nucleus in detail, taking into account density variations and tracking in the nuclear field. The Binary model is based on hadron collisions with nucleons, giving resonances that decay according to their quantum numbers. The Bertini model is based on the cascade code reported in [39] and hadron collisions are assumed to proceed according to free-space partial cross sections corrected for nuclear field effects and final-state distributions measured for the incident-particle types. For the figures the Bertini model was chosen as representative of the cascade models since it was found to perform better than the Binary model.
At higher energies, instead, two parton string models, the quark-gluon string (QGS) model [36,40] and the Fritiof (FTF) model [41] are used, in addition to a high-energy parametrized model (HEP) derived from the high-energy part of the GHEISHA code used inside GEANT3 [42]. The parametrized models of GEANT4 (HEP and LEP) are intended to be fast, but conserve energy and momentum on average and not event by event. In the figures a low-energy version of the FTF model ("FTFB") is shown for the data up to 8 GeV/c, and the higher energy version ("FTFP") to compare with the 12 GeV/c data. The MARS code system [35] uses as basic model an inclusive approach to multiparticle production. Above 5 GeV phenomenological particle production models are used: while below 5 GeV a cascade-exciton model [43] combined with the Fermi breakup model, the coalescence model, an evaporation model and a multifragmentation extension are used instead.
The comparison between data and models reveals sizable differences. Discrepancies up to a factor of three are seen. One should note that the models have more difficulties with the incoming pion data than the incoming proton data presented previously [12].
Let us first examine the comparison for the carbon target. The FTFB model provides a fair description of the π + -C data. There are some detailed differences, e.g. the absolute level of the first angular bin in the π − production in the π + -C data at 3 GeV/c. This model also provides a good description of the π − production in the π − -C data, although it has difficulties with the backward direction, especially for the low-energy data. It has much more difficulty with the π + production in the π − -C data. The other models (MARS and Bertini) are both based on a cascade description at this energy and predict features in the low-energy data which are not visible in the data. At 5 GeV/c the MARS prediction is very different from its 3 GeV/c prediction. This is not justified by the data, which look very similar for all beam energies. The Bertini model is not very successful for any of the π-C data sets, while MARS gives a reasonable description in the mid-energy range and for some of the highest energy data sets. Especially in the more forward bins the Bertini model underestimates the production of pions. At 12 GeV/c the FTFP and QGSP models give a better description for the opposite-charge pion production, while MARS does better for same-charge pions.
In the tantalum data one observes again the structure of the cascade based models at 3 GeV/c, which is not representative of the data. The FTFB model does not have this feature and gives a better description at this energy. All models have considerable difficulties with the Ta data from 5 GeV/c onward, with large differences among the models. The π − production in the π − -Ta data is better described from 5 GeV/c onward by FTFB and MARS.
In general, the Bertini model produces a too isotropic pion production. This may be due to the lack of an explicit diffractive process. On the contrary, the Fritiof based models do contain these processes and describe the forward production better. At low energies MARS uses a cascade model and has similar problems as the Bertini model. One notes that with the Ta target for 5 GeV/c and 8 GeV/c, the predictions for the π + and π − beams are identical. A similar problem was found with the Pb target for 3 GeV/c and 5 GeV/c. This appears to be a technical problem with the model code and the authors have been informed [46].
The overall impression is that all investigated models need to be considerably improved before they can give a reliable description. Especially, the choice to change description between 3 GeV/c and 5 GeV/c does not find support in the data.

VI. SUMMARY AND CONCLUSIONS
An analysis of the production of pions at large angles with respect to the beam direction for incoming charged pions of 3 GeV/c, 5 GeV/c, 8 GeV/c, 8.9 GeV/c (Be only), 12 GeV/c and 12.9 GeV/c (Al only) beam momentum impinging on thin (5% interaction length) beryllium, carbon, aluminium, copper, tin, tantalum and lead targets is described. The secondary pion yield is measured in a large angular and momentum range and double-differential FIG. 20: Dependence on the atomic number A of the pion production yields in π ± -Be, π ± -Al, π ± -C, π ± -Cu, π ± -Sn, π ± -Ta, π ± -Pb interactions averaged over the forward angular region (0.350 rad ≤ θ < 1.550 rad) and momentum region 100 MeV/c ≤ p < 700 MeV/c. The top-left panel (a) refers to π + production in π + beams, the top-right panel (b) to π − production in π + beams, while to bottom-left (c) and bottom-right (d) panels refer respectively to π + and π − production in π − beams.
FIG. 21: Comparison of HARP double-differential π ± production cross sections for π + -C at 3 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.

FIG. 22:
Comparison of HARP double-differential π ± production cross sections for π + -C at 5 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 23: Comparison of HARP double-differential π ± production cross sections for π + -C at 8 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 24: Comparison of HARP double-differential π ± production cross sections for π + -C at 12 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the QGSP model, the black solid line FTFP and the dashed line MARS.
FIG. 25: Comparison of HARP double-differential π ± production cross sections for π + -Ta at 3 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 26: Comparison of HARP double-differential π ± production cross sections for π + -Ta at 5 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 27: Comparison of HARP double-differential π ± production cross sections for π + -Ta at 8 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 28: Comparison of HARP double-differential π ± production cross sections for π + -Ta at 12 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the QGSP model, the black solid line FTFP and the dashed line MARS.
FIG. 29: Comparison of HARP double-differential π ± production cross sections for π − -C at 3 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 30: Comparison of HARP double-differential π ± production cross sections for π − -C at 5 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 31: Comparison of HARP double-differential π ± production cross sections for π − -C at 8 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFB and the dashed line MARS.
FIG. 32: Comparison of HARP double-differential π ± production cross sections for π − -C at 12 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see text for details). The left (right) panel shows π − (π + ) production. The gray line shows the Bertini model, the black solid line FTFP and the dashed line MARS. cross-sections are obtained. Results on the dependence of pion production on the atomic number A and incoming beam momentum are also presented. The comparisons of the π − /π + production ratios by pion beams of opposite polarity and proton beams show interesting features as a function of incoming beam momentum, target nucleus, and momentum of the secondary pions.
The use of a single detector for a range of beam momenta makes it possible to measure the dependence of the pion yield on the secondary particle momentum and emission angle θ with high precision. The A dependence of the crosssection can be studied, using data from a single experiment. Some hadronic production models (from GEANT4 and MARS) describing this energy range have been compared with our new results. The description is not yet satisfactory for heavy nuclei, while in some cases the carbon data are reasonably well reproduced. These new HARP data together with the incoming proton results published earlier [12] provide a unique possibility for validating and tuning hadron production models.

VII. ACKNOWLEDGEMENTS
Appendix A: Cross-section data TABLE VII: HARP results for the double-differential π + production cross-section in the laboratory system, d 2 σ π + /(dpdθ) for π + -Be interactions. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.  HARP results for the double-differential π − production cross-section in the laboratory system, d 2 σ π − /(dpdθ) for π + -Be interactions. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.  HARP results for the double-differential π + production cross-section in the laboratory system, d 2 σ π + /(dpdθ) for π − -Be interactions. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.
θ min θmax p min pmax  HARP results for the double-differential π − production cross-section in the laboratory system, d 2 σ π − /(dpdθ) for π − -Be interactions. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.  HARP results for the double-differential π + production cross-section in the laboratory system, d 2 σ π + /(dpdθ) for π + -C interactions. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.
θ min θmax p min pmax  HARP results for the double-differential π − production cross-section in the laboratory system, d 2 σ π − /(dpdθ) for π + -C interactions. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.