Feasibility studies for quarkonium production at a fixed-target experiment using the LHC proton and lead beams (AFTER@LHC)

Used in the fixed-target mode, the multi-TeV LHC proton and lead beams allow for studies of heavy-flavour hadroproduction with unprecedented precision at backward rapidities - far negative Feyman-x - using conventional detection techniques. At the nominal LHC energies, quarkonia can be studies in detail in p+p, p+d and p+A collisions at sqrt(s_NN) ~ 115 GeV as well as in Pb+p and Pb+A collisions at sqrt(s_NN) ~ 72 GeV with luminosities roughly equivalent to that of the collider mode, i.e. up to 20 fb-1 yr-1 in p+p and p+d collisions, up to 0.6 fb-1 yr-1 in p+A collisions and up to 10 nb-1 yr-1 in Pb+A collisions. In this paper, we assess the feasibility of such studies by performing fast simulations using the performance of a LHCb-like detector.


I. INTRODUCTION
Since its start-up, the large hadron collider (LHC) -the most energetic hadron collider ever built so far-has already made the demonstration of its outstanding capabilities. These can greatly be complemented by the addition of a fixed-target physics program. Its multi-TeV beams indeed allow one to study p + p, p + d and p + A collisions at a center-of-mass energy √ s NN 115 GeV as well as Pb + p and Pb + A collisions at √ s NN 72 GeV, with the high precision typical of the fixed-target mode. In this context the proposal of a fixed target experiment at the LHC [1], referred to AFTER@LHC -A Fixed Target Experiment -, has been promoted [1] in order to complement the existing collider experiments such as the Relativistic Heavy Ion collider (RHIC) or the future Electron-Ion Collider (EIC) project in a similar energy range. The idea underlying the AFTER@LHC proposal is a multi-purpose detector allowing for the study of a multitude of probes. Various technological ways to perform fixed-target experiment at the LHC exist. On the one hand, the beam can be extracted by means of a bent crystal. This technology [2,3] is currently developed as a smart beam-collimation solution and is studied by the UA9/LUA9 collaboration respectively at SPS and LHC. A bent crystal installed in the halo of the LHC beam would deflect the particles of the halo onto a target, with a flux of 5 ×10 8 proton/s without any impact on the LHC performances [3][4][5].
On the other hand, the LHC beam can go through an internal-gas-target system in an existing (or new) LHC experiment. Such a system is already tested at low gas pressure by the LHCb collaboration to monitor the luminosity of the beam [6][7][8]. Data were taken at a center-of-mass energy of √ s NN = 87 (54) GeV with p + Ne (Pb + Ne) collisions during pilot runs in 2012 and 2013. Although this system, called SMOG, was tested during only few hours in a row, no decrease of the LHC performances was observed.
In the bent-crystal case, the luminosity achievable with AF-TER@LHC would surpass that of RHIC by 3 order of magni-tudes [1]. We have reported in Table (I) the instantaneous and yearly integrated luminosities expected with the proton and Pb beams on various target species of various thicknesses, for the bent-crystal as well as internal-gas-target options. Integrated luminosities as large as 20 fb −1 can be delivered during a oneyear run of p + H collisions with a bent crystal. Besides, it is worth mentioning that both technologies allow one to polarise the target, which is an important requirement to lead an extensive spin physics programme [1,9].
Overall, thanks to the large luminosity expected, AF-TER@LHC would become a quarkonium, prompt photon and heavy-flavour observatory [1,10] in p + p and p + A collisions where, by instrumenting the target-rapidity region, gluon and heavy-quark distributions of the proton, the neutron and the nuclei can be accessed at large x and even at x larger than unity in the nuclear case [11]. In addition, the fixed-target mode allows for single-target-spin-asymmetry measurements over the full backward rapidity domain up to x F −1 [12,13]. In addition, the versatility in the target choices offer a unique opportunity to study the nuclear matter versus the hot and dense matter formed in heavy-ion collisions which can be studied during the one-month lead run. In the latter case, modern detection technology (such as high granularity calorimeter) should allow for extensive studies of quarkonium excited states, from the ψ(2S ) to the χ c and χ b resonances thanks to the boost of the fixed-target mode [14].
In this paper, we report on a feasibility study of quarkonium production at a fixed-target experiment using LHC beams. In section II, we outline the simulation framework which was used. In section III, we describe how a fast simulation of a detector response has been implemented, following a LHCb-like detector setup. In section IV, we present the charmonium and bottomonium family studies performed with the p + H simulations at √ s = 115 GeV. In section V, we present multiplicity studies in p + A and A + p collisions as well as the expected nuclear modification factors for J/ψ and Υ in p + Pb collisions at √ s NN = 115 GeV. Finally in section VI some prospects for

II. SIMULATION INPUTS
In order to get the most realistic minimum bias simulations at AFTER@LHC energy for quarkonium studies in the dimuon decay channels, we have simulated the quarkonium signal and all the background sources separately to have under control the transverse momentum and rapidity input distributions as well as the normalisation of the different sources.
The simulation has been performed for p + p collisions at √ s = 115 GeV. On the one hand, the quarkonium signal and the correlated background (Drell-Yan, cc, bb) were simulated with HELAC-Onia [15] which produces outputs following the format of Les Houches Event Files [16]. The outputs were then processed with Pythia (Pythia 8.185 [17]) to perform the hadronisation, the initial/final-state radiations and the decay of the resonances. On the other hand, the uncorrelated background was obtained from minimum bias p + p collisions generated with Pythia.
The relative normalisation of the signal and background sources was performed according to the production cross section of the process (taking into account initial phase space cuts, if any). Values of the cross section and the number of simulated events N sim -not to be confused with the expected events for a specific luminosity-are reported in Table (II). The cross section values are integrated over rapidity and p T .
The parameters κ, λ, n and p T were determined by fitting the differential cross section d 2 σ/d p T dy to the experimental data. The dedicated codes used to perform the fit and to generate unweighted events for quarkonium production have been implemented in HELAC-Onia [15] and we used MSTW2008NLO PDF set [22] provided in LHAPDF5 [23] and the factorisation scale µ F = M 2 Q + p 2 T . In order to constrain the non-trivial energy dependence of quarkonium production, we used the differential measurements of charmonium production performed by the PHENIX collaboration at RHIC, in p + p collisions at √ s = 200 GeV [20] to predict the corresponding yields at √ s = 115 GeV. Given the lack of such measurements for Υ at RHIC, we performed a combined fit to CDF [24], ATLAS [25], CMS [26] and LHCb [21,27] data on Υ production. The values of the fitted parameters are listed in Table (III). For illustration, the comparison between fits and the selected experimental data is shown in Fig. 1.   [20] for charmonium production and to CDF [24], ATLAS [25], CMS [26] and LHCb [21,27] data for bottomonium production. We have fixed n = 2 and p T = 4.5 (13.5) GeV/c for charmonium (bottomonium) production. The number of fitted data points is also reported.
In order to increase the statistics of the simulated data sample, the decay of the quarkonium in Pythia is forced into the dimuon decay channel. The simulated yields are then weighted by the cross section for this process multiplied by the Branching Ratio (BR).

Open charm
Open-charm production was simulated with the process gg → cc in HELAC-Onia. In order to avoid the huge theoretical uncertainties in the state-of-the-art perturbative calculations, open charm yields at √ s = 115 GeV are also computed in a data-driven way following the method described in the previous section. Similarly, the matrix element of gg → cc is determined using Eq. (1). The parameters are obtained from a fit to the p T -differential cc cross section measured by the STAR experiment [28] in p+ p collisions at √ s = 200 GeV (see Fig. 3). We obtained κ = 0.437, λ = 3.04 and p T = 2.86 GeV/c when n = 2 by using CTEQ6L1 [29] and by fixing the c quark mass to m c = 1.5 GeV/c 2 and the factorisation scale The χ 2 of the fit is equal to 4.39 with 10 experimental data points. The tuned result is shown in Fig. 2. The evolution of the cross section with the energy down to √ s = 115 GeV is then given by HELAC-Onia. A comparison between fit and the STAR data [28] in p + p collisions at √ s = 200 GeV, for cc production.
After embedding the Les Houches Event File into Pythia, muons from the underlying Pythia event can be produced on top of muons from the initial cc pair. The combination of those additional muons with a muon from the initial cc pair is not included in our definition of open charm correlated background. We have however checked that this contribution is negligible. In order to increase statistics, the D 0 ,D 0 , D +/− and D +/− s were forced to decay into muons and only those decay muons were considered as correlated background. s . This fraction is obtained from Pythia and found to be 95%.

Open beauty
The theoretical uncertainty on open beauty production is relatively smaller than the one on open charm production.
We therefore calculated open-beauty-production yields with a Leading Order (LO) matrix element and which was normalised to the Next-To-Leading-Order (NLO) K factor. The NLO cross section with the same setup was calculated by MadGraph5 aMC@NLO [30]. We used CTEQ6L1 (CTEQ6M) for the LO (NLO) calculation. The K factor is found to be 1.83. The renormalisation and factorisation scale is

Drell-Yan
Drell-Yan (DY) correlated background was simulated with the process qq → γ /Z → µ + µ − at LO where qq is a pair of same flavour light quarks. The LO calculation was done with the CTEQ6L1 pdf set and the renormalisation and factorisation scale was set to µ R = µ F = p T . In order to have enough statistics in the J/ψ and ψ(2S ) mass window, a phase space cut requesting that the invariant mass of the dimuons (M) is greater than 2.5 GeV/c 2 was applied. For the simulation of the DY background under the Υ family peaks, a phase space cut M > 7 GeV/c 2 was applied. The DY cross section obtained with HELAC-Onia at √ s = 38.8 GeV is compared to the existing E866 data at the same energy [31]. A K factor 1.2 is needed to match the data and therefore it was also applied at √ s = 115 GeV. Such a K factor is known to approximately account for the higher-order QCD corrections.

B. Uncorrelated background
The uncorrelated background was obtained from a minimum bias Pythia p + p simulation at √ s = 115 GeV using the process SoftQCD:nonDiffractive with the MRSTMCal.LHgrid LHAPDF (6.1.4) set [33]. By comparing our simulation of open charm with a low statistic pure minimum bias Pythia simulation, we have checked that the contribution of dimuons originating from a muon from charm/beauty and a muon from π/K is negligible. The dominant source of uncorrelated opposite-sign muon pairs is the simultaneous semi-muonic decay of uncorrelated π and/or K. In order to avoid possible double counting of signal and correlated background processes, the following hard processes have been switched off from the minimum bias simulations: HardQCD:hardccbar, HardQCD:hardbbbar, WeakSingleBoson:ffbar2gmZ 1 , Charmonium:all and Bottomonium:all, . 1 in order to avoid Drell-Yan pair production.

III. FAST SIMULATION OF THE RESPONSE OF A LHCb-LIKE DETECTOR
The HELAC-Onia and Pythia generators provide the opposite-sign muon pairs from quarkonia decays, correlated and uncorrelated backgrounds sources, as defined in the previous section. In order to account for the detector resolution and the particle identification capabilities of a given detector and to investigate the feasibility of the quarkonium studies in p + p collisions at √ s 115 GeV, the detector response needs to be simulated. For this purpose we have chosen a detector setup similar to the LHCb detector [34]. The forward detector is very well suited as a fixed-target experiment setup as well, with a good tracking and particle identification capabilities.
According to LHCb analysis cuts, muons in our simulations are required to have their transverse momentum satisfying p T > 0.7 GeV/c [35] and their pseudo-rapidity in the laboratory frame satisfying 2 < η < 5. The η cut range corresponds to the LHCb detector coverage. Since the momentum resolution reported by LHCb is δp/p ∼ 0.4 (0.6)% for a momentum of 3 (100) GeV/c [36] we consider a momentum resolution of δp/p = 0.5 %. The single µ identification efficiency is taken to be P = 98%, which is an average efficiency obtained by LHCb for muons coming from J/ψ decays, for p > 3 GeV/c and p T > 0.8 GeV/c [36]. These cuts and the abovementioned detector response on the muons are applied to simulate the quarkonium states and all the background sources.
In the case of uncorrelated background, as discussed in section II, most of the µ originate from π +/− or K +/− decays. If a π or K decays to a µ before 12 m along the z axis, the µ is rejected by the tracking system and it is not considered in the simulation. 12 m corresponds to the distance for uncorrelated muon background. The efficiency takes into account the identification efficiency of the prompt muons, and the π and K misidentification probability , P MID (π → µ) and P MID (K → µ).
The pair efficiency is extracted in each kinematic phasespace point and is shown as a function of the dimuon invariant mass, transverse momentum and rapidity in Fig. 5. This efficiency is used to correct dimuon spectra obtained with the uncorrelated background Pythia simulations.
In this section, we show results on the quarkonium production studies in the dimuon decay channels, with the dominant background sources. Simulations have been performed for a 7 TeV proton beam on a hydrogen target (p + p), which gives √ s = 115 GeV. We consider an integrated luminosity of 10 fb −1 which is expected to be obtained after half of a LHC year with the crystal mode, as described in section I and Table (I). A. Background studies These simulations allow to quantify the background sources in the quarkonium studies in the dimuon decay channel, which could potentially make the quarkonium signal extraction more difficult or even prevent from obtaining a clear signal. In particular, this may be critical for the excited states. We present here simulations of invariant mass of opposite-sign muon pairs, µ + µ − , from the quarkonia and from the domi- In the J/ψ and ψ(2S ) invariant mass window, the dominant background source is from uncorrelated µ + µ − pairs, mostly from π +/− and K +/− decays. Contributions from Drell-Yan and bb continuum are very small. In the case of Υ(nS ) states the Drell-Yan contribution is the dominant one. Under the Υ(nS ) peak, a contribution from the cc continuum is negligible, and is not considered here.
The significance (sig = S / √ (S + B), where S is the number of signal counts and B the number of background counts, in the invariant mass range M Q ± 3 σ Q ) and the signal to background ratio (S /B) of each quarkonium state are given in the following: • sig J/ψ = 134. 6  Transverse momentum and rapidity distributions for the quarkonium signals and for each background source were also studied. As an example, p T and y distributions in the J/ψ mass range, 3.063 < M µ + µ − < 3.129 GeV/c 2 (corresponding to M J/ψ ± 3 σ J/ψ ) are shown in Fig. 7. It is visible that the distributions for the J/ψ and different backgrounds differ. In more backward or forward rapidity regions, the signal to background ratio increases. This can also be seen in Fig. 8, where the dimuon invariant mass distributions in J/ψ and ψ(2S ) mass  window are shown in three rapidity ranges. In terms of transverse momentum, one can obtain a very clean signal when going to higher p T . Above ∼ 4 GeV/c, the uncorrelated background starts to vanish. Since cc, bb and Drell-Yan simulations are LO simulations, the p T spectra of these correlated background sources are not shown here.

B. Quarkonium simulations
We have also studied the p T and rapidity coverage reach of the quarkonium signals. The transverse momentum distributions are shown in Fig. 9 (a), for J/ψ, ψ(2S ), Υ(1S ), Υ(2S ) and Υ(3S ), from the top to the bottom distribution. Similarly, Fig. 9 (b) shows the rapidity distribution for each quarkonium state. With an integrated luminosity of 10 fb −1 the quarkonium studies can be carried out in a wide rapidity and p T range. It should be possible to study Υ(nS ) signals up to p T 10 GeV/c, and J/ψ and ψ(2S ) could be studied even up to p T 15 GeV/c. All the quarkonium states can be measured down to p T = 0 GeV/c. This study is limited by the rapidity range of 2 < y < 5, in the laboratory frame, due to the pseudorapidity cuts on the decay µ. The red x-axis on the top of Fig. 9 (b) denotes the rapidity in the center-of-mass frame. The rapidity shift for a 7 TeV proton beam on a fixed target is -4.8, i.e. y CM = 0 → y lab = 4.8. J/ψ and ψ(2S ) signals can be studied in the whole mentioned rapidity range, while the lowest rapidity reach for Υ(nS ) is ∼ 2.5 -3. In proton-nucleus collisions, the high track multiplicity may induce a high detector occupancy and lead to a reduction of the detector capabilities. Since LHCb has successfully measured the J/ψ and Υ production in p+Pb collisions at √ s NN = 5 TeV [38,39], one would expect a good capability of such detector under similar particle multiplicity environment. In the following, the charged particle multiplicity has been generated with the EPOS generator [40,41] in different configurations: p + Pb collisions at √ s NN = 5 TeV in collider mode, p + Pb collisions at √ s NN = 115 GeV and Pb + H collisions at √ s NN = 72 GeV in fixed-target mode. The chargedparticle multiplicity is dominated by the π multiplicity. By comparing these three distributions as a function of the pseudorapidity of the particle in the laboratory frame as shown in Fig. 10, one can conclude that the charged particle multiplicity in a fixed-target mode never exceeds the one obtained in p + Pb collisions at √ s NN = 5 TeV in the collider mode in the full pseudorapidity range: a detector with the LHCb capabilities will be able to run in such conditions. To illustrate the potential offered by AFTER@LHC in p + Pb collisions at √ s NN = 115 GeV, we have evaluated, in this section, the impact of the nuclear modification of the gluon densities in nucleons within large nucleus -generically referred to as gluon shadowing-and its uncertainty as encoded in the nuclear PDF set EPS09. For that, we have used the probabilistic Glauber Monte-Carlo framework, JIN [42,43], which allows us to encode different mechanisms for the partonic production and to interface these production processes with different cold nuclear matter effects, such as the aforementioned shadowing, in order to get the production cross sections for proton-nucleus and nucleus-nucleus collisions. JIN also straightforwardly computes any nuclear modification factor for minimum bias collisions or in specific centrality classes. In the case of proton-nucleus (p + A) collisions, it is the ratio of the yield per inelastic collision in p + A collisions to the yield in pp collisions at the same energy multiplied by the average number of binary collisions in a typical p + p collision, N coll : In the presence of a net nuclear effect, R pA is defined such that it differs from unity. In the simplest case of minimum bias collisions, one should have As in [44], we have used the central curve of EPS09 as well as four specific extreme curves (minimal/maximal shadowing, minimal/maximal EMC effect), which reproduce the envelope of the gluon nPDF uncertainty encoded in EPS09 LO [45].
In addition to the modification of the partonic densities, quarkonium production in p + A collisions can be affected by other effects, for instance by the nuclear absorption which depends much on the nature of the object traversing the nuclear medium. If the meson is already formed, it may be affected more than a smaller pre-resonant pair. To discuss such an effect, it is useful to introduce the concept of the formation time, t f , based on the Heisenberg uncertainty principle and the time -in the rest frame of the meson-to discriminate between two S states, for instance the J/ψ and the ψ(2S ). In fact, one finds [44,46] that such a time is similar for the charmonium and bottomonium states and is on the order of 0.3-0.4 fm. Obviously, this time has to be boosted in the frame where the nuclear matter sits. For t f smaller than the nucleus radius, the quarkonium is formed before escaping it. In the fixed-target mode with a proton beam and a nuclear target, the boost factor is simply γ(y lab ) = cosh(y lab ). We therefore obtain t f as in Table (IV). One sees that looking at quarkonium production in p + Pb collisions at different backward rapidities allows one to look at quarkonia traversing the nuclear matter at very different stages of their evolution. This effect could theoretically be studied by giving an ad-hoc rapidity dependence to the effective absorption cross section, σ effective abs . This is left for a future study since, here, we wish to consider only the nPDF effects and the expected statistics. Other effects to be considered are the coherent energy loss [47] (expected to grow in the forward region) and the rescattering by comovers [48] (expected to grow with the multiplicity along the J/ψ direction).
Since we wish to assess the descriminating power of possible data to be taken with AFTER@LHC, we attribute to the EPS09 central values statistical uncertainties which directly follow from the differential yields repectively expected in p+p and p+Pb collisions. For that we take an integrated luminosity of 10 fb −1 for the p + p runs and 100 pb −1 for the p + Pb runs, in accordance with the luminosities discussed above (see Table (I)). As this stage, we do not consider additional systematical uncertainties. This simplifying assumption could be lifted in a more detail study which would also take into account a   possible detector acceptance (and related efficiencies) as done in the previous section. In particular, we do not expect that the rapidity region for y CMS > 1.5 would be easily accessible.
In Fig. 11, we show the rapidity dependence of R p+Pb for Υ and its p T dependence near y = 0. The million of Υ to be collected per year allows for the measurement of a R p+Pb with a much better precision than the gluon nPDF, nearly up to x → 1. In addition, one notes that the nuclear modification factor is certainly measurable up to p T 10 GeV/c.
In Fig. 12, we also show the rapidity dependence of R p+Pb for J/ψ and its p T dependence near y = 0. In both cases, the luminosity to be taken in a year at AFTER@LHC yields to statistical uncertainties which are largely negligible as compared to the nPDF uncertainties -the statistical uncertainties are not even visible on Fig. 12. We except this to hold also for the ψ(2S ) although its yields are down by a factor of 100.
As aforementioned, the nPDFs do not account for all the expected nuclear matter effects. However, it is clear that combining the measurements of Υ, J/ψ and ψ(2S ) for −3 < y CMS < 0 (as a LHCb-like detector would do) will allow one to pin down the existence of a possible gluon EMC and antishadowing effect. We also stress that the complications induced by a rapidity dependence of σ effective abs could be avoided by the parallel measurement of R p+Pb for non-prompt J/ψ which can only be sensitive to the energy loss since the b quark decay (weakly) into the J/ψ, way outside the nucleus. Fig. 13 shows that the trend is similar than for Υ. Measuring the p T dependence of R p+Pb for prompt J/ψ and Υ should also avoid the sensitivity on formation time effects.

VI. PROSPECTS OF Pb+A MEASUREMENTS AT
√ s = 72 GEV The charged particle multiplicity has been generated with the EPOS generator [40,41] in different configurations: Pb + Pb at √ s NN = 5.5 TeV in collider mode, Pb + Ar, Pb + Xe and Pb + Pb at √ s NN = 72 GeV in fixed-target mode. By comparing these three distributions in the pseudorapidity of the particle in the laboratory frame as shown in Fig. 14, one can conclude that the charged particle multiplicity in a fixed-target mode never exceeds the one obtained in Pb + Pb collisions at √ s NN = 5.5 TeV obtained in a collider mode in the full pseudorapidity range: a detector with the ALICE MFT+Muon detector [49] capability will be able to run in such conditions. Details studies are needed to evaluate up to which multiplicity a detector such as LHCb would be able to take good quality data.

VII. CONCLUSION
In summary, we have shown that in a fixed target mode with an integrated luminosity of 10 fb −1 , using 7 TeV LHC proton beam on a hydrogen target, and with a detector setup and performances similar to the LHCb detector, quarkonium studies in the dimuon decay channel can be performed over a wide transverse momentum range and rapidity in the center of mass from ∼ −2.8 for J/ψ and ψ(2S ), and ∼ −2 for Υ states, up to ∼ 0. We have performed simulations of the dominant background sources contributing to the µ + µ − invariant mass spectrum. The uncorrelated background was obtained using Pythia generator and dimuons from correlated background sources: cc, bb and Drell-Yan, were simulated using both HELAC-Onia and Pythia generators. The estimated background level allows for J/ψ, ψ(2S ), Υ(1S ), Υ(2S ) and Υ(3S ) measurements in the dimuon decay channel with good signal to background ratios.
These simulations set the stage for further ones including, on the one hand, the detection of photon from P wave or η c decay or from the production of a J/ψ + γ pair, whose studies at low transverse momentum can provide important insight on the gluon transverse dynamics [50][51][52][53] and, on the other hand, the large combinatorial background typical of p + A and A + A collisions in which the study of excited quarkonium at AF-TER@LHC energies is of paramount importance [1,10]. We note that the Delphes [54] framework seems particularly well suited to account for the photon detectability in such prospective studies.
Along our investigations, we have also noted that main source of dimuons around the Υ(nS ) masses is from the Drell-Yan process (see Fig. 6 (right)). This gives us great confidence that the corresponding cross section can easily be extracted in this mass region in p + p collisions, a fortiori with a vertex detector allowing for tagging the heavy-flavour muons. We therefore consider that the single-spin asymetries for Drell-Yan pair production can indeed be extracted using a light polarised target. Motivations for such studies are discussed in [12,13,55]. Quarkonium polarisation measurement are of course also possible given the large statistical samples.
As regards the case of p + A collisions, we have had a first look at the charged particle multiplicities as a function of the laboratory pseudo-rapidity. We have found out that, for all the possible fixed target modes, p+Pb, Pb+H, Pb+Pb, these are smaller than the ones reached in the collider modes where the LHCb was used (p+Pb and Pb+p at 5 TeV). We therefore believe that a detector with similar characteristics as compared to LHCb can very well be used in the fixed-target mode 2 .
In view of the above, we have evaluated the impact -and its uncertainty-on the nuclear modification of the gluon densities on prompt and non-prompt J/ψ and Υ in form of R p+Pb . We have found that the measurements at backward rapidities allows one to search for the gluon antishadowing, the gluon EMC effect and even the Fermi motion effect on the gluons with unheard of statistical precisions. The statistics are large enough to perform such measurement with the ψ(2S ) and probably also with Υ(2S ) and Υ(3S ) allowing for thorough investigations of formation time effect of the meson propagating in the nuclear matter. Overall, our results confirm the great potential of AFTER@LHC for heavy-quark and quarkonium physics.