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

M. Apollonio,1,* A. Artamonov,2,† A. Bagulya,3 G. Barr,4 A. Blondel,5 F. Bobisut,6,7 M. Bogomilov,8 M. Bonesini,9,‡ C. Booth,10 S. Borghi,5,§ S. Bunyatov,11 J. Burguet-Castell,12 M. G. Catanesi,13 A. Cervera-Villanueva,12 P. Chimenti,1 L. Coney,14,‖ E. Di Capua,15 U. Dore,16 J. Dumarchez,17 R. Edgecock,18 M. Ellis,18,¶ F. Ferri,9 U. Gastaldi,19 S. Giani,2 G. Giannini,1 D. Gibin,6,7 S. Gilardoni,2 P. Gorbunov,2,† C. Gößling,20 J. J. Gómez-Cadenas,12 A. Grant,2 J. S. Graulich,21,** G. Grégoire,21 V. Grichine,3 A. Grossheim,2,†† A. Guglielmi,6 L. Howlett,10 A. Ivanchenko,2,‡‡ V. Ivanchenko,2,§§ A. Kayis-Topaksu,2,‖‖ M. Kirsanov,22 D. Kolev,8 A. Krasnoperov,11 J. Martı́n-Albo,12 C. Meurer,23 M. Mezzetto,6 G. B. Mills,24,‖ M. C. Morone,5,¶¶ P. Novella,12 D. Orestano,25,26 V. Palladino,27 J. Panman,2 I. Papadopoulos,2 F. Pastore,25,26 S. Piperov,28 N. Polukhina,3 B. Popov,11,*** G. Prior,5,§ E. Radicioni,13 D. Schmitz,14,‖ R. Schroeter,5 V. Serdiouk,11 G. Skoro,10 M. Sorel,12 E. Tcherniaev,2 P. Temnikov,28 V. Tereschenko,11 A. Tonazzo,25,26 L. Tortora,25 R. Tsenov,8 I. Tsukerman,2,† G. Vidal-Sitjes,15,* C. Wiebusch,2,††† and P. Zucchelli2,‡‡‡ (HARP Collaboration) 1Università degli Studi e Sezione INFN, Trieste, Italy 2CERN, Geneva, Switzerland 3P. N. Lebedev Institute of Physics (FIAN), Russian Academy of Sciences, Moscow, Russia 4Nuclear and Astrophysics Laboratory, University of Oxford, United Kingdom 5Section de Physique, Université de Genève, Switzerland 6Sezione INFN, Padova, Italy 7Universitá degli Studi, Padova, Italy 8Faculty of Physics, St. Kliment Ohridski University, Sofia, Bulgaria 9Sezione INFN Milano Bicocca, Milano, Italy 10Department of Physics, University of Sheffield, United Kingdom 11Joint Institute for Nuclear Research, JINR Dubna, Russia 12Instituto de Fı́sica Corpuscular, IFIC, CSIC and Universidad de Valencia, Spain 13Sezione INFN, Bari, Italy 14Columbia University, New York, New York, USA 15Università degli Studi e Sezione INFN, Ferrara, Italy 16Università “La Sapienza” e Sezione INFN Roma I, Rome, Italy 17LPNHE, Universités de Paris VI et VII, Paris, France 18Rutherford Appleton Laboratory, Chilton, Didcot, United Kingdom 19Laboratori Nazionali di Legnaro dell’ INFN, Legnaro, Italy 20Institut für Physik, Universität Dortmund, Germany 21Institut de Physique Nucléaire, UCL, Louvain-la Neuve, Belgium 22Institute for Nuclear Research, Moscow, Russia 23Institut für Physik, Forschungszentrum Karlsruhe, Germany 24Los Alamos National Laboratory, Los Alamos, New Mexico, USA 25Sezione INFN, Rome, Italy 26Universitá Roma Tre, Rome, Italy 27Università “Federico II” e Sezione INFN, Napoli, Italy 28Institute for Nuclear Research and Nuclear Energy, Academy of Sciences, Sofia, Bulgaria


I. INTRODUCTION
From the HAdRon Production (HARP) experiment, final measurements are presented of the double-differential cross section, d2 σ /dp d , for π ± forward production (in the range of momentum 0.5 p 8.0 GeV/c and angle 0.025 θ 0.25 rad) by incident protons of 3, 5, 8, 8.9 (Be only), 12, and 12.9 (Al only) GeV/c momenta impinging on thin solid beryllium, carbon, aluminum, copper, tin, tantalum, or lead targets of 5% nuclear interaction length.Our results on the forward production of π + in p-Al interactions at 12.9 GeV/c [1] and p-Be interactions at 8.9 GeV/c [2], useful for understanding the K2K [3], MiniBoone, and SciBooNE neutrino beams, have been reanalyzed now with a consistent binning.In addition, results at 12 GeV/c on thin carbon targets and on cryogenic N 2 and O 2 targets, relevant for a precise calculation of atmospheric neutrino fluxes and improved modeling of extensive air showers, previously reported in Refs.[4,5], are included for completeness.
The HARP experiment [6] at the CERN Proton Synchrotron (PS) was designed to measure hadron yields from a large range of solid and cryogenic nuclear targets and for incident particle momenta from 1.5 to 15 GeV/c.This corresponds to a proton momentum region of great interest for neutrino beams and far from coverage by earlier dedicated hadroproduction experiments [7,8].The main motivations were to measure pion yields to allow substantially improved calculations of the atmospheric neutrino flux [9] to be made, to measure particle yields as input for the flux calculation of accelerator neutrino experiments [10], to help design the targetry for a future neutrino factory [11] and to improve the reliability of extensive air shower simulations, by reducing the uncertainties of hadronic interaction models in the low energy range [12].
By 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 dataset, taken with many different targets, full particle identification, and large detector acceptance.These data, together with those published in Ref. [13], are also § Now at CERN.MiniBooNE Collaboration.¶ Now at FNAL, Batavia, Illinois, USA.
‡ ‡ On leave from Novosibirsk University, Russia.
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 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 and 17.5 GeV/c with protons on Be, Cu, and Au targets [17].
In this study, data were taken in the T9 beam of the CERN PS.

II. EXPERIMENTAL APPARATUS
The HARP experiment uses a large-acceptance 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 powers of the TOFW below 3 GeV/c and the Cherenkov detector 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 greater 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 ranges 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 .The convention for the coordinate system is shown in the lower-right corner.The three most downstream (unlabeled) drift chamber modules are only partly equipped with electronics and are not used for tracking.The detector covers a total length of 13.5 m along the beam axis and has a maximum width of 6.5 m perpendicular to the beam.
A set of four multiwire 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 a signal in both beam Cherenkov detectors (BCA and BCB).At the lowest beam energy, 3 GeV/c, the BTOF is used to reject pions that 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, contaminations from all other particle species is negligible.Only events with a single reconstructed beam track in the four MWPCs, 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 resulted in a measured efficiency of >99.8%.Particle identification was 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 in which 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 of 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.[2,4,5,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-offlight measurement while maintaining a high reconstruction efficiency.The method to combine the information from these PID detectors is described in Ref. [23].In the kinematic range FIG. 4. Same as Fig. 3, but for p-C differential cross sections. of the current analysis, the pion identification efficiency is about 98%, while the background from misidentified 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,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 datasets.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 summarized in Table I.
FIG. 6. Same as Fig. 3, but for p-Cu differential cross sections.FIG. 7. Same as Fig. 3, but for p-Sn differential cross sections.

A. Cross-section calculation
The double-differential cross section is calculated as follows: where d 2 σ α dpd (p i , θ j ) 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.
N α (p i , θ j ) 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.
FIG. 8. Same as Fig. 3, but for p-Ta differential cross sections.

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 particle number to take into account effects such as reconstruction efficiency, acceptance, absorption, pion decay, tertiary production, PID efficiency and PID misidentification rate.
N pot , and A is the number of target nuclei per unit area; 2 N pot is the number of protons on target; and p i and j are the bin sizes in momentum and solid angle, respectively. 3 We do not make a correction for the attenuation of the proton beam in the target, so strictly speaking the cross sections are valid for λ I = 5% targets.
The calculation of the correction matrix M cor p i θ j αp i θ j α is done with the unfolding method introduced by D'Agostini 2 A-atomic mass, N A -Avogadro number, ρ-target density, and t-target thickness. 3 . [26]. 4 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 reinteractions.The simulation is independent of the beam particle type, because it only generates for 4 The unfolding method tries to put in correspondence the vector of measured observables (such as particle momentum, polar angle, and particle type) x meas with the vector of true values x true using a migration matrix: x meas = M migr × x true .The goal of the method is to compute a transformation (correction matrix) to obtain the expected values for x true 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 neighboring 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.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 multiparticle hadron-production model.A similar difference is seen between the single-particle 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 magnitudes 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 momentum.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 the Cherenkov threshold (∼3 GeV/c), the TOFW separates pions and protons with negligible background and an efficiency of ≈98% for pions.Above the Cherenkov threshold, the efficiency for pions is greater than 99% with only 1.5% of the protons misidentified as pions.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 reinteract 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 FIG.14. Same as Fig. 13, but for p-Be at 5 GeV/c.FIG. 15.Same as Fig. 13, but for p-Be at 8 GeV/c.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 Monte Carlo 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 aluminum and carbon targets at lower momenta (3 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 using the upstream drift chamber module (which FIG.16.Same as Fig. 13, but for p-Be at 12 GeV/c.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 day-by-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 incurred while obtaining the unfolding matrix from the data.The latter component increases the direct error by a factor of 2. The procedure is outlined in Refs.[4,28].FIG. 17. Same as Fig. 13, but for p-Ta at 3 GeV/c.Different types of sources induce systematic errors for the analysis described here: track yield corrections (∼5%), particle identification (∼0.1%), and momentum and angular reconstruction (∼1%). 5The strategy to calculate these systematic errors and the different methods used for their evaluation are described in Ref. [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 Ref. [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 FIG.18. Same as Fig. 13, but for p-Ta at 5 GeV/c.by a covariance matrix.The sum of these matrices describes the total systematic error, as explained in Ref. [4].
The experimental uncertainties are shown for a typical light target (Be) in Fig. 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 same order.

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 Figs.3-9 for solid targets from Be to Pb.The error bars 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. 6The results of this analysis are also fully tabulated and reported in Ref. [29].
An overall fit with a Sanford-Wang parametrization [30] has been done by using all solid targets and available beam momenta data (p 5 GeV/c), see Sec.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 π − and π + yields, averaged over two angular regions (0.05 θ < 0.15 and 0.15 θ < 0.25 rad) and four momentum regions (0.5 p < 1.5, 1.5 p < 2.5, 2.5 p < 3.5, and 3.5 p < 4.5 GeV/c), are shown 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 and 12 GeV/c with the p-Al, p-Cu, and p-Sn showing an intermediate behavior.
The dependence of the averaged pion yields on the atomic number A is instead shown in Fig. 11.The π − and π + yields, averaged over two angular regions (0.05 θ < 0.15 and 0.15 θ < 0.25 rad) and four momentum regions (0.5 p < 1.5, 1.5 p < 2.5, 2.5 p < 3.5, and 3.5 p < 4.5 GeV/c), are shown for four different beam momenta.One observes a smooth behavior of the averaged yields.The A dependence is slightly different for π − and π + production, the latter saturating earlier toward higher A, especially at lower beam momenta.

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 [30].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 σ/(dp d ) 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 target (typically Be) to another nuclear target A, a correction factor corr = (A/A Be ) α (4) was introduced, with α = α 0 + α 1 × x F + α 2 × x 2 F , where x F is the Feynman x, see Refs.[10,31] for details.A nineparameter (11-parameter) fit was done over 24 (25) π − (π + ) datasets, 7 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 = 13 030 (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 double-differential 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 what can be accommodated within the Sanford-Wang formula, where this ratio is given by g(θ ) = exp[−c 6 θ (p − p c )], with p c ≡ c 7 p beam cos c 8 θ .
The overall fit may be used as a fast approximation of HARP data valid within a factor of 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 11 (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, Fig. 12 reports the comparison, at 8 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 ad hoc Sanford-Wang parametrizations, using only the relevant HARP datasets, have been published in Refs.[1,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 Ref. [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 because of 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 initial 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 [32], 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-10 GeV/c, has to be done.
At intermediate energies (up to 5-10 GeV/c), GEANT4 uses two types of intranuclear cascade models: the Bertini model [33,34] (valid up to ∼10 GeV/c) and the binary model [35] (valid up to ∼3 GeV/c).Both models treat the target nucleus in detail, taking into account density variations and FIG.20.Same as Fig. 13, but for p-Ta at 12 GeV/c.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 Ref. [36], 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, two parton string models-the quarkgluon string (QGS) model [33,37] and the Fritiof (FTP) model [37]-are used instead, in addition to a high-energy parametrized (HEP) model derived from the high-energy part of the GHEISHA code used inside GEANT3 [38].The parametrized models of GEANT4 (HEP and LEP) are intended to be fast, but they 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."In high-energy calorimetry, the two most commonly used are the QGSP physics list-based on the QGS model, the pre-compound nucleus model, and some of the low-energy parametrized (LEP) model8 -and the LHEP physics list [39] based on the parametrized LEP and HEP models.The MARS code system [32] uses instead as its basic model an inclusive approach multiparticle production originated by R. Feynman.Above 3 GeV, phenomenological particle production models are used; below 5 GeV, a cascade-exciton model [40] combined with the Fermi breakup model, the coalescence model, an evaporation model, and a multifragmentation extension are used instead.
In Figs.13-20, data have been compared against Monte Carlo predictions, using beryllium and tantalum as examples of a light and a heavy target.
Computing the χ 2 between models and data themselves, where a systematic error going from 0% to 50% has 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 using 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 Ref. [41].

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 p π 8 GeV/c and 0.025 θ π 0.25 rad from the collisions of protons of 3, 5, 8, and 12 GeV/c on beryllium, aluminum, carbon, copper, tin, and tantalum targets of 5% λ I .In addition, results at 8.9 (Be only), 12.9 (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 of 2-3 of the quoted errors.
The pion yield averaged over different momentum and angular ranges increases 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 toward higher A, especially at lower beam momenta.
Comparisons with GEANT4 and MARS generators are presented.

FIG. 1 .
FIG.1.Schematic layout of the HARP detector.The convention for the coordinate system is shown in the lower-right corner.The three most downstream (unlabeled) drift chamber modules are only partly equipped with electronics and are not used for tracking.The detector covers a total length of 13.5 m along the beam axis and has a maximum width of 6.5 m perpendicular to the beam.

FIG. 3 .
FIG.3.p-Be differential cross sections for π − and π + 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. 12 .
FIG. 12.Comparison of HARP production data at 8 and 12 GeV/c with the Sanford-Wang global fit.

FIG. 13 .
FIG.13.Comparison of HARP double-differential π ± cross sections for p-Be at 3 GeV/c with GEANT4 and MARS MC predictions, using several generator models (see key in figure, and text for details).

TABLE I .
Total number of events and tracks used in the various nuclear 5% λ I solid target data sets and the number of accepted protons on target as calculated from the prescaled incident beam triggers.Numbers are for incident protons in units of 10 3 events.

TABLE II .
Summary of the systematic uncertainties affecting the computed π + double-differential cross sections and the integrated cross-section measurements for p-C interactions at 12 GeV/c.Entries are weighted bin by bin with the pion production yields.

TABLE IV .
Correlation coefficients among the Sanford-Wang fit parameters, obtained by fitting the data for π + production.

TABLE V .
Correlation coefficients among the Sanford-Wang fit parameters, obtained by fitting the data for π − production.

TABLE VII .
Normalization factors data simulation.