Rivet analyses
Bulk properties of identified hadrons in the Au Au medium from the RHIC beam-energy scan
Experiment: STAR (RHIC)
Inspire ID: 1510593
Status: UNVALIDATED
Authors: - Johannes Jahan - Gabriela Pokropska - Maria Stefaniak
References: - Phys.Rev. C96 (2017) no.4, 044904 - DOI: 10.1103/PhysRevC.96.044904 - arXiv: 1701.07065
Beams: 1000791970 1000791970
Beam energies: (758.5, 758.5); (1132.8, 1132.8); (1930.6, 1930.6); (2659.5, 2659.5); (3841.5, 3841.5)GeV
Run details: - Minimum bias AuAu events at various collision energies.
Results of measurements of bulk properties of the matter in Au+Au collisions at energies of the Beam Energy Scan program, using pions, kaons and protons identified : $\frac{dN}{dy}$ spectra, ⟨pT⟩ spectra and particles ratios are shown.
Source
code:STAR_2017_I1510593.cc
// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Analyses/RHICCommon.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/HepMCHeavyIon.hh"
#include "Rivet/Projections/ImpactParameterProjection.hh"
#include "Rivet/Projections/SingleValueProjection.hh"
#include "Rivet/Tools/Percentile.hh"
namespace Rivet {
/// pT distributions, ratios and production yields of hadrons in STAR
class STAR_2017_I1510593 : public Analysis {
public:
RIVET_DEFAULT_ANALYSIS_CTOR(STAR_2017_I1510593);
/// Book histograms and initialise projections before the run
void init() {
// Initialise and register projections
declareCentrality(STAR_BES_Centrality(), "STAR_BES_CALIB", "CMULT", "CMULT");
// The observed particles.
declare(ChargedFinalState(Cuts::abseta < 0.5 && Cuts::absrap < 0.1 && Cuts::pT > 0.2), "CFS");
// Access the HepMC heavy ion info
declare(HepMCHeavyIon(), "HepMC");
/// Booking npart histograms
size_t ih = 0;
for (double eval : allowedEnergies()) {
double enuc = eval / 197.;
const string en = toString(round(enuc));
if (isCompatibleWithSqrtS(eval)) {
_sqs = en;
_sqsEdge = enuc;
// Energy bins for Fig. 25
if (ih == 0)
_sqs25 = "0";
else if (ih == 4)
_sqs25 = "1";
}
for (size_t ic = 0; ic < centAxis.numBins() + 1; ++ic) {
const string cent = toString(round(centAxis.max(ic)));
dualbook(en + "dpT_Pi" + cent, 2 + ih * 6, 1, ic + 1);
dualbook(en + "dpT_Piplus" + cent, 3 + ih * 6, 1, ic + 1);
dualbook(en + "dpT_Kaon" + cent, 4 + ih * 6, 1, ic + 1);
dualbook(en + "dpT_Kaonplus" + cent, 5 + ih * 6, 1, ic + 1);
dualbook(en + "dpT_AntiProton" + cent, 6 + ih * 6, 1, ic + 1);
dualbook(en + "dpT_Proton" + cent, 7 + ih * 6, 1, ic + 1);
}
tribook(en + "npart_PiMinus", 32 + ih, 1, 1);
tribook(en + "npart_PiPlus", 32 + ih, 1, 2);
tribook(en + "npart_KaMinus", 32 + ih, 1, 3);
tribook(en + "npart_KaPlus", 32 + ih, 1, 4);
tribook(en + "npart_AntiProton", 32 + ih, 1, 5);
tribook(en + "npart_Proton", 32 + ih, 1, 6);
book(_p[en + "npart_Piratio"], 42 + ih, 1, 1);
book(_p[en + "npart_Karatio"], 42 + ih, 1, 2);
book(_p[en + "npart_Pratio"], 42 + ih, 1, 3);
book(_p[en + "npart_KaPi"], 47 + ih, 1, 1);
book(_p[en + "npart_AntiPPi"], 47 + ih, 1, 2);
book(_p[en + "npart_KaPiplus"], 47 + ih, 1, 3);
book(_p[en + "npart_PPiplus"], 47 + ih, 1, 4);
++ih;
}
raiseBeamErrorIf(_sqs.empty());
book(_h["npart_PiMinus"], 52, 1, 1);
book(_h["npart_PiPlus"], 52, 1, 2);
book(_h["npart_KaMinus"], 52, 1, 3);
book(_h["npart_KaPlus"], 52, 1, 4);
book(_h["npart_AntiProton"], 52, 1, 5);
book(_h["npart_Proton"], 52, 1, 6);
book(_p["mt_PiPlus"], 53, 1, 1);
book(_p["mt_PiMinus"], 53, 1, 2);
book(_p["mt_KaPlus"], 53, 1, 3);
book(_p["mt_KaMinus"], 53, 1, 4);
book(_p["mt_Proton"], 53, 1, 5);
book(_p["mt_AntiProton"], 53, 1, 6);
book(_p["Piratio"], 54, 1, 1);
book(_p["Karatio"], 54, 1, 2);
book(_p["Pratio"], 54, 1, 3);
book(_p["KaPiplus"], 55, 1, 1);
book(_p["KaPiminus"], 55, 1, 2);
book(_p["yields0"], 56, 1, 1);
book(_p["yields1"], 56, 1, 2);
book(_p["ratios0"], 57, 1, 1);
book(_p["ratios1"], 57, 1, 2);
}
void dualbook(const string& tag, unsigned int d, unsigned int x, unsigned int y) {
book(_c[tag], "TMP/" + mkAxisCode(d, x, y));
book(_h[tag], d, x, y);
}
void tribook(const string& tag, unsigned int d, unsigned int x, unsigned int y) {
book(_c[tag], "TMP/" + mkAxisCode(d, x, y));
book(_h[tag], d, x, y);
book(_p[tag], d + 5, x, y);
}
/// Perform the per-event analysis
void analyze(const Event& event) {
const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
// Require at least two charged particles for the analysis to
// make sense. No further triggers are described in the paper.
const Particles& particles = cfs.particles();
if (particles.size() < 2) return;
/// Determine the centrality
const CentralityProjection& cent = apply<CentralityProjection>(event, "CMULT");
/// Determine the centrality bin
const double c = int(centAxis.index(cent()));
const string centbin = toString(round(centAxis.max(c)));
/// Determine the impact parameter
const HepMCHeavyIon& hi = apply<HepMCHeavyIon>(event, "HepMC");
const double Npart = hi.Npart_targ();
// The following vector contains the counters for all particles used in
// Fig. 25. In the right order : pi+, pi-, K+, K-, p, Antip, Lambda,
// AntiLambda, Xi, AntiXi
map<int, double> nparts, nparts25;
/// Loop over all charged particles of the CFS
for (const Particle& p : cfs.particles()) {
double pT = p.pT() / GeV;
double mass = p.mass() / GeV;
double mTm = sqrt(pT * pT + mass * mass) - mass;
if (p.absrap() < 0.1) {
const PdgId id = p.pid();
switch (id) {
case 211:
if (c < 80) {
_h[_sqs + "dpT_Piplus" + centbin]->fill(pT, 1. / pT);
_h[_sqs + "npart_PiPlus"]->fill(Npart, 1. / (0.2 * 0.5 * Npart));
_p[_sqs + "npart_PiPlus"]->fill(Npart, pT, 5);
}
if (c < 5) {
++nparts25[PID::PIPLUS];
_h["npart_PiPlus"]->fill(_sqsEdge, 1.0 / (0.2 * 0.5 * Npart));
_p["mt_PiPlus"]->fill(_sqsEdge, mTm);
}
++nparts[PID::PIPLUS];
break;
case -211:
if (c < 80) {
_h[_sqs + "dpT_Pi" + centbin]->fill(pT, 1.0 / pT);
_h[_sqs + "npart_PiMinus"]->fill(Npart, 1.0 / (0.2 * 0.5 * Npart));
_p[_sqs + "npart_PiMinus"]->fill(Npart, pT, 5);
}
if (c < 5) {
++nparts25[PID::PIMINUS];
_h["npart_PiMinus"]->fill(_sqsEdge, 1.0 / (0.2 * 0.5 * Npart));
_p["mt_PiMinus"]->fill(_sqsEdge, mTm);
}
++nparts[PID::PIMINUS];
break;
case 321:
if (c < 80) {
_h[_sqs + "dpT_Kaonplus" + centbin]->fill(pT, 1.0 / pT);
_h[_sqs + "npart_KaPlus"]->fill(Npart, 1.0 / (0.2 * 0.5 * Npart));
_p[_sqs + "npart_KaPlus"]->fill(Npart, pT, 5);
}
if (c < 5) {
++nparts25[PID::KPLUS];
_h["npart_KaPlus"]->fill(_sqsEdge, 1.0 / (0.2 * 0.5 * Npart));
_p["mt_KaPlus"]->fill(_sqsEdge, mTm);
}
++nparts[PID::KPLUS];
break;
case -321:
if (c < 80) {
_h[_sqs + "dpT_Kaon" + centbin]->fill(pT, 1.0 / pT);
_h[_sqs + "npart_KaMinus"]->fill(Npart, 1.0 / (0.2 * 0.5 * Npart));
_p[_sqs + "npart_KaMinus"]->fill(Npart, pT, 5);
}
if (c < 5) {
++nparts25[PID::KMINUS];
_h["npart_KaMinus"]->fill(_sqsEdge, 1.0 / (0.2 * 0.5 * Npart));
_p["mt_KaMinus"]->fill(_sqsEdge, mTm);
}
++nparts[PID::KMINUS];
break;
case 2212:
if (c < 80) {
_h[_sqs + "dpT_Proton" + centbin]->fill(pT, 1.0 / pT);
_h[_sqs + "npart_Proton"]->fill(Npart, 1.0 / (0.2 * 0.5 * Npart));
_p[_sqs + "npart_Proton"]->fill(Npart, pT, 5);
}
if (c < 5) {
++nparts25[PID::PROTON];
_h["npart_Proton"]->fill(_sqsEdge, 1.0 / (0.2 * 0.5 * Npart));
_p["mt_Proton"]->fill(_sqsEdge, mTm);
}
++nparts25[PID::PROTON];
break;
case -2212:
if (c < 80) {
_h[_sqs + "dpT_AntiProton" + centbin]->fill(pT, 1.0 / pT);
_h[_sqs + "npart_AntiProton"]->fill(Npart, 1.0 / (0.2 * 0.5 * Npart));
_p[_sqs + "npart_AntiProton"]->fill(Npart, pT, 5);
}
if (c < 5) {
++nparts25[PID::ANTIPROTON];
_h["npart_AntiProton"]->fill(_sqsEdge, 1.0 / (0.2 * 0.5 * Npart));
_p["mt_AntiProton"]->fill(_sqsEdge, mTm);
}
++nparts[PID::ANTIPROTON];
break;
case 3122:
if (c < 5) ++nparts25[PID::LAMBDA];
break;
case -3122:
if (c < 5) ++nparts25[-PID::LAMBDA];
break;
case 3312:
if (c < 5) ++nparts25[PID::XIPLUS];
break;
case -3312:
if (c < 5) ++nparts25[PID::XIMINUS];
break;
}
}
}
/// Particle Ratios
if (!isZero(nparts[PID::PIPLUS])) {
_p[_sqs + "npart_Piratio"]->fill(Npart, nparts[PID::PIMINUS] / nparts[PID::PIPLUS], 5);
_p[_sqs + "npart_KaPiplus"]->fill(Npart, nparts[PID::KPLUS] / nparts[PID::PIPLUS], 5);
_p[_sqs + "npart_PPiplus"]->fill(Npart, nparts[PID::PROTON] / nparts[PID::PIPLUS], 5);
}
if (!isZero(nparts[PID::PIMINUS])) {
_p[_sqs + "npart_KaPi"]->fill(Npart, nparts[PID::KMINUS] / nparts[PID::PIMINUS], 5);
_p[_sqs + "npart_AntiPPi"]->fill(Npart, nparts[PID::ANTIPROTON] / nparts[PID::PIMINUS], 5);
}
if (!isZero(nparts[PID::KPLUS])) {
_p[_sqs + "npart_Karatio"]->fill(Npart, nparts[PID::KMINUS] / nparts[PID::KPLUS], 5);
}
if (!isZero(nparts[PID::PROTON])) {
_p[_sqs + "npart_Pratio"]->fill(Npart, nparts[PID::ANTIPROTON] / nparts[PID::PROTON], 5);
}
/// Particle Yields
if (!_sqs25.empty()) {
if (!isZero(nparts25[PID::PIPLUS])) {
_p["yields" + _sqs25]->fill(1., nparts25[PID::PIMINUS], 5);
_p["ratios" + _sqs25]->fill(1., nparts25[PID::PIMINUS] / nparts25[PID::PIPLUS], 5);
}
if (!isZero(nparts25[PID::KPLUS])) {
_p["yields" + _sqs25]->fill(3., nparts25[PID::KPLUS], 5);
_p["ratios" + _sqs25]->fill(2., nparts25[PID::KMINUS] / nparts25[PID::KPLUS], 5);
}
if (!isZero(nparts25[PID::KMINUS])) {
_p["yields" + _sqs25]->fill(4., nparts25[PID::KMINUS], 5);
}
if (!isZero(nparts25[PID::PROTON])) {
_p["yields" + _sqs25]->fill(5., nparts25[PID::PROTON], 5);
_p["ratios" + _sqs25]->fill(3., nparts25[PID::ANTIPROTON] / nparts25[PID::PROTON], 5);
}
if (!isZero(nparts25[PID::ANTIPROTON])) {
_p["yields" + _sqs25]->fill(6., nparts25[PID::ANTIPROTON], 5);
}
if (!isZero(nparts25[PID::LAMBDA])) {
_p["yields" + _sqs25]->fill(7., nparts25[PID::LAMBDA], 5);
_p["ratios" + _sqs25]->fill(4., nparts25[-PID::LAMBDA] / nparts25[PID::LAMBDA], 5);
}
if (!isZero(nparts25[-PID::LAMBDA])) {
_p["yields" + _sqs25]->fill(8., nparts25[-PID::LAMBDA], 5);
}
if (!isZero(PID::XIPLUS)) {
_p["yields" + _sqs25]->fill(9., nparts25[PID::XIPLUS], 5);
_p["ratios" + _sqs25]->fill(5., nparts25[PID::XIMINUS] / nparts25[PID::XIPLUS], 5);
}
if (!isZero(PID::XIMINUS)) {
_p["yields" + _sqs25]->fill(10., nparts25[PID::XIMINUS], 5);
}
if (!isZero(nparts25[PID::PIMINUS])) {
_p["yields" + _sqs25]->fill(2., nparts25[PID::KMINUS], 5);
_p["ratios" + _sqs25]->fill(6., nparts25[PID::KMINUS] / nparts25[PID::PIMINUS], 5);
_p["yields" + _sqs25]->fill(7., nparts25[PID::ANTIPROTON], 5);
_p["ratios" + _sqs25]->fill(7., nparts25[PID::ANTIPROTON] / nparts25[PID::PIMINUS], 5);
_p["ratios" + _sqs25]->fill(8., nparts25[PID::LAMBDA] / nparts25[PID::PIMINUS], 5);
_p["ratios" + _sqs25]->fill(9., nparts25[PID::XIMINUS] / nparts25[PID::PIMINUS], 5);
}
}
if (!isZero(nparts[PID::PIPLUS])) {
_p["Piratio"]->fill(_sqsEdge, nparts[PID::PIMINUS] / nparts[PID::PIPLUS], 5);
_p["KaPiplus"]->fill(_sqsEdge, nparts[PID::KPLUS] / nparts[PID::PIPLUS], 5);
}
if (!isZero(nparts[PID::PIMINUS])) {
_p["KaPiminus"]->fill(_sqsEdge, nparts[PID::KMINUS] / nparts[PID::PIMINUS], 5);
}
if (!isZero(nparts[PID::KPLUS])) {
_p["Karatio"]->fill(_sqsEdge, nparts[PID::KMINUS] / nparts[PID::KPLUS], 5);
}
if (!isZero(nparts[PID::PROTON])) {
_p["Pratio"]->fill(_sqsEdge, nparts[PID::ANTIPROTON] / nparts[PID::PROTON], 5);
}
/// Sum the weight of the event
if (c < 80) {
_c[_sqs + "dpT_Pi" + centbin]->fill();
_c[_sqs + "dpT_Piplus" + centbin]->fill();
_c[_sqs + "dpT_Kaon" + centbin]->fill();
_c[_sqs + "dpT_Kaonplus" + centbin]->fill();
_c[_sqs + "dpT_Proton" + centbin]->fill();
_c[_sqs + "dpT_AntiProton" + centbin]->fill();
_c[_sqs + "npart_PiPlus"]->fill();
_c[_sqs + "npart_PiMinus"]->fill();
_c[_sqs + "npart_KaPlus"]->fill();
_c[_sqs + "npart_KaMinus"]->fill();
_c[_sqs + "npart_Proton"]->fill();
_c[_sqs + "npart_AntiProton"]->fill();
}
}
/// Normalise histograms etc., after the run
void finalize() {
/// Normalisation
for (const auto& item : _c) {
if (!item.second->sumW()) continue;
scale(_h[item.first], 1.0 / (TWOPI * 0.2 * item.second->sumW()));
}
}
private:
/// @name Histograms
/// @{
map<string, Profile1DPtr> _p;
map<string, Histo1DPtr> _h;
map<string, CounterPtr> _c;
/// @}
/// @name Variables
/// @{
string _sqs = "", _sqs25 = "";
double _sqsEdge;
/// @}
/// @name Bin edges
/// @{
YODA::Axis<double> centAxis{5., 10., 20., 30., 40., 50., 60., 70., 80.};
/// @}
};
RIVET_DECLARE_PLUGIN(STAR_2017_I1510593);
}