1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/LeptonFinder.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/ChargedFinalState.hh"
  9#include "Rivet/Projections/MissingMomentum.hh"
 10#include "Rivet/Projections/InvisibleFinalState.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @brief Colour reconnection in ttbar dilepton events
 16  class ATLAS_2022_I2152933 : public Analysis {
 17  public:
 18
 19    /// Constructor
 20    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2022_I2152933);
 21
 22
 23    /// Book histograms and initialise projections before the run
 24    void init() {
 25
 26      // Initialise and register projections
 27
 28    	// Photons
 29    	PromptFinalState photons(Cuts::abspid == PID::PHOTON);
 30
 31    	// Muons
 32    	Cut lepton_cuts = Cuts::abseta < 2.5 && Cuts::pT > 25*GeV;
 33    	PromptFinalState bare_mu(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
 34    	LeptonFinder all_dressed_mu(bare_mu, photons, 0.1, lepton_cuts);
 35    	declare(all_dressed_mu, "muons");
 36
 37    	// Electrons
 38    	PromptFinalState bare_el(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
 39    	LeptonFinder all_dressed_el(bare_el, photons, 0.1, lepton_cuts);
 40    	declare(all_dressed_el, "electrons");
 41
 42    	//Jet forming
 43    	const InvisibleFinalState neutrinos(OnlyPrompt::YES, TauDecaysAs::PROMPT);
 44
 45    	VetoedFinalState vfs(FinalState(Cuts::abseta < 5.0));
 46    	vfs.addVetoOnThisFinalState(all_dressed_el);
 47    	vfs.addVetoOnThisFinalState(all_dressed_mu);
 48    	vfs.addVetoOnThisFinalState(neutrinos);
 49    	FastJets jetfs(vfs, JetAlg::ANTIKT, 0.4, JetMuons::ALL, JetInvisibles::ALL);
 50    	declare(jetfs, "jets");
 51
 52      // FinalState charged particles
 53      ChargedFinalState tracks(Cuts::pT > 0.5*GeV && Cuts::abseta < 2.5);
 54      declare(tracks, "tracks");
 55
 56      declare(MissingMomentum(), "ETmiss");
 57
 58      // Book histograms
 59      book(_h["nch"], 7,1,1);
 60      book(_h["sumPt"], 8,1,1);
 61
 62      book(_h_multi, { 0, 19.5, 39.5, 59.5, 79.5, 101 });
 63      for (auto& b : _h_multi->bins()) {
 64        book(b, 8+b.index(), 1, 1);
 65      }
 66
 67    }
 68
 69
 70    /// Perform the per-event analysis
 71    void analyze(const Event& event) {
 72      // Retrieve dressed leptons, sorted by pT
 73      DressedLeptons electrons = apply<LeptonFinder>(event, "electrons").dressedLeptons();
 74      DressedLeptons muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
 75      // Retrieve clustered jets, sorted by pT, with a minimum pT cut
 76      Jets jets = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
 77
 78      // Retrieve charge final state particles, sorted by pT
 79      const Particles& tracks = apply<ChargedFinalState>(event, "tracks").particlesByPt();
 80
 81      // OVERLAP REMOVAL
 82      idiscardIfAnyDeltaRLess(muons, jets, 0.4);
 83      idiscardIfAnyDeltaRLess(electrons, jets, 0.4);
 84
 85      // Select jets ghost-associated to B-hadrons with a certain fiducial selection
 86      Jets bjets = select(jets, [](const Jet& jet) {
 87        return  jet.bTagged(Cuts::pT > 5*GeV);
 88      });
 89
 90      // Veto event if there are not exactly one electron and one muon
 91      if ( electrons.size() != 1 || muons.size() != 1)  vetoEvent;
 92
 93      // Veto event if the selected electron and muon are not OS
 94      if (electrons[0].charge() == muons[0].charge())  vetoEvent;
 95
 96      // Veto event if there are no 2 jets or 3 jets
 97      if (jets.size() < 2 || jets.size() > 3)  vetoEvent;
 98
 99      // Veto event if the number of b-jets not exactly 2
100      if (bjets.size() != 2)  vetoEvent;
101
102      // Veto event with dilepton mass < 15 GeV
103      FourMomentum ll = electrons[0].momentum() + muons[0].momentum();
104      if (ll.mass() <= 15*GeV ) vetoEvent;
105
106      // Remove tracks inside jets and leptons: within 0.4 delatR from a jet and 0.01 from a lepton
107
108      Particles myTracks;
109      for (const Particle& p : tracks) {
110        bool isJetAssoc = false;
111        for (unsigned int i = 0; i < jets.size(); ++i) {
112          if(deltaR(jets[i], p, PSEUDORAPIDITY) < 0.4) {
113          	isJetAssoc= true;
114          	break;
115          }
116        }
117        if (isJetAssoc) continue;
118        if(deltaR(muons[0], p, PSEUDORAPIDITY) < 0.01) continue;
119        if(deltaR(electrons[0], p, PSEUDORAPIDITY) < 0.01) continue;
120        if(p.pid() == 0 || p.isStable() != 1 || p.charge() == 0) continue;
121        myTracks.push_back(p);
122      }
123
124      // Veto event if the number of selected charged particles higher than 100
125      if (myTracks.size() > 100)  vetoEvent;
126
127      // Fill nch, sumpt, and sumpt vs. nch
128      const double nch = min(myTracks.size(), (size_t)99); // put overflow in the last bin
129      _h["nch"]->fill(nch);
130
131      const double sumPt = min(sum(myTracks, Kin::pT, 0.0), 119*GeV);
132      _h["sumPt"]->fill(sumPt/GeV);
133
134      const double sumPt2d = min(sumPt, 99*GeV);
135      _h_multi->fill(nch, sumPt2d/GeV);
136    }
137
138
139    /// Normalise histograms etc., after the run
140    void finalize() {
141
142      // Normalize to unity
143      normalize(_h);
144      normalizeGroup(_h_multi, 1.0, false);
145      divByGroupWidth(_h_multi);
146    }
147
148
149    /// @name Histograms
150    /// @{
151    map<string, Histo1DPtr> _h;
152    map<string, CounterPtr> _c;
153    Histo1DGroupPtr _h_multi;
154
155    /// @}
156
157  };
158
159
160  RIVET_DECLARE_PLUGIN(ATLAS_2022_I2152933);
161
162}