1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5using namespace std;
  6
  7namespace Rivet {
  8
  9
 10  /// Measurement of the NSD charged particle multiplicity at 0.9, 2.36, and 7 TeV
 11  class CMS_2011_I879315 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2011_I879315);
 15
 16
 17    void init() {
 18      ChargedFinalState cfs((Cuts::etaIn(-2.4, 2.4)));
 19      declare(cfs, "CFS");
 20
 21      // eta bins
 22      _etabins.push_back(0.5);
 23      _etabins.push_back(1.0);
 24      _etabins.push_back(1.5);
 25      _etabins.push_back(2.0);
 26      _etabins.push_back(2.4) ;
 27
 28      if (isCompatibleWithSqrtS(900*GeV)) {
 29        for (size_t ietabin=0; ietabin < _etabins.size(); ietabin++) {
 30          _h_dNch_dn.push_back( Histo1DPtr() );
 31          book( _h_dNch_dn.back(), 2 + ietabin, 1, 1);
 32        }
 33        book(_h_dNch_dn_pt500_eta24 ,20, 1, 1);
 34        book(_h_dmpt_dNch_eta24 ,23, 1, 1);
 35      }
 36
 37      if (isCompatibleWithSqrtS(2360*GeV)) {
 38        for (size_t ietabin=0; ietabin < _etabins.size(); ietabin++) {
 39          _h_dNch_dn.push_back( Histo1DPtr() );
 40          book(_h_dNch_dn.back(), 7 + ietabin, 1, 1);
 41        }
 42        book(_h_dNch_dn_pt500_eta24 ,21, 1, 1);
 43        book(_h_dmpt_dNch_eta24 ,24, 1, 1);
 44      }
 45
 46      if (isCompatibleWithSqrtS(7000*GeV)) {
 47        for (size_t ietabin=0; ietabin < _etabins.size(); ietabin++) {
 48          _h_dNch_dn.push_back( Histo1DPtr() );
 49          book(_h_dNch_dn.back(), 12 + ietabin, 1, 1);
 50        }
 51        book(_h_dNch_dn_pt500_eta24 ,22, 1, 1);
 52        book(_h_dmpt_dNch_eta24 ,25, 1, 1);
 53      }
 54    }
 55
 56
 57    void analyze(const Event& event) {
 58
 59      // Get the charged particles
 60      const ChargedFinalState& charged = apply<ChargedFinalState>(event, "CFS");
 61
 62      // Resetting the multiplicity for the event to 0;
 63      vector<int> _nch_in_Evt;
 64      vector<int> _nch_in_Evt_pt500;
 65      _nch_in_Evt.assign(_etabins.size(), 0);
 66      _nch_in_Evt_pt500.assign(_etabins.size(), 0);
 67      double sumpt = 0;
 68
 69      // Loop over particles in event
 70      for (const Particle& p : charged.particles()) {
 71        // Selecting only charged hadrons
 72        if (! PID::isHadron(p.pid())) continue;
 73
 74        double pT = p.pT();
 75        double eta = p.eta();
 76        sumpt += pT;
 77        for (size_t ietabin = _etabins.size(); ietabin > 0; --ietabin) {
 78          if (fabs(eta) > _etabins[ietabin-1]) break;
 79          ++_nch_in_Evt[ietabin-1];
 80          if (pT > 0.5/GeV) ++_nch_in_Evt_pt500[ietabin-1];
 81        }
 82      }
 83
 84      // Filling multiplicity-dependent histogramms
 85      for (size_t ietabin = 0; ietabin < _etabins.size(); ietabin++) {
 86        _h_dNch_dn[ietabin]->fill(_nch_in_Evt[ietabin]);
 87      }
 88
 89      // Do only if eta bins are the needed ones
 90      if (_etabins[4] == 2.4 && _etabins[0] == 0.5) {
 91        if (_nch_in_Evt[4] != 0) {
 92          _h_dmpt_dNch_eta24->fill(_nch_in_Evt[4], sumpt/GeV / _nch_in_Evt[4]);
 93        }
 94        _h_dNch_dn_pt500_eta24->fill(_nch_in_Evt_pt500[4]);
 95      } else {
 96        MSG_WARNING("You changed the number of eta bins, but forgot to propagate it everywhere !!");
 97      }
 98    }
 99
100
101    void finalize() {
102      for (size_t ietabin = 0; ietabin < _etabins.size(); ietabin++){
103        normalize(_h_dNch_dn[ietabin]);
104      }
105      normalize(_h_dNch_dn_pt500_eta24);
106    }
107
108
109  private:
110
111    /// @{
112    vector<Histo1DPtr> _h_dNch_dn;
113    Histo1DPtr _h_dNch_dn_pt500_eta24;
114    Profile1DPtr _h_dmpt_dNch_eta24;
115    /// @}
116
117    vector<double> _etabins;
118
119  };
120
121
122
123  RIVET_DECLARE_ALIASED_PLUGIN(CMS_2011_I879315, CMS_2011_S8884919);
124
125}