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/LeptonFinder.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      /// @todo Clean this up with something in-place
 41      fastjet::contrib::PseudoXConePlugin* plugin_xcone = new fastjet::contrib::PseudoXConePlugin(2, 1.2, 2.0);
 42      declare(FastJets(fs_jets, plugin_xcone), "FatJets");
 43
 44      // Partonic tops for decay channel definition
 45      declare(PartonicTops(TopDecay::E_MU, PromptEMuFromTau::NO), "LeptonicTops");
 46      declare(PartonicTops(TopDecay::HADRONIC), "HadronicTops");
 47
 48      // Book histograms
 49      book(_hist_mass, "d01-x01-y01");
 50      book(_hist_mass_norm, "d02-x01-y01");
 51    }
 52
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56
 57      // Decay mode check
 58      const Particles& leptonicTops = apply<PartonicTops>(event, "LeptonicTops").particlesByPt();
 59      const Particles& hadronicTops = apply<PartonicTops>(event, "HadronicTops").particlesByPt();
 60      if (leptonicTops.size() != 1 || hadronicTops.size() != 1) vetoEvent;
 61
 62
 63      // Get prompt leptons
 64      const PromptFinalState& prompt_leptons = apply<PromptFinalState>(event, "PromptLeptons");
 65      const Particles& leptons = prompt_leptons.particles();
 66      if (leptons.empty()) vetoEvent;
 67
 68      // Select leading lepton
 69      Particle lepton;
 70      for (const Particle& l : leptons){
 71        if (l.pT() > lepton.pT()) lepton = l;
 72      }
 73      if (lepton.pT() < 60*GeV) vetoEvent;
 74
 75      // Get the fat jets
 76      const Jets& fatjets = apply<FastJets>(event, "FatJets").jets();
 77
 78      // Get index of hadronic jet by distance to lepton
 79      int ihad = 0;
 80      int ilep = 1;
 81
 82      double dR0 = deltaR(lepton, fatjets.at(0));
 83      double dR1 = deltaR(lepton, fatjets.at(1));
 84
 85      if (dR0 < dR1) {
 86        ihad = 1;
 87        ilep = 0;
 88      }
 89
 90      // Get jet constituents
 91      const Particles& phad = fatjets.at(ihad).particles();
 92      const Particles& plep = fatjets.at(ilep).particles();
 93
 94      // Cluster subjets
 95      FinalState fs_dummy;
 96      /// @todo Avoid rebuilding the plugin for every event
 97      fastjet::JetDefinition::Plugin* plugin_subhad = new fastjet::contrib::PseudoXConePlugin(3, 0.4, 2.0);
 98      fastjet::contrib::PseudoXConePlugin* plugin_sublep = new fastjet::contrib::PseudoXConePlugin(3, 0.4, 2.0);
 99      FastJets hadsubcluster(fs_dummy, plugin_subhad);
100      FastJets lepsubcluster(fs_dummy, plugin_sublep);
101      hadsubcluster.calc(phad);
102      lepsubcluster.calc(plep);
103
104      Jets subjets_had = hadsubcluster.jets();
105      Jets subjets_lep = lepsubcluster.jets();
106
107      // Subtract the lepton four vector from closest subjet if dR<0.4
108      Jets subjets_had_clean;
109      double dRmin_had = 0.4;
110      unsigned int i_dRmin_had = 0;
111      bool found_match_had = false;
112      for (unsigned int i=0; i<subjets_had.size(); i++) {
113        double dR = deltaR(subjets_had[i], lepton);
114        if (dR < dRmin_had) {
115          dRmin_had = dR;
116          i_dRmin_had = i;
117          found_match_had = true;
118        }
119      }
120      for (unsigned int i=0; i<subjets_had.size(); i++) {
121        Jet subjet = subjets_had[i];
122        if (found_match_had && i == i_dRmin_had) {
123	  subjet = Jet(subjets_had[i].momentum()-lepton.momentum(), subjets_had[i].particles(), subjets_had[i].tags());
124	}
125        subjets_had_clean.push_back(subjet);
126      }
127      std::sort(subjets_had_clean.begin(), subjets_had_clean.end(), cmpMomByPt);
128
129      // do the same for lep jets
130      Jets subjets_lep_clean;
131      double dRmin_lep = 0.4;
132      unsigned int i_dRmin_lep = 0;
133      bool found_match_lep = false;
134      for (unsigned int i=0; i<subjets_lep.size(); i++) {
135        double dR = deltaR(subjets_lep[i], lepton);
136        if (dR < dRmin_lep) {
137          dRmin_lep = dR;
138          i_dRmin_lep = i;
139          found_match_lep = true;
140        }
141      }
142      for (unsigned int i=0; i<subjets_lep.size(); i++) {
143        Jet subjet = subjets_lep[i];
144        if (found_match_lep && i == i_dRmin_lep) {
145	  subjet = Jet(subjets_lep[i].momentum()-lepton.momentum(), subjets_lep[i].particles(), subjets_lep[i].tags());
146	}
147        subjets_lep_clean.push_back(subjet);
148      }
149      std::sort(subjets_lep_clean.begin(), subjets_lep_clean.end(), cmpMomByPt);
150
151      // Subjet cuts
152      if (subjets_had_clean.size() != 3) vetoEvent;
153      if (subjets_lep_clean.size() != 3) vetoEvent;
154      for (Jet jet : subjets_had_clean){
155	if (jet.pT() < 30*GeV) vetoEvent;
156	if (jet.abseta() > 2.5) vetoEvent;
157      }
158
159      // Combine subjets to final jets
160      FourMomentum hadjet;
161      for (Jet subjet : subjets_had_clean) {
162        if (subjet.abseta() < 2.5) hadjet += subjet.momentum();
163      }
164      FourMomentum lepjet;
165      for (Jet subjet : subjets_lep_clean) {
166        if (subjet.abseta() < 2.5) lepjet += subjet.momentum();
167      }
168
169      // Jet pT cuts
170      if (hadjet.pT() < 400*GeV) vetoEvent;
171      if (lepjet.pT() < 10*GeV) vetoEvent;
172
173      // m(hadjet) > m(lepjet+lepton)
174      FourMomentum secondJetLepton = lepjet + lepton.momentum();
175      if (hadjet.mass() < secondJetLepton.mass()) vetoEvent;
176
177      // Fill histograms
178      _hist_mass->fill(hadjet.mass()/GeV);
179      _hist_mass_norm->fill(hadjet.mass()/GeV);
180    }
181    
182
183    /// Normalise and scale histograms
184    void finalize() {
185      const double sf = crossSection() / femtobarn / sumOfWeights();
186      scale(_hist_mass, sf);
187      normalize(_hist_mass_norm, 1.0, false);
188    }
189
190    /// @}
191
192
193  private:
194
195    // Histograms
196    Histo1DPtr _hist_mass, _hist_mass_norm;
197
198  };
199
200
201
202  RIVET_DECLARE_PLUGIN(CMS_2019_I1764472);
203
204}