1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PartonicTops.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8#include "Rivet/Projections/DressedLeptons.hh"
  9#include "Rivet/Projections/ChargedLeptons.hh"
 10
 11#include "fastjet/contrib/Nsubjettiness.hh"
 12#include "fastjet/contrib/XConePlugin.hh"
 13
 14namespace Rivet {
 15
 16
 17  /// @brief Measurement of the jet mass for boosted top quarks at 13 TeV
 18  class CMS_2019_I1764472 : public Analysis {
 19  public:
 20
 21    /// Constructor
 22    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2019_I1764472);
 23
 24
 25    /// @name Analysis methods
 26    //@{
 27
 28    void init() {
 29
 30      // Prompt leptons
 31      ChargedLeptons charged_leptons;
 32      PromptFinalState prompt_leptons(charged_leptons);
 33      declare(prompt_leptons, "PromptLeptons");
 34
 35      // Final state particles for jet clustering
 36      VetoedFinalState fs_jets;
 37      fs_jets.vetoNeutrinos();
 38
 39      // First XCone jet clustering step
 40      fastjet::contrib::PseudoXConePlugin* plugin_xcone = new fastjet::contrib::PseudoXConePlugin(2, 1.2, 2.0);
 41      declare(FastJets(fs_jets, plugin_xcone), "FatJets");
 42
 43      // Partonic tops for decay channel definition
 44      declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicTops");
 45      declare(PartonicTops(PartonicTops::DecayMode::HADRONIC), "HadronicTops");
 46
 47      // Book histograms
 48      book(_hist_mass, "d01-x01-y01");
 49      book(_hist_mass_norm, "d02-x01-y01");
 50    }
 51
 52
 53    /// Perform the per-event analysis
 54    void analyze(const Event& event) {
 55
 56      // Decay mode check
 57      const Particles& leptonicTops = apply<PartonicTops>(event, "LeptonicTops").particlesByPt();
 58      const Particles& hadronicTops = apply<PartonicTops>(event, "HadronicTops").particlesByPt();
 59      if (leptonicTops.size() != 1 || hadronicTops.size() != 1) vetoEvent;
 60
 61
 62      // Get prompt leptons
 63      const PromptFinalState& prompt_leptons = apply<PromptFinalState>(event, "PromptLeptons");
 64      const Particles & leptons = prompt_leptons.particles();
 65      if(leptons.empty()) vetoEvent;
 66
 67      // Select leading lepton
 68      Particle lepton;
 69      for(const Particle& l : leptons){
 70        if(l.pT() > lepton.pT()) lepton = l;
 71      }
 72      if(lepton.pT() < 60*GeV) vetoEvent;
 73
 74      // Get the fat jets
 75      const Jets& fatjets = applyProjection<FastJets>(event, "FatJets").jets();
 76
 77      // Get index of hadronic jet by distance to lepton
 78      int ihad = 0;
 79      int ilep = 1;
 80
 81      double dR0 = deltaR(lepton, fatjets.at(0));
 82      double dR1 = deltaR(lepton, fatjets.at(1));
 83
 84      if(dR0 < dR1){
 85        ihad = 1;
 86        ilep = 0;
 87      }
 88
 89      // Get jet constituents
 90      const Particles & phad = fatjets.at(ihad).particles();
 91      const Particles & plep = fatjets.at(ilep).particles();
 92
 93      // Cluster subjets
 94      FinalState fs_dummy;
 95      fastjet::JetDefinition::Plugin* plugin_subhad = new fastjet::contrib::PseudoXConePlugin(3, 0.4, 2.0);
 96      fastjet::contrib::PseudoXConePlugin* plugin_sublep = new fastjet::contrib::PseudoXConePlugin(3, 0.4, 2.0);
 97      FastJets hadsubcluster(fs_dummy, plugin_subhad);
 98      FastJets lepsubcluster(fs_dummy, plugin_sublep);
 99      hadsubcluster.calc(phad);
100      lepsubcluster.calc(plep);
101
102      Jets subjets_had = hadsubcluster.jets();
103      Jets subjets_lep = lepsubcluster.jets();
104
105      // Subtract the lepton four vector from closest subjet if dR<0.4
106      Jets subjets_had_clean;
107      double dRmin_had = 0.4;
108      unsigned int i_dRmin_had = 0;
109      bool found_match_had = false;
110      for(unsigned int i=0; i<subjets_had.size(); i++){
111        double dR = deltaR(subjets_had[i], lepton);
112        if(dR < dRmin_had){
113          dRmin_had = dR;
114          i_dRmin_had = i;
115          found_match_had = true;
116        }
117      }
118      for(unsigned int i=0; i<subjets_had.size(); i++){
119        Jet subjet = subjets_had[i];
120        if(found_match_had && i == i_dRmin_had) subjet = Jet(subjets_had[i].momentum()-lepton.momentum(), subjets_had[i].particles(), subjets_had[i].tags());
121        subjets_had_clean.push_back(subjet);
122      }
123      std::sort(subjets_had_clean.begin(), subjets_had_clean.end(), cmpMomByPt);
124
125      // do the same for lep jets
126      Jets subjets_lep_clean;
127      double dRmin_lep = 0.4;
128      unsigned int i_dRmin_lep = 0;
129      bool found_match_lep = false;
130      for(unsigned int i=0; i<subjets_lep.size(); i++){
131        double dR = deltaR(subjets_lep[i], lepton);
132        if(dR < dRmin_lep){
133          dRmin_lep = dR;
134          i_dRmin_lep = i;
135          found_match_lep = true;
136        }
137      }
138      for(unsigned int i=0; i<subjets_lep.size(); i++){
139        Jet subjet = subjets_lep[i];
140        if(found_match_lep && i == i_dRmin_lep) subjet = Jet(subjets_lep[i].momentum()-lepton.momentum(), subjets_lep[i].particles(), subjets_lep[i].tags());
141        subjets_lep_clean.push_back(subjet);
142      }
143      std::sort(subjets_lep_clean.begin(), subjets_lep_clean.end(), cmpMomByPt);
144
145      // Subjet cuts
146      if(subjets_had_clean.size() != 3) vetoEvent;
147      if(subjets_lep_clean.size() != 3) vetoEvent;
148      for (Jet jet : subjets_had_clean){
149          if(jet.pT() < 30*GeV) vetoEvent;
150          if(jet.abseta() > 2.5) vetoEvent;
151      }
152
153      // Combine subjets to final jets
154      FourMomentum hadjet;
155      for(Jet subjet : subjets_had_clean){
156        if(subjet.abseta() < 2.5) hadjet += subjet.momentum();
157      }
158      FourMomentum lepjet;
159      for(Jet subjet : subjets_lep_clean){
160        if(subjet.abseta() < 2.5) lepjet += subjet.momentum();
161      }
162
163      // Jet pT cuts
164      if(hadjet.pT() < 400*GeV) vetoEvent;
165      if(lepjet.pT() < 10*GeV) vetoEvent;
166
167      // m(hadjet) > m(lepjet+lepton)
168      FourMomentum secondJetLepton = lepjet + lepton.momentum();
169      if(hadjet.mass() < secondJetLepton.mass()) vetoEvent;
170
171      // Fill histograms
172      _hist_mass->fill(hadjet.mass()/GeV);
173      _hist_mass_norm->fill(hadjet.mass()/GeV);
174
175    }
176
177    /// Normalise and scale histograms
178    void finalize() {
179      const double sf = crossSection() / femtobarn / sumOfWeights();
180      scale(_hist_mass, sf);
181      normalize(_hist_mass_norm, 1.0, false);
182    }
183
184    //@}
185
186
187  private:
188
189    // Histograms
190    Histo1DPtr _hist_mass, _hist_mass_norm;
191
192  };
193
194
195
196  RIVET_DECLARE_PLUGIN(CMS_2019_I1764472);
197
198}