Responses to PLB Referees

 Responses to Reviewers' comments:

We would like to thank the referees for the insightful and constructive comments. We discuss below our detailed replies to your questions and the corresponding explanation of changes to the manuscript. But before we go into the replies to the comments, we want to make the referees aware of changes to the results that were prompted via our studies of the systematic effects on the yield extraction.  Since this paper deals with cross sections and with nuclear-modification factors, both of which involve obtaining the yields of the Upsilon states, this change affects all the results in the paper. We therefore wanted to discuss this change first. Please note that the magnitude of these effects do not change the overall message of the paper.

We wanted to alert the referees up-front about this important change before we proceeded into the detailed responses.  This study was indirectly prompted by one of the questions from Referee 2 regarding systematic effects from yield extraction.

In the process of investigating the systematic difference between extracting the upsilon yield through simultaneous fitting compared to background subtraction as requested by the referee, we also studied the effects of chi^2 fits (specifically of Modified Least-Squares fits) compared to maximum-likelihood fits. We used chi^2 fits in our original submission. We were aware that extracting yields using a chi^2 fit introduces a bias (e.g see Glen Cowan's "Statistical Data Analysis", Sec. 7.4). The size of the bias is proportional to the value of the chi^2 of the fit.  In the case of the Modified Least-squares fit, when fitting a histogram including the total yield as a fit parameter, the yield will on average be lower than the true yield by an amount equal to chi^2.  The relative bias, i.e. the size of the bias divided by the extracted yield, goes to zero in the large yield limit, which is why for cases with large statistics this effect can be negligible.  We had attempted to mitigate the effects of this bias by using the integral of the data, since this removes the bias completely in the signal-only case.  But a bias remains in the case where there are both signal and background present. For our case, the yield extracted from the fit for the background is also biased toward lower values, and since we used this background estimate to subtract from the integral of the data in the extraction of the Upsilon yields, these biased the Upsilon yields towards higher values.   Through simulation studies, where we include signal and 3 background components as in our analysis, we were able to quantify these effects. Given that in some cases the biases could be of order 10-20%, the fits needed to be redone in order to remove the bias. The solution is straightforward since the extraction of yields using a maximum-likelihood fit is unbiased.  We have studied the difference of a modified-least squares fit and a maximum-likelihood fit and confirmed that the yield extraction in the latter method is essentially unbiased. We therefore have redone all the fits to extract the Upsilon yields via maximum-likelihood fits. The revised results are now quoted in the paper. The overall message of the paper is not affected by these changes.

We proceed next to answer the specific points raised by the reviewers.

Reviewer #1: This paper reports results on Y production in pp, dAu, and AuAu
collisions at top RHIC energy. It contains original and important
results and clearly qualifies for publication in PLB. However, there
are many aspects of the paper which need attention and/or improvement
prior to publication. They are detailed below:

1. the introduction is carelessly written. For example, the value
quoted for the pseudo-critical temperature near mu = 0 of 173 MeV is
taken from an old publication in 2003. Recent lattice results from the
Wuppertal-Budapest group (PoS LATTICE2013 (2013) 155) and the Hot QCD
Collaboration (Phys.Rev. D81 (2010) 054504) imply much lower T values
near 150 MeV and are far superior in terms of lattice sizes and spacing.

There are certainly newer results, which we now cite in the paper. However, we note that the
the results from the Hot QCD collaboration (Phys.Rev. D81 (2010) 054504) do not imply
much lower T values.  In that paper, in section IV "Deconfinement and Chiral aspects of the QCD transition", when discussing the deconfinement transition temperature range the authors write:
"...we have seen that the energy density shows a rapid rise in the temperature interval T = 170–200. MeV. This is usually interpreted to be due to deconfinement, i.e., liberation of many new degrees of freedom". 
Therefore, this does not indicate T values near 150 MeV. In addition, they also mention this range when discussing their results for the renormalized Polyakov loop, which
is the parameter most closely related to the deconfiment transition, being that it is the exact order parameter in the pure
gauge case: 
"The renormalized Polyakov loop rises in the temperature interval T = 170–200 MeV where we also see the rapid increase of the energy density."
Therefore, the results from the Hot QCD collaboration do not imply T values near 150 MeV.

In addition, in reference 9 of the Wuppertal-Budapest group (JHEP 1009 (2010) 073 arXv:1005.3508), which is a paper comparing the various results for Tc between the Wuppertal-Budapest and HotQCD groups, again the results for the renormalized Polyakov loop (figure 7, right) indicate a broad transition region in the region T=160-200 MeV.  They do have a table discussing values of Tc of about 147 MeV, but that is for the chiral transition, which is not the most relevant one for quarkonium suppression.   
When they look at the trace anomaly (e-3p)/T^4, they see 154 MeV for the Tc value.  They in addtion make the point that the transition is a broad crossover, which is something we also say in our paper.  The fact that the transition is a broad crossover leads to differences in the estimates of the pseudo-critical temperatures depending on which observable is used.  As an example, in the caption of Table 2, where they give the values of Tc for many observables, they mention that the Bielefeld-Brookhaven-Columbia-Riken Collaboration obtained Tc=192. They also note "It is more informative to look at the complete T dependence of observables, than
just at the definition-dependent characteristic poins of them." So given the above, we will modify the paper to give a range of temperatures, 150-190, and cite the papers from the
Wuppertal-Budapest and HotQCD collaborations.  

also the discussion on whether charmonium or bottomonium 'is a cleaner
probe..' does not get to any of the real issues, such as the complex
feeding pattern in the Y sector and the crucial question of whether Y
mesons reach equilibrium in the hot fireball as required to interpret
the apparent sequential melting pattern in terms of 'break-up'
temperatures. 
 

 
The issues we discuss, in our opinion, are real issues.  We discuss co-mover absorption and the interplay between suppression and statistical recombination of uncorrelated charm pairs. These have been a topic of intense interest in the charmonium case for over a decade.  We certainly agree that these are not the only issues, but in this paper we aim to present the result of our measurement, so we cannot give a detailed review of all issues. However, the aim was to point out that for the bottomonium case, the expected contributions to either suppression or enhancement from both of these mechanisms are much smaller, and hence studying Upsilons is cleaner. The reviewer brings up the importaint issue of feed-down that affects the bottomonium as it does the charmonium sector.  We have added a few sentences regarding feed-down.  Regarding bottomonium, the feed-down contributions to the Upsilon states are not measured at RHIC energies yet. It is therefore assumed in the models used by Strickland, Rapp, etc. that the fraction of directly-produced Upsilon(1S) is ~51%, as measured in pp collisions at high pT at Tevatron energy. The original paper motivating the quarkonium sequential suppression by Digal, Petreczky, and Satz discussed feed-down as part of the impetus for looking for suppression of the Upsilon(1S). Given the ~51% direct Upsilon fraction, an R_AA of the Upsilon(1S) as low as ~0.51 would not necessarily imply suppression of the direct 1S, but could be due solely to suppression of the excited states. We have added text about this point in the paper, in discussing the R_AA(1S) result.

The reviewer also mentions that there is a crucial question as to whether the Upsilon mesons reach equilibrium with the fireball as a requirement to interpret the sequential melting pattern. We respectfully disagree with the referee in this matter. The Upsilon is by definition not in equilibrium. The only requirement of course is that the medium is deconfined. In lattice QCD studies only the medium is thermalized; the potential between the heavy quarks is screened independent of whether the Upsilon is in equilibrium or not. We discussed this issue with lattice expert Peter Petreczky who confirmed our view.

 

furthermore, statistical recombination is not a 'complication' but a
direct measure of deconfinement. And the smallness of off-diagonal
terms in the recombination matrix does not imply absence of
recombination as the diagonal terms can be substantial.

We agree that statistical recombination is an indication of deconfinement, but from the experimental side, it has made the interpretation of the results more complicated, because one needs to take into account the interplay of suppression and recombination.  Because this effect is negligible for the bottomonium states even at LHC energies, the quantitative interpretation of the experimental results is less complicated. It is in this sense that the word "complication" is meant.
 

Also the newest results on p-Pb collisions from the LHC are entirely
ignored, see, e.g., arXiv:1308.6726.

We are aware of the quarkonium pPb results from LHC, however there is not a way to make a direct comparison to LHC results, because there are no pPb results on the nuclear modification factor of Upsilons.  The results from ALICE in the reference given above are for the J/psi meson, and they are also for the forward-backward kinematic range.  There are also results from CMS (arXiv:1312.6300) for Upsilon mesons at midrapidity, but these are in the form of ratios of the yield of the excited states to the ground states in a given system (pp, pPb, PbPb), and of double ratios, i.e. excited-to-ground-state ratios in pPb divided by excited-to-ground-state ratios in pp.  These give us relative suppression of the excited states, whereas our results are for absolute suppression, and are therefore not directly comparable. The only quantitative comparison to the CMS data we can make is to estimate a double ratio for the excited states. The double ratio we find is consistent with the result from pPb from CMS, but it is also consistent with 1, i.e. no suppression of the excited states relative to the ground state in pPb compared to pp (We find the double ratio = 0.72 +/- 0.37). We will make a comment about this in the paper, and cite the CMS pPb result.  But the advantage of the results we are presenting is that we have fully corrected nuclear modification factors, which convey more information than relative suppression. 
 

2. section on experimental methods

no detail is given concerning the crucial momentum resolution but it
is stated at the end of this section that cuts were adjusted for
different systems such that 'tracking and electron id would be the
same across the 3 data sets'. On the other hand, already in Fig. 1 we
see a strong dependence of the mass resolution on the system even for
low multiplicities as in pp and p-Pb. The effect must be much stronger
in Pb-Pb as is indeed visible in Fig. 4. Especially in view of the
importance of resolution for the separation of excited Y states this
referee has to be convinced that the systematic errors are under
control for momentum and pid measurements as a function of
multiplicity. Also how the systematic errors for the separation of Y'
and Y'' from Y are determined as a function of multiplicity needs to
be demonstrated explicitely.

We agree that the mass resolution is very important for the results of the paper, and need to be discussed.  We added text to discuss the Upsilon mass resolution, and how it was studied as a function of TPC multiplicity (we focus on mass resolution, but this is directly related to the momentum resolution of the electron tracks used to reconstruct the Upsilon).  We studied the mass resolution using both simulations and data-driven methods.  Regarding the momentum resolution and the difference of the mass width seen in the pp compared to the dAu and AuAu plots, the majority of the difference between the pp lineshape and the dAu/AuAu lineshapes comes from a miscalibration in the TPC which was corrected in the pp dataset via a reproduction of the raw data, but due to time constraints was not corrected in the dAu and AuAu datasets. With the distortion correction, the pp mass resolution is found to be 1.3%.  If there were no distortions in Au+Au, we find in peripheral events a mass resolution for the Upsilon(1S) of 1.7%, which widens to 2.0% for central events, based on simulations. In order to ensure that we had this resolution effect under quantitative control, in addition to studying it via embedding simulated tracks into real collision events, we also studied the difference in the reconstructed mass of every reconstructed di-electron pair between the corrected and uncorrected pp datasets on an event-by-event basis. We were able to determine the additional mass smearing introduced by this TPC distortion, and this data-driven knowledge was used in the determination of the line shapes in dAu and in Au+Au.  The additional smearing introduced by the TPC distortion resulted in a mass resolution for the Upsilon(1S) of 2.7%, widening to 2.9% for central events.  For d+Au, the mass resolution is also found to be 2.7%, consistent with the peripheral events in AuAu. Finally, given the importance of the resolution for the separation of the states, we also used one additional data-driven method to check the resolution. We used the data from Au+Au and performed a chi^2 scan varying the mass-width parameter of the Upsilon line shape to see if this would give the same results as those found from the pp event-by-event mass-smearing data-driven study.  The results were consistent with each other, giving us confidence that the mass resolution is under control.  We used the shape of the chi^2 vs. resolution-parameter derived from the data to assign an uncertainty to our mass-resolution parameter knowledge, and used this to estimate a systematic uncertainty on the yields.
We have added a few sentences to the text towards the end of the Experimental  Methods section to give the relevant mass resolution numbers, and to make it clear that the mass resolution is different for pp compared to both dAu and AuAu due to the TPC distortion.  We also added some sentences to discuss this in the description of the invariant mass figures.
The systematic uncertainties due to our knowledge of the line-shapes, which are directly related to the separation between the Upsilon(1S) and the excited states, are also listed in the systematic uncertainty table that we added to the paper. The rows listing the uncertainty due to the line shape, including these mass resolution effects and the uncertainty in the knowledge of these TPC distortions, are given in the table.  Finally, we also added a systematic uncertainty to the extraction of the Upsilon(1S) yield based on the purity of our mass cut.  This is affected not only by the knowledge of the line shapes, as discussed above, but also by the possible suppression of the excited states.  We estimated this uncertainty by comparing the case of no suppression to that of complete suppression of the excited states, and recalculated the Upsilon(1S) purity for each.  The systematic uncertainty table has a row giveing the value of this uncertainty on the Upsilon(1S) yield, and we also added text to clarify it in the paper.


3.  Fig. 2b  

even at y = 0 the difference between data and models is less than 2
sigma, taking uncertainties due to the pp reference into account and I
don't believe it makes sense to argue about effects beyond shadowing
and initial state parton energy loss in these data. 

With the new fit results, only the point at y=0 is different from the models, and while the difference is now of order ~3sigma, the two other points do not show any deviation from the models. The text has been changed removing the sentence mentioning effects beyond shadowing.  Also, thanks to this comment, we realized that we had plotted the full pp cross section systematic uncertainty on the figure, but for the purposes of R_dAu, some of these systematic uncertainties cancel.  Therefore, the band illustrating the systematics due to the pp reference should have been smaller, and this has now been corrected.
 

4. in Fig. 3 the size of the systematic errors should be indicated.

The problem is that none of the data points from E772 had systematic uncertainties, so we cannot include them.  We have indicated the size of our systematics in the plots. 
 

5.  in Fig. 5 it is demonstrated that the observed suppression near
midrapidity is independent of system size (N_part). This could imply
that the higher Y states are already disappeared in dAu
collisions. This is mentioned briefly in the conclusion, but could be
stressed more. 

The statement we made in the paper is complementary to the one suggested by the referee. We state at the end of the paper that we cannot claim that the suppression in AuAu is unambiguously from color deconfinement in AuAu given the suppression in dAu.  The reason we stated it like this is that the expectation was for only a minimal amount of suppression in dAu, but our results are a call for caution and for considering other hypotheses.  The referee's comment is a call to consider a specific hypothesis: that the higher Upsilon states already dissappear in dAu collisions. This is a hypothesis that is not accounted for yet in any model. While we are not advocating any particular hypothesis for dAu suppression, we can add a sentence phrased as suggested here, just to stress that our data elicit new thinking about Upsilon suppression in dAu.


6.  At LHC energy, the anisotropic model of Strickland reproduces well
the centrality dependence of R_AA but not the rapidity dependence,
see, e.g. the final session of the recent hard probes meeting in South
Africa.

The predictions from Strickland shown in Fig. 4 are rapidity-dependent. Our data agree fairly well in both rapidity ranges. The rapidity range examined by the slides shown in Strickland's Hard Probes talk (|y|<4.0) is much wider than the ranges we examine (|y|<1.0). Looking at the models, we do not expect to see much variation at all in the range we examine which is consistent with our observation. To see this variation, we would need to examine a much wider rapidity range which is not the focus of this paper.  (To constrain models at larger rapidities, it will be interesting to see the PHENIX results near |y|~2, which should be submitted for publication soon.)
 

7.  The presentation in Fig. 6 on the quantitative evaluation of
different model assumptions compared to data depends again strongly on
the size of the systematic errors, see the comment in section 2.

We've added a table and a new plot summarizing efficiencies and systematics.

 

----------------------------------------------------------------------------------------------------------------------------------------------------------------------

Reviewer #2: I have read the manuscript PLB-D-13-01645 submitted to me for review.
The authors present a detail analysis on the suppression of Y production in d+Au and Au+Au collisions at sqrt(s_NN)=200 GeV using the STAR detector at RHIC. The article is very well written and deserves publication. However, I would like to suggest considering the following remarks to improve the understandability of the article:

1. Page 1, column 1, paragraph 1: The now accepted value for the critical temperature (chiral transition) is Tc = 150 - 160 MeV (depending on the exact definition of the observables). Reference 3 is outdated and should be replaced by more recent publications, i.e. arXiv:1005.3508 [hep-lat]

See response to first reviewer's comment #1. References have been updated. Furthermore, the Tc noted here (the chiral transition) is not the relevant phase transition for quarkonium suppression. The more relevant one is the deconfinement transition (which is somewhat broad as currently noted in the text).


2. Page 2, column 2, paragraph 1: Please quantify the corrections due to the trigger bias w.r.t. the event centrality. Same for the tracking efficiency as a function of N_part. How does acceptance times efficiency for detecting Y as a function of rapidity and N_part looks like?

Added a new figure summarizing the total efficiency as a function of N_part and rapidity.


3. Page 3, column 1, paragraph 1: statement "some information will be lost" is too general! What are the systematic uncertainties arising from the different methods (same-event like-sign CB, fit to the CB) of the combinatorial background subtraction? What is the signal significance, in particular in the d+Au measurement? How does the signal looks like after CB and physical background subtraction? Systematic errors should be clearly mentioned.

Many thanks to the referee for mentioning this. In the process of investigating the systematic difference between extracting the upsilon yield through simultaneous fitting compared to background subtraction, we also investigated the effects of chi^2 fits compared to likelihood fits. We used chi^2 fits in our original submission. We were aware that extracting yields using a chi^2 fit introduces a bias. We attempted to minimize this bias by extracting the yields using the integral of the data. However, we still needed to subtract the backgrounds and the bakground yields came from the fits. Through investigation, we have demonstrated that chi^2 fits systematically underestimate the background yield, leading to an overestimate of the upsilon signals.

We have studied these effects through various MC simulations in order to extract the biases. The likelihood fits have negligible biases. Furthermore, and to get back to the original question posed by the referee, in these simulations we obtained the variance of our results when doing simultaneous fits when compared to background-subtracted fits. We found a reduction in the variance when using simultaneous fits which was our original impetus. We also found no systematic effect in the expectation values of the yields obtained by the two different fitting methods. However, given the reduction in the variance of the extracted yield (i.e. in their error) in the simultaneous fit, we favor this method since it introduces a smaller uncertainty. We have redone all of our fits using the likelihood method and we corrected for any extraction biases seen through simulation.

Regarding signal significance, in all cases we see significant signals in d+Au. This can be infered by examining Fig. 3a and comparing the size of the statistical+fit error bars to the measured value of the cross section. This ratio is a good indicator of the statistical significance of our signal. For example, the dAu signal at |y|<0.5 has a significance of 11.7/3.2 = 3.7 sigma.

Since the referee also asked about systematic uncertainties, we have added a full table covering all measured sources of systematic uncertainty and added additional comments about the main sources in the text. 

4. Fig 1: It would be easier for reader if the range of the y axis would be the same in Fig 1a and Fig 1b. Why is the explanation of the grey curves in the figure discussed in this complicated way, to my understanding the gray band simply shows the pp yield scaled by the number of binary collisions? If so, the label could read simply pp*<N_coll>.

The axes in Figs. 1a and 1b now match. We've relabeled the gray band.
 

5. Fig 1a: From where the line shape for pp comes from? It seems NOT to fit experimental data, i.e. all data points around 9 GeV/c^2 and below. Is it then evident to take as a cross section the integral of the data points?

The line shape in pp comes from simulations embedded in real data. Below 9 GeV, the lineshape threads between high and low datapoints. It cannot fit exactly to all of them without introducing wiggles in the function.
 

6. Page 3, column 2, paragraph 1: How was the measured Y(1S+2S+3S) yield transformed to cross section?

The cross section was calculated by correcting for EID efficiencies, triggering efficiencies, and acceptance to get a corrected yield. We then divided by the integrated luminosity to get a cross section.
 


7. Page 3, column 2, paragraph 3 (wording):  "Hence, averaging between forward and reverse rapidities is not warranted as it is in
p+p." -->  "Hence, averaging between forward and backward rapidities is not justified as it is in p+p." sounds more understandable.

Since the two words are very closely related (Merriam-Webster includes "justification" as one of the definitions of "warrant"), this is more a matter of style. The authors prefer the word "warranted."


8. Page 4: Try to arrange the placement of Figs such that there will not be a single line of the text within one column.

Done

9. Fig 2: also here Fig a and Fig b could be presented with the same range on the Y axis, e.g. from -3 to 3.

Done.


10. Fig 2a : what is shown here is Y(1S+2S+3S), moreover PHENIX results on Y -> mu+mu- are shown in the same plot, that is why the figure label should be changed, i.e. Y->e+e- should be replaced by Y(1S+2S+3S)

We changed the label to Y->l+l- to represent leptons.

11. Page 4, column 2, paragraph 2: <N_coll> (not <N_bin>) is commonly used as notation for the number of binary collisions. Sigma_AA is sigma^tot_AA (same for pp). It is important to indicate in the text the values for the total inelastic cross sections in pp, dA and AA and <N_coll> used to calculate R_AA.

Now using <N_coll>. Inelastic cross sections are provided inline.


 

12. Page 4, column 2, paragraph 3: In view of the discussion would it be helpful to also show R_AuAu vs. Rapidity?

We have addede new plots to the paper, and given the length considerations, and that this plot doesn't really add any new information beyond the existing tables and figures, for the sake of space, we would prefer to leave this plot out.


13. Page 6, column 1, paragraph 1: Which function has been used to fit the CB - exponential? Again, what are the systematic uncertainty arising from the different methods (same-event like-sign CB, fit to the CB) of the combinatorial background subtraction. See also comment 4. concerning the label.

The function used to model the CB is now discussed in the text. Systematics from the fit methods are summarized in Tab. I.

14. Page 6, column 1, paragraph 2: The statement "Similar suppression is found by CMS in PbPb collisions (37)" should be moved to the paragraph 4 where the authors discuss Y(1S) suppression. Actually, for the same value of N_part=325 R_AuAu=0.54+-0.09 as for R_PbPb=0.45

Done.

15. Page 6, column 1, paragraph 4: How did the authors derived: R_AA(1S+2S+3S) = R_AA(1S)*0.69?

We calculate this number by relating the two nuclear modification factors. For the (1S+2S+3S) case, this needs the ratio of the yields of (1S+2S+3S) in AA to the same yield in pp.  Since the R_AA(1S) is the ratio of the yield of the 1S state in AA to that in pp, one can take this out as a common factor in the R_AA(1S+2S+3S), obtaining the relation R_AA(1S+2S+3S)=R_AA(1S) * (1+ N_AA(2S+3S)/N_AA(1S))/(1+N_pp(2S+3S)/N_pp(1S)), where N_AA refers to the yield obtained in AA collisions and N_pp refers to the yield obtained in pp collisions.  This equation makes no assumptions.  When one takes the additional hypothesis that the yield in AA of the excited states is zero, the factor becomes 1/(1+N_pp(2S+3S)/N_pp(1S)). So with the ratio of excited states to ground state in pp collisions, one can find the multiplicative factor.  We calculated this ratio in two manners: first, by using the PDG branching ratios together with NLO pQCD calculations for the upsilon production cross sections (from Ref 21 by Frawley, Ullrich, and Vogt), and, second, by using the measured 2S/1S and 3S/1S ratios. For example, these ratios have been measured at sqrt{s}=38.8 GeV and also at sqrt(s)=2.76 TeV from CMS, and are relatively independent of sqrt(s), or even whether the collision system is pp or pA. In the first case where we used the pQCD cross section and PDF branching ratios, we get 0.69 for the multiplicative factor. In the second case where we only used measured ratios, we get 0.72 +/- 0.02.  The difference between using the low-sqrt(s)-pA data or the CMS pp data at 2.76 TeV is 0.01, which is smaller than the statistical error of the CMS data.
We had aimed to keep the text brief, since we were mindful of the space constraints, but given this question, we have added a few more sentences and references to clarify the R_AA(1S+2S+3S)=R_AA(1S)*0.7 statement, and also reduced our significant figures, quoting only a factor of 0.7.
 

16. Page 6, column 2, paragraph 2: What are the uncertainties on Drell-Yan and bbbar cross sections and how does it influence the significance of the signal.

Various normalizations are used in the fit. This is accounted for in the correllation

17. http://arxiv.org/pdf/1109.3891.pdf reports on the first measurement of the Y nuclear modification factor with STAR. It is probably worth to mention this work in the ms.

We are certainly aware of the proceeding mentioned here, which showed preliminary results for these analyses. The author of the proceedings was a member of the institute where the primary analsys shown in this paper was done. The reason we omit the citation to this and to other proceedings where the preliminary results have been shown is that it is a policy of the STAR collaboration to not cite our own proceedings showing preliminary data. This is partly with the goal to make it clear that the final results presented in a given paper, which have gone through the full collaboration review and through the external peer-review process, are the ones that should be referenced once they are available.


18. The R_AA of J/psi (p_T > 5 GeV), Y(1S) and an upper limit on the R_AA (2S+3S) was obtained in STAR. I would like to suggest to add a plot showing R_AA as a function of binding energy as a summary figure (also as a key figure to the long discussion on the extraction of the upper limit on R_AA(2S+3S)).

Good idea. We added this figure towards the end of the paper.


In summary, this ms. contains very interesting results and I propose publication in Phys. Letter B after the authors have taken care of the remarks above.

We thank the referee for her/his comments and remarks, which have helped improve the paper.  We hope that we have addressed the issues raised, and adequately answered the questions posed, and look forward to the publication of the paper.