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