Rivet analyses referenceALICE_2018_I1669819Strange meson production at 5 TeVExperiment: ALICE (LHC) Inspire ID: 1669819 Status: VALIDATED Authors:
Beam energies: (2510.0, 2510.0); (522080.0, 522080.0) GeV Run details:
We report measurements of the production of prompt $D^0$, $D^+$, $D^{∗+}$ and $D^+_\text{s}$ mesons in Pb-Pb collisions at the centre-of-mass energy per nucleon-nucleon pair $\sqrt{s_\text{NN}}$ = 5.02 TeV, in the centrality classes 0-10%, 30-50% and 60-80%. The $D$-meson production yields are measured at mid-rapidity (|y|<0.5) as a function of transverse momentum ($p_\text{T}$). The $p_\text{T}$ intervals covered in central collisions are $1 < p_\text{T} < 50$ GeV/c for $D^0$, $2 < p_\text{T} < 50$ GeV/c for $D^+$, $3 < p_\text{T} < 50$ GeV/c for $D^{*+}$, and $4 < p_\text{T} < 16$ GeV/c for $D^+_\text{s}$ mesons. The nuclear modification factors ($R_\text{AA}$) for non-strange $D$ mesons ($D^0$, $D^+$, $D^{∗+}$) show minimum values of about 0.2 for $p_\text{T}$ = 6-10 GeV/c in the most central collisions and are compatible within uncertainties with those measured at $\sqrt{s_\text{NN}}$ = 2.76 TeV. For $D^+_\text{s}$ mesons, the values of $R_\text{AA}$ are larger than those of non-strange $D$ mesons, but compatible within uncertainties. In central collisions the average $R_\text{AA}$ of non-strange $D$ mesons is compatible with that of charged particles for $p_\text{T} > 8$ GeV/c, while it is larger at lower $p_\text{T}$. The nuclear modification factors for strange and non-strange $D$ mesons are also compared to theoretical models with different implementations of in-medium energy loss. Source code: ALICE_2018_I1669819.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Tools/AliceCommon.hh"
5#include "Rivet/Projections/AliceCommon.hh"
6#include "Rivet/Projections/UnstableParticles.hh"
7
8namespace Rivet {
9
10
11 /// @brief Strange meson production at 5 TeV
12 class ALICE_2018_I1669819 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2018_I1669819);
17
18 void mkAverage(const string& avgname, const vector<string>& estnames) {
19 for (auto& b : _e[avgname]->bins()) {
20 double wtotal = 0., vtotal = 0., etotal = 0.;
21 for (const string& ename : estnames) {
22 const auto& est = _e[ename]->binAt(b.xMid());
23 if (!_e[ename]->isVisible(est.index())) continue;
24 const double w = 1.0 / sqr( est.relErrAvg() );
25 wtotal += w;
26 vtotal += est.val() * w;
27 etotal += sqr( est.errAvg() * w );
28 }
29 b.set(vtotal / wtotal, sqrt(etotal) / wtotal);
30 }
31 }
32
33 /// @name Analysis methods
34 ///@{
35
36 /// Book histograms and initialise projections before the run
37 void init() {
38
39 // Initialise and register projections
40
41 declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M", "V0M");
42
43 const ALICE::PrimaryParticles app(Cuts::abseta < 0.8 && Cuts::pT > 1*GeV && Cuts::abscharge > 0);
44 declare(app, "app");
45
46 const UnstableParticles ufsD0(Cuts::absrap < 0.5 && Cuts::pT > 1*GeV && Cuts::abspid == PID::D0);
47 declare(ufsD0, "ufsD0");
48
49 const UnstableParticles ufsDplus(Cuts::absrap < 0.5 && Cuts::pT > 2*GeV && Cuts::abspid == PID::DPLUS);
50 declare(ufsDplus, "ufsDplus");
51
52 const UnstableParticles ufsDstar(Cuts::absrap < 0.5 && Cuts::pT > 1*GeV && Cuts::abspid == PID::DSTARPLUS);
53 declare(ufsDstar, "ufsDstar");
54
55 const UnstableParticles ufsDs(Cuts::absrap < 0.5 && Cuts::pT > 2*GeV && Cuts::abspid == PID::DSPLUS);
56 declare(ufsDs, "ufsDs");
57
58 book(_c["sow_pp5TeV"], "_sow_pp5TeV");
59 book(_c["sow_pp2TeV"], "_sow_pp2TeV");
60 book(_c["sow_PbPb2TeV"], "_sow_PbPb2TeV");
61
62 unsigned int pt_idx = 1, part_idx = 13, beam_idx = 25;
63 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
64 book(_c["sow_PbPb5TeV_"+cen], "_sow_PbPb5TeV_"+cen);
65 for (const string& part : vector<string>{ "D0", "Dplus", "Dstar", "Ds" }) {
66
67 book(_h[part+"Pt_"+cen], pt_idx++, 1, 1);
68
69 if (part != "D0") {
70 const string ref1name = mkAxisCode(part_idx++, 1, 1);
71 const YODA::Estimate1D ref1 = refData(ref1name);
72 string rname(part+"_D0"+cen);
73 book(_h["num_"+rname], "_num_"+rname, ref1);
74 book(_h["den_"+rname], "_den_"+rname, ref1);
75 book(_e[rname], ref1name);
76 if (part == "Ds") {
77 const string ref2name = mkAxisCode(part_idx++, 1, 1);
78 const YODA::Estimate1D ref2 = refData(ref2name);
79 rname = part+"_Dplus"+cen;
80 book(_h["num_"+rname], "_num_"+rname, ref2);
81 book(_h["den_"+rname], "_den_"+rname, ref2);
82 book(_e[rname], ref2name);
83 }
84 }
85 const string brefname = mkAxisCode(beam_idx++, 1, 1);
86 const YODA::Estimate1D bref = refData(brefname);
87 string rname(part+"PbPb_pp"+cen);
88 book(_h["num_"+rname], "_num_"+rname, bref);
89 book(_h["den_"+rname], "_den_"+rname, bref);
90 book(_e[rname], brefname);
91
92 if (part == "Ds") {
93 book(_e["average"+cen], beam_idx++, 1, 1);
94 }
95 }
96 }
97
98 //Table 40 is PbPb 2.76 TeV
99 for (const string& part : vector<string>{ "D0", "Dplus", "Dstar" }) {
100 vector<double> bins;
101 if (part != "D0") bins = {3., 4., 5., 6., 8., 12., 16., 24., 36.};
102 else bins = {1., 2., 3., 4., 5., 6., 8., 12., 16., 24.};
103 book(_h["num_"+part+"PbPb_pp2TeV"], "_"+part+"_PbPb", bins);
104 book(_h["den_"+part+"PbPb_pp2TeV"], "_"+part+"_pp", bins);
105 book(_e[part+"PbPb_pp2TeV"], "_"+part+"_PbPb_pp", bins);
106 }
107 book(_e["average2TeV"], 40, 1, 1);
108
109 unsigned int avg_idx = 41;
110 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
111 const string refname = mkAxisCode(avg_idx++, 1, 1);
112 const Estimate1D& ref = refData(refname);
113 book(_h["num_mult_PbPb_pp"+cen], "_num_mult_PbPb_pp_"+cen, ref);
114 book(_h["den_mult_PbPb_pp"+cen], "_den_mult_PbPb_pp_"+cen, ref);
115 book(_e["mult_PbPb_pp"+cen], "_mult_PbPb_pp_"+cen);
116 book(_e["avgDmult"+cen], refname);
117 }
118
119 }
120
121
122 /// Perform the per-event analysis
123 void analyze(const Event& event) {
124
125 const ParticlePair& beam = beams();
126 string CollSystem = "Empty";
127 const double NN = 208;
128
129 if (beam.first.pid() == PID::LEAD && beam.second.pid() == PID::LEAD) {
130 CollSystem = "PBPB";
131 if(fuzzyEquals(sqrtS()/GeV, 2760*NN, 1E-3)) CollSystem += "2TeV";
132 else if(fuzzyEquals(sqrtS()/GeV, 5020*NN, 1E-3)) CollSystem += "5TeV";
133 }
134 if (beam.first.pid() == PID::PROTON && beam.second.pid() == PID::PROTON) {
135 CollSystem = "PP";
136 if(fuzzyEquals(sqrtS()/GeV, 2760, 1E-3)) CollSystem += "2TeV";
137 else if(fuzzyEquals(sqrtS()/GeV, 5020, 1E-3)) CollSystem += "5TeV";
138 }
139
140 const Particles& particlesD0 = apply<UnstableParticles>(event,"ufsD0").particles();
141 const Particles& particlesDplus = apply<UnstableParticles>(event,"ufsDplus").particles();
142 const Particles& particlesDstar = apply<UnstableParticles>(event,"ufsDstar").particles();
143 const Particles& particlesDs = apply<UnstableParticles>(event,"ufsDs").particles();
144 const Particles& chParticles = apply<ALICE::PrimaryParticles>(event, "app").particles();
145
146 if (CollSystem == "PP5TeV") {
147
148 _c["sow_pp5TeV"]->fill();
149
150 for (const Particle& p : particlesD0) {
151 if (p.fromBottom()) continue;
152 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
153 _h["den_D0PbPb_pp"+cen]->fill(p.pT()/GeV);
154 }
155 }
156
157 for (const Particle& p : particlesDplus) {
158 if (p.fromBottom()) continue;
159 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
160 _h["den_DplusPbPb_pp"+cen]->fill(p.pT()/GeV);
161 }
162 }
163
164 for (const Particle& p : particlesDstar) {
165 if(p.fromBottom()) continue;
166 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
167 _h["den_DstarPbPb_pp"+cen]->fill(p.pT()/GeV);
168 }
169 }
170
171 for (const Particle& p : particlesDs) {
172 if (p.fromBottom()) continue;
173 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
174 _h["den_DsPbPb_pp"+cen]->fill(p.pT()/GeV);
175 }
176 }
177
178 for (const Particle& p : chParticles) {
179 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
180 _h["den_mult_PbPb_pp"+cen]->fill(p.pT()/GeV);
181 }
182 }
183
184 }
185
186 if (CollSystem == "PP2TeV") {
187
188 _c["sow_pp2TeV"]->fill();
189
190 for (const Particle& p : particlesD0) {
191 if (p.fromBottom()) continue;
192 _h["den_D0PbPb_pp2TeV"]->fill(p.pT()/GeV);
193 }
194
195 for (const Particle& p : particlesDplus) {
196 if (p.fromBottom()) continue;
197 _h["den_DplusPbPb_pp2TeV"]->fill(p.pT()/GeV);
198 }
199
200 for (const Particle& p : particlesDstar) {
201 if (p.fromBottom()) continue;
202 _h["den_DstarPbPb_pp2TeV"]->fill(p.pT()/GeV);
203 }
204
205 }
206
207
208 // The centrality projection.
209 const CentralityProjection& centProj = apply<CentralityProjection>(event,"V0M");
210 const double cent = centProj();
211 if(cent >= 80.) vetoEvent;
212
213 if (CollSystem == "PBPB5TeV") {
214
215 string cen("");
216 if (cent < 10.) cen = "00-10";
217 else if(cent >= 30. && cent < 50.) cen = "30-50";
218 else if (cent >= 60. && cent < 80.) cen = "60_80";
219
220 if (cen != "") {
221
222 _c["sow_PbPb5TeV"+cen]->fill();
223
224 for (const Particle& p : particlesD0) {
225 if (p.fromBottom()) continue;
226 _h["D0Pt"+cen]->fill(p.pT()/GeV);
227 _h["den_Dplus_D0"+cen]->fill(p.pT()/GeV);
228 _h["den_Dstar_D0"+cen]->fill(p.pT()/GeV);
229 _h["den_Ds_D0"+cen]->fill(p.pT()/GeV);
230 _h["num_D0PbPb_pp"+cen]->fill(p.pT()/GeV);
231 }
232
233 for (const Particle& p : particlesDplus) {
234 if (p.fromBottom()) continue;
235 _h["DplusPt_"+cen]->fill(p.pT()/GeV);
236 _h["num_Dplus_D0"+cen]->fill(p.pT()/GeV);
237 _h["den_Ds_Dplus"+cen]->fill(p.pT()/GeV);
238 _h["num_DplusPbPb_pp"+cen]->fill(p.pT()/GeV);
239 }
240
241 for (const Particle& p : particlesDstar) {
242 if (p.fromBottom()) continue;
243 _h["DstarPt_"+cen]->fill(p.pT()/GeV);
244 _h["num_Dstar_D0"+cen]->fill(p.pT()/GeV);
245 _h["num_DstarPbPb_pp"+cen]->fill(p.pT()/GeV);
246 }
247
248 for (const Particle& p : particlesDs) {
249 if (p.fromBottom()) continue;
250 _h["DsPt_"+cen]->fill(p.pT()/GeV);
251 _h["num_Ds_D0"+cen]->fill(p.pT()/GeV);
252 _h["num_Ds_Dplus"+cen]->fill(p.pT()/GeV);
253 _h["num_DsPbPb_pp"+cen]->fill(p.pT()/GeV);
254 }
255
256 for (const Particle& p : chParticles){
257 _h["num_mult_PbPb_pp"+cen]->fill(p.pT()/GeV);
258 }
259 }
260
261 }
262 else if (CollSystem == "PBPB2TeV") {
263 if (cent < 10.) {
264 _c["sow_PbPb2TeV"]->fill();
265 for (const Particle& p : particlesD0) {
266 if (p.fromBottom()) continue;
267 _h["num_D0PbPb_pp2TeV"]->fill(p.pT()/GeV);
268 }
269 for (const Particle& p : particlesDplus) {
270 if (p.fromBottom()) continue;
271 _h["num_DplusPbPb_pp2TeV"]->fill(p.pT()/GeV);
272 }
273 for (const Particle& p : particlesDstar) {
274 if (p.fromBottom()) continue;
275 _h["num_DstarPbPb_pp2TeV"]->fill(p.pT()/GeV);
276 }
277 }
278 }
279
280 }
281
282
283 /// Normalise histograms etc., after the run
284 void finalize() {
285
286 for (auto& item : _h) {
287 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80" }) {
288 if (item.first.substr(0, 4) == "den_" && item.first.find("_pp") != string::npos) {
289 scale(item.second, _n[cen]/_c["sow_pp5TeV"]->sumW());
290 continue;
291 }
292 if (item.first.find("2TeV") != string::npos) {
293 const double sf = (item.first.find("num_") != string::npos)? 1.0 : _n["00-10"];
294 scale(item.second, sf / _c["sow_pp2TeV"]->sumW());
295 continue;
296 }
297 if (item.first.find("mult") != string::npos) {
298 if (item.first.find("num_") != string::npos) {
299 scale(item.second, 1./_c["sow_PbPb5TeV_"+cen]->sumW());
300 }
301 else {
302 scale(item.second, _n[cen] / _c["sow_pp5TeV"]->sumW());
303 }
304 continue;
305 }
306 if (item.first.find(cen) != string::npos) {
307 const double sf = (item.first.find("PbPb") != string::npos)? 1.0 : 0.5;
308 scale(item.second, sf / _c["sow_PbPb5TeV_"+cen]->sumW());
309 }
310 }
311 }
312
313 for (auto& item : _e) {
314 if (item.first.find("_") == string::npos) continue;
315 divide(_h["num_"+item.first], _h["den_"+item.first], item.second);
316 }
317
318 for (const string& cen : vector<string>{ "00-10", "30-50", "60-80", "2TeV" }) {
319 mkAverage("average"+cen, { "D0PbPb_pp"+cen, "DplusPbPb_pp"+cen, "DstarPbPb_pp"+cen });
320 if (cen == "2TeV") continue;
321 divide(_e["average"+cen], _e["mult_PbPb_pp"+cen], _e["avgDmult"+cen]);
322 }
323
324 }
325
326 ///@}
327
328
329 /// @name Histograms
330 ///@{
331
332 map<string, Histo1DPtr> _h;
333 map<string, CounterPtr> _c;
334 map<string, Estimate1DPtr> _e;
335 map<string, double> _n{ {"00-10", 1572.}, { "30-50", 264.8 }, { "60-80", 28.31} };
336
337 ///@}
338
339
340 };
341
342
343 RIVET_DECLARE_PLUGIN(ALICE_2018_I1669819);
344
345}
|