1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/CentralityProjection.hh"
  4#include "Rivet/Projections/HepMCHeavyIon.hh"
  5#include "Rivet/Analyses/AliceCommon.hh"
  6#include "Rivet/Tools/Cuts.hh"
  7
  8namespace Rivet {
  9
 10
 11  class ALICE_2014_I1243865 : public Analysis {
 12
 13  // @brief Multi-strange baryon production at mid-rapidity in 2.76 TeV Pb--Pb collisions
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2014_I1243865);
 18
 19
 20    // Book histograms and projections etc.
 21    void init() {
 22      // Particles of interest.
 23      declare(ALICE::PrimaryParticles(Cuts::absrap < 0.5),"CFS");
 24
 25      // The event trigger.
 26      declare(ALICE::V0AndTrigger(), "V0-AND");
 27
 28      // The centrality projection.
 29      declareCentrality(ALICE::V0MMultiplicity(),
 30           "ALICE_2015_PBPBCentrality", "V0M", "V0M");
 31
 32      // Access the HepMC heavy ion info
 33      declare(HepMCHeavyIon(), "HepMC");
 34
 35      // Xi Baryons.
 36      size_t ixi = 0;
 37      _histPtXi.resize(5);
 38      _histPtXi_bar.resize(5);
 39      for (string str : {"d01-","d02-","d03-","d04-","d05-"}){
 40        book(_histPtXi[ixi], str+"x01-y01");
 41        book(_histPtXi_bar[ixi], str+"x01-y02");
 42        ixi += 1;
 43      }
 44
 45      // Omega Baryons.
 46      size_t iom = 0;
 47      _histPtOmega.resize(5);
 48      _histPtOmega_bar.resize(5);
 49      for (string str : {"d06-","d07-","d08-","d09-","d10-"}){
 50        book(_histPtOmega[iom], str+"x01-y01");
 51        book(_histPtOmega_bar[iom], str+"x01-y02");
 52        iom += 1;
 53      }
 54
 55      // Sum of weights for the centrality intervals.
 56      sow.resize(_histPtOmega.size());
 57      for (int i = 0, N = _histPtOmega.size(); i < N; ++i) {
 58        book(sow[i], "sow_" + toString(i));
 59      }
 60
 61      book(_histXitoPi, "d14-x01-y01");
 62      book(_histOmegatoPi, "d14-x01-y02");
 63    }
 64
 65
 66    void analyze(const Event& event) {
 67      // Event trigger.
 68      if (!apply<ALICE::V0AndTrigger>(event, "V0-AND")() ) vetoEvent;
 69
 70      // Centrality.
 71      const CentralityProjection& cent = apply<CentralityProjection>(event,"V0M");
 72      const double c = cent();
 73
 74      int centralityclass = -1;
 75      if(c > 0. && c <= 10) centralityclass = 0;
 76      if(c > 10. && c <= 20) centralityclass = 1;
 77      if(c > 20. && c <= 40) centralityclass = 2;
 78      if(c > 40. && c <= 60) centralityclass = 3;
 79      if(c > 60. && c <= 80) centralityclass = 4;
 80      if (centralityclass == -1) vetoEvent;
 81      // Fill sum of weights
 82      sow[centralityclass]->fill();
 83      int nPions = 0;
 84      int nXi = 0;
 85      int nOmega = 0;
 86      for (const Particle& p : apply<ALICE::PrimaryParticles>(event,"CFS").particles()) {
 87        const double pT = p.pT() / GeV;
 88        switch (p.pid()){
 89	  case 211:
 90	    nPions++;
 91	    break;
 92           case 3312:
 93	     _histPtXi[centralityclass]->fill(pT);
 94	     nXi++;
 95	     break;
 96	   case -3312:
 97	     _histPtXi_bar[centralityclass]->fill(pT);
 98	     nXi++;
 99	     break;
100	   case 3334:
101	     _histPtOmega[centralityclass]->fill(pT);
102	     nOmega++;
103	   break;
104	     case -3334:
105	     _histPtOmega_bar[centralityclass]->fill(pT);
106	     nOmega++;
107	   break;
108         }
109       }
110       // Extract Npart form GenEvent. TODO: Unclear how to do
111       // this in HepMC3
112      const HepMCHeavyIon & hi = apply<HepMCHeavyIon>(event, "HepMC");
113       if ( nPions != 0){
114         const double npart = hi.Npart_proj() + hi.Npart_targ();
115         if (nXi != 0)
116           _histXitoPi->fill(npart, double(nXi) / double(nPions));
117         if (nOmega != 0)
118           _histOmegatoPi->fill(npart, double(nOmega) / double(nPions));
119        }
120     }
121
122    void finalize() {
123      for (int i = 0, N = _histPtOmega.size(); i < N; ++i) {
124        const double s = 1./sow[i]->sumW();
125        _histPtXi[i]->scaleW(s);
126        _histPtXi_bar[i]->scaleW(s);
127        _histPtOmega[i]->scaleW(s);
128        _histPtOmega_bar[i]->scaleW(s);
129      }
130   }
131
132
133  private:
134
135    vector<Histo1DPtr> _histPtXi;
136    vector<Histo1DPtr> _histPtXi_bar;
137    vector<Histo1DPtr> _histPtOmega;
138    vector<Histo1DPtr> _histPtOmega_bar;
139    vector<CounterPtr> sow;
140    Profile1DPtr _histXitoPi;
141    Profile1DPtr _histOmegatoPi;
142
143  };
144
145
146  RIVET_DECLARE_PLUGIN(ALICE_2014_I1243865);
147
148}