1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/PromptFinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/InvMassFinalState.hh"
  7#include "Rivet/Projections/VisibleFinalState.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9#include "Rivet/Projections/ZFinder.hh"
 10
 11namespace Rivet {
 12
 13
 14  class CMS_2015_I1310737 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2015_I1310737);
 19
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      FinalState fs; ///< @todo No cuts?
 25      VisibleFinalState visfs(fs);
 26      // Prompt leptons only
 27      PromptFinalState pfs(fs);
 28
 29      ZFinder zeeFinder(pfs, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV, PID::ELECTRON, 71.0*GeV, 111.0*GeV);
 30      declare(zeeFinder, "ZeeFinder");
 31
 32      ZFinder zmumuFinder(pfs, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV, PID::MUON, 71.0*GeV, 111.0*GeV);
 33      declare(zmumuFinder, "ZmumuFinder");
 34
 35      VetoedFinalState jetConstits(visfs);
 36      jetConstits.addVetoOnThisFinalState(zeeFinder);
 37      jetConstits.addVetoOnThisFinalState(zmumuFinder);
 38
 39      FastJets akt05Jets(jetConstits, FastJets::ANTIKT, 0.5);
 40      declare(akt05Jets, "AntiKt05Jets");
 41
 42
 43      book(_h_excmult_jets_tot ,1, 1, 1);
 44      book(_h_incmult_jets_tot ,2, 1, 1);
 45      book(_h_leading_jet_pt_tot ,3, 1, 1);
 46      book(_h_second_jet_pt_tot ,4, 1, 1);
 47      book(_h_third_jet_pt_tot ,5, 1, 1);
 48      book(_h_fourth_jet_pt_tot ,6, 1, 1);
 49      book(_h_leading_jet_eta_tot ,7, 1, 1);
 50      book(_h_second_jet_eta_tot ,8, 1, 1);
 51      book(_h_third_jet_eta_tot ,9, 1, 1);
 52      book(_h_fourth_jet_eta_tot ,10, 1, 1);
 53      book(_h_ht1_tot ,11, 1, 1);
 54      book(_h_ht2_tot ,12, 1, 1);
 55      book(_h_ht3_tot ,13, 1, 1);
 56      book(_h_ht4_tot ,14, 1, 1);
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {;
 62
 63      const ZFinder& zeeFS = apply<ZFinder>(event, "ZeeFinder");
 64      const ZFinder& zmumuFS = apply<ZFinder>(event, "ZmumuFinder");
 65
 66      const Particles& zees = zeeFS.bosons();
 67      const Particles& zmumus = zmumuFS.bosons();
 68
 69      // We did not find exactly one Z. No good.
 70      if (zees.size() + zmumus.size() != 1) {
 71        MSG_DEBUG("Did not find exactly one good Z candidate");
 72        vetoEvent;
 73      }
 74
 75      // Find the (dressed!) leptons
 76      const Particles& dressedLeptons = zees.size() ? zeeFS.constituents() : zmumuFS.constituents();
 77
 78      // Cluster jets
 79      // NB. Veto has already been applied on leptons and photons used for dressing
 80      const FastJets& fj = apply<FastJets>(event, "AntiKt05Jets");
 81      const Jets& jets = fj.jetsByPt(Cuts::abseta < 2.4 && Cuts::pT > 30*GeV);
 82
 83      // Perform lepton-jet overlap and HT calculation
 84      double ht = 0;
 85      Jets goodjets;
 86      for (const Jet& j : jets) {
 87        // Decide if this jet is "good", i.e. isolated from the leptons
 88        /// @todo Nice use-case for any() and a C++11 lambda
 89        bool overlap = false;
 90        for (const Particle& l : dressedLeptons) {
 91          if (deltaR(j, l) < 0.5) {
 92            overlap = true;
 93            break;
 94          }
 95        }
 96
 97        // Fill HT and good-jets collection
 98        if (overlap) continue;
 99        goodjets.push_back(j);
100        ht += j.pT();
101      }
102
103      // We don't care about events with no isolated jets
104      if (goodjets.empty()) {
105        MSG_DEBUG("No jets in event");
106        vetoEvent;
107      }
108
109
110      /////////////////
111
112
113      // Weight to be used for histo filling
114      const double w = 0.5 * 1.0;
115
116      // Fill jet number integral histograms
117      _h_excmult_jets_tot->fill(goodjets.size(), w);
118      /// @todo Could be better computed by toIntegral transform on exclusive histo
119      for (size_t iJet = 1; iJet <= goodjets.size(); iJet++ )
120        _h_incmult_jets_tot->fill(iJet, w);
121
122      // Fill leading jet histograms
123      const Jet& j1 = goodjets[0];
124      _h_leading_jet_pt_tot->fill(j1.pT()/GeV, w);
125      _h_leading_jet_eta_tot->fill(j1.abseta(), w);
126      _h_ht1_tot->fill(ht/GeV, w);
127
128      // Fill 2nd jet histograms
129      if (goodjets.size() < 2) return;
130      const Jet& j2 = goodjets[1];
131      _h_second_jet_pt_tot->fill(j2.pT()/GeV, w);
132      _h_second_jet_eta_tot->fill(j2.abseta(), w);
133      _h_ht2_tot->fill(ht/GeV, w);
134
135      // Fill 3rd jet histograms
136      if (goodjets.size() < 3) return;
137      const Jet& j3 = goodjets[2];
138      _h_third_jet_pt_tot->fill(j3.pT()/GeV, w);
139      _h_third_jet_eta_tot->fill(j3.abseta(), w);
140      _h_ht3_tot->fill(ht/GeV, w);
141
142      // Fill 4th jet histograms
143      if (goodjets.size() < 4) return;
144      const Jet& j4 = goodjets[3];
145      _h_fourth_jet_pt_tot->fill(j4.pT()/GeV, w);
146      _h_fourth_jet_eta_tot->fill(j4.abseta(), w);
147      _h_ht4_tot->fill(ht/GeV, w);
148    }
149
150
151    /// Normalise histograms etc., after the run
152    void finalize() {
153
154      const double norm = (sumOfWeights() != 0) ? crossSection()/sumOfWeights() : 1.0;
155
156      scale(_h_excmult_jets_tot, norm );
157      scale(_h_incmult_jets_tot, norm );
158      scale(_h_leading_jet_pt_tot, norm );
159      scale(_h_second_jet_pt_tot, norm );
160      scale(_h_third_jet_pt_tot, norm );
161      scale(_h_fourth_jet_pt_tot, norm );
162      scale(_h_leading_jet_eta_tot, norm );
163      scale(_h_second_jet_eta_tot, norm );
164      scale(_h_third_jet_eta_tot, norm );
165      scale(_h_fourth_jet_eta_tot, norm );
166      scale(_h_ht1_tot, norm );
167      scale(_h_ht2_tot, norm );
168      scale(_h_ht3_tot, norm );
169      scale(_h_ht4_tot, norm );
170    }
171
172
173  private:
174
175    /// @name Histograms
176
177    Histo1DPtr _h_excmult_jets_tot,  _h_incmult_jets_tot;
178    Histo1DPtr _h_leading_jet_pt_tot, _h_second_jet_pt_tot, _h_third_jet_pt_tot, _h_fourth_jet_pt_tot;
179    Histo1DPtr _h_leading_jet_eta_tot, _h_second_jet_eta_tot, _h_third_jet_eta_tot, _h_fourth_jet_eta_tot;
180    Histo1DPtr _h_ht1_tot, _h_ht2_tot, _h_ht3_tot, _h_ht4_tot;
181
182  };
183
184
185  RIVET_DECLARE_PLUGIN(CMS_2015_I1310737);
186
187
188}