1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FastJets.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/LeptonFinder.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/UnstableParticles.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief ttbar + gamma at 8 TeV
 13  class ATLAS_2017_I1604029 : public Analysis {
 14  public:
 15
 16    // Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2017_I1604029);
 18
 19    // Book histograms and initialise projections before the run
 20    void init() {
 21
 22      const FinalState fs;
 23
 24      // signal photons
 25      PromptFinalState prompt_ph(Cuts::abspid == PID::PHOTON && Cuts::pT > 15*GeV && Cuts::abseta < 2.37);
 26      declare(prompt_ph, "photons");
 27
 28      // bare leptons
 29      Cut base_cuts = (Cuts::abseta < 2.7) && (Cuts::pT > 10*GeV);
 30      IdentifiedFinalState bare_leps(base_cuts);
 31      bare_leps.acceptIdPair(PID::MUON);
 32      bare_leps.acceptIdPair(PID::ELECTRON);
 33      declare(bare_leps, "bare_leptons");
 34
 35      // dressed leptons
 36      Cut dressed_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
 37      PromptFinalState prompt_mu(base_cuts && Cuts::abspid == PID::MUON);
 38      PromptFinalState prompt_el(base_cuts && Cuts::abspid == PID::ELECTRON);
 39      IdentifiedFinalState all_photons(fs, PID::PHOTON);
 40      LeptonFinder elecs(prompt_el, all_photons, 0.1, dressed_cuts);
 41      declare(elecs, "elecs");
 42      LeptonFinder muons(prompt_mu, all_photons, 0.1, dressed_cuts);
 43      declare(muons, "muons");
 44
 45      // auxiliary projections for 'single-lepton ttbar filter'
 46      PromptFinalState prompt_lep(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON);
 47      declare(prompt_lep, "prompt_leps");
 48      declare(UnstableParticles(), "ufs");
 49
 50      // jets
 51      FastJets jets(fs, JetAlg::ANTIKT, 0.4, JetMuons::NONE, JetInvisibles::NONE);
 52      declare(jets, "jets");
 53
 54
 55      // BOOK HISTOGRAMS
 56      book(_h["pt"],  2,1,1);
 57      book(_h["eta"], 3,1,1);
 58
 59    }
 60
 61
 62    // Perform the per-event analysis
 63    void analyze(const Event& event) {
 64
 65      // analysis extrapolated to 1-lepton-plus-jets channel, where "lepton" cannot be a tau
 66      // (i.e. contribution from dileptonic ttbar where one of the leptons is outside
 67      // the detector acceptance has been subtracted as a background)
 68      if (apply<PromptFinalState>(event, "prompt_leps").particles().size() != 1)  vetoEvent;
 69      for (const auto& p : apply<UnstableParticles>(event, "ufs").particles()) {
 70        if (p.fromPromptTau())  vetoEvent;
 71      }
 72
 73      // photon selection
 74      Particles photons = apply<PromptFinalState>(event, "photons").particlesByPt();
 75      Particles bare_leps  = apply<IdentifiedFinalState>(event, "bare_leptons").particles();
 76      for (const Particle& lep : bare_leps)
 77        idiscard(photons, deltaRLess(lep, 0.1));
 78      if (photons.size() != 1)  vetoEvent;
 79      const Particle& photon = photons[0];
 80
 81      // jet selection
 82      Jets jets = apply<JetFinder>(event, "jets").jetsByPt(Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
 83
 84      // lepton selection
 85      const DressedLeptons& elecs = apply<LeptonFinder>(event, "elecs").dressedLeptons();
 86      const DressedLeptons& all_muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
 87
 88      // jet photon/electron overlap removal
 89      for (const DressedLepton& e : elecs)
 90        idiscard(jets, deltaRLess(e, 0.2, RAPIDITY));
 91      for (const Particle& ph : photons)
 92        idiscard(jets, deltaRLess(ph, 0.1, RAPIDITY));
 93	    if (jets.size() < 4)  vetoEvent;
 94
 95      // photon-jet minimum deltaR
 96      double mindR_phjet = 999.;
 97      for (Jet jet : jets) {
 98        const double dR_phjet = deltaR(photon, jet);
 99        if (dR_phjet < mindR_phjet) mindR_phjet = dR_phjet;
100      }
101      if (mindR_phjet < 0.5)  vetoEvent;
102
103      // muon jet overlap removal
104      DressedLeptons muons;
105      for (const DressedLepton& mu : all_muons) {
106        bool overlaps = false;
107        for (const Jet& jet : jets) {
108          if (deltaR(mu, jet) < 0.4) {
109            overlaps = true;
110            break;
111          }
112        }
113        if (overlaps) continue;
114        muons.push_back(mu);
115      }
116
117      // one electron XOR one muon
118      bool isEl = elecs.size() == 1 && muons.size() == 0;
119      bool isMu = muons.size() == 1 && elecs.size() == 0;
120      if (!isEl && !isMu)  vetoEvent;
121
122      // photon-lepton deltaR
123      double mindR_phlep = deltaR(photon, isEl? elecs[0] : muons[0]);
124      if (mindR_phlep < 0.7)  vetoEvent;
125
126      // b-tagging
127      Jets bjets;
128      for (const Jet& jet : jets) {
129        if (jet.bTagged(Cuts::pT > 5*GeV)) bjets +=jet;
130      }
131      if (bjets.empty())  vetoEvent;
132
133      _h["pt"]->fill(photon.pT()/GeV);
134      _h["eta"]->fill(photon.abseta());
135    }
136
137
138    // Normalise histograms etc., after the run
139    void finalize() {
140      const double normto(crossSection() / femtobarn / sumOfWeights());
141      for (auto &hist : _h) {  scale(hist.second, normto);  }
142    }
143
144
145  private:
146
147    map<string, Histo1DPtr> _h;
148
149  };
150
151  RIVET_DECLARE_PLUGIN(ATLAS_2017_I1604029);
152}