1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/CentralityProjection.hh"
  5#include "Rivet/Projections/AliceCommon.hh"
  6#include "Rivet/Tools/AliceCommon.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// @brief Strangeness enhancement in pp 7 TeV by ALICE.
 12  class ALICE_2016_I1471838 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2016_I1471838);
 17
 18    int profileIndex(vector<double> cBins, double c) {
 19      int index = 100;
 20      if (c > 0 && c <= cBins[0]) return cBins.size() - 1;
 21      for (size_t i = 0; i < cBins.size() - 1; ++i) {
 22        if (c > cBins[i] && c <= cBins[i + 1]) {
 23	  index = i;
 24	  break;
 25	}
 26      }
 27      return max(0, int(cBins.size() - index - 2));
 28    }
 29
 30    /// Book histograms and initialise projections before the run
 31    void init() {
 32      // Centrality projection.
 33      declareCentrality(ALICE::V0MMultiplicity(),
 34        "ALICE_2015_PPCentrality","V0M","V0M");
 35      // Central primary particles
 36      declare(ChargedFinalState(Cuts::abseta < 1.0),"PP");
 37      declare(ALICE::PrimaryParticles(Cuts::absrap < 0.5),"PPy");
 38      centralityBins = {1.,5.,10.,15.,20., 30., 40., 50., 70., 100.};
 39      centralityBinsOmega = {5.,15.,30.,50.,100.};
 40      // Book histograms
 41      for (int i = 0; i < 10; ++i) {
 42        book(K0SpT[centralityBins[i]], i+1,1,1);
 43        book(LambdapT[centralityBins[i]], i+11,1,1);
 44        book(XipT[centralityBins[i]], i+21,1,1);
 45        book(sow[centralityBins[i]], "sow_" + toString(i));
 46      }
 47      for (int i = 0; i < 5; ++i) {
 48        book(OmegapT[centralityBinsOmega[i]], i+31,1,1);
 49        book(sowOmega[centralityBinsOmega[i]], "sowO_" + toString(i));
 50      }
 51      book(piYield, 40,1,1);
 52      book(pYield, 41,1,1);
 53      book(kYield, 42,1,1);
 54      book(lambdaYield, 43,1,1);
 55      book(xiYield, 44,1,1);
 56      book(omegaYield, 45,1,1);
 57      book(piRebinned, "/piRebinned", refData(45,1,1));
 58
 59      // Make the ratios
 60      book(kpi, 36, 1, 1, true);
 61      book(ppi, 47, 1, 1, true);
 62      book(lpi, 37, 1, 1, true);
 63      book(xpi, 38, 1, 1, true);
 64      book(opi, 39, 1, 1, true);
 65      book(lk, 46, 1, 1, true);
 66    }
 67
 68
 69    /// Perform the per-event analysis
 70    void analyze(const Event& event) {
 71      if (apply<ChargedFinalState>(event,"PP").particles().size() < 1) vetoEvent;
 72      const ALICE::PrimaryParticles& prim = apply<ALICE::PrimaryParticles>(event,"PPy");
 73      const CentralityProjection& cent = apply<CentralityProjection>(event,"V0M");
 74      double c  = cent();
 75      // Find the correct histograms
 76      auto kptItr = K0SpT.upper_bound(c);
 77      if (kptItr == K0SpT.end()) return;
 78      auto lptItr = LambdapT.upper_bound(c);
 79      if (lptItr == LambdapT.end()) return;
 80      auto xptItr = XipT.upper_bound(c);
 81      if (xptItr == XipT.end()) return;
 82      auto optItr = OmegapT.upper_bound(c);
 83      if (optItr == OmegapT.end()) return;
 84      // Fill the sow.
 85      auto sowItr = sow.upper_bound(c);
 86      if (sowItr == sow.end()) return;
 87      auto sowOmegaItr = sowOmega.upper_bound(c);
 88      if (sowOmegaItr == sowOmega.end()) return;
 89      sowItr->second->fill();
 90      sowOmegaItr->second->fill();
 91      // Fill the pt histograms and count yields.
 92      int npi = 0, npr = 0, nk = 0;
 93      int nla = 0, nxi = 0, nom = 0;
 94      for (auto p : prim.particles()) {
 95        const double pT = p.pT();
 96	const int pid = abs(p.pid());
 97	if (pid == 211) ++npi;
 98	else if (pid == 2212) ++npr;
 99	else if (pid == 310) {
100	  kptItr->second->fill(pT);
101	  ++nk;
102	}
103	else if (pid == 3122) {
104	  lptItr->second->fill(pT);
105	  ++nla;
106	}
107	else if (pid == 3312) {
108	  xptItr->second->fill(pT);
109	  ++nxi;
110	}
111	else if (pid == 3334) {
112	  optItr->second->fill(pT);
113	  ++nom;
114	}
115      }
116      // Fill the profiles of yields.
117      int index = profileIndex(centralityBins,c);
118      piYield->fillBin(index, double(npi));
119      pYield->fillBin(index, double(npr));
120      kYield->fillBin(index, double(nk));
121      lambdaYield->fillBin(index, double(nla));
122      xiYield->fillBin(index, double(nxi));
123      index = profileIndex(centralityBinsOmega, c);
124      omegaYield->fillBin(index, double(nom));
125      piRebinned->fillBin(index,double(npi));
126    }
127
128
129    /// Normalise histograms etc., after the run
130    void finalize() {
131      // Normalize the spectra
132      for (int i = 0; i < 10; ++i) {
133        K0SpT[centralityBins[i]]->scaleW(1./sow[centralityBins[i]]->sumW());
134        XipT[centralityBins[i]]->scaleW(1./sow[centralityBins[i]]->sumW());
135        LambdapT[centralityBins[i]]->scaleW(1./sow[centralityBins[i]]->sumW());
136      }
137      for (int i = 0; i < 5; ++i) {
138        OmegapT[centralityBinsOmega[i]]->scaleW(1./sowOmega[centralityBinsOmega[i]]->sumW());
139      }
140
141      divide(kYield, piYield, kpi);
142      kpi->scaleY(2.);
143      divide(pYield, piYield, ppi);
144      divide(lambdaYield, piYield, lpi);
145      divide(xiYield, piYield, xpi);
146      divide(omegaYield, piRebinned, opi);
147      divide(lambdaYield, kYield, lk);
148      lk->scaleY(0.5);
149    }
150
151    //@}
152
153
154    /// @name Histograms
155    //@{
156    // Histograms ordered in centrality classes
157    vector<double> centralityBins;
158    vector<double> centralityBinsOmega;
159
160    // pT spectra
161    map<double, Histo1DPtr> K0SpT;
162    map<double, Histo1DPtr> LambdapT;
163    map<double, Histo1DPtr> XipT;
164    map<double, Histo1DPtr> OmegapT;
165    map<double, CounterPtr> sow;
166    map<double, CounterPtr> sowOmega;
167
168    // Total yields
169    Profile1DPtr piYield;
170    Profile1DPtr pYield;
171    Profile1DPtr kYield;
172    Profile1DPtr lambdaYield;
173    Profile1DPtr xiYield;
174    Profile1DPtr omegaYield;
175    Profile1DPtr piRebinned;
176
177    // Ratios
178    Scatter2DPtr kpi;
179    Scatter2DPtr ppi;
180    Scatter2DPtr lpi;
181    Scatter2DPtr xpi;
182    Scatter2DPtr opi;
183    Scatter2DPtr lk;
184    //@}
185  };
186
187
188  // The hook for the plugin system
189  RIVET_DECLARE_PLUGIN(ALICE_2016_I1471838);
190
191
192}