Forward production of charged pions with incident protons on nuclear targets at the CERN PS

Measurements of the double-differential charged pion production cross-section in the range of momentum 0.5 GeV/c<p<8.0 GeV/c and angle 0.025 rad<theta<0.25 rad in collisions of protons on beryllium, carbon, nitrogen, oxygen, aluminium, copper, tin, tantalum and lead are presented. The data were taken with the large acceptance HARP detector in the T9 beam line of the CERN PS. Incident particles were identified by an elaborate system of beam detectors. The data were taken with thin targets of 5% of a nuclear interaction length. The tracking and identification of the produced particles was performed using the forward system of the HARP experiment. Results are obtained for the double-differential cross section mainly at four incident proton beam momenta (3 GeV/c, 5 GeV/c, 8 GeV/c and 12 GeV/c). Measurements are compared with the GEANT4 and MARS Monte Carlo generators. A global parametrization is provided as an approximation of all the collected datasets which can serve as a tool for quick yields estimates.

Covering an extended range of solid targets in the same experiment, it is 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. Pion production data at low momenta (≤ 25 GeV/c) are extremely scarce and HARP is the first experiment to provide a large data set, taken with many different targets, full particle identification and large detector acceptance. These data, together with the ones published in reference [13], are also of great interest for corrections of secondary interactions in detector studies in particle physics experiments.
Existing data are mainly at fixed production angles with Be targets [14] and are affected usually by large uncertainties. Only two experiments [15,16] at 19.2 GeV/c and 24 GeV/c were performed with an extended set of nuclear targets. The E910 experiment at BNL has recently published data at 12.3 GeV/c and 17.5 GeV/c with protons on Be, Cu and Au targets [17].
Data were taken in the T9 beam of the CERN PS.

II. EXPERIMENTAL APPARATUS
The HARP experiment makes use of a largeacceptance spectrometer consisting of a forward and a large-angle detection system. The HARP detector is shown in Fig. 1 and is fully described in Ref. [18]. The forward spectrometer -based on five modules of large area drift chambers (NDC1-5) [19] and a dipole magnet complemented by a set of detectors for particle identification (PID): a time-of-flight wall (TOFW) [20], a large Cherenkov detector (CHE) and an electromagnetic calorimeter (ECAL) -covers polar angles up to 250 mrad, which is well matched to the angular range of interest for the measurement of hadron production to calculate the properties of conventional neutrino beams. As the most downstream detection element visible in the figure, the 1.4 m wide beam muon identifier (BMI) is used to measure the muon contamination in the beam.
The discrimination power of TOFW time-of-flight below 3 GeV/c and the Cherenkov dectector above 3 GeV/c are combined to provide powerful separation of forward pions and protons. The calorimeter is used only for separating pions and electrons when characterizing the response of the other detectors. The resulting pion identification efficiency is bigger than 96% (98%) for pion momenta larger than 0.5 (3.0) GeV/c, with purity around 99.5%.
The large-angle spectrometer -based on a Time Projection Chamber (TPC) and Resistive Plate Chambers (RPCs) located inside a solenoidal magnet -has a large acceptance in the momentum and angular range for the pions relevant to the production of the muons in a neutrino factory. For the analysis described here only the forward spectrometer and the beam instrumentation are used.
The HARP experiment took data in 2001 and 2002. The momentum definition of the T9 beam is known with a precision of the order of 1% [21].
The target is placed inside the inner field cage of the TPC, in an assembly that can be moved in and out of the solenoid magnet. The solid targets used for the measurements reported here have a cylindrical shape with a nominal diameter of about 30 mm. Their thickness along the beam is equivalent to about 5% of an interaction length λ I [41].
A set of four multi-wire proportional chambers (MWPCs) measures the position and direction of the incoming beam particles with an accuracy of ≈1 mm in position and ≈0.2 mrad in angle per projection. A beam time-of-flight system (BTOF) measures the time difference of particles over a 21.4 m path-length. It is made of two identical scintillation hodoscopes, TOFA and TOFB (originally built for the NA52 experiment [22]), which, together with a small target-defining trigger counter, also used for the trigger, provide particle identification at low energies. This provides separation of pions, kaons and protons up to 5 GeV/c and determines the initial time at the interaction vertex (t 0 ). The timing resolution of the combined BTOF system is about 70 ps. A system of two N 2 -filled Cherenkov detectors (BCA and BCB) is used to tag electrons at low energies and pions at higher energies. The electron and pion tagging efficiency is found to be close to 100%. A set of trigger detectors completes the beam instrumentation.
The positive beam used in this analysis contains mainly positrons, pions and protons, with a small admixture of kaons and deuterons and heavier ions. The proton fraction in the incoming beam goes from ∼ 35% at 3 GeV/c to ∼ 92% at 12 GeV/c. The selection of beam protons is done requiring a pulse height consistent with the absence of signal in both beam Cherenkov detectors (BCA and BCB). At the lowest beam energy, 3 GeV/c, the BTOF is used to reject pions which at that momentum do not produce Cherenkov light. At 5 GeV/c, the π/p separation is obtained by the BTOF system and only one of the Cherenkov (usually BCB), while the other is used to tag e ± . Ions are rejected by the BTOF system at all momenta. In the 12.9 GeV/c beam the two Cherenkov counters were operated with different pressures to make it possible to measure the kaons separately from the protons and pions. After tagging the beam contamination from all particle species is negligible. Only events with a single reconstructed beam track in the four MW-PCs, good timing measurements in BTOF and no signal in the beam halo counters are accepted.
A downstream trigger in the forward scintillator trigger plane was required to record the event, accepting only tracks with a trajectory outside the central hole (60 mm diameter) which allows beam particles to pass. Accepting only tracks with a trajectory outside the central hole, the efficiency is measured to be >99.8%. Particle identification is done later, at the analysis stage, via the downstream PID detectors.
The length of the accelerator spill is 400 ms with a typical intensity of 15 000 beam particles per spill, and less for the lower momentum settings. The average number of events recorded by the data acquisition ranges from 300 to 350 per spill.
The absolute normalization of the cross-section was performed using 'incident-proton' triggers. These are triggers where the same selection on the beam particle was applied but no selection on the interaction was performed. The rate of this trigger was down-scaled by a factor 64. These unbiased events are used to determine N pot in the cross-section formula (1), see later.

III. DATA ANALYSIS
The analysis procedure is similar to the one reported in Refs. [4], [5], [2] and [13] where with respect to our first paper on forward pion production in p-Al interactions [1], a number of improvements to the analysis techniques and detector simulation had been made. Tracks are reconstructed in the drift chambers downstream of the magnet. Only tracks with at least one hit in the TOFW are accepted for the analysis. Hits are searched for in the Cherenkov detector consistent with these tracks to complete the particle identification. Secondary track selection criteria are optimized to ensure the quality of momentum reconstruction and a clean time-of-flight measurement while maintaining a high reconstruction efficiency. The method to combine the information from these PID detectors is described in [23]. In the kinematic range of the current analysis the pion identification efficiency is about 98%, while the background from mis-identified protons is well below 1%. The momentum reconstruction is performed by means of a special implementation [24] of the Kalman filter [25] which uses the position of the target as a constraint. For the final stage of the analysis the same unfolding technique as in Refs. [4], [5] and [13] has been applied.
The background induced by interactions of beam particles in the materials outside the target is measured by taking data without a target in the target holder ("empty target data"). These data are subject to the same event and track selection criteria as the standard data sets. To take into account this background the number of particles of the observed particle type in the "empty target data" are subtracted bin-by-bin (momentum and angular bins) from the number of particles of the same type. The uncertainty induced by this method is labeled "empty target subtraction". The collected event statistics on the different solid targets is summarised in Table I.

A. Cross-section calculation
The double-differential cross-section is calculated as follows: is the cross-section in mb/(GeV/c sr) for the produced particle type α (p, π + or π − ) for each true momentum and angle bin (p i , θ j ) covered in this analysis; is the number of reconstructed particles of type α ′ in bins of reconstructed momentum p ′ i and angle θ ′ j in the raw data, after subtraction of empty target data (due to beam protons interacting in material other than the nuclear target). These particles must satisfy the event, track and PID selection criteria.
• M cor pθαp ′ θ ′ α ′ is the correction matrix which accounts for the finite efficiency and resolution of the detector. It unfolds the true variables p i , θ j , α from the reconstructed variables p ′ i , θ ′ j , α ′ and corrects the observed We do not make a correction for the attenuation of the proton beam in the target, so that strictly speaking the cross-sections are valid for λ I = 5% targets.
The calculation of the correction matrix M cor piθj αp ′ i θ ′ j α ′ is done with the unfolding method introduced by D'Agostini [26] [44]. This method has been used in the recent HARP publications [4,5,13] and it is also applied in the analysis described here.
The Monte Carlo simulation of the HARP setup is based on GEANT4 [27]. The detector materials are accurately described in this simulation as well as the relevant features of the detector response and the digitization process. All relevant physics processes are considered, including multiple scattering, energy loss, absorption and re-interactions. The simulation is independent of the beam particle type because it only generates for each event exactly one secondary particle of a specific particle type inside the target material and propagates it through the complete detector. A small difference (at the few percent level) is observed between the efficiency calculated for events simulated with the single-particle Monte Carlo and with a simulation using a multi-particle hadron-production model. A similar difference is seen between the singleparticle Monte Carlo and the efficiencies measured directly from the data. A momentum-dependent correction factor determined using the efficiency measured with the data is applied to take this into account. The track reconstruction used in this analysis and the simulation are identical to the ones used for the π + production in p-Be collisions [2]. A detailed description of the corrections and their magnitude can be found there.
The reconstruction efficiency (inside the geometrical acceptance) is larger than 95% above 1.5 GeV/c and drops to 80% at 0.5 GeV/c. The requirement of a match with a TOFW hit has an efficiency between 90% and 95% independent of mo-mentum. The electron veto rejects about 1% of the pions and protons below 3 GeV/c with a remaining background of less than 0.5%. Below Cherenkov threshold ( ∼ 3 GeV/c) the TOFW separates pions and protons with negligible background and an efficiency of ≈98% for pions. Above Cherenkov threshold the efficiency for pions is greater than 99% with only 1.5% of the protons mis-identified as a pion. The kaon background in the pion spectra is smaller than 1%.
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, thus giving rise to tracks in the dipole spectrometer. In such cases also the additional tracks are entered into the migration matrix thereby taking into account the combined effect of the generated particle and any secondaries it creates. The absorption correction is on average 20%, approximately independent of momentum. Uncertainties in the absorption of secondaries in the dipole spectrometer material are taken into account by a variation of 10% of this effect in the simulation. The effect of pion decay is treated in the same way as the absorption and is 20% at 500 MeV/c and negligible at 3 GeV/c. The uncertainty in the production of background due to tertiary particles is larger. The average correction is ≈10% and up to 20% at 1 GeV/c. The correction includes reinteractions in the detector material (mainly carbon) as well as a small component coming from reinteractions in the target. The subtraction may be computed by MonteCarlo simulations and as most of the encountered material is carbon, the check of the inelastic interactions of low-energy protons or pions in carbon is essential. The validity of the generators used in the simulation was checked by the analysis of our data with incoming protons and charged pions on aluminium and carbon targets at lower momenta (3 GeV/c and 5 GeV/c). A 30% uncertainty on the secondary production was assumed.
The average empty-target subtraction amounts to ≈20%.
Owing to the redundancy of the tracking system downstream of the target the detection efficiency is very robust under the usual variations of the detector performance during the long data taking periods. Since the momentum is reconstructed without making use of the upstream drift chamber module (which is more sensitive in its performance to the beam intensity) the reconstruction efficiency is uniquely determined by the downstream system. No variation of the overall efficiency has been observed. The performance of the TOFW and CHE system have been monitored to be constant for the data taking periods used in this analysis. The calibration of the detectors was performed on a dayby-day basis.

B. Error estimation
The total statistical error takes into account the direct error propagation of the statistical errors in the raw data and the statistical error which is incurred while obtaining the unfolding matrix from the data. The latter component increases the direct error by a factor two. The procedure is outlined in references [4], [28].
Different types of sources induce systematic errors for the analysis described here: track yield corrections (∼ 5%), particle identification (∼ 0.1%), momentum and angular reconstruction (∼ 1%) [45]. The strategy to calculate these systematic errors and the different methods used for their evaluation are described in [4]. An additional source of error is due to misidentified secondary kaons, which are not considered in the particle identification method used for this analysis and are subtracted on the basis of a Monte Carlo simulation, as in [4]. No explicit correction is made for pions coming from decays of other particles created in the target, as they give a very small contribution according to the selection criteria applied in the analysis. As a result of these systematic error studies, each error source can be represented by a covariance matrix. The sum of these matrices describes the total systematic error, as explained in [4].
The experimental uncertainties are shown for a typical light target (Be) in Figure 2 for π + at some typical incident beam momenta and in Table II for another target (C) at the incident beam momentum 12 GeV/c. They are very similar for π − and at the other beam energies. Going from lighter (Be, C) to heavier targets (Ta, Pb) the corrections for π 0 (conversion) and absorption/tertiaries increase.
The overall normalization has an uncertainty of ∼ 2% and is mainly due to the uncertainty in the efficiency that beam protons counted in the normalization actually hit the target, with smaller components from the target density and the beam particle counting procedure.
On average the total integrated systematic error is around 5 − 6%, with a differential bin to bin systematic error of the order of 10 − 11%, to be compared with a statistical integrated (bin-to-bin differential) error of ∼ 2 − 3% (∼ 10 − 13%). Systematic and statistical errors are roughly of the

IV. 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 9 for solid 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. 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 overall normalization error (< 2%) is not included in the error bars [46]. The results of this analysis are also fully tabulated and reported in appendix A.
An overall fit with a Sanford-Wang parametrization [29] has been done by using all solid targets and available beam momenta data (p ≥ 5 GeV/c), see Section IV A for details, and is shown as a solid line on all figures.
The dependence of the averaged pion yields on the incident beam momentum is shown in Fig. 10. The π − yields, averaged over two angular regions (0.05 rad ≤ θ < 0.15 rad and 0.15 rad ≤ θ < 0.25 rad) and four momentum regions (0.5 GeV/c ≤ p < 1.5 GeV/c, 1.5 GeV/c ≤ p < 2.5 GeV/c, 2.5 GeV/c ≤ p < 3.5 GeV/c and 3.5 GeV/c ≤ p < 4.5 GeV/c), are shown in the left panel and the π + data averaged over the same regions in the right panel, for four different beam momenta. Whereas the beam energy dependence of the yields in the p-Be, p-C data differs clearly from the dependence in the p-Ta, p-Pb data one can observe that the p-Al, p-Cu and p-Sn data display a smooth transition between them. The dependence in the p-Be, p-C data is much more flat with a saturation of the yield between 8 GeV/c and 12 GeV/c with the p-Al, p-Cu and p-Sn showing an intermediate behaviour.
The dependence of the averaged pion yields on the atomic number A is instead shown in Fig. 11. The π − yields, averaged over two angular regions (0.05 rad ≤ θ < 0.15 rad and 0.15 rad ≤ θ < 0.25 rad) and four momentum regions (0.5 GeV/c ≤ p < 1.5 GeV/c, 1.5 GeV/c ≤ p < 2.5 GeV/c, 2.5 GeV/c ≤ p < 3.5 GeV/c and 3.5 GeV/c ≤ p < 4.5 GeV/c), are shown in the left panel and the π + data averaged over the same regions 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 latter saturating earlier towards higher A, especially at lower beam momenta.
FIG. 3: p-Be differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 4: p-C differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 5: p-Al differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 6: p-Cu differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 7: p-Sn differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 8: p-Ta differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 9: p-Pb differential cross sections: left panel π − production, right panel π + production. The curves represent the global parametrization as described in the text. In the top right corner of each plot the covered angular range is shown in mrad.
FIG. 10: The dependence on the beam momentum of the π − (top panel) and π + (bottom panel) production yields in p-Be, p-C, p-Al, p-Cu, p-Sn, p-Ta, p-Pb interactions averaged over two forward angular regions (0.05 rad ≤ θ < 0.15 rad and 0.15 rad ≤ θ < 0.25 rad) and four momentum regions (0.5 GeV/c ≤ p < 1.5 GeV/c, 1.5 GeV/c ≤ p < 2.5 GeV/c, 2.5 GeV/c ≤ p < 3.5 GeV/c and 3.5 GeV/c ≤ p < 4.5 GeV/c), for the four different incoming beam energies. Data points for different target nuclei and equal momenta are slightly shifted horizontally with respect to each other to increase the visibility.

A. Pion production data parametrization
At low energies, it is common to use the empirical data parametrization for pion production in proton-nucleus interactions originally developed by Sanford and Wang [29]. This parametrization has the functional form: where: X denotes any system of other particles in the final state, p beam is the proton beam momentum in GeV/c, p and θ are the π ± momentum and angle in units of GeV/c and radians, respectively, d 2 σ/(dpdΩ) is expressed in units of mb/(GeV/c)/sr and the parameters c 1 , . . . , c 8 are obtained from fits to the meson production data.
The parameter c 1 is an overall normalization factor, the four parameters c 2 , c 3 , c 4 , c 5 describe the momentum distribution of the secondary pions in the forward direction, and the three parameters c 6 , c 7 , c 8 describe the corrections to the pion momentum distribution for pion production angles that are different from zero. The π ± production data on solid targets reported here have been fitted simultaneously to the empirical Sanford-Wang formula. In the χ 2 minimization, the full error matrix was used.
To go from a baseline nuclear targets (typically Be) to another nuclear target A, a correction factor where x F is the Feynman x, see references [10] and [30] for details. A 9-parameter (11-parameter) fit was done over 24 (25) π − (π + ) datasets [47], corresponding to 1440 (1472) experimental points. The goodness-of-fit of the Sanford-Wang parametrization hypothesis for the HARP results can be assessed by considering the best-fit χ 2 value of χ 2 min = 13030 (8061) for 1431 (1461) degrees of freedom for the π − (π + ) production, indicating a very poor fit quality. In particular, inspection of the HARP inclusive pion production doubledifferential cross-section, and resulting Sanford-Wang parametrization, points to a description of the ratio g(θ) of the pion momentum distribution at θ = 0 with respect to the θ = 0 pion momentum distribution that is more complicated than TABLE III: Sanford-Wang parameters and errors obtained by fitting the π + (π − ) datasets. The errors refer to the 68.27% confidence level for eleven parameters (∆χ 2 = 12.6) for π + and nine parameters (∆χ 2 = 10.4) for π − .
what can be accommodated within the Sanford-Wang formula, where this ratio is given by The overall fit may be used as a fast approximation of HARP data valid within a factor 2-3 of the quoted experimental errors. The best-fit values of the parameters are reported in Table III together with their errors. The fit parameter errors are estimated by requiring ∆χ 2 ≡ χ 2 −χ 2 min = 12.6 (10.4), corresponding to the 68.27% confidence level region for eleven (nine) variable parameters. Significant correlations among fit parameters are found, as shown by the correlation matrix given in Tables IV and V.
To show the trend of the Sanford-Wang global fit of all HARP datasets, figure 12 reports the comparison, at 8 GeV/c and 12 GeV/c, between pion production data and the above parametrization.
For the 8.9 GeV/c MiniBooNE/SciBooNE beamline and the 12.9 GeV/c K2K beamline two adhoc Sanford-Wang parametrizations, using only the relevant HARP datasets, have been published in references [1] and [2]. Given the poor description of HARP pion production data in terms of the original Sanford-Wang parametrization, one extra parameter to better describe the angular dependence was introduced in the fit reported in reference [2]. In the global fit presented here, this extra parameter is not used for simplicity, as it was found that it did not improve the fit quality in a significant way.
As a final remark, we stress again that due to the poor fit quality our global fit may be just considered as a simple way to summarize an extended set    of data (∼ 1000 experimental data points) using a formula with about ten parameters. This may be useful in the intial phase of an experiment design or a Monte Carlo validation.

B. Comparison with Monte Carlo generators.
In the following we will show only some comparisons with two widely available Monte Carlo generators: GEANT4 version 7.1 [27] and MARS version 15.07 [31], using different models. The comparison will be shown for a limited set of plots and only for the Be and Ta targets, as examples of a light and a heavy target. In both generators, no single model is applicable to all energies and a transition between low energy models and high energy models, at about 5 GeV/c-10 GeV/c, has to be done.
At intermediate energies (up to 5 GeV/c-10 GeV/c), GEANT4 uses two types of intranuclear cascade models: the Bertini model [32,33] (valid up to ∼10 GeV/c) and the Binary model [34] (valid up to ∼3 GeV/c). 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 [35] 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, two parton string models, the quark-gluon string (QGS) model [32,36] and the Fritiof (FTP) model [36] are used, in addition to a High Energy Parametrized model (HEP) derived from the high energy part of the GHEISHA code used inside GEANT3 [37]. 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.
A realistic GEANT4 simulation is built by combining models and physics processes into what is called a "physics list". Computing the χ 2 between models and data themselves, where a systematic error going from 0% to 50% have been added to simulation results we obtain the results shown in Table VI. Normalization factors, between data and Monte Carlo simulations, are shown instead in Table VII.
Over the full energy range covered by the HARP experiment, the best comparison is obtained with the MARS Monte Carlo. This may be owing to the fact that MARS is using different models in different energy regions, equivalent to use a collection of models as implemented in the "physics lists" of GEANT4.
The full set of HARP data, taken with targets spanning the full periodic table of elements, with small total errors and full coverage of the solid angle in a single detector may help the validation of models used in hadronic simulations in the difficult energy range between 3 and 15 GeV/c of incident momentum, as done in reference [40].

V. SUMMARY AND CONCLUSIONS
In this paper we report our final results on doubledifferential cross sections for the production of positive and negative pions in the kinematic range 0.5 GeV/c≤ p π ≤ 8 GeV/c and 0.025 rad ≤ θ π ≤ 0.25 rad from the collisions of protons of 3, 5, 8 and 12 GeV/c on beryllium, aluminium, carbon, copper, tin, tantalum targets of 5% λ I . In addition results at 8.9 GeV/c (Be only), 12.9 GeV/c (Al only) and 12 GeV/c (N 2 and O 2 ) are reported.
A parametrization, inspired by the Sanford-Wang formula, of all our datasets at all energies is also presented. This may be used as a fast approximation to HARP data, valid within a factor 2-3 of the quoted errors.
The pion yield averaged over different momentum and angular ranges increase smoothly with the atomic number A of the target and with the energy of the incoming proton beam. The A-dependence is slightly different for π − and π + production, the latter saturating earlier towards higher A, especially at lower beam momenta.
Comparisons with GEANT4 and MARS generators are presented.

VI. ACKNOWLEDGMENTS
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 [44] The unfolding method tries to put in correspondence the vector of measured observables (such as particle momentum, polar angle and particle type) xmeas with the vector of true values xtrue using a migration matrix: xmeas = Mmigr × xtrue. The goal of the method is to compute a transformation (correction matrix) to obtain the expected values for xtrue from the measured ones. The most simple and obvious solution, based on simple matrix inversion M −1 migr , is usually unstable and is dominated by large variances and strong negative correlations between neighbouring bins. In the method of D' Agostini, the correction matrix M UFO tries to connect the measurement space (effects) with the space of the true values (causes) using an iterative Bayesian approach, based on Monte Carlo simulations to estimate the probability for a given effect to be produced by a certain cause.
[45] The quoted error in parenthesis refers to fractional error of the integrated cross-section (δ π int (%)) in the kinematic range covered by the HARP experiment [46] This makes it possible to calculate e.g. integrated particle ratios taking these normalization errors into account only when applicable, i.e. when different beams are compared [47] In the fit procedure only data with p beam ≥ 5 GeV/c were used, as the inclusion of the 3 GeV/c data gave problems in the convergence. The c7 parameter was fixed to zero in the π − fit, as the results were found insensitive to this parameter in an extended range around zero [48] Also this model, at low energy, has its root in the GHEISHA code inside GEANT3.

APPENDIX A: CROSS-SECTION DATA
The following tables report the measured differential cross-section for positive and negative pion production in interactions of 3, 5, 8 and 12 GeV/c momentum protons on different types of nuclear targets. The data are presented in the kinematic range of 0.5 GeV/c ≤ p π ≤ 8 GeV/c and 0.05 rad ≤ θ π ≤ 0.25 rad. The overall normalization uncertainty (≤ 2%) is not included in the reported errors.
Results at higher incoming beam momenta, from 8 GeV/c to 12 GeV/c, are also presented extending to a lower value of the polar angle θ (0.025 rad), using a finer binning.
The cross-section values in some bins of the following tables are omitted and replaced by the symbol * ± * , due to instabilities in the unfolding procedure coming from the low statistics available. This appears especially evident for the 3 GeV/c data where the total negative final state pion statistics is often well below one hundred events. HARP results for the double-differential π + production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-Be interactions at 3, 5, 8, 8.9, 12 GeV/c. 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/(sr GeV/c)) 3 GeV/c 5 GeV/c 8 GeV/c 8.  HARP results for the double-differential π − production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-Be interactions at 3, 5, 8, 8.9, 12 GeV/c. 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 p-C interactions at 3, 5, 8, 12 GeV/c. 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/(sr GeV/c)) 3 GeV/c  HARP results for the double-differential π ± production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p − N2 interactions at 12 GeV/c. 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/(sr GeV/c)) 3 GeV/c  HARP results for the double-differential π − production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-Cu interactions at 3, 5, 8, 12 GeV/c. 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 p-Sn interactions at 3, 5, 8, 12 GeV/c. 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/(sr GeV/c)) 3 GeV/c  HARP results for the double-differential π + production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-Be interactions at 8,8.9,12 GeV/c. 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/(sr GeV/c)) 8 GeV/c 8.  HARP results for the double-differential π + and π − production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-C interactions at 8 and 12 GeV/c. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. A finer angular binning than in the previous set of tables is used. 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/(sr GeV/c)) 8 GeV/c 12 GeV/c 12.  HARP results for the double-differential π + and π − production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-Cu interactions at 8 and 12 GeV/c. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. A finer angular binning than in the previous set of tables is used. 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Ω) d 2 σ π − /(dpdΩ) (rad) (rad) (GeV/c) (GeV/c) (barn/(sr GeV/c)) (barn/(sr GeV/c)) 8  HARP results for the double-differential π + and π − production cross-section in the laboratory system, d 2 σ π /(dpdΩ), for p-Sn interactions at 8 and 12 GeV/c. Each row refers to a different (p min ≤ p < pmax, θ min ≤ θ < θmax) bin, where p and θ are the pion momentum and polar angle, respectively. A finer angular binning than in the previous set of tables is used. The central value as well as the square-root of the diagonal elements of the covariance matrix are given.