Kinematic distributions and nuclear effects of $J/\psi$ production in 920 GeV fixed-target proton-nucleus collisions

Measurements of the kinematic distributions of $J/\psi$ mesons produced in $p-$C, $p-$Ti and $p-$W collisions at $\sqrt{s}=41.6 \mathrm{GeV}$ in the Feynman-$x$ region $-0.34<x_{F}<0.14$ and for transverse momentum up to $p_T = 5.4 \mathrm{GeV}/c$ are presented. The $x_F$ and $p_T$ dependencies of the nuclear suppression parameter, $\alpha$, are also given. The results are based on $2.4 \cdot 10^{5}$ $J/\psi$ mesons in both the $e^+ e^-$ and $\mu^{+}\mu^{-}$ decay channels. The data have been collected by the HERA-B experiment at the HERA proton ring of the DESY laboratory. The measurement explores the negative region of $x_{F}$ for the first time. The average value of $\alpha$ in the measured $x_{F}$ region is $0.981 \pm 0.015$. The data suggest that the strong nuclear suppression of $J/\psi$ production previously observed at high $x_F$ turns into an enhancement at negative $x_F$.


Introduction
The DESY experiment HERA-B has measured inclusive J/ψ production in proton-carbon, protontitanium and proton-tungsten collisions at a centerof-mass energy √ s = 41. 6 GeV. The results are based on a sample of about 2.4 · 10 5 J/ψ mesons reconstructed in both dilepton decay channels. A measurement of the atomic mass number dependence of J/ψ production is derived from a comparison of the different samples. The atomic number dependence of inclusive particle production is often characterized by a power law: σ pA = σ pN · A α where σ pN is the proton-nucleon cross section and σ pA is corresponding proton-nucleus cross section for a target of atomic mass number A. Previous measurements by E866 at Fermilab [1,2] at √ s = 38.8 GeV and NA50 at CERN [3] at lower energy indicate that α ∼ 0.94 − 0.95 at x F ∼ 0 and decreases to ∼ 0.65 as x F approaches unity [2]. The results presented here provide a first measurement of nuclear effects in charmonium production extending into the negative part of the Feynman-x spectrum, −0.34 < x F < 0.14 1 .
An understanding of the basic mechanisms responsible for the suppression of charmonium production in proton-nucleus collisions relative to proton-nucleon collisions is a prerequisite for the identification of possible signals of new physics in high-energy heavy-ion data. Interpretations of the existing proton-nucleus data at positive x F rely on a delicate balance of several processes: nuclear absorption, shadowing of parton densities, energy loss, interactions with comovers, hadronization of intrinsic cc components of the scattering nucleons, and so on. Ad hoc combinations of such elementary mechanisms, considered within various theoretical frameworks [4,5,6] are able to qualitatively reproduce the observed strong increase of J/ψ suppression as x F approaches unity [2].
However the presently available data give little guidance in the largely unexplored negative-x F region where other mechanisms such as formation-time effects can influence the effective nuclear path length of produced states [7,8,9,10] and can potentially lead to a decidedly different behavior of α. In contrast to J/ψ production in the positive x F region, at negative x F the produced cc pair preferentially evolves into a charmonium state before leaving the nucleus and nuclear effects influencing the J/ψ itself become important. Especially in this region, different models and approaches lead to contrasting predictions, for example arising from differing assump-tions on energy loss of beam partons or produced cc pairs. By extending the current experimental knowledge of the nuclear dependence of J/ψ production towards negative x F , the measurement described here provides new constraints for possible explanations of the observed nuclear modification pattern. The wide range of transverse momenta (up to 5.4 GeV/c) of the present measurement also permits a complementary measurement of the p T -broadening effect obtained at lower energies.
The paper is divided into four main sections: an overview of the apparatus, trigger and data samples (Sect. 2), a description of the methods used for the selection and counting of J/ψs (Sect. 3), the measurements of the kinematic distributions (Sect. 4) and the measurements of the nuclear dependence (Sect. 5), followed by concluding remarks (Sect. 6).

Apparatus, data taking and
Monte Carlo simulation HERA-B was a fixed-target experiment which studied particles produced in interactions of 920 GeV/c protons with wire targets positioned in the halo of the proton beam. The apparatus [11], shown in Fig. 1, was a forward spectrometer with acceptance ranging from 15 to 220 mrad and from 15 to 160 mrad in the bending (xz) and vertical (yz) planes, respectively. Because of this large acceptance and the fact that the J/ψ decay tracks were measured before the muon detector (MUON) and the electromagnetic calorimeter (ECAL), HERA-B was the first fixed-target experiment with significant coverage in the region of negative x F with an accessible range of x F ∈ [−0.34, 0.14].
The target system [12] consisted of eight wires which were grouped into two stations separated by 4 cm along the beam line. Each wire could be individually steered in the beam halo by a servo system in order to maintain a constant interaction rate. A total of five wires, differing in shape (round or rectangular), dimensions (between 50 µm and 500 µm) and material (C, Ti and W) were used. Depending on running conditions and run-plan, either a single wire or a pair of wires was active for any given data-taking run. The interaction rate was maintained in the range of 2 to 6 MHz, depending on beam conditions and target.
The vertex detector (VDS) [13] comprised eight planar stations of double-sided silicon micro-strip modules, seven of which were mounted in Roman pots built into the vacuum vessel and operated at a minimum distance of 10 mm from the beam. A track traversed typically three stations, yielding twelve mea- surement points in four stereo views. Vertex resolutions of 450 and 50 µm in the beam direction and in the transverse plane, respectively, were achieved. The eighth station was on a fixed mount immediately following the exit window of the main vertex system, 2 m downstream from the target.
The momenta of charged tracks were measured from their bending through a vertical magnetic field of integral 2.13 Tm. The main tracker was located between 2 to 13 m downstream of the target, with one station preceding the magnet, four stations immediately after the magnet (Pattern Chambers, PCs) and another two (Trigger Chambers, TCs) after the Ring Imaging Cherenkov detector (RICH). Each station contained an inner part [14] (made of micro-strip gas chambers and covering angles less than 20 mrad) which was not included in the trigger system and therefore does not play a role in the analysis presented here. The outer part (OTR) [15] was composed of honeycomb drift chambers, with wire pitches of 5 mm closer to the beam pipe and 10 mm elsewhere. The final momentum resolution for muons was found to be σ p /p[%] = (1.61 ± 0.02) + (0.0051 ± 0.0006) · p [GeV/c] [15].
The identification of the J/ψ in its dilepton decay modes as well as the first stage of the trigger system relied mainly on the signals provided by the ECAL [17] and MUON [18] systems. The ECAL was a sampling calorimeter using shashlik technology with Pb and W absorbers sandwiched between scintillator layers. It was divided into three sections (Inner, Middle and Outer) with cell widths of 2.2, 5.5 and 11 cm, respectively, to roughly equalize occupancies. The inner section used W absorbers, while the middle and outer sections used Pb. The design was optimized for good electron/photon energy resolution and for electron-hadron discrimination. The final energy resolution reached by the detector can be written in the form , with A = 0.206, 0.118 and 0.108 and B = 0.012, 0.014 and 0.014, for the Inner, Middle and Outer sections, respectively. The spatial resolution ranged from 1 to 10 mm depending on the calorimeter section and on the energy of the particle [17].
The MUON system consisted of four tracking stations interleaved with iron or concrete absorbers. As in the main tracker, two different technologies were used: gas pixel chambers in the innermost region and conventional tube chambers in the outer part. For the last two stations, not only the anode wires of the tubes but also the segmented cathodes were read out. Signals from the cathode pads were also given as inputs to the trigger.
The RICH [16] detector relied on a C 4 F 10 radiator and was used extensively in other analyses for π/K/p separation. It also played a small role in the present analysis as a means to reject backgrounds in the dilepton analysis, particularly from kaon decays.
The trigger system selected both e + e − and µ + µ − signatures and was organized into three levels: a pretrigger [17,19], a First Level Trigger (FLT [20]) and a software-based Second Level Trigger (SLT [21]). The pretrigger used signals from the ECAL and MUON detectors and required the presence of at least two reconstructed ECAL clusters with transverse energy above 1.1 GeV or the presence of two muon candidates, defined as coincidences of projective pads in the last two MUON layers. Starting from pretrigger seeds, the FLT attempted to find tracks in a subset of the OTR tracking layers and required that at least one of the seeds resulted in a reconstructed track. Starting again from the pretrigger seeds, the SLT searched for tracks inside regions-of-interest generated by the pretriggers using all OTR layers and continued the tracking through the VDS. Finally, at least two fully reconstructed tracks consistent with the hypothesis of a common vertex were required. Events passing the SLT were transferred to a computer farm which provided full online reconstruction of a fraction of the events for data quality monitoring. The global trigger suppression factor, 5 × 10 4 , resulted in an event archival rate of about 100 Hz.
A total of 160 million dilepton triggered events were recorded between October 2002 and February 2003, together with an approximately fixed 10 Hz rate of minimum bias events which were used for monitoring and luminosity determination. The event samples were distributed between three target materials: carbon (64%), tungsten (32%) and a small fraction with titanium (4%).
A full Monte Carlo (MC) simulation is used to determine the triggering (except for FLT, see below), reconstruction and selection efficiencies. In view of the range of physics topics addressed by the experiment (pA inelastic interactions, meson decays and heavy flavor physics), the MC generator is built as a combination of two standard tools: Pythia 5.7 [22] for heavy flavor (b or c) quark production in pN interactions and subsequent hadronization and Fritiof 7.02 [23] for light quark production, secondary interactions in detector materials and generic pA inelastic interactions. The production of the J/ψ is simulated by first generating the basic hard process pN → ccX and subsequent cc hadronization with Pythia and then giving the remaining energy and momentum (X) of the interaction to Fritiof for generation of further interactions inside the hit nucleus. The generated particles are then given as input to the Geant 3.21 based package [24] for full simulation of active (instrumented) and inactive (support structure) elements of the detector and for the digitization of the electronic signals.
To describe the kinematic characteristics of the produced J/ψs as accurately as possible, an x F , p T , and decay-angle dependent weight is assigned to each event and used in the subsequent analysis to force the simulated J/ψ production and decay distributions to agree with the corresponding measured distributions for each target material. The weights are determined by an iterative procedure in which computed corrections are based on comparisons of MC event distributions after reconstruction and selection cuts with the corresponding distributions from data.
The FLT efficiency is derived from an efficiency map which is determined from the data itself. Since the SLT result is completely independent of the FLT, and since the FLT triggered on only one of the two lepton tracks, the efficiency map can be determined by a so-called tag-and-probe method. Using the efficiency map, each event is assigned a weight which multiplies the kinematic weight discussed above.
To accurately reproduce the actual working configuration at the time of data-taking and properly account for time variations of working conditions, the full data taking period is divided into five calibration periods of similar lengths, each matched by a corresponding simulation sample for which the efficiencies of the individual detector cells are evaluated and given as input to the MC. The MC samples are reconstructed and analyzed with the same methods and software packages used for the analysis of the real data.

J/ψ selection and counting
The electron and muon candidates are selected with common track-and vertex-selection criteria while channel-specific methods are applied for lepton identification and the treatment of the J/ψ signal.
All tracks passing the SLT are initially considered as lepton candidates. The track reconstruction procedure consists of finding straight segments in the VDS and PCs independently, matching them to each other and also to segments in the TCs. A full, iterative fit of found tracks is then performed. To reject incorrectly reconstructed or ghost tracks, loose cuts on the minimum number of hits in the VDS and OTR, as well as on the χ 2 probability of the track fit are applied by the event reconstruction algorithms. The detector acceptance and trigger requirements effectively limit the momentum and transverse momentum ranges of lepton tracks to 5 < p < 200 GeV/c and 0.7 < p T < 5.0 GeV/c, respectively. For each event, among all possible pairs of oppositely charged lepton candidates consistent with a common vertex (χ 2 probability greater than 1%), only the pair with the best particle identification (see below) for both leptons is accepted.

Dimuon channel
Since muons are the only particles having a significant probability of penetrating through the absorbers of the MUON detector, only minimal selection cuts are needed to obtain a clean sample (see also [25]). The background of muons from pion and kaon decays is reduced by imposing tighter cuts on the quality of the track fit and on the matching of track segments in the VDS, OTR and MUON. Doing so rejects the typical "broken trajectories" produced by decays into leptons. Contamination from kaons is further reduced by discarding tracks with high values of the corresponding RICH likelihood. After all selection cuts, the background under the J/ψ signal is reduced by a factor of 2.5 with respect to the triggered data with a loss of about 11% of the signal. Fig. 2(a) shows the resulting dimuon mass spectrum together with the result of a fit to a sum of three functions [26], which model the J/ψ and ψ(2S) signals and the exponential background. The J/ψ and ψ(2S) signals are each modeled as a superposition of three Gaussians with a common mean which takes into account track resolution and effects of Molière scatter-Channel C Ti W Total µ + µ − 94800 8060 48100 152000 e + e − 57700 4280 25300 87200 Table 1: The numbers of reconstructed J/ψs after all selection cuts in the dimuon and dielectron channels, and for different target materials.
ing and a function representing the radiative tail due to the decay J/ψ → µ + µ − γ [26]. The background is described by an exponential of a second-order polynomial. The fitted position and width of the J/ψ peak are 3.0930 ± 0.0002 GeV/c 2 and 40 ± 1 MeV/c 2 , respectively.

Dielectron channel
The J/ψ selection in the dielectron channel is affected by major background contributions from charged hadrons which produce energetic ECAL clusters and by overlapping photon and charged-hadron energy deposits in the ECAL. For this reason, the electron identification requirements were the subject of careful optimization studies which resulted in much more stringent selection cuts than for the muon sample.
A cut on the transverse energy of the ECAL cluster (E T > 1.15 GeV) is applied in order to mask different threshold cuts applied at the pretrigger level for the various acquisition periods.
The reconstructed momentum vectors of electrons and positrons are corrected for energy loss from bremsstrahlung (BR) emission in the materials in front of the magnet. For each electron track, an attempt to identify an ECAL cluster due to a BR photon is made by looking for an energy deposition in coincidence with the extrapolation of the associated VDS track segment to the ECAL. Any recovered energy (about 18% of the initial electron energy on average) is then added to the momentum measured by the tracking system. Since BR is a clear signature for an electron track, it is also exploited to obtain substantial background reduction which is essential for accurate J/ψ counting. The baseline results of the analysis are obtained by requiring that at least one lepton of a decaying J/ψ has an associated BR cluster. This requirement reduces the signal by about 30% and suppresses the background by more than a factor of two. Alternative requirements (no BR requirement, only one, or both electrons emitting BR) lead to very different background shapes and amounts. The differences are exploited for systematic studies on the stability of signal counting and on the correctness of Adjustments to the measured momenta of electron tracks were applied to compensate for differences in multiple scattering between electrons and muons since the track-fitter had been calibrated for muons. For this purpose, a correction map, determined from the shift of the J/ψ peak position in different kinematic regions was used.
Additional selection cuts are applied to further improve the significance of the dielectron signal. A particularly discriminating variable is the ratio E/p, where E is the energy of an ECAL cluster and p is the momentum of the associated track. The E/p distribution for electrons has a Gaussian shape with a mean value close to 1 and width varying between 6.4% and 7.4% depending on calorimeter sector. Values of E/p much lower than 1 are mainly due to particles, mostly hadrons, which release only part of their energy in the calorimeter. Further selection variables used in the analysis are the distances ∆x and ∆yalong the x and y directions -between the reconstructed cluster and the track position extrapolated to the ECAL. The ∆x and ∆y distributions for electrons are, apart from a small tail, well described by Gaussians centered at zero with widths between 0.2 and 1.0 cm depending on calorimeter sector. Cuts on these quantities lead to a significant reduction of the contamination from hadrons and random clustertrack matches which are characterized by significantly wider distributions. The selection of the candidate electron-positron pairs is further refined by putting an upper bound on the distance of closest approach (∆b) between the two accepted tracks near the vertex.
All the requirements described above have been simultaneously optimized by maximizing the significance S/ √ S + B of the J/ψ signal (S) -taken from the MC (scaled to the number of J/ψ in the data) -with respect to the background (B) -evaluated from the data. The optimal ranges for the different cut variables depend on the number (one or two) of BR clusters associated to the electron-positron pair. When both electron and positron have a BR cluster correlated to the track, where the cluster position was determined by hierarchical clustering [17], the event is already rather cleanly reconstructed and only one additional request (for each lepton candidate), (E/p − 1)/σ E/p > −3.6, is applied. When only one of the two possible BR clusters is found, the accepted ranges are determined for each lepton as −3.6 < (E/p − 1)/σ E/p < 5.4, |∆x|/σ ∆x < 7.0, |∆y|/σ ∆y < 3.3 and ∆b < 370 µm, respectively.
The combined selection cuts increase the S/B ratio of the J/ψ by about a factor of 10 with respect to triggered events, and have an overall efficiency of (45 ± 4)% as evaluated using the data and verified with the simulation -the rather large uncertainty is due to the difficulty of counting the signal when no cuts are applied. As can be seen in Fig. 2(b), the significance of the optimized J/ψ signal, although less than that of the muon channel shown in Fig. 2(a), is nonetheless such that the electron sample significantly enhances the statistical precision of the final results. The method adopted for counting the signal uses a Gaussian shape for the right part of the peak and a Breit-Wigner shape for the left part to account for the sizable asymmetry of the signal caused by missing BR energy and the contribution of the radiative decay J/ψ → e + e − γ. The background is parametrized with a Gaussian at lower mass values and an exponential at higher mass, joined together such that the resulting curve is smooth. The position and width of the J/ψ peak as determined from the fit are 3.110 ± 0.001 GeV/c 2 and 72 ± 1 MeV/c 2 , respectively.
The yields of selected J/ψ candidates for the two decay channels and for each target material are listed in Table 1. Fig. 3 shows (for the carbon data) a comparison between data and MC of the distributions of reconstructed J/ψs as a function of the kinematic variables p T and x F .

Kinematic distributions 4.1 Results
The present analysis adopts the degrees of freedom p T , x F and Φ (azimuthal production angle) for the description of the J/ψ production kinematics. Singlevariable distributions are obtained according to the formula (here written e.g. for where ∆N rec J/ψ (x F ) is the fraction of J/ψs reconstructed in a given x F interval of width ∆x F and ǫ J/ψ (x F ) is the global (trigger, reconstruction and selection) efficiency for that interval integrated over all other kinematic variables (including the J/ψ decay degrees of freedom) using the tuned MC. In each case, the signal is also integrated over all other kinematic variables. All distributions are normalized to unit area 2 .
The final p T and x F distributions for the three  Table 2: J/ψ p T distributions (dN/dp T , normalized to their integrals over the measured range) for three target materials with statistical and systematic uncertainties. • The impact of selection and optimization requirements is evaluated by changing the cuts on momentum and transverse momentum of muon and electron candidates with respect to the intrinsic thresholds of the trigger selection, and by scanning systematically the values of all cut variables used for the optimization of the dimuon and dielectron signals (including, for the latter, different BR requirements).
• The uncertainty associated to the signal counting method has been estimated from the variation of the results obtained with the adoption of modified background and signal functions. Special attention is given to the background evaluation of the dielectron channel: as an alternative to the fit with an assumed background function, a background shape constructed by mixing real events (combining each track with one of opposite-charge from a different event) has been used in an unbinned maximum-likelihood fit of the invariant mass spectrum. A further crosscheck is represented by the comparison between the efficiency-corrected J/ψ yield obtained with different BR requirements (and therefore different background shapes). When the BR requirement is removed, the evaluation of the number of J/ψs becomes less stable due to increased  Table 3: J/ψ x F distributions (dN/dx F , normalized to their integrals over the measured range) for the three target materials with statistical and systematic uncertainties.
background, but the variation of the efficiencycorrected yield with respect to the standard selection is estimated to be lower than 5%. • The systematic uncertainties also account for the stability of the results when specific acquisition periods and conditions are selected. A large subsample of the collected events was produced on two target wires of different materials operated simultaneously. The comparison of these data with those acquired with a single target provides an indication of the extent to which the measurements are affected by variations of the experimental conditions. • The results are sensitive to the shape of the J/ψ decay angular distribution assumed in the MC generator. The hypothesis -made by previous experiments -that the J/ψ is produced in an unpolarized state is not supported by the HERA-B data [28]. There is, moreover, an indication that the polarization increases in magnitude with decreasing p T , while no significant x F dependence is found. Since a longitudinally polarized J/ψ is detected more efficiently due to the lower probability that its decay leptons escape detection by passing through the uninstrumented region near the beam, the kinematic dependence of the polarization assumed in the MC influences the shape of the efficiency-corrected p T and/or x F spectra. The systematic stability tests therefore include a variety of different assumptions for polarization (including longitudinal, p T -dependent polarization -also considering the possibility of an A-dependent polarization -and the absence of polarization).
The distribution of the azimuthal production angle Φ has been evaluated as a systematic check of the uniformity of the MC description of the geometrical acceptance. Fig. 6 shows the result obtained when combining the full data reconstructed in both decay channels: the points are consistent, within the statistical uncertainties, with the expected flat distribution.
The dimuon and dielectron results (which, as discussed below, are found to be compatible) are averaged such that correlations in their systematic uncertainties are taken into account. Such correlations have been estimated by maintaining (when possible) a parallelism among the two channels when evaluating the effect of each single systematic test. The assumed Channel Param. ∆x F -0.0024 ± 0.0016 ± 0.0022 -0.0052 ± 0.0051 ± 0.0015 -0.0096 ± 0.0024 ± 0.0012 γ 1.723 ± 0.036 ± 0.011 1.48 ± 0.14 ± 0.01 1.820 ± 0.063 ± 0.018 Table 4: Parameter values obtained from the fit of the kinematic distributions of each of the three target samples to the functions described in the text (Eq.s 2 and 3). The first of the given uncertainty ranges is statistical and the second is systematic. polarization hypothesis is found to be the dominating source of uncertainty in the final results -especially for the lower part of the p T distribution -as well as the most important cause of correlation between the two analyses. Relative to this uncertainty, the signal selection cuts are in general responsible for negligible systematic variations, except for the most positive part of the x F spectrum, where acceptance corrections increase dramatically due to the low-angle detector acceptance cut-off near the beam.
The results are well represented -and can therefore be described -by the following interpolating functions which are further motivated in the following The results of the present analysis (black filled circles with total and statistical uncertainties) are compared to previous measurements performed with different beam energies at Fermilab [1] and at the SPS [3]. The data are fitted with linear functions. section: The parameters p 2 T , β, w xF (width at half maximum), ∆x F (shift of the center of the distribution with respect to x F = 0) and γ are left free in the fit of the distributions. The resulting values are listed in Table 4 for each target material. The systematic uncertainties of the parameters have been determined from the maximum variation (divided by √ 12) of the results obtained by re-fitting the distributions after each of the stability tests described above. The normalized χ 2 obtained from the fits of the combined dimuon-dielectron results are, respectively, 1.5, 0.8 and 2.1 for the p T distributions of carbon, titanium and tungsten, and 1.1, 1.5 and 1.2 for the three x F distributions (only the statistical uncertainties are taken into account). Table 4 are the results obtained separately in the dimuon and dielectron channels with the respective systematic uncertainties. The good agreement between the results of the two analyses within the statistical uncertainties confirms that the and at the SPS [3]. The data are fitted with linear functions.

Also included in
channel-specific issues of particle identification and counting are not responsible for large systematic variations in the shapes of the distributions.

Parameterization and interpretation of kinematic distributions
Among the parameters adopted for the description of the data, the width of the p T distribution ( p 2 T ), the position of the maximum of the x F distribution (∆x F ) and, possibly but less significantly, its width (w xF ) show a trend with the mass number A.
It is well known [29] that the average p 2 T of particles produced in nuclear collisions increases with the mass of the target nucleus. This trend is confirmed by the HERA-B data with high significance. The "p Tbroadening" effect is commonly explained as a consequence of multiple elastic scattering of the incoming beam parton in the surrounding nucleus before the hard scattering process takes place. The measured increase of p 2 T with A is shown in Fig. 7 together with the results of experiments at lower energies. The variable of the abscissa, A 1/3 − 1, is approximately proportional to the radius of the target nucleus i.e. to the average path length of the parton inside the nucleus with the shift of −1 such that the magnitude of the effect is measured with respect to A = 1. All measurements are actually consistent with the parameterization to which they are fitted in the plot. As summarized in Fig. 8(a), the results for the average p 2 T extrapolated to proton-nucleon interactions ( p 2 T pp ) are compatible with a linear growth with the square of the centerof-mass production energy s. On the other hand, ρ is approximately independent of center-of-mass energy as can be seen in Fig. 8(b).
Furthermore, HERA-B observes a difference in shape between the x F distributions of the J/ψ for different target nuclei consisting of an increasing displacement of the center of the distribution towards negative values. As shown in Fig. 9, J/ψs are produced in tungsten with an x F distribution which 2.030 ± 0.014 ± 0.014 ρ (GeV 2 /c 2 ) 0.0852 ± 0.0053 ± 0.0043 β pp 6.97 ± 0.22 ± 0.23 β ′ 0.261 ± 0.087 ± 0.070 χ 2 /N DoF 106.5/71 x F distribution, τ free κ -0.00211 ± 0.00044 ± 0.00022 w pp xF 0.1435 ± 0.0017 ± 0.0046 τ 0.00277 ± 0.00066 ± 0.00092 γ pp 1.735 ± 0.029 ± 0.014 1.725 ± 0.0028 ± 0.0015 χ 2 /N DoF 91.5/61 Table 5: Results of the global fits of p T and x F distributions described in the text. The first of the given uncertainty ranges is statistical and the second is systematic.
has equal or slightly greater width with respect to those produced in carbon and tends to be asymmetrically centered at a lower value. This behavior is also supported by a fit of the E789 gold data [1] (−0.035 < x F < 0.135, E b = 800 GeV) with Eq. 3. The fitted width of 0.11 ± 0.01 is lower than our value for tungsten suggesting not only that the maximum is shifted but that the shape becomes asymmetric. As a possible interpretation, the effect may be attributed to the energy loss undergone by the incident parton and/or the produced state in their path through the nucleus, causing a reduction of the average x F of the J/ψ and a possible additional smearing of the momentum distribution. This hypothesis motivates the choice of representing the data also in this case as a function of A 1/3 − 1. The points in Fig. 9 (a) and b) are fitted respectively with: where κ, w pp xF and τ are free parameters.
To obtain the best description of the dependence of the p T and x F spectra on the target nucleus, a simultaneous fit of the three (C, Ti and W) distributions according to the functions given in Eq. 2 and 3 has been done. According to the hypothesis that energy loss is responsible for the observed nuclear dependence of the shape of the kinematic distributions, the fit has been constrained by imposing the relations from Eqs. 4, 5 and 6, and, moreover, where β pp , β ′ and γ pp are additional parameters of the fit. In Table 5 the results of this procedure are summarized. The fit of the x F distributions has been performed in two variants, with the parameter τ left free or fixed to zero -therefore assuming in the latter case that w xF is independent of A. The resulting best-fit curves (with τ left free) are the interpolating lines plotted in Figs. 4 and 5. The fit results indicate a significant nuclear dependence not only of the p T distribution (parameter ρ), but also of the x F distribution: there is a significance of 7σ for κ = 0 when τ is fixed to zero, which changes to 4 and 3σ, respectively, for κ and τ (with a strong anti-correlation between the two) when both are left free.

Nuclear dependence of J/ψ production
The Glauber Model [30] suggests that the dependence of the J/ψ production cross section on atomic mass number (A) can be approximated by a power law: where σ pN is the proton-nucleon cross section and α, the "suppression" parameter, characterizes the nuclear dependence. Pure hard scattering in the absence of any nuclear effects would correspond to α equal to unity. A suppression of J/ψ production would lead to α < 1 while an enhancement (anti-screening effect) would be signaled by α > 1. Usually, α is described and measured as a function of x F and p T (see for example [2,3,31]).
Eq. 9 is generally used to describe data and predictions independently of particular mechanisms of nuclear modification. In general, however, α may depend on A and thus depend on the targets used to make the measurement. For the measurement presented here, α is evaluated by comparing the J/ψ yields from two different targets: carbon and tungsten 3 as a function of x F and p T .

The α measurement
Using Eq. 9, the nuclear suppression parameter α can be extracted from a measurement of the ratio of crosssections for carbon (C) and tungsten (W) targets: The measurement of the cross section ratio requires a measurement of the ratio of the integrated luminosities of the carbon and tungsten target samples. For the HERA-B setup, this can be done using data samples where two different targets are operated simultaneously (double-target runs) since most of the systematic uncertainties cancel and an absolute luminosity measurement can be avoided. On the other hand, for studies of the dependence of α on the kinematic variables, greater statistical precision can be obtained by also using the single-target runs. The HERA-B measurement of α thus consists of two submeasurements: a measurement of the average value of α, α , over the full visible kinematic range and a measurement of α − α as a function of x F and p T . The shape distributions are then corrected using α to produce an absolute measurement of the distribution of α over the measured range.
More specifically, α is evaluated using double target runs based on the formula: where N X (X=C,W) denotes the total number of reconstructed J/ψ mesons originating from the corresponding target wire, L X is the luminosity, ǫ X is the overall detection efficiency (see Sec. 2) and the event yields are derived as discussed in Sec. 3. The measurement of the luminosity ratios is described in Sec. 5.2.
The dependence of α on x F and p T is obtained from the full carbon and tungsten target data samples (double-and single-target runs). The full carbon sample is roughly twice the size of the double-target subsample while the full tungsten sample is 10% larger than the double-target subsample. The shape of the differential distributions are given by (here written e.g. for x F ): where σ J/ψ = σ(pA → J/ψ + X) is the total visible J/ψ cross section.
The measurement of nuclear effects can then be derived from the distributions of Eq. 12 using (here written e. g. for x F ):

Luminosity ratios
The luminosity acquired on target X, where X is either C (carbon) or W (tungsten), can be expressed as: where N X is the total number of inelastic interactions occurring on target X during the measurement, σ inel X is the total inelastic cross section, N BX is the corresponding total number of filled bunch crossings (BX) and λ X is the average number of interactions per filled BX. The total cross sections σ inel X for each target material together with some details on this topic can be found in [32]. The luminosity ratio R L needed in Eq. 11 is then given by: Assuming the interaction probability on target X follows a Poisson distribution, λ X can be calculated from the observed number of events with at least one interaction (N (≥ 1) obs ) using: where ǫ inel X is the probability to observe a single interaction.
The determination of λ X relies on random-trigger events which were accumulated together with the dilepton-trigger events used for J/ψ counting. Five methods which differ by the event characteristics used to define the presence of an interaction are used to count events. All methods rely on tracks found in the vertex detector. To maintain high efficiency, the requirements imposed are minimal but sufficient to also keep the probability of incorrect target wire assignment at a low level. The methods are based on the following five criteria: For all events, 1. ≥ 1 primary vertex on wire X where the primary vertex is formed from tracks measured both in the VDS and OTR ("long tracks"), 2. ≥ 2 tracks (including long tracks and tracks seen only in the VDS) with impact parameter ≤ 3σ of wire X and ≥ 5σ from the other wire, where σ is the impact parameter measurement uncertainty, and, using the subset of events with no vertex found on the other wire, 3. ≥ 1 primary vertex on wire X using all tracks, 4. ≥ 1 primary vertex on wire X using long tracks only, 5. ≥ 2 long tracks within ≤ 3σ from wire X.
Both as a cross check and as an estimate of the systematic uncertainty on the efficiency, all counting methods are checked in parallel. For the final luminosity ratio determination, the average of the five determinations is used and the rms spread of the five is factored into the systematic uncertainty. Furthermore, Eq. 16 assumes that the interaction probabilities for each wire follow a Poisson distribution. However, the individual bunch fillings are often uneven and the interaction rate varies in time by typically 20%. To quantify this influence, an alternative luminosity calculation is performed in which the detailed bunch filling structure and the interaction rate distribution on each wire are taken into account. The differences between the resulting ratios and those computed directly from Eq. 15 are negligible compared to other systematic uncertainties.
To minimize the dependence of the efficiency estimate on MC, the efficiency of each of the above methods is calibrated by comparing the luminosity estimate found using it with that found by the methods described in [32]. These latter methods rely on very simple criteria to identify events with interactions, such as a minimal number of hits in the RICH detector (typically twenty, compared to thirty hits expected for a fully accepted fast charged particle) or a small (1 GeV) energy deposit in the ECAL and are estimated to be sensitive to roughly 95% of the total non-diffractive cross section.
Each target of a two-wire configuration is calibrated separately using single-wire data runs taken nearby in time to the run being calibrated. A "ghost" wire is introduced at the location of the other wire of the configuration. Thus the MC is not relied on to model tracking in the VDS or vertex finding, but only to estimate the efficiency of simpler and more robust event counting techniques described in [32]. Overall efficiencies in the range of 60 -80 % are found, depending on method and wire. The efficiencycalibration method based on real data is also used to evaluate the probability that interactions are as-signed to the wrong wire. This probability is method and configuration dependent and is never more than 0.4%.
The average relative systematic uncertainties on the luminosity ratios due to interaction counting and MC calibration [32] are 1.3% and 3.2%, respectively, giving an overall scale uncertainty of 3.4% on R L . Depending on wire configuration between 0.6 and 1.2 million events were used for the determination of λ C and λ W , thus the statistical uncertainty on the luminosity ratio is negligible.

Results
Based on Eq. 11, an average suppression value of α = 0.981 ± 0.004 stat. ± 0.016 sys. (17) in the visible range of x F is obtained. As explained in Sect. 2, the target system of HERA-B consisted of eight different wires grouped in two stations. The data were taken with four different two-wire configurations which were analyzed separately and averaged to obtain the final value. Also, the electron and the muon decay channels represent two statistically independent measurements. The average value is calculated as a weighted mean with weights being the squared quadratic sum of statistical and luminosity ratio uncertainties. The luminosity ratios are determined for each wire configuration separately as discussed in Sec. 5.2 with a contribution to the systematic uncertainty on α of 3.4 %/ ln(A W /A C ) = 1.24 %. A systematic effect of 1.1% due to time variatons of detector performance and imprecisions in detector or trigger simulations was estimated from the variations of α among the four samples and two decay channels. The total systematic uncertainty on α is thus 1.66 %. The statistical precision of the α measurement contributes an uncorrelated uncertainty 0.4 %.
The values of α for individual x F and p T bins are given in Table 6 and shown in Fig. 10 for p T and in Fig. 11 for x F . They will be further discussed in Sec. 5.4. The values and systematic uncertainty estimates were derived according to the procedure described in Sec. 4.1. The error bars on the figures show both statistical and total contributions. The systematic uncertainties in the estimate of α are largely uncorrelated with those from the α − α measurement. The final systematic uncertainty estimate is the quadratic sum of the two. The systematic uncertainty is substantially correlated from bin to bin.

Discussion
The results of the α measurement as functions of p T and x F are given in Table 6. The distributions of the α(p T ) and α(x F ) are presented in Figs. 10 and 11 where they are compared with measurements performed by other fixed target experiments: E866 [2] (E p = 800 GeV) and NA50 [3] (E p = 450 GeV). As already seen (e.g. in Fig. 7), the measured p T dependence of the nuclear modification effects is very similar for HERA-B and E866. In Fig. 11 the E866 and HERA-B measurements are seen to be compatible within statistical and systematic uncertainties in the overlap region. The NA50 results are based on lower energy collisions and are systematically below both HERA-B and E866. At lower values of x F , the HERA-B α(x F ) measurement indicates a reversal of the suppression trend seen at high x F : the strong suppression established by previous measurements at high x F turns into a slight tendency towards enhancement in the negative x F region.
The dependence of J/ψ production in hadronnucleus interactions on x F has been modeled by Vogt [4,5]. Nuclear effects caused by final-state absorption, interactions with co-movers, shadowing of parton distributions, energy loss and intrinsic charm quark components are described separately and integrated into the model. It is further assumed that the cc pair is subject to more severe energy losses if produced in a color octet state. Four curves from this model which differ in their descriptions of nuclear Parton Density Functions (nPDF) and energy loss are shown in Fig. 11. All calculations shown here were done for the center-of-mass energy of HERA-B at √ s = 41.6 GeV. The nPDF distributions of Eskola,  [3] (empty triangles). The curves were calculated by Vogt [4,5] based on three different nuclear parton distribution functions: (EPS [37], EKS [33,39] and HKN [38]) and two models of initial state energy-loss: GM [7] and BH [8]. For all approaches, energy loss, intrinsic charm and shadowing are taken into account.
Kolhinen and Salgado (EKS) [33] describe the scale dependence of the ratios of nPDFs of a proton inside a nucleus to those of a free proton within the framework of lowest order leading-twist DGLAP evolution [34] by evolving the initial PDF from the CTEQ4L [35] and leading order GRV [36] parameterizations. An improved leading-order DGLAP analysis of nPDFs including next to leading order calculations has been published recently by Eskola, Paukkunen and Salgado (EPS) [37]. In another approach by Hirai, Kumano and Nagai (HKN) [38], nuclear structure function ratios F A 2 /F A ′ 2 and Drell-Yan cross section ratios are analyzed to obtain nPDFs. The HKN analysis shows weak anti-shadowing at negative x F . Initial state energy loss as described by Gavin and Milana (GM) [7] and modified by Brodsky and Hoyer (BH) [8] is based on a multiple scattering approach that essentially depletes the projectile parton momentum fraction as the parton moves through the nucleus. Both quarks and gluons can scatter elastically and therefore lose energy prior to the hard process resulting in an effective reduction of J/ψ production for The measurement of HERA-B shows that α increases with decreasing x F and suggests enhanced J/ψ production for x F < −0.1. The HERA-B data favors the nPDFs of EPS and HKN over EKS. The BH description of energy loss is clearly ruled out. None of the variants of the Vogt model give a satisfactory description of both HERA-B and E866 data. For example, while the HKN curve is compatible with most of the HERA-B data points at negative x F , it lies significantly above the E866 points and furthermore fails to adequately describe RHIC data [37].
mesons, and therefore has a large interaction cross section. As it propagates through the nucleus it interacts and loses energy. Ultimately, the observed J/ψ mesons are projected out of the energy-depleted colorless state. The measurements of HERA-B are qualitatively compatible with the calculations described in BK [6].

Conclusions
HERA-B has performed the first determination of the nuclear dependence of J/ψ production kinematics at negative x F in proton-nucleus collisions. The analyzed data samples were obtained in collisions of protons from the 920 GeV HERA-proton beam with carbon, tungsten and titanium targets. The J/ψ mesons are observed in both dimuon and dielectron decay channels. The comparison of results from the two channels affords some additional control over systematic uncertainties arising from triggering and selection procedures.
The measurement covers the kinematic range −0.34 < x F < 0.14 and p T < 5.4 GeV/c. The measured dN/dp T distribution is seen to become broader with increasing atomic mass number as has already been observed by experiments at lower center-of-mass energies [1,3]. The data indicates that the dN/dx F distribution also tends to become broader and that its center moves towards negative x F values with increasing A.
The dependences of the nuclear suppression parameter, α, on p T and x F are also presented. The α parameter is seen to increase with increasing p T in agreement with data from E866 [2]. In the x F region of overlap of the two experiments, the two α measurements are mutually consistent. As x F decreases, α increases and becomes greater than 1 below x F ≈ −0.15, although the data remains compatible with a value of 1 to within 2 σ. Thus instead of the strong suppression observed at high positive x F , HERA-B measured no suppression or a possible enhancement of J/ψ production at negative x F . Hard-scattering based models [4,5] have difficulty simultaneously accommodating the HERA-B, E866 and RHIC measurements, while the reggeon-inspired model of Boreskov and Kaidalov [6] is in qualitative agreement with the data.