Comparison of large-angle production of charged pions with incident protons on cylindrical long and short targets

The HARP collaboration has presented measurements of the double-differential pi+/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 with proton beams hitting thin nuclear targets. In many applications the extrapolation to long targets is necessary. In this paper the analysis of data taken with long (one interaction length) solid cylindrical targets made of carbon, tantalum and lead is presented. The data were taken with the large acceptance HARP detector in the T9 beam line of the CERN PS. The secondary pions were produced by beams of protons with momenta 5 GeV/c, 8 GeV/c and 12 GeV/c. The tracking and identification of the produced particles were performed using a small-radius cylindrical time projection chamber (TPC) placed inside a solenoidal magnet. Incident protons were identified by an elaborate system of beam detectors. Results are obtained for the double-differential yields per target nucleon d2 sigma / dp dtheta. The measurements are compared with predictions of the MARS and GEANT4 Monte Carlo simulations.


I. INTRODUCTION
The main objectives of the HARP experiment [1] are to measure charged pion production yields to help design the proton driver of a future 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. The HARP experiment 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 and a large-angle spectrometer constructed in a solenoidal magnet.
A large amount of data collected by the HARP experiment with thin (5% of nuclear interaction length) and cryogenic targets have already been analyzed and published [12,13,14,15,16,17,18,19,20,21], covering all the physics subjects discussed above. Cross-sections for some of these data-sets based on the same raw data have been published by a different group [22]. We disagree with the analysis of Ref. [22] (see Ref. [23] for differences in detector calibration).
In this paper effective measurements of the double-differential production cross-section, d 2 σ π /dpdθ for π ± production valid for solid cylindrical long targets (one interaction length) made of carbon, tantalum and lead are presented. The secondary pions were produced by beams of protons in a momentum range from 5 GeV/c to 12 GeV/c impinging perpendicularly on one of the flat surfaces of the target rods. In earlier papers the measurements of the doubledifferential production cross-sections were presented for data taken using protons hitting thin beryllium, carbon, aluminium, copper, tin, tantalum and lead thin targets of 5% nuclear interaction length. Final results can be found in Ref. [17] for the proton beam data. The results presented in this paper provide a means to check the ability of hadron production models to simulate pion production with realistic extended targets by comparing short [43] and long target data taken with the same experiment. The choice of target materials covers an often used low A material, carbon, and examples of large A targets like tantalum and lead. We have limited ourselves to the momentum range of the incoming beam where the most statistics is available (5 GeV/c-12 GeV/c).
Effective measurements of particle production using extended targets are not unambiguously defined. In this case the results have to be understood as follows. The absorption and re-interactions of the beam proton in the target are not corrected for, thus the measurements encompass the full behaviour of the beam particle in the one interaction length (λ I ) of the target. This holds e.g. when the proton is (elastically) scattered or when protons or pions are produced in a relatively small forward cone (up to ≈200 mrad). As long as particles are produced within this forward cone they are regarded as part of the beam. Particles emitted outside this cone are regarded as "products". The effect of absorption and interactions of these secondary pions and protons emitted outside this cone are corrected for as losses of the produced pions or as background in the case of production of tertiaries. Thus, the analysis procedure aims to correct for these effects such that the effective target is transparent for the secondary "product" pions and one λ I long for the "beam particles". Some ambiguities still persist, and for the precise, quantitative interpretation of the results it is important to take into account the consequences of this procedure. One should note that the diameter of the targets is 30.3 mm both for the long and short targets. The relevant target parameters are listed in Table I. Nevertheless, we believe that this definition is the most useful one to provide input for the simulation of the production targets of super-beams and neutrino factories. Another procedure could be envisaged if one would know the target characteristics for a given installation exactly. In that case measurements on a replica target could be made and the particles emerging from the target could be registered without bothering what happened inside the target. Since these measurements are aiming at design studies for future installations, the latter approach cannot be followed.

Target material C Ta Pb
Length (cm) 38.01 11.14 17.05 Density (g/cm 3 ) 1.72 16.67 11.34 Data were taken in the T9 beam of the CERN PS. Contrary to the short-target data, no interaction trigger was applied while taking data with long targets. All good beam triggers were recorded, regardless of the activity in the downstream trigger detectors. The collected statistics for the different nuclear targets are reported in Table II. The analysis proceeds by selecting tracks in the Time Projection Chamber (TPC) in events with an incident proton. 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 yields. Otherwise, the analysis follows the same methods as used for the determination of π ± production by protons on short targets. These analysis methods are documented in Ref. [19] with improvements described in Ref. [17] and will be only briefly outlined here.
The long-target data were taken with relatively high beam intensities (typically 5000 particles per spill). Under these circumstances a space charge due to ions built up inside the drift volume of the TPC. These charges were responsible for a time-dependent ("dynamic") distortion of the track images measured by the TPC. These effects are more severe for the long-target data than for the short-target data due to the high interaction rate. A first set of results on pion production at large angles by protons using short-target data had been published [19,20,21] based on the analysis of the data limited to the beginning of each accelerator spill for this reason. It was no longer necessary to limit the analysis of the short-target data set to the events taken in the beginning of each spill after corrections had been developed [24]. These corrections were fully applied in the analysis of Ref. [17]. However, in the case of the long-target data the distortion effects are larger, and the corrections could only be reliably applied to a part of the 400 ms spill, so that typically 30%-50% of the statistics are available for the analysis reported in this paper. More details will be given in Section II A.

II. EXPERIMENTAL APPARATUS AND DATA SELECTION
The HARP detector is shown in Fig. 1 and is described in detail in Ref. [25]. The forward spectrometer, mainly used in the analysis for the conventional neutrino beams and atmospheric neutrino flux, comprises a dipole magnet, large planar drift chambers (NDC) [26], a time-of-flight wall (TOFW) [27], 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 [28]. A set of resistive plate chambers (RPC) forms a barrel inside the solenoid around the TPC to measure the arrival time of the secondary particles [29]. Charged particle identification (PID) can be achieved 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 provide an independent calibration of the PID based on dE/dx.  In addition to the data taken with the solid targets of 5% and 100% λ I runs were also taken with an empty target holder to check backgrounds. Data taken with a liquid hydrogen target at 3 GeV/c, 5 GeV/c and 8 GeV/c incident-beam momentum together with cosmic-ray data were used to provide an absolute calibration of the efficiency, momentum scale and resolution of the detector.
The momentum of the T9 beam is known with a precision of the order of 1% [30]. The beam profiles projected to the target position as measured by the beam MWPCs are shown for the three different beam momenta in Fig. 2. The absolute normalization of the number of incident protons is performed by just counting the incoming beam particles with the same trigger as used for the analysis of the secondary particles. This is possible since no selection on the interaction is performed in the trigger for the data sets used in the present analysis. The dimensions and mass of the solid targets were carefully measured. The purity of the target materials exceeded 99.9%. The uncertainties in thickness and density of the targets are well below 1%.
Beam instrumentation provides identification of the incoming particle, the determination of the time when it hits the target, and the impact point and direction of the beam particle on the target. It is based on a set of four multi-wire proportional chambers (MWPC) to measure the position and direction of the incoming beam particles and time-offlight (TOF) detectors and two N 2 -filled Cherenkov counters to identify incoming particles. Several trigger detectors are installed to select events with an interaction and to define the normalization. The beam of positive particles used for this measurement consists mainly of positrons, pions and protons, with 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. 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. 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 using a dedicated combination of the pressure setting of the two Cherenkov counters and are found to contribute less than 0.5%, and hence are negligible in the proton beam sample. Electrons radiate in the Cherenkov counters and would be counted as pions. More details on the beam particle selection can be found in Refs. [25] and [12,13,14].
The length of the accelerator spill is 400 ms with a typical intensity of 5000 beam particles per spill. The average number of events recorded by the data acquisition ranges up to 450 per spill for the data taken with long targets. The analysis proceeds by first selecting a beam proton not accompanied by other beam tracks. After the event selection the sample of tracks to be used for analysis is defined. Table II shows the number of events and the number of π ± selected in the analysis. The large difference between the first and second set of rows ("Total DAQ events" and "Accepted beam protons") is due to the relatively large fraction of pions in the beam and to the number of calibration triggers used. The entry "Fraction of triggers used" shows the part of the data for which distortions in the TPC could be calibrated reliably as explained below.

A. Distortions in the TPC
Besides the usual need for calibration of the detector, a number of hardware shortfalls, discovered mainly after the end of data-taking, had to be overcome to use the TPC data reliably in the analysis. The TPC is affected by a relatively large number of dead or noisy pads and static and dynamic distortions of the reconstructed trajectories. The corrections applied to the measurements are identical to the ones used for our analysis of short-target data taken with incoming protons and a description can be found in Ref. [17,24] [44].
The size of the corrections for dynamic distortions grows as function of the time within each accelerator spill and for each data set the part of the spill which can be reliably corrected is checked. The fraction of events usable for the analysis is typically 30%-50%, but varies for the different data sets (see Table II). The presence of a possible residual momentum bias in the TPC measurement due to the dynamic distortions was investigated using a large set of calibration methods. A dedicated paper [23] addresses this point and shows that our estimation of the knowledge of the absolute momentum scale is better than 3.5%. Due to the large event rate in the data taken with the long targets, the dynamic distortions are more severe than for the short-target data. It is possible to correct only a relatively small fraction of the data reliably, i.e. the first ≈30% of the Ta and Pb data and the first half of the C data. The track impact distance[45] with respect to the trajectory of the incoming beam particle, d ′ 0 , is a very sensitive probe to measure the distortion strength [24]. As an example of the effect of the distortions and the quality of the corrections the distributions of d ′ 0 in the 12 GeV/c C and Pb data are shown in Fig. 3, as examples of a better and worse situation, respectively. Figure 4 shows as a further check the distribution in q/p T , where q is the measured charge of the particle and p T its transverse momentum. Tracks have been divided into six groups depending on the number N evt of their event in the spill. The six groups correspond to 50n < N evt ≤ 50(n + 1) (for n ranging from zero to five), thus displaying the distribution of tracks from early events in the spill separately from tracks in events measured later in the spill. To make the absolute normalization meaningful, the distributions have been scaled to an equal number of incident-beam particles compared to the first group of 50 events. In the left panel, no dynamic distortion corrections have been applied and a clear difference of the distributions is visible especially for Pb. One should note that the momentum measurement as well as the efficiency is modified by the distortions. The right panel shows the distributions after the corrections. The distributions are more equal, although especially for the Pb data still important differences are observed at the end of the spill. To understand the asymmetry of positively and negatively charged tracks, one should keep in mind that no particle identification was performed. Thus both protons and pions contribute to the positives while the π − 's are the only component of the negative particles. Since the statistical errors in this analysis are smaller than the systematic errors a conservative approach was chosen, and only the first part of the spill where the dynamic distortion corrections could be applied was used. GeV/c C data in the top panels and for 12 GeV/c Pb data in the bottom panels. The left panels (a: p-C; c: p-Pb) show uncorrected data and in the right panels (b: p-C; d: p-Pb) dynamic distortion corrections have been applied. After the "default" correction for the static distortions (equal for each setting) a small residual effect at the beginning of the spill is visible at Nevt = 0 (left panel). This is due to the fact that the inner and outer field cages are powered with individual HV supplies. A setting-by-setting correction compatible with the reproducibility of the power supplies is applied for the data of the right panel together with the dynamic distortion correction. The value of d ′ 0 at Nevt = 0 in the right panel has a small negative value as expected from the fact that the energy-loss is not described in the track-model used in the fit. The difference observed in the results for the two charges shows that the model can correct the distortion in the 12 GeV/c C data to within about ±2 mm. In the bottom left panel (12 GeV/c Pb data) one observes a steep (almost linear) behaviour from Nevt ≈ 50 events, while the data early in the spill are practically not affected. In this case the distortion can be corrected up to about ±3 mm, as measured by d ′ 0 . However, beyond Nevt ≈ 150 events other benchmarks show that momentum biases are not kept under control.
FIG. 4: Distribution in q/pT for the 12 GeV/c C data (top panels) and for the 12 GeV/c Pb data (bottom panels), where q is the measured charge of the particle and pT its transverse momentum. The six curves show six regions in event number in spill (each in groups of 50 events in spill). Groups are labelled with the last event number accepted in the group, e.g. "50" stands for the group with event number from 1 to 50. The six groups are normalized to the same number of incoming beam particles, taking the first group as reference. Left panels (a: p-C; c: p-Pb): without dynamic distortion corrections; right panels (b: p-C; d: p-Pb): with dynamic distortion corrections. In the top left panel only the first three groups of 50 events in spill are equivalent, while in the top right panel the groups are nearly indistinguishable. In the bottom left panel (Pb) a very large difference between the groups of 50 events in spill are observed (only the first two groups are equivalent). The very large loss of tracks at high event numbers is due to the fact that particles no longer point back to the incident-beam particle and are rejected by this selection criterion. In the (bottom) right panel the first five groups are nearly indistinguishable.
FIG. 5: The average momentum as a function of Nevt observed for protons selected within a high dE/dx region for the p-Pb data at 12 GeV/c. The angle of the particles is restricted in a range with sin θ ≈ 0.9. The data are corrected for dynamic distortions and stay stable within 3% up to about 150 Nevt. The corrections become too large to be corrected reliably beyond 150 events (30% of the spill). The shaded band shows a ±3% variation.
observed. This is explained by the relatively high interaction rate using the long target and high multiplicity in p-Pb interactions without sufficient reduction of the beam intensity during data taking.
It cannot be taken for granted that no residual momentum bias is incurred when corrections have to be applied corresponding to d ′ 0 values larger than 15 mm. Therefore, the results of the corrections have to be checked using a benchmark which ensures good momentum reconstruction. A direct test of the effect of the correction on the measurement of momentum is shown in Fig. 5 for the worst case (12 GeV/c Pb). A sample of relatively high momentum protons was selected setting a fixed window with relatively high values of dE/dx in the TPC. This high dE/dx window ensured that the particles were correctly identified as protons and simultaneously selected a momentum band. It is visible that the momentum measurement starts to show deviations beyond event number 150 inside spills. Whereas the other estimators do not reveal a problem, the momentum estimator using dE/dx reveals a deviation beyond the accepted systematic error of ±3%. This explains the relatively small part of the full spill used in the analysis.

III. DATA ANALYSIS
Only a short outline of the data analysis procedure is presented here, for further details see Refs. [17,19]. The double-differential yieldper target nucleon 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.
The matrix M −1 ijαi ′ j ′ α ′ corrects for the efficiency and the resolution of the detector. It unfolds the true variables ijα 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 method used to correct for the various effects is described in more detail in Ref. [19].
In order 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 [33], 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. [19]. 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). As explained above, we do not make a correction for the attenuation of the beam in the target, so that the yields are valid for a λ I = 100% target. The result is normalized to the number of incident protons 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 has uncertainties smaller than 2% for all beam momentum settings.
The background due to interactions of the primary pions outside the target (called '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 are used for the subtraction, while the difference between π + and π − production is used to estimate the systematic error. In the region below 125 MeV/c a large fraction of the electrons can be unambiguously identified. These tracks are used as a relative normalization between data and simulation. An additional systematic error of 10% is assigned to the normalization of the π 0 subtraction using the identified electrons and positrons. 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. In the low momentum and large angle region the corrections for tertiary particles amount to 10%-15%.
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 overall systematic errors will be described in Section IV.

IV. EXPERIMENTAL RESULTS
The measured double-differential yields per target nucleon 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  6 to 8 for the three long targets studied here. 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 correlation of the statistical errors (introduced by the unfolding procedure) are typically smaller than 20% for adjacent momentum bins and smaller for adjacent angular bins. The correlations of III: Experimental uncertainties for the analysis of the data taken with carbon and tantalum targets in the 5 GeV/c, 8 GeV/c and 12 GeV/c beams. The numbers represent the uncertainty in percent of the cross-section integrated over the angle and momentum region indicated. The systematic errors for the Pb data are very similar to the Ta data.   the systematic errors are larger, typically 90% for adjacent bins. The correlation between systematic errors is shown in Table V for selected bins in the carbon data taken with the 12 GeV/c beam. The overall scale error (< 2%) is not shown. The double-differential results of this analysis are also tabulated in Appendix A. The integrated π − /π + ratio in the forward direction is displayed in Fig. 9 as a function of the secondary momentum. The ratios are very similar to the ones observed in the data with the short targets [17]. In the part of the momentum range shown in most bins more π + 's are produced than π − 's. The π − /π + ratio is larger for higher incoming beam momenta than for lower momenta and drops with increasing secondary momentum. 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 had already been observed in the short-target data. 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 [34]. A plausible explanation has been put forward in Ref. [35], where it was shown that this effect is due to an asymmetry in the production of ∆ resonances given the large neutron excess in these heavy nuclei.
The experimental uncertainties are summarized in Table III for the C and Ta targets. The systematic error breakdown is shown in Table IV for the C and Ta data taken with the 12 GeV/c beam. The numbers for the other momenta are very similar. The errors for the Pb target are very similar to the ones for the Ta target. The relative sizes of the different systematic error sources are very similar for π − and π + (only π + is shown) and for the different beam energies. Going from lighter (C) to heavier targets (Ta, Pb) the corrections for π 0 (conversion, concentrated at low secondary momentum), absorption, tertiaries and target pointing uncertainties are bigger. The discussion and figures shown in [17] give a reliable indication of the momentum and angular dependence of the systematic error components.
One observes that the statistical error is small compared to the systematic errors. 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 [46] and hence also contributes to the statistical error. This procedure almost doubles the statistical error, but avoids an important systematic error which would otherwise be introduced by assuming a cross-section model calculating 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 [19]. This error dominates at large momenta where the resolution is worse and in the bins with the largest angles where the distributions are steep. At low momentum in the forward direction the uncertainty in the subtraction of the electron and positron background due to π 0 production is large (∼ 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 absolute knowledge of the angular and the momentum resolution. The correction for tertiaries is relatively large at low momenta and large angles. Its uncertainty is responsible for a systematic error of ∼ 3 − 5% in these regions. Compared to the data taken with the shorter targets the systematic errors related to the efficiency of the cut requiring pions to point to the beam particle trajectory are much larger, especially for the Ta and Pb targets. This is not surprising given the large energy loss of these particles, especially in the forward direction when traversing a large part of the target.
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 particles counted in the normalization actually hit the target, with smaller components from the target density and beam particle counting procedure.
As a cross-check of the experimental procedures a comparison of the π + production yield measured with the first 0.1 λ I of the long carbon target and the full 0.05 λ I of the short carbon target was made. The carbon target is optimal for this cross-check owing to its lower density which favours the ratio of the resolution of the z measurement of the origin of the secondary particles and the interaction length. Nevertheless, an additional systematic error is introduced by selecting the outgoing particles from the partial target. The downstream boundary of the selection has to be made in a region with high particle population. The 8 GeV/c carbon data were chosen for their relatively high statistics. The choice to select the first 10% rather than the first 5% of the long target was made as a compromise between the reduced statistics, the additional systematic error due to the track selection near the boundary and the aim to remain as close as possible to the 5% λ I target. Due to the limited statistics once only 10% of the particles are selected, the comparison is only meaningful for π + production and for the first eight out of nine angular bins. The result is shown in Fig. 10. The ratio is compatible with unity within the relatively large uncertainties.

A. Comparisons with short target data
Our final results, obtained using the long targets can be compared with the previously presented data on short targets [17]. We stress here that the data sets using the two sets of targets have been analysed with methods which are as equal as possible. The results are shown in Figs. 11-19. The error bars are an estimate of the uncertainties of the ratios, taking into account the correlation of the common errors, such as momentum scale and assumptions on background subtraction.
The ratio of the long over short-target data can reveal the effect of the "degrading" of the incoming beam particles in the target. If all interacting protons would be absorbed in each interaction, the ratio would be roughly 0.65 [47]. If, on the other hand, the beam particles would continue unchanged, the ratio would be unity. One observes that in most cases the value of the ratio stays between the two limits mentioned above (0.65-1.0). It is closer to unity for the C data, especially at 12 GeV/c. In some part of the phase-space the C data displays a ratio above one. The Ta and Pb data are closer to 0.65. The behaviour as function of target material is perhaps understood by the fact that the p-C pion production cross-section at large angles is nearly independent of incoming momentum between 5 GeV/c and 12 GeV/c while the production cross-sections increase with beam momentum in this momentum range for heavier nuclei, as shown in Ref. [17]. Thus for p-C forward scattered protons with lower momentum are almost as effective in producing pions as the original beam proton. This is not the case for the heavier nuclei.
In Figs. 11-19 we also show comparisons of the ratio of the long-target data and the corresponding λ I = 5% targets with a selected list of Monte Carlo generators of GEANT4 [33] and the MARS [36] model. When a range of incident energies has to be simulated, GEANT4 offers the possibility to define a "physics list" of models allowing the user to describe the different energies with optimized models. The physics list defined for the present comparisons uses an intra-nuclear cascade model (the "Bertini model" [37,38]) for the lower energies. 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. At higher energies, instead, a parton string model (QGSP) [40] is used. The MARS code system [36] uses as basic model an inclusive approach multi-particle production originated by R. Feynman. Above 5 GeV phenomenological particle production models are used, while below 5 GeV a cascade-exciton model [41] combined with the Fermi break-up model, the coalescence model, an evaporation model and a multi-fragmentation extension are used instead.
An extra curve is shown for the C and Ta data representing the fraction of pions produced by the original ("first generation") beam proton compared to the overall number of pions produced by the beam. More explicitly, we repeat here that all particles produced by the primary beam proton within a forward cone with half-angle 0.2 rad are regarded as "beam particles".
FIG. 6: Double-differential yields per target nucleon for π + production (left) and π − production (right) in p-C interactions as a function of momentum displayed in different angular bins (shown in mrad in the panels). For better visibility the cross-sections have been scaled by a factor 0.5 (0.25) for the data taken at 8 GeV/c (5 GeV/c). The error bars represent the combination of statistical and systematic uncertainties.
FIG. 7: Double-differential yields per target nucleon for π + production (left) and π − production (right) in p-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. 8: Double-differential yields per target nucleon for π + production (left) and π − production (right) in p-Pb 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. 9: The ratio of the differential yields for π − and π + production in p-C (left panel), p-Ta (middle panel) and p-Pb (right panel) interactions as a function of the secondary momentum integrated over the forward angular region.
FIG. 10: Ratio of the double-differential π + production yields per target nucleon for p-C at 8 GeV/c measured with the first 0.1 λI of the 100% λI target and the 5% λI target for eight angular bins (shown in mrad).
The comparison, between data and models is good for the heavier nuclei, but discrepancies for carbon are visible for the comparison with MARS. Here MARS predicts a lower ratio than observed experimentally, while GEANT4 provides a more accurate description of the data. From an inspection of the line representing the ratio of pions produced by "beam particles" and all produced pions as calculated by MARS it appears that in the simulation the large majority of secondary pions are generated by the primary protons of the beam.
The ratio increases for all target materials with increasing beam momentum. This trend can be investigated by looking at the distribution of the the origin of pion tracks as a function of the longitudinal coordinate z. Figure 20 shows these distributions as a function of normalized z, z N , which runs from -1 to +1 over the length of the target. We define z as the value of the longitudinal coordinate of the point of closest approach of the secondary pion track and the beam particle [48]. The slopes of the z N distributions are steeper for the lower beam momenta. One also observes that there is a shower effect for the higher beam momenta, i.e. the maximum is not at the very beginning of the target indicating that some of the secondary particles are still efficient in producing high angle pions. Since the resolution is constant in absolute space coordinates, the resolution tail is longer for the higher density targets in the variable z N .

V. SUMMARY AND CONCLUSIONS
An analysis of the production of pions at large angles with respect to the beam direction for incoming protons of 5 GeV/c, 8 GeV/c and 12 GeV/c beam momentum impinging on long (100% interaction length) carbon, tantalum and lead targets is described. The secondary pion yield is measured in a large angular and momentum range and effective double-differential cross-sections are obtained. The present measurements are important to understand the simulations needed to design realistic pion production targets, e.g. for a neutrino factory [42]. The choice has been made to correct the absorption and re-interaction of secondary pions to make the results more universally usable for different target geometries. Thus the effect of the long target on the primary beam is highlighted.
The data are compared with corresponding results obtained with thin (5% λ I ) targets using the same detector and analysis techniques. It is observed that the effective attenuation of the incoming particle beam is larger for heavier targets than for carbon. This effect may have its origin in the stronger energy dependence of the production cross-section for the high-A targets.
The hadronic production models used for the comparison of the ratio compare well with the data for tantalum and lead. The GEANT4 model provides a good description for carbon, while MARS does not describe the relatively high ratio observed in the carbon data.

VI. ACKNOWLEDGEMENTS
We gratefully acknowledge the help and support of the PS beam staff and of the numerous technical collaborators who contributed to the detector design, construction, commissioning and operation. In particular, we would like to thank G. Barichello  Comparison of the ratio of double-differential π + (left) and π − (right) production yields per target nucleon for p-C at 5 GeV/c taken with the 100% λI and 5% λI target (circles) with MC predictions. The black curve represents the MARS prediction and the grey curve the GEANT4 simulation, while the dotted line shows the ratio of pions produced by a "first generation" beam proton and all pions as calculated by MARS. FIG. 13: Comparison of the ratio of double-differential π + (left) and π − (right) production yields per target nucleon for p-C at 12 GeV/c taken with the 100% λI and 5% λI target (circles) with MC predictions. The black curve represents the MARS prediction and the grey curve the GEANT4 simulation, while the dotted line shows the ratio of pions produced by a "first generation" beam proton and all pions as calculated by MARS.
FIG. 14: Comparison of the ratio of double-differential π + (left) and π − (right) production yields per target nucleon for p-Ta at 5 GeV/c taken with the 100% λI and 5% λI target (circles) with MC predictions. The black curve represents the MARS prediction and the grey curve the GEANT4 simulation, while the dotted line shows the ratio of pions produced by a "first generation" beam proton and all pions as calculated by MARS. FIG. 15: Comparison of the ratio of double-differential π + (left) and π − (right) production yields per target nucleon for p-Ta at 8 GeV/c taken with the 100% λI and 5% λI target (circles) with MC predictions. The black curve represents the MARS prediction and the grey curve the GEANT4 simulation, while the dotted line shows the ratio of pions produced by a "first generation" beam proton and all pions as calculated by MARS. FIG. 16: Comparison of the ratio of double-differential π + (left) and π − (right) production yields per target nucleon for p-Ta at 12 GeV/c taken with the 100% λI and 5% λI target (circles) with MC predictions. The black curve represents the MARS prediction and the grey curve the GEANT4 simulation, while the dotted line shows the ratio of pions produced by a "first generation" beam proton and all pions as calculated by MARS. [43] In previous papers we used the term "thin" instead.
[44] The results of Ref. [22] have been analysed by another group with a treatment of the dynamic distortions [31] with which we disagree.
[45] The d ′ 0 sign indicates if the helix encircles the beam particle trajectory (positive sign) or not (negative sign).
[46] The migration matrix is calculated without prior knowledge of the cross-sections, while the unfolding procedure determined the unfolding matrix from the migration matrix and the distributions found in the data.
[47] This follows from the ratio of the integral of the exponential degrading of the beam for 1λI and 0.05λI targets,  VI: HARP results for the double-differential π + production yield per target nucleon in the laboratory system, d 2 σ π + /(dpdθ) for p-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.  HARP results for the double-differential π − production yield per target nucleon in the laboratory system, d 2 σ π − /(dpdθ) for p-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 d 2 σ π − /(dpdθ) (rad) (rad) (GeV/c) (GeV/c) (barn/(GeV/c · rad)) GeV/c GeV/c GeV/c 0.  HARP results for the double-differential π + production yield per target nucleon in the laboratory system, d 2 σ π + /(dpdθ) for p-Ta 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 d 2 σ π + /(dpdθ) (rad) (rad) (GeV/c) (GeV/c) (barn/(GeV/c · rad)) GeV/c GeV/c GeV/c 0. 35   HARP results for the double-differential π − production yield per target nucleon in the laboratory system, d 2 σ π − /(dpdθ) for p-Ta 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 d 2 σ π − /(dpdθ) (rad) (rad) (GeV/c) (GeV/c) (barn/(GeV/c · rad)) GeV/c GeV/c GeV/c 0. 35   HARP results for the double-differential π + production yield per target nucleon in the laboratory system, d 2 σ π + /(dpdθ) for p-Pb 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.