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/VetoedFinalState.hh"
  7
  8namespace Rivet {
  9
 10
 11  /// Forward energy flow in MB and dijet events at 0.9 and 7 TeV
 12  class CMS_2011_I930319 : public Analysis {
 13  public:
 14
 15    /// Constructor
 16    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2011_I930319);
 17
 18
 19    void init() {
 20      const FinalState fs(Cuts::abseta < 6.0);
 21      declare(fs, "FS");
 22      declare(FastJets(fs, JetAlg::ANTIKT, 0.5), "Jets");
 23
 24      VetoedFinalState fsv(fs);
 25      fsv.vetoNeutrinos();
 26      fsv.addVetoPair(PID::MUON);
 27      declare(fsv, "fsv");
 28
 29      // For the MB ND selection
 30      const ChargedFinalState fschrgd(Cuts::abseta < 6.0);
 31      declare(fschrgd, "fschrgd");
 32      VetoedFinalState fschrgdv(fschrgd);
 33      fschrgdv.vetoNeutrinos();
 34      declare(fschrgdv, "fschrgdv");
 35
 36      for (double eVal : allowedEnergies()) {
 37        const string en = toString(int(eVal));
 38        if (isCompatibleWithSqrtS(eVal))  _sqs = en;
 39        size_t offset = (en == "900"s)? 0 : 2;
 40        book(_h[en+"mb"],    1+offset, 1, 1); // energy flow in MB
 41        book(_h[en+"dijet"], 2+offset, 1, 1); // energy flow in dijet events
 42        book(_c[en+"mb"],    "/tmp/weightMB"+en);
 43        book(_c[en+"dijet"], "/tmp/weightDijet"+en);
 44      }
 45      if (_sqs == "" && !merging()) {
 46        throw BeamError("Invalid beam energy for " + name() + "\n");
 47      }
 48
 49    }
 50
 51
 52    void analyze(const Event& event) {
 53      // Skip if the event is empty
 54      const FinalState& fsv = apply<FinalState>(event, "fsv");
 55      if (fsv.empty()) vetoEvent;
 56
 57      // Veto diffractive topologies according to defined hadron level
 58      double count_chrg_forward = 0;
 59      double count_chrg_backward = 0;
 60      const FinalState& fschrgdv = apply<FinalState>(event, "fschrgdv");
 61      for (const Particle& p : fschrgdv.particles()) {
 62        if (3.9 < p.eta() && p.eta() < 4.4)     ++count_chrg_forward;
 63        if (-4.4 < p.eta() && p.eta() < -3.9)  ++count_chrg_backward;
 64      }
 65      if (count_chrg_forward == 0 || count_chrg_backward == 0) vetoEvent;
 66      /// @todo "Diffractive" veto should really also veto dijet events?
 67
 68
 69      // MINIMUM BIAS EVENTS
 70      _c[_sqs+"mb"]->fill();
 71      for (const Particle& p: fsv.particles()) {
 72        _h[_sqs+"mb"]->fill(p.abseta(), p.E()/GeV);
 73      }
 74
 75
 76      // DIJET EVENTS
 77      const FastJets& jetpro = apply<FastJets>(event, "Jets");
 78      const double PTCUT = (_sqs=="900"s)? 8*GeV : 20*GeV;
 79      const Jets jets = jetpro.jetsByPt(Cuts::pT > PTCUT);
 80      if (jets.size() >= 2) {
 81        // eta cut for the central jets
 82        if (jets[0].abseta() < 2.5 && jets[1].abseta() < 2.5) {
 83          // Back to back condition of the jets
 84          const double diffphi = deltaPhi(jets[1].phi(), jets[0].phi());
 85          if (diffphi-PI < 1.0) {
 86            _c[_sqs+"dijet"]->fill();
 87            for (const Particle& p: fsv.particles()) {
 88              _h[_sqs+"dijet"]->fill(p.abseta(), p.E()/GeV);
 89            }
 90          }
 91        }
 92      }
 93
 94    }
 95
 96
 97    void finalize() {
 98      for (auto& item : _h) {
 99        scale(item.second, 0.5/_c[item.first]->sumW());
100      }
101    }
102
103
104  private:
105
106    /// @{
107    map<string,Histo1DPtr> _h;
108    map<string,CounterPtr> _c;
109
110    string _sqs = "";
111    /// @}
112
113  };
114
115
116
117  RIVET_DECLARE_ALIASED_PLUGIN(CMS_2011_I930319, CMS_2011_S9215166);
118
119}