1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/VetoedFinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Projections/PromptFinalState.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  // ATLAS $tt+b(b)$ at 8~TeV
 13  class ATLAS_2015_I1390114 : public Analysis {
 14  public:
 15
 16    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1390114);
 17
 18
 19    void init() {
 20
 21      // Eta ranges
 22      Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT >= 1.0*MeV;
 23      Cut eta_lep  = Cuts::abseta < 2.5;
 24
 25      // All final state particles
 26      FinalState fs(eta_full);
 27
 28      // Get photons to dress leptons
 29      IdentifiedFinalState photons(fs);
 30      photons.acceptIdPair(PID::PHOTON);
 31
 32      // Projection to find the electrons
 33      IdentifiedFinalState el_id(fs);
 34      el_id.acceptIdPair(PID::ELECTRON);
 35      PromptFinalState electrons(el_id);
 36      electrons.acceptTauDecays(true);
 37      declare(electrons, "electrons");
 38      LeptonFinder dressedelectrons(electrons, photons, 0.1, eta_lep && Cuts::pT > 25*GeV);
 39      declare(dressedelectrons, "dressedelectrons");
 40      LeptonFinder ewdressedelectrons(electrons, photons, 0.1, eta_full);
 41
 42      // Projection to find the muons
 43      IdentifiedFinalState mu_id(fs);
 44      mu_id.acceptIdPair(PID::MUON);
 45      PromptFinalState muons(mu_id);
 46      muons.acceptTauDecays(true);
 47      declare(muons, "muons");
 48      LeptonFinder dressedmuons(muons, photons, 0.1, eta_lep && Cuts::pT > 25*GeV);
 49      declare(dressedmuons, "dressedmuons");
 50      LeptonFinder ewdressedmuons(muons, photons, 0.1, eta_full);
 51
 52      // Projection to find neutrinos and produce MET
 53      IdentifiedFinalState nu_id;
 54      nu_id.acceptNeutrinos();
 55      PromptFinalState neutrinos(nu_id);
 56      neutrinos.acceptTauDecays(true);
 57
 58      // Jet clustering.
 59      VetoedFinalState vfs;
 60      vfs.addVetoOnThisFinalState(ewdressedelectrons);
 61      vfs.addVetoOnThisFinalState(ewdressedmuons);
 62      vfs.addVetoOnThisFinalState(neutrinos);
 63      FastJets jets(vfs,JetAlg::ANTIKT, 0.4);
 64      jets.useInvisibles();
 65      declare(jets, "jets");
 66
 67      book(_histo ,1,1,1);
 68      book(_ratio, 2,1,1);
 69      book(_aux   ,"_aux");
 70    }
 71
 72
 73    void analyze(const Event& event) {
 74
 75      // Get the selected objects, using the projections.
 76      DressedLeptons electrons = apply<LeptonFinder>(event, "dressedelectrons").dressedLeptons();
 77      DressedLeptons muons = apply<LeptonFinder>(event, "dressedmuons").dressedLeptons();
 78      if (electrons.empty() && muons.empty())  vetoEvent;
 79
 80      Jets jets = apply<FastJets>(event, "jets").jets((Cuts::pT>20*GeV) && (Cuts::abseta < 2.5), cmpMomByPt);
 81      Jets bjets;
 82
 83      // Check overlap of jets/leptons.
 84      for (Jet jet : jets) {
 85        // if dR(el,jet) < 0.4 skip the event
 86        for (DressedLepton el : electrons) {
 87          if (deltaR(jet, el) < 0.4)  vetoEvent;
 88        }
 89        // if dR(mu,jet) < 0.4 skip the event
 90        for (DressedLepton mu : muons) {
 91          if (deltaR(jet, mu) < 0.4)  vetoEvent;
 92        }
 93        // Count the number of b-tags
 94        // We have to check that the ghost-matched B hadrons have pT > 5 GeV
 95        // By default jet.bTags() returns all B hadrons without cuts
 96        bool btagged = jet.bTags(Cuts::pT >= 5*GeV).size();
 97        if (btagged)  bjets += jet;
 98      }
 99
100      // Evaluate basic event selection
101      bool pass_1lep = (electrons.size() == 1 && muons.empty()) || (muons.size() == 1 && electrons.empty());
102      if (pass_1lep && bjets.size() >= 3 && jets.size() >= 5)  _histo->fill(1);
103
104      if (muons.size() == 1 && electrons.size() == 1 && bjets.size() >= 3) {
105        if (muons[0].charge() * electrons[0].charge() < 0.0)  _histo->fill(2);
106      }
107
108      DressedLepton *lep1 = nullptr, *lep2 = nullptr;
109      bool zveto = false;
110      if (electrons.size() == 2 && muons.empty()) {
111        lep1 = &electrons[0];
112        lep2 = &electrons[1];
113      }
114      else if (muons.size() == 2 && electrons.empty()) {
115        lep1 = &muons[0];
116        lep2 = &muons[1];
117      }
118      else if (electrons.size() == 1 && muons.size() == 1) {
119        lep1 = &electrons[0];
120        lep2 = &muons[0];
121        zveto = true;
122      }
123      if (lep1 && lep2 && lep1->charge() * lep2->charge() < 0.0 && bjets.size() >= 2) {
124        double mass = (lep1->momentum() + lep2->momentum()).mass();
125        bool pass_2lep = mass > 15*GeV && (zveto || !(mass > 81*GeV && mass < 101*GeV));
126        if ( pass_2lep && bjets.size() >= 4) {
127          _histo->fill(3);
128          _histo->fill(4);
129        }
130        if ( pass_2lep && jets.size() >= 4) {
131          _aux->fill();
132        }
133      }
134    }
135
136
137    void finalize() {
138      const double sf(crossSection()/sumOfWeights()/femtobarn);
139      scale(_histo, sf);
140      scale(_aux,  sf);
141
142
143      const double  n = _histo->bin(4).sumW();
144      const double dN = _histo->bin(4).sumW2();
145      const double  d = _aux->sumW();
146      const double dD = _aux->sumW2();
147      const double  r = safediv(n, d);
148      double e = sqrt( safediv(r * (1 - r), d) );
149      if ( _aux->effNumEntries() != _aux->numEntries() ) {
150        // use F. James's approximation for weighted events:
151        e = sqrt( safediv((1 - 2 * r) * dN + r * r * dD, d * d) );
152      }
153      _ratio->bin(1).set(100.0 * r, 100.0 * e); // convert into percentage
154    }
155
156
157  private:
158
159    Histo1DPtr _histo;
160    CounterPtr _aux;
161    Estimate1DPtr _ratio;
162
163  };
164
165
166  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1390114);
167
168}