rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ALICE_2021_I1946131

Prompt charm-meson production at 5 TeV
Experiment: ALICE (LHC)
Inspire ID: 1946131
Status: VALIDATED
Authors:
  • Yonne Lourens
References: Beams: p+ p+, 1000822080 1000822080
Beam energies: (2510.0, 2510.0); (522080.0, 522080.0) GeV
Run details:
  • Minimum bias events

The production of prompt $D^0$, $D^+$, and $D^{*+}$ mesons was measured at midrapidity ($|y| < $0.5) in Pb-Pb collisions at the centre-of-mass energy per nucleon-nucleon pair $\sqrt{s}_\text{NN} = 5.0$ TeV with the ALICE detector at the LHC. The $D$ mesons were reconstructed via their hadronic decay channels and their production yields were measured in central (0-10%) and semicentral (30-50%) collisions. The measurement was performed up to a transverse momentum ($p_\text{T}$) of 36 or 50 GeV/c depending on the $D$ meson species and the centrality interval. For the first time in Pb-Pb collisions at the LHC, the yield of $D^0$ mesons was measured down to $p_\text{T} =$ 0, which allowed a model-independent determination of the $p_\text{T}$-integrated yield per unit of rapidity ($\text{d}N/\text{d}y$). A maximum suppression by a factor 5 and 2.5 was observed with the nuclear modification factor ($R_\text{AA}$) of prompt $D$ mesons at $p_\text{T} =$ 6-8 GeV/c for the 0-10% and 30-50% centrality classes, respectively. The $D$-meson $R_\text{AA}$ is compared with that of charged pions, charged hadrons, and $J/\psi$ mesons as well as with theoretical predictions. The analysis of the agreement between the measured $R_\text{AA}$, elliptic ($v_2$) and triangular ($v_3$) flow, and the model predictions allowed us to constrain the charm spatial diffusion coefficient $D_\text{s}$. Furthermore the comparison of $R_\text{AA}$ and $v_2$ with different implementations of the same models provides an important insight into the role of radiative energy loss as well as charm quark recombination in the hadronisation mechanisms.

Source code: ALICE_2021_I1946131.cc
  1// -*- C++ -*-
  2#include <iostream>
  3#include "Rivet/Analysis.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Tools/AliceCommon.hh"
  7#include "Rivet/Projections/AliceCommon.hh"
  8#include "Rivet/Projections/UnstableParticles.hh"
  9#include "Rivet/Projections/CentralityProjection.hh"
 10#include "Rivet/Analysis.hh"
 11#include "Rivet/Projections/UnstableParticles.hh"
 12#include "Rivet/Tools/Cuts.hh"
 13#include "Rivet/Projections/HepMCHeavyIon.hh"
 14
 15namespace Rivet {
 16
 17  /// @brief Prompt strange-meson production at 5 TeV
 18  class ALICE_2021_I1946131 : public Analysis {
 19  public:
 20
 21    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2021_I1946131);
 22
 23    void mkAverage(const string& avgname, const vector<string>& estnames) {
 24      for (auto& b : _e[avgname]->bins()) {
 25        double wtotal = 0., vtotal = 0., etotal = 0.;
 26        for (const string& ename : estnames) {
 27          const auto& est = _e[ename]->binAt(b.xMid());
 28          if (!_e[ename]->isVisible(est.index()))  continue;
 29          const double w = 1.0 / sqr( est.relErrAvg() );
 30          wtotal += w;
 31          vtotal += est.val() * w;
 32          etotal += sqr( est.errAvg() * w );
 33        }
 34        b.set(vtotal / wtotal, sqrt(etotal) / wtotal);
 35      }
 36    }
 37
 38    void init() {
 39
 40      declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_CENT_PBPB", "V0M","V0M");
 41
 42      const UnstableParticles ufsD0(Cuts::absrap < 0.5 && Cuts::pT > 0.*GeV && Cuts::abspid == 421);
 43      declare(ufsD0, "ufsD0");
 44
 45      const UnstableParticles ufsDplus(Cuts::absrap < 0.5 && Cuts::pT > 2.*GeV && Cuts::abspid == 411);
 46      declare(ufsDplus, "ufsDplus");
 47
 48      const UnstableParticles ufsDstar(Cuts::absrap < 0.5 && Cuts::pT > 2.*GeV && Cuts::abspid == 413);
 49      declare(ufsDstar, "ufsDstar");
 50
 51      book(_c["sow_pp5TeV"], "_sow_pp5TeV");
 52
 53      size_t pt_idx = 1, part_idx = 7, beam_idx = 11, avg_idx = 17;
 54      for (const string& part : vector<string>{ "D0", "Dplus", "Dstar" }) {
 55        for (const string& cent : vector<string>{ "00-10", "30-50" }) {
 56          book(_h[part+"Pt_"+cent], pt_idx++, 1, 1);
 57          if (part != "D0") {
 58            const string refname = mkAxisCode(part_idx++, 1, 1);
 59            const YODA::Estimate1D ref = refData(refname);
 60            string rname(part+"_D0"+cent);
 61            book(_h["num_"+rname], "_num_"+rname, ref);
 62            book(_h["den_"+rname], "_den_"+rname, ref);
 63            book(_e[rname], refname);
 64          }
 65          else {
 66            book(_c["sow_PbPb5TeV_"+cent], "_sow_PbPb5TeV_"+cent);
 67            book(_e["average"+cent], avg_idx++, 1, 1);
 68          }
 69          size_t offset = cent == "00-10"? 0 : 3;
 70          const string brefname = mkAxisCode(beam_idx + offset, 1, 1);
 71          const YODA::Estimate1D bref = refData(brefname);
 72          string rname(part+"PbPb_pp"+cent);
 73          book(_h["num_"+rname], "_num_"+rname, bref);
 74          book(_h["den_"+rname], "_den_"+rname, bref);
 75          book(_e[rname], brefname);
 76        }
 77        ++beam_idx;
 78      }
 79
 80      const string RAAname = mkAxisCode(19, 1, 1);
 81      const Estimate1D& RAAref = refData(RAAname);
 82      book(_h["num_RAAPbPb_pp"], "_num_RAA", RAAref);
 83      book(_h["den_RAAPbPb_pp"], "_den_RAA", RAAref);
 84      book(_e["RAAPbPb_pp"], RAAname);
 85
 86    }
 87
 88    void analyze(const Event& event) {
 89
 90      const ParticlePair& beam = beams();
 91      string CollSystem = "Empty";
 92      const double NN = 208;
 93
 94      if (beam.first.pid() == PID::LEAD && beam.second.pid() == PID::LEAD) {
 95        CollSystem = "PBPB";
 96        if (fuzzyEquals(sqrtS()/GeV, 5020*NN, 1E-3)) CollSystem += "5TeV";
 97      }
 98      if (beam.first.pid() == PID::PROTON && beam.second.pid() == PID::PROTON) {
 99        CollSystem = "PP";
100        if (fuzzyEquals(sqrtS()/GeV, 5020, 1E-3)) CollSystem += "5TeV";
101      }
102
103      const Particles& particlesD0 = apply<UnstableParticles>(event,"ufsD0").particles();
104      const Particles& particlesDplus = apply<UnstableParticles>(event,"ufsDplus").particles();
105      const Particles& particlesDstar = apply<UnstableParticles>(event,"ufsDstar").particles();
106
107      if (CollSystem == "PP5TeV") {
108        _c["sow_pp5TeV"]->fill();
109
110        for (const Particle& p : particlesD0) {
111          if (p.fromBottom()) continue;
112          for (const string& cent : vector<string>{ "00-10", "30-50" }) {
113            _h["den_D0PbPb_pp"+cent]->fill(p.pT()/GeV);
114            _h["den_RAAPbPb_pp"]->fill(cent == "00-10"? 5. : 40.);
115          }
116        }
117
118        for (const Particle& p : particlesDplus) {
119          if (p.fromBottom()) continue;
120          for (const string& cent : vector<string>{ "00-10", "30-50" }) {
121            _h["den_DplusPbPb_pp"+cent]->fill(p.pT()/GeV);
122          }
123        }
124
125        for (const Particle& p : particlesDstar) {
126          if(p.fromBottom()) continue;
127          for (const string& cent : vector<string>{ "00-10", "30-50" }) {
128            _h["den_DstarPbPb_pp"+cent]->fill(p.pT()/GeV);
129          }
130        }
131      }
132
133      const CentralityProjection& centProj = apply<CentralityProjection>(event, "V0M");
134
135      const double cent_val = centProj();
136
137      if (cent_val >= 50.) vetoEvent;
138
139      if (CollSystem == "PBPB5TeV") {
140        string cent("");
141        if (cent_val < 10.)  cent = "00-10";
142        else if(cent_val >= 30. && cent_val < 50.)  cent = "30-50";
143
144        if (cent == "")  vetoEvent;
145        _c["sow_PbPb5TeV_"+cent]->fill();
146        for (const Particle& p : particlesD0) {
147          if (p.fromBottom()) continue;
148          _h["D0Pt_"+cent]->fill(p.pT()/GeV);
149          _h["den_Dplus_D0"+cent]->fill(p.pT()/GeV);
150          _h["den_Dstar_D0"+cent]->fill(p.pT()/GeV);
151          _h["num_D0PbPb_pp"+cent]->fill(p.pT()/GeV);
152          _h["num_RAAPbPb_pp"]->fill(cent == "00-10"? 5. : 40.);
153        }
154
155        for (const Particle& p : particlesDplus) {
156          if (p.fromBottom()) continue;
157          _h["DplusPt_"+cent]->fill(p.pT()/GeV);
158          _h["num_Dplus_D0"+cent]->fill(p.pT()/GeV);
159          _h["num_DplusPbPb_pp"+cent]->fill(p.pT()/GeV);
160        }
161
162        for (const Particle& p : particlesDstar) {
163          if (p.fromBottom()) continue;
164          _h["DstarPt_"+cent]->fill(p.pT()/GeV);
165          _h["num_Dstar_D0"+cent]->fill(p.pT()/GeV);
166          _h["num_DstarPbPb_pp"+cent]->fill(p.pT()/GeV);
167        }
168      }
169    }
170
171    void finalize() {
172
173      for (auto& item : _h) {
174        if (item.first.find("_D0") == string::npos)  continue;
175        if (item.first.substr(0, 4) == "den_" && item.first.find("_pp") != string::npos) {
176          if (item.first.find("RAA") != string::npos) {
177            item.second->bin(1).scaleW( _n["00-10"]/_c["sow_pp5TeV"]->sumW() );
178            item.second->bin(3).scaleW( _n["30-50"]/_c["sow_pp5TeV"]->sumW() );
179            continue;
180          }
181          const string cent(item.first.substr(item.first.length()-5));
182          scale(item.second, _n[cent]/_c["sow_pp5TeV"]->sumW());
183          continue;
184        }
185        string sfname = "sow_PbPb5TeV_";
186        if (item.first.find("00-10") != string::npos)  sfname += "00-10";
187        else                                           sfname += "30-50";
188        const double sf = (item.first.find("PbPb") != string::npos)? 1.0 : 0.5;
189        if (item.first.find("RAA") != string::npos) {
190          item.second->bin(1).scaleW( sf/_c["sow_pp5TeV00-10"]->sumW() );
191          item.second->bin(3).scaleW( sf/_c["sow_pp5TeV30-50"]->sumW() );
192          continue;
193        }
194        scale(item.second, sf / _c[sfname]->sumW());
195      }
196
197      for (auto& item : _e) {
198        if (item.first.find("_") == string::npos)  continue;
199        divide(_h["num_"+item.first], _h["den_"+item.first], item.second);
200      }
201
202      for (const string& cen : vector<string>{ "00-10", "30-50" }) {
203        mkAverage("average"+cen, { "D0PbPb_pp"+cen, "DplusPbPb_pp"+cen, "DstarPbPb_pp"+cen });
204      }
205
206    }
207
208    map<string, Histo1DPtr> _h;
209    map<string, CounterPtr> _c;
210    map<string, Estimate1DPtr> _e;
211    map<string, double> _n{ {"00-10", 1572.}, {"30-50", 264.8} };
212
213  };
214
215  RIVET_DECLARE_PLUGIN(ALICE_2021_I1946131);
216
217}