Rivet analyses

title: STAR_2017_I1510593

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, $\langle p_\mathrm{T} \rangle$ spectra and particles ratios are shown.

Source code:STAR_2017_I1510593.cc

```c++ // -- 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); } ```