1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/Beam.hh"
  7#include "Rivet/Projections/VetoedFinalState.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// Jet and underlying event properties as a function of particle multiplicity
 13  class CMS_2013_I1261026 : public Analysis {
 14  public:
 15
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1261026);
 17
 18    void init() {
 19      const ChargedFinalState cfs(Cuts::abseta < 2.4 && Cuts::pT > 0.25*GeV);
 20      declare(cfs, "CFS250");
 21
 22      FastJets jetpro(cfs, JetAlg::ANTIKT, 0.5);
 23      declare(jetpro, "Jets");
 24
 25      // For min bias trigger
 26      const ChargedFinalState cfsBSCplus(Cuts::etaIn(3.23, 4.65) && Cuts::pT > 500*MeV);
 27      declare(cfsBSCplus, "cfsBSCplus");
 28
 29      const ChargedFinalState cfsBSCminus(Cuts::etaIn(-4.65, -3.23) && Cuts::pT > 500*MeV);
 30      declare(cfsBSCminus, "cfsBSCminus");
 31
 32      // Histograms:
 33      book(_h_AllTrkMeanPt            ,1, 1, 1);
 34      book(_h_SoftTrkMeanPt           ,2, 1, 1);
 35      book(_h_IntrajetTrkMeanPt       ,3, 1, 1);
 36      book(_h_IntrajetLeaderTrkMeanPt ,4, 1, 1);
 37      book(_h_MeanJetPt               ,5, 1, 1);
 38      book(_h_JetRate5GeV             ,6, 1, 1);
 39      book(_h_JetRate30GeV            ,7, 1, 1);
 40
 41      for (int ihist = 0; ihist < 5; ++ihist) {
 42        book(_h_JetSpectrum[ihist] ,ihist+8, 1, 1);
 43        book(_h_JetStruct[ihist]   ,ihist+13, 1, 1);
 44      }
 45    }
 46
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      // MinBias trigger
 51      const ChargedFinalState& cfsBSCplus = apply<ChargedFinalState>(event, "cfsBSCplus");
 52      if (cfsBSCplus.empty()) vetoEvent;
 53      const ChargedFinalState& cfsBSCminus = apply<ChargedFinalState>(event, "cfsBSCminus");
 54      if (cfsBSCminus.empty()) vetoEvent;
 55
 56      const ChargedFinalState& cfsp = apply<ChargedFinalState>(event, "CFS250");
 57      if (cfsp.empty()) vetoEvent;
 58
 59      const FastJets& jetpro = apply<FastJets>(event, "Jets");
 60      const Jets& jets = jetpro.jetsByPt(Cuts::pT > 5*GeV);
 61
 62      const int mult = cfsp.size();
 63
 64      int multbin[6] = { 10, 30, 50, 80, 110, 140 };
 65      for (int ibin = 0; ibin < 5; ++ibin) {
 66        if (mult > multbin[ibin] && mult <= multbin[ibin + 1]) {
 67          eventDecomp(event, mult, ibin);
 68          unsigned int jetCounter5GeV(0), jetCounter30GeV(0);
 69          for (size_t ijets = 0; ijets < jets.size(); ++ijets) {
 70            if (jets[ijets].abseta() < 1.9) {
 71              _h_JetSpectrum[ibin]->fill(jets[ijets].pT()/GeV);
 72              _h_MeanJetPt->fill(mult, jets[ijets].pT()/GeV);
 73              if (jets[ijets].pT() > 5*GeV)   ++jetCounter5GeV;
 74              if (jets[ijets].pT() > 30*GeV)  ++jetCounter30GeV;
 75            }
 76          }
 77          _h_JetRate5GeV->fill( mult,  jetCounter5GeV);
 78          _h_JetRate30GeV->fill(mult, jetCounter30GeV);
 79        }
 80      }
 81    }
 82
 83
 84    /// Normalise histograms etc., after the run
 85    void finalize() {
 86      for (size_t i = 0; i < 5; ++i) {
 87        normalize(_h_JetSpectrum[i],  4.0);
 88        normalize(_h_JetStruct[i]  , 0.08);
 89      }
 90
 91    }
 92
 93    void eventDecomp(const Event& event, int mult, size_t ibin) {
 94
 95      struct TrkInJet { double pt; double eta; double phi; double R; };
 96      TrkInJet jetConstituents[100][100]; //1-st index - the number of the jet, 2-nd index - track in the jet
 97      TrkInJet jetsEv[100];
 98      size_t j[100];
 99      size_t jCount = 0;
100
101      for (size_t i = 0; i < 100; ++i) {
102        j[i] = 0;
103        jetsEv[i].pt = 0;
104        jetsEv[i].eta = 0;
105        jetsEv[i].phi = 0;
106        for (size_t k = 0; k < 100; ++k) {
107          jetConstituents[i][k].pt = 0;
108          jetConstituents[i][k].phi = 0;
109          jetConstituents[i][k].eta = 0;
110          jetConstituents[i][k].R = 0;
111        }
112      }
113
114      const FastJets& jetpro = apply<FastJets>(event, "Jets");
115      const Jets& jets = jetpro.jetsByPt(Cuts::pT > 5*GeV);
116
117      // Start event decomp
118
119      for (size_t ijets = 0; ijets < jets.size(); ++ijets) {
120        jetsEv[ijets].pt = jets[ijets].pT();
121        jetsEv[ijets].eta = jets[ijets].eta();
122        jetsEv[ijets].phi = jets[ijets].phi();
123        jCount += 1;
124      }
125
126      const ChargedFinalState& cfsp = apply<ChargedFinalState>(event, "CFS250");
127      for (const Particle& p : cfsp.particles()) {
128        _h_AllTrkMeanPt->fill(mult, p.pT()/GeV);
129        int flag = 0;
130        for (size_t i = 0; i < jCount; ++i) {
131          const double delta_phi = deltaPhi(jetsEv[i].phi, p.phi());
132          const double delta_eta = jetsEv[i].eta - p.eta();
133          const double R = sqrt(delta_phi * delta_phi + delta_eta * delta_eta);
134          if (R <= 0.5) {
135            flag++;
136            jetConstituents[i][j[i]].pt = p.pT();
137            jetConstituents[i][j[i]].R = R;
138            j[i]++;
139          }
140        }
141        if (flag == 0) _h_SoftTrkMeanPt->fill(mult, p.pT()/GeV);
142      }
143
144      for (size_t i = 0; i < jCount; ++i) {
145        double ptInjetLeader = 0;
146        if (!inRange(jetsEv[i].eta, -1.9, 1.9)) continue; // only fully reconstructed jets for internal jet studies
147        for (size_t k = 0; k < j[i]; ++k) {
148          _h_IntrajetTrkMeanPt->fill(mult, jetConstituents[i][k].pt);
149          _h_JetStruct[ibin]->fill(jetConstituents[i][k].R, jetConstituents[i][k].pt/jetsEv[i].pt);
150          if (ptInjetLeader < jetConstituents[i][k].pt) ptInjetLeader = jetConstituents[i][k].pt;
151        }
152        if (ptInjetLeader != 0) _h_IntrajetLeaderTrkMeanPt->fill(mult, ptInjetLeader);
153      }
154
155    }
156
157
158  private:
159
160    Profile1DPtr _h_AllTrkMeanPt, _h_SoftTrkMeanPt;
161    Profile1DPtr _h_IntrajetTrkMeanPt, _h_IntrajetLeaderTrkMeanPt;
162    Profile1DPtr _h_MeanJetPt;
163    Profile1DPtr _h_JetRate5GeV, _h_JetRate30GeV;
164
165    array<Histo1DPtr,5> _h_JetSpectrum;
166    array<Histo1DPtr,5> _h_JetStruct;
167
168  };
169
170  RIVET_DECLARE_PLUGIN(CMS_2013_I1261026);
171
172}