1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/IdentifiedFinalState.hh"
  7#include "Rivet/Projections/PromptFinalState.hh"
  8#include "Rivet/Projections/LeptonFinder.hh"
  9#include "Rivet/Projections/FastJets.hh"
 10
 11namespace Rivet {
 12
 13
 14  class ATLAS_2015_I1397637 : public Analysis {
 15  public:
 16
 17    /// Constructor
 18    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1397637);
 19
 20
 21    /// Book projections and histograms
 22    void init() {
 23
 24      // Base final state definition
 25      const FinalState fs(Cuts::abseta < 4.5);
 26
 27      // Neutrinos for MET
 28      IdentifiedFinalState nu_id;
 29      nu_id.acceptNeutrinos();
 30      PromptFinalState neutrinos(nu_id);
 31      neutrinos.acceptTauDecays(true);
 32      declare(neutrinos, "neutrinos");
 33
 34      // Get photons used to dress leptons
 35      IdentifiedFinalState photons(fs);
 36      photons.acceptId(PID::PHOTON);
 37
 38      // Use all bare muons as input to the DressedMuons projection
 39      IdentifiedFinalState mu_id(fs);
 40      mu_id.acceptIdPair(PID::MUON);
 41      PromptFinalState bare_mu(mu_id);
 42      bare_mu.acceptTauDecays(true);
 43
 44      // Use all bare electrons as input to the DressedElectrons projection
 45      IdentifiedFinalState el_id(fs);
 46      el_id.acceptIdPair(PID::ELECTRON);
 47      PromptFinalState bare_el(el_id);
 48      bare_el.acceptTauDecays(true);
 49
 50      // Use all bare leptons including taus for single-lepton filter
 51      IdentifiedFinalState lep_id(fs);
 52      lep_id.acceptIdPair(PID::MUON);
 53      lep_id.acceptIdPair(PID::ELECTRON);
 54      PromptFinalState bare_lep(lep_id);
 55      declare(bare_lep, "bare_lep");
 56
 57      // Tau finding
 58      /// @todo Use TauFinder
 59      UnstableParticles ufs;
 60      IdentifiedFinalState tau_id(ufs);
 61      tau_id.acceptIdPair(PID::TAU);
 62      PromptFinalState bare_tau(tau_id);
 63      declare(bare_tau, "bare_tau");
 64
 65      // Muons and electrons must have |eta| < 2.5
 66      Cut eta_ranges = Cuts::abseta < 2.5;
 67
 68      // Get dressed muons and the good muons (pt>25GeV)
 69      LeptonFinder all_dressed_mu(bare_mu, photons, 0.1, eta_ranges);
 70      LeptonFinder dressed_mu(bare_mu, photons, 0.1, eta_ranges && Cuts::pT > 25*GeV);
 71      declare(dressed_mu, "muons");
 72
 73      // Get dressed electrons and the good electrons (pt>25GeV)
 74      LeptonFinder all_dressed_el(bare_el, photons, 0.1, eta_ranges);
 75      LeptonFinder dressed_el(bare_el, photons, 0.1, eta_ranges && Cuts::pT > 25*GeV);
 76      declare(dressed_el, "electrons");
 77
 78      // Jet clustering
 79      VetoedFinalState vfs(fs);
 80      vfs.addVetoOnThisFinalState(all_dressed_el);
 81      vfs.addVetoOnThisFinalState(all_dressed_mu);
 82      vfs.addVetoOnThisFinalState(neutrinos);
 83
 84      // Small-R jets
 85      /// @todo Use extra constructor args
 86      FastJets jets(vfs, JetAlg::ANTIKT, 0.4);
 87      jets.useInvisibles(JetInvisibles::ALL);
 88      jets.useMuons(JetMuons::DECAY);
 89      declare(jets, "jets");
 90
 91      // Large-R jets
 92      /// @todo Use extra constructor args
 93      FastJets large_jets(vfs, JetAlg::ANTIKT, 1.0);
 94      large_jets.useInvisibles(JetInvisibles::ALL);
 95      large_jets.useMuons(JetMuons::DECAY);
 96      declare(large_jets, "fat_jets");
 97
 98
 99      /// Book histogram
100      book(_h_pttop ,1, 1, 1);
101    }
102
103
104    /// Perform the per-event analysis
105    void analyze(const Event& event) {
106
107      // Single lepton filter on bare leptons with no cuts
108      const Particles& bare_lep = apply<PromptFinalState>(event, "bare_lep").particles();
109      const Particles& bare_tau = apply<PromptFinalState>(event, "bare_tau").particles();
110      if (bare_lep.size() + bare_tau.size() != 1) vetoEvent;
111
112      // Electrons and muons
113      const DressedLeptons& electrons = apply<LeptonFinder>(event, "electrons").dressedLeptons();
114      const DressedLeptons& muons = apply<LeptonFinder>(event, "muons").dressedLeptons();
115      if (electrons.size() + muons.size() != 1) vetoEvent;
116      const DressedLepton& lepton = muons.empty() ? electrons[0] : muons[0];
117
118      // Get the neutrinos from the event record (they have pT > 0.0 and |eta| < 4.5 at this stage
119      const Particles& neutrinos = apply<PromptFinalState>(event, "neutrinos").particlesByPt();
120      FourMomentum met;
121      for (const Particle& nu : neutrinos) met += nu.momentum();
122      if (met.pT() < 20*GeV) vetoEvent;
123
124
125      // Thin jets and trimmed fat jets
126      /// @todo Use Rivet built-in FJ trimming support
127      const Jets& jets  = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
128      const PseudoJets& fat_pjets = apply<FastJets>(event, "fat_jets").pseudojetsByPt();
129      const double Rfilt = 0.3, ptFrac_min = 0.05; ///< @todo Need to be careful about the units for the pT cut passed to FJ?
130      PseudoJets trimmed_fat_pjets;
131      fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, Rfilt), fastjet::SelectorPtFractionMin(ptFrac_min));
132      for (const PseudoJet& pjet : fat_pjets) trimmed_fat_pjets += trimmer(pjet);
133      trimmed_fat_pjets = fastjet::sorted_by_pt(trimmed_fat_pjets);
134
135      // Jet reclustering
136      // Use a kT cluster sequence to recluster the trimmed jets so that a d12 can then be obtained from the reclustered jet
137      vector<double> splittingScales;
138      for (const PseudoJet& tpjet : trimmed_fat_pjets) {
139        const PseudoJets tpjet_constits = tpjet.constituents();
140        const fastjet::ClusterSequence kt_cs(tpjet_constits, fastjet::JetDefinition(fastjet::kt_algorithm, 1.5, fastjet::E_scheme, fastjet::Best));
141        const PseudoJets kt_jets = kt_cs.inclusive_jets();
142        const double d12 = 1.5 * sqrt(kt_jets[0].exclusive_subdmerge(1));
143        splittingScales += d12;
144      }
145      Jets trimmed_fat_jets;
146      for (size_t i = 0; i < trimmed_fat_pjets.size(); ++i) {
147        const Jet tj = trimmed_fat_pjets[i];
148        if (tj.mass()          <= 100*GeV) continue;
149        if (tj.pT()            <= 300*GeV) continue;
150        if (splittingScales[i] <=  40*GeV) continue;
151        if (tj.abseta()        >=     2.0) continue;
152        trimmed_fat_jets += tj;
153      }
154      if (trimmed_fat_jets.empty()) vetoEvent;
155
156      // Jet b-tagging
157      Jets bjets, non_bjets;
158      for (const Jet& jet : jets)
159        (jet.bTagged() ? bjets : non_bjets) += jet;
160      if (bjets.empty()) vetoEvent;
161
162
163      // Boosted selection: lepton/jet overlap
164      const double transmass = sqrt( 2 * lepton.pT() * met.pT() * (1 - cos(deltaPhi(lepton, met))) );
165      if (transmass + met.pt() <= 60*GeV) vetoEvent;
166      int lepJetIndex = -1;
167      for (size_t i = 0; i < jets.size(); ++i) {
168        const Jet& jet = jets[i];
169        if (deltaR(jet, lepton) < 1.5) {
170          lepJetIndex = i;
171          break;
172        }
173      }
174      if (lepJetIndex < 0) vetoEvent;
175      const Jet& ljet = jets[lepJetIndex];
176
177      // Boosted selection: lepton-jet/fat-jet matching
178      int fatJetIndex = -1;
179      for (size_t j = 0; j < trimmed_fat_jets.size(); ++j) {
180        const Jet& fjet = trimmed_fat_jets[j];
181        const double dR_fatjet = deltaR(ljet, fjet);
182        const double dPhi_fatjet = deltaPhi(lepton, fjet);
183        if (dR_fatjet > 1.5 && dPhi_fatjet > 2.3) {
184          fatJetIndex = j;
185          break;
186        }
187      }
188      if (fatJetIndex < 0) vetoEvent;
189      const Jet& fjet = trimmed_fat_jets[fatJetIndex];
190
191      // Boosted selection: b-tag matching
192      const bool lepbtag = ljet.bTagged();
193      bool hadbtag = false;
194      for (const Jet& bjet : bjets) {
195        hadbtag |= (deltaR(fjet, bjet) < 1.0);
196      }
197
198      // Fill histo if selection passed
199      if (hadbtag || lepbtag) _h_pttop->fill(fjet.pT()/GeV);
200    }
201
202
203    /// Normalise histograms etc., after the run
204    void finalize() {
205      scale(_h_pttop, crossSection()/femtobarn / sumOfWeights());
206    }
207
208
209  private:
210
211    Histo1DPtr _h_pttop;
212
213  };
214
215
216  RIVET_DECLARE_PLUGIN(ATLAS_2015_I1397637);
217
218}