1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/SingleValueProjection.hh"
  4#include "Rivet/Projections/ImpactParameterProjection.hh"
  5#include "Rivet/Projections/FinalState.hh"
  6#include "Rivet/Projections/UnstableParticles.hh"
  7#include "Rivet/Projections/ChargedFinalState.hh"
  8#include "Rivet/Analyses/RHICCommon.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief Brahms pT spectra for id particles (pi+, pi-, K+, K-)
 14  //  in small bins of rapidity, 5% central collisions.
 15  //  System: AuAu @ 200GeV/nn.
 16  class BRAHMS_2004_I647076 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(BRAHMS_2004_I647076);
 21
 22
 23    /// @name Analysis methods
 24    /// @{
 25
 26    /// Book histograms and initialise projections before the run
 27    void init() {
 28
 29      // Initialise and register projections
 30      // Centrality Projection.
 31      declareCentrality(BRAHMSCentrality(), "BRAHMS_2004_AUAUCentrality","mult","BCEN");
 32      // TODO: Feed down correction is unclear.
 33      declare(FinalState(Cuts::rap < 4 && Cuts::rap > -0.1 && Cuts::pT > 100*MeV), "FS");
 34      // The measured rapidity intervals for pions.
 35      rapIntervalsPi = {{-0.1,0.},{0.,0.1},{0.4,0.6},{0.6,0.8},{0.8,1.0},
 36        {1.0,1.2},{1.2,1.4},{2.1,2.3},{2.4,2.6},{3.0,3.1},{3.1,3.2},{3.2,3.3},
 37        {3.3,3.4},{3.4,3.66}};
 38      // The measured rapidity intervals for kaons.
 39      rapIntervalsK = {{-0.1,0.},{0.,0.1},{0.4,0.6},{0.6,0.8},{0.8,1.0},
 40        {1.0,1.2},{2.0,2.2},{2.3,2.5},{2.9,3.0},{3.0,3.1},{3.1,3.2},{3.2,3.4}};
 41      // Book histograms
 42      piPlus.resize(rapIntervalsPi.size());
 43      piMinus.resize(rapIntervalsPi.size());
 44      for (int i = 0, N = rapIntervalsPi.size(); i < N; ++i) {
 45        book(piPlus[i], 1, 1, 1 + i);
 46        book(piMinus[i], 1, 1, 15 + i);
 47      }
 48
 49      kPlus.resize(rapIntervalsK.size());
 50      kMinus.resize(rapIntervalsK.size());
 51      for (int i = 0, N = rapIntervalsK.size(); i < N; ++i) {
 52        book(kPlus[i], 2, 1, 1 + i);
 53        book(kMinus[i], 2, 1, 13 + i);
 54      }
 55      // Counter for accepted sum of weights (centrality cut).
 56      book(centSow, "centSow");
 57
 58    }
 59
 60
 61    /// Perform the per-event analysis
 62    void analyze(const Event& event) {
 63      // Reject all non-central events. The paper does not speak of
 64      // any other event trigger, which in any case should matter
 65      // little for central events.
 66      if(apply<CentralityProjection>(event,"BCEN")() > 5.0) return;
 67      // Keep track of sum of weights.
 68      centSow->fill();
 69      const FinalState& fs = apply<FinalState>(event,"FS");
 70      // Loop over particles.
 71      for (const auto& p : fs.particles()) {
 72        const double y = p.rapidity();
 73	const double pT = p.pT();
 74	const int id = p.pid();
 75	// First pions.
 76	if (abs(id) == 211) {
 77          // Protect against decaying K0S and Lambda
 78	  if (p.hasAncestorWith(Cuts::pid == 310) || p.hasAncestorWith(Cuts::pid == -310) ||
 79	    p.hasAncestorWith(Cuts::pid == 3122) || p.hasAncestorWith(Cuts::pid == 3122)) continue;
 80	  for (int i = 0, N = rapIntervalsPi.size(); i < N; ++i) {
 81	    if (y > rapIntervalsPi[i].first && y <= rapIntervalsPi[i].second) {
 82	      const double dy = rapIntervalsPi[i].second - rapIntervalsPi[i].first;
 83	      const double nWeight = 1.0 / ( 2.*M_PI*pT*dy);
 84	      if (id == 211) piPlus[i]->fill(pT, nWeight);
 85	      else piMinus[i]->fill(pT, nWeight);
 86	      break;
 87	    }
 88	  }
 89	}
 90	// Then kaons.
 91	else if (abs(id) == 321) {
 92	  for (int i = 0, N = rapIntervalsK.size(); i < N; ++i) {
 93	    if (y > rapIntervalsK[i].first && y <= rapIntervalsK[i].second) {
 94	      const double dy = rapIntervalsK[i].second - rapIntervalsK[i].first;
 95	      const double nWeight = 1.0 / ( 2.*M_PI*pT*dy);
 96	      if (id == 321) kPlus[i]->fill(pT, nWeight);
 97	      else kMinus[i]->fill(pT, nWeight);
 98	      break;
 99	    }
100	  }
101	}
102      }
103    }
104
105
106    /// Normalise histograms etc., after the run
107    void finalize() {
108      // Normalize all histograms to per-event yields.
109      for (int i = 0, N = rapIntervalsPi.size(); i < N; ++i) {
110        piPlus[i]->scaleW(1./centSow->sumW());
111        piMinus[i]->scaleW(1./centSow->sumW());
112      }
113      for (int i = 0, N = rapIntervalsK.size(); i < N; ++i) {
114        kPlus[i]->scaleW(1./centSow->sumW());
115        kMinus[i]->scaleW(1./centSow->sumW());
116      }
117
118    }
119
120    /// @}
121
122    // The rapidity intervals.
123    vector<pair<double, double> > rapIntervalsPi;
124    vector<pair<double, double> > rapIntervalsK;
125
126    /// @name Histograms
127    /// @{
128    vector<Histo1DPtr> piPlus;
129    vector<Histo1DPtr> piMinus;
130    vector<Histo1DPtr> kPlus;
131    vector<Histo1DPtr> kMinus;
132    CounterPtr centSow;
133    /// @}
134  };
135
136  RIVET_DECLARE_PLUGIN(BRAHMS_2004_I647076);
137}