rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ALICE_2018_I1669819

Charm meson production at 5 TeV
Experiment: ALICE (LHC)
Inspire ID: 1669819
Status: VALIDATED
Authors:
  • Antonio Silva
  • Yonne Lourens
  • Christian Gutschow
References: Beams: p+ p+, 1000822080 1000822080
Beam energies: (2510.0, 2510.0); (522080.0, 522080.0) GeV
Run details:
  • Minimum bias events

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_CENT_PBPB", "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}