1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/FinalState.hh"
 4#include "Rivet/Projections/ChargedFinalState.hh"
 5#include "Rivet/Projections/SingleValueProjection.hh"
 6#include "Rivet/Projections/ImpactParameterProjection.hh"
 7#include "Rivet/Tools/Percentile.hh"
 8#include "Rivet/Analyses/RHICCommon.hh"
 9
10namespace Rivet {
11
12
13  /// @brief Third harmonic of azimuthal correlations in Au+Au collisions at different COM energies
14  class STAR_2016_I1414638 : public Analysis {
15  public:
16
17    /// Constructor
18    STAR_2016_I1414638 ()
19      : Analysis("STAR_2016_I1414638") {}
20
21
22    void init() {
23      /// Projections
24      declareCentrality(STAR_BES_Centrality(), "STAR_BES_CALIB", "CMULT", "CMULT");
25      declare(ChargedFinalState(Cuts::abseta < 1.0 && Cuts::pT > 0.2*GeV), "CFS");
26
27      // Histograms
28      int energy = -1;
29      for (int i = 0, N = energies.size(); i < N; ++i) {
30        if (isCompatibleWithSqrtS(197.*energies[i]*GeV, 1e-1)) energy = i;
31      }
32      if (energy == -1) MSG_ERROR("Incompatible beam energy!");
33      for (int i = 0; i < 9; ++i) book(h_v32[centralityBins[i]], 1 + i + 9 * energy, 1, 1);
34    }
35
36
37    /// Perform the per-event analysis
38    void analyze(const Event& event) {
39      const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
40      // Require at least two charged particles for the analysis to make sense.
41      // No further triggers are described in the paper.
42      const Particles& particles = cfs.particles();
43      if (particles.size() < 2) return;
44      // The centrality projection
45      const CentralityProjection& cent = apply<CentralityProjection>(event,"CMULT");
46      const double c = cent();
47      // Find the correct histogram to fill.
48      auto hItr = h_v32.upper_bound(c);
49      if (hItr == h_v32.end()) return;
50      for (int i = 0, N = particles.size(); i < N; ++i){
51        for (int j = i + 1; j < N; ++j) {
52          const double eta1 = particles[i].eta();
53          const double eta2 = particles[j].eta();
54          if (eta1 * eta2 < 0) {
55            const double deltaPhi = abs(particles[i].phi() - particles[j].phi());
56            // Fill profile with v_2(2)^2 from eq. (1) in the paper.
57            hItr->second->fill(abs(eta1 - eta2), cos(3.*deltaPhi));
58          }
59        }
60      }
61    }
62
63    /// Normalise histograms etc., after the run
64    // void finalize() {}
65
66    /// @}
67
68
69    /// @name Bin edges
70    /// @{
71    /// The centrality bins
72    const doubles centralityBins = {5., 10., 20., 30., 40., 50., 60., 70., 80.};
73    /// The analysis energies
74    const doubles energies = {7.7, 11.5, 14.5, 19.6, 27.0, 39.0, 62.4, 200.0};
75    /// @}
76
77
78    /// The histograms
79    map<double, Profile1DPtr> h_v32;
80
81  };
82
83
84  RIVET_DECLARE_PLUGIN(STAR_2016_I1414638);
85
86}