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  /// Ratio of energy in -6.6 < eta < -5.2 for events with a charged-particle jet
 13  class CMS_2013_I1218372 : public Analysis {
 14  public:
 15
 16    /// Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1218372);
 18
 19
 20    void init() {
 21
 22      // gives the range of eta and min pT for the final state from which I get the jets
 23      FastJets jetpro (ChargedFinalState(Cuts::abseta < 2.5 && Cuts::pT >= 0.3*GeV), JetAlg::ANTIKT, 0.5);
 24      declare(jetpro, "Jets");
 25
 26      // skip Neutrinos and Muons
 27      VetoedFinalState fsv(FinalState(Cuts::etaIn(-7.0, -4.0)));
 28      fsv.vetoNeutrinos();
 29      fsv.addVetoPairId(PID::MUON);
 30      declare(fsv, "fsv");
 31
 32      // for the hadron level selection
 33      VetoedFinalState sfsv;
 34      sfsv.vetoNeutrinos();
 35      sfsv.addVetoPairId(PID::MUON);
 36      declare(sfsv, "sfsv");
 37
 38      size_t id = 0;
 39      for (double eVal : allowedEnergies()) {
 40        const string en = toString(int(eVal));
 41        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
 42        ++id;
 43
 44        //counters
 45        book(_c[en+"pass"], "_pass"+en);
 46        book(_c[en+"incl"], "_incl"+en);
 47
 48        // Temporary histograms to fill the energy flow for leading jet events.
 49        // Ratios are calculated in finalize().
 50        book(_e[en+"ratio"], id, 1, 1);
 51        book(_h[en+"num"], "TMP/eflow_jet"+en,  refData(id, 1, 1)); // Leading jet energy flow in pt
 52        book(_h[en+"den"], "TMP/number_jet"+en, refData(id, 1, 1)); // Number of events in pt
 53      }
 54      if (_sqs == "" && !merging()) {
 55        throw BeamError("Invalid beam energy for " + name() + "\n");
 56      }
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62
 63      // Skip if the event is empty
 64      const FinalState& fsv = apply<FinalState>(event, "fsv");
 65      if (fsv.empty()) vetoEvent;
 66
 67      // ====================== Minimum Bias selection
 68
 69      const FinalState& sfsv = apply<FinalState>(event, "sfsv");
 70      Particles parts = sfsv.particles(cmpMomByRap);
 71      if (parts.empty()) vetoEvent;
 72
 73      // find dymax
 74      double dymax = 0;
 75      int gap_pos  = -1;
 76      for (size_t i = 0; i < parts.size()-1; ++i) {
 77        double dy = parts[i+1].rap() - parts[i].rap();
 78        if (dy > dymax) {
 79          dymax = dy;
 80          gap_pos = i;
 81        }
 82      }
 83
 84      // calculate mx2 and my2
 85      FourMomentum xmom;
 86      for (int i=0; i<=gap_pos; ++i) {
 87        xmom += parts[i].mom();
 88      }
 89      double mx2 = xmom.mass2();
 90      if (mx2<0) vetoEvent;
 91
 92      FourMomentum ymom;
 93      for (size_t i=gap_pos+1; i<parts.size(); ++i) {
 94        ymom += parts[i].mom();
 95      }
 96      double my2 = ymom.mass2();
 97      if (my2<0) vetoEvent;
 98
 99      // calculate xix and xiy and xidd
100      double xix  = mx2 / sqr(sqrtS());
101      double xiy  = my2 / sqr(sqrtS());
102      double xidd = mx2*my2 / sqr(sqrtS()*0.938*GeV);
103
104      // combine the selection: xi cuts
105      bool passedHadronCuts = false;
106      if (_sqs=="900"s  && (xix > 0.1  || xiy > 0.4 || xidd > 0.5))  passedHadronCuts = true;
107      if (_sqs=="2760"s && (xix > 0.07 || xiy > 0.2 || xidd > 0.5))  passedHadronCuts = true;
108      if (_sqs=="7000"s && (xix > 0.04 || xiy > 0.1 || xidd > 0.5))  passedHadronCuts = true;
109      if (!passedHadronCuts) vetoEvent;
110
111      //  ============================== MINIMUM BIAS EVENTS
112
113      // loop over particles to calculate the energy
114      _c[_sqs+"pass"]->fill();
115
116      for (const Particle& p : fsv.particles()) {
117        if (-5.2 > p.eta() && p.eta() > -6.6) _c[_sqs+"incl"]->fill(p.E()/GeV);
118      }
119
120      //  ============================== JET EVENTS
121
122      const FastJets& jetpro = apply<FastJets>(event, "Jets");
123      const Jets& jets = jetpro.jetsByPt(Cuts::pT > 1.0*GeV);
124      if (jets.size() < 1) vetoEvent;
125
126      if (jets[0].abseta() < 2.0) {
127        _h[_sqs+"den"]->fill(jets[0].pT()/GeV);
128
129        // energy flow
130        for (const Particle& p : fsv.particles()) {
131          if (p.eta() > -6.6 && p.eta() < -5.2) {  // ask for the CASTOR region
132            _h[_sqs+"num"]->fill(jets[0].pT()/GeV, p.E()/GeV);
133          }
134        }
135      }
136
137    }// analysis
138
139    void finalize() {
140      for (double eVal : allowedEnergies()) {
141        const string en = toString(int(eVal));
142        if (_c[en+"incl"]->sumW()) scale(_h[en+"num"], *_c[en+"pass"] / *_c[en+"incl"]);
143        divide(_h[en+"num"], _h[en+"den"], _e[en+"ratio"]);
144      }
145    }
146
147  private:
148    // counters
149    map<string,CounterPtr> _c;
150    map<string,Histo1DPtr> _h;
151    map<string,Estimate1DPtr> _e;
152
153    string _sqs = "";
154  };
155
156
157  RIVET_DECLARE_PLUGIN(CMS_2013_I1218372);
158
159}