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