1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/LeptonFinder.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/InvisibleFinalState.hh"
  8#include "Rivet/Projections/VetoedFinalState.hh"
  9
 10namespace Rivet {
 11
 12
 13  /// @brief ttbar + gamma at 13 TeV
 14  class ATLAS_2018_I1707015 : public Analysis {
 15  public:
 16
 17
 18    // Constructor
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1707015);
 20
 21
 22    // Book histograms and initialise projections before the run
 23    void init() {
 24
 25      // Set default running mode to 3 (all)
 26      _mode = 3;
 27
 28      // Get running mode
 29      if ( getOption("LMODE") == "SINGLE" )   _mode = 1;
 30      if ( getOption("LMODE") == "DILEPTON" ) _mode = 2;
 31      if ( getOption("LMODE") == "ALL" )      _mode = 3;
 32
 33      // All final state particles
 34      const FinalState fs;
 35
 36      // Charged particles for signal photon isolation
 37      ChargedFinalState cfs(fs);
 38      declare(cfs, "CFS");
 39
 40      // Signal photons
 41      PromptFinalState photons(Cuts::abspid == PID::PHOTON && Cuts::pT > 20*GeV && Cuts::abseta < 2.37, TauDecaysAs::PROMPT);
 42      declare(photons, "Photons");
 43
 44      // Leptons
 45      PromptFinalState leptons(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 46
 47      // Dress the leptons
 48      FinalState dressPhotons(Cuts::abspid == PID::PHOTON);
 49      Cut lepCuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV);
 50      LeptonFinder dressedLeptons(leptons, dressPhotons, 0.1, lepCuts);
 51      declare(dressedLeptons, "Leptons");
 52
 53      // Jet alg input
 54      VetoedFinalState vfs(fs);
 55
 56      // Remove prompt invisibles from jet input
 57      vfs.addVetoOnThisFinalState(InvisibleFinalState(OnlyPrompt::YES, TauDecaysAs::PROMPT));
 58
 59      // Remove prompt dressed muons (muons + associated photons) from jet input
 60      PromptFinalState muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 61      LeptonFinder dressedmuons(muons, dressPhotons, 0.1);
 62      vfs.addVetoOnThisFinalState(dressedmuons);
 63
 64      // Jet clustering
 65      FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
 66      declare(jets, "Jets");
 67
 68      // Book histograms
 69      if ( _mode == 1 or _mode == 3 ){  // Single lepton channel
 70        book(_h["sl_ph_pt"],   3,1,1);  // photon pT
 71        book(_h["sl_ph_eta"],  4,1,1);  // photon eta
 72        book(_h["sl_ph_l_dR"], 5,1,1);  // photon-lepton dR
 73      }
 74      if ( _mode == 2 or _mode == 3 ){  // Dilepton channel
 75        book(_h["dl_ph_pt"],   6,1,1);  // photon pT
 76        book(_h["dl_ph_eta"],  7,1,1);  // photon eta
 77        book(_h["dl_ph_l_dR"], 8,1,1);  // min photon-lepton dR
 78        book(_h["dl_l_dEta"],  9,1,1);  // lepton-lepton dEta
 79        book(_h["dl_l_dPhi"],  10,1,1); // lepton-lepton dPhi
 80      }
 81    }
 82
 83
 84    void analyze(const Event& event) {
 85
 86      // Fetch objects
 87      const DressedLeptons& leptons = apply<LeptonFinder>(event, "Leptons").dressedLeptons();
 88      const Particles& photons = apply<PromptFinalState>(event, "Photons").particles();
 89      ChargedFinalState charged = apply<ChargedFinalState>(event, "CFS");
 90      Jets jets = apply<JetFinder>(event, "Jets").jetsByPt(Cuts::abseta < 2.5 && Cuts::pT > 25*GeV);
 91
 92      // Immediate veto on events without one good photon
 93      if (photons.size() != 1) vetoEvent;
 94      const Particle& photon = photons[0];
 95
 96      // Veto event if photon too close to a lepton
 97      for (const DressedLepton& lep : leptons) {
 98        if ( deltaR(lep, photon) < 1.0 ) vetoEvent;
 99      }
100
101      // Overlap removel of jets near leptons
102      idiscardIfAnyDeltaRLess(jets, leptons, 0.4);
103
104      // Overlap removel of jets near isolated photon
105      const double conePt = sum(charged.particles(deltaRLess(photon, 0.3)), Kin::pT, 0.0);
106      if ( conePt / photon.pT() < 0.1 ) {
107        idiscard(jets, deltaRLess(photon, 0.4));
108      }
109
110      // Veto event if photon too close to good jets
111      for (const Jet& jet : jets) {
112        if ( deltaR(jet, photon) < 0.4 ) vetoEvent;
113      }
114
115      // Require at least one bjet
116      unsigned int nbjets = 0;
117      for (const Jet& jet : jets) {
118        if ( jet.bTagged(Cuts::pT > 5*GeV) ) nbjets += 1;
119      }
120      if ( nbjets == 0 ) vetoEvent;
121
122      // Fill histos in single lepton channel
123      if ( _mode == 1 || _mode == 3 ) {
124        if ( jets.size() >= 4 && leptons.size() == 1 ) {
125          _h["sl_ph_pt"]->fill(   photon.pT()/GeV );
126          _h["sl_ph_eta"]->fill(  photon.abseta() );
127          _h["sl_ph_l_dR"]->fill( deltaR(leptons[0], photon) );
128        }
129      }
130
131      // Fill histos in dilepton channel
132      if ( _mode == 2 || _mode == 3 ) {
133        if ( jets.size() >= 2 && leptons.size() == 2 ) {
134          double deltaRNew;
135          double deltaRMin = 999.0;
136          for (const DressedLepton& lep : leptons) {
137            deltaRNew = deltaR(lep, photon);
138            if ( deltaRNew < deltaRMin ) deltaRMin = deltaRNew;
139          }
140          _h["dl_ph_pt"]->fill( photon.pT()/GeV );
141          _h["dl_ph_eta"]->fill( photon.abseta() );
142          _h["dl_ph_l_dR"]->fill( deltaRMin );
143          _h["dl_l_dEta"]->fill( deltaEta(leptons[0], leptons[1]) );
144          _h["dl_l_dPhi"]->fill( deltaPhi(leptons[0], leptons[1]) );
145        }
146      }
147    }
148
149
150    void finalize() {
151      // Normalise histograms after the run
152      normalize(_h, 1.0, true);
153    }
154
155
156  private:
157    int _mode;
158    map<string, Histo1DPtr> _h;
159
160  };
161
162  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1707015);
163}