rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2019_I1750330

Semileptonic ttbar at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1750330
Status: VALIDATED
Authors:
  • Federica Fabbri
  • Francesco La Ruffa
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • pp -> non-allhadronic ttbar production at 13 TeV

Single- and double-differential cross-section measurements are presented for the production of top-quark pairs, in the lepton + jets channel at particle and parton level. Two topologies, resolved and boosted, are considered and the results are presented as a function of several kinematic variables characterising the top and the system and jet multiplicities. The study was performed using data from $pp$ collisions at centre-of-mass energy of 13 TeV collected in 2015 and 2016 by the ATLAS detector at the CERN Large Hadron Collider (LHC), corresponding to an integrated luminosity of 36 fb$^{-1}$. Due to the large $t\bar{t}$ cross-section at the LHC, such measurements allow a detailed study of the properties of top-quark production and decay, enabling precision tests of several Monte Carlo generators and fixed-order Standard Model predictions. Overall, there is good agreement between the theoretical predictions and the data.

Source code: ATLAS_2019_I1750330.cc
  1// -*- C++ -*
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/InvisibleFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8#include "Rivet/Projections/MissingMomentum.hh"
  9#include "Rivet/Projections/FastJets.hh"
 10#include "Rivet/Tools/BinnedHistogram.hh"
 11
 12
 13namespace Rivet {
 14
 15  /// @brief Semileptonic ttbar at 13 TeV
 16  class ATLAS_2019_I1750330 : public Analysis {
 17  public:
 18
 19    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2019_I1750330);
 20
 21    void init() {
 22
 23      _doBoosted = true, _doResolved = true;
 24      if ( getOption("TYPE") == "BOOSTED" )  _doResolved = false;
 25      else if ( getOption("TYPE") == "RESOLVED" )  _doBoosted = false;
 26
 27      Cut eta_full =  (Cuts::abseta < 5.0);
 28      Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 27*GeV);
 29      const FinalState fs(eta_full);
 30
 31      FinalState all_photons(fs, Cuts::abspid == PID::PHOTON);
 32
 33      PromptFinalState photons(all_photons, false);
 34      declare(photons, "photons");
 35
 36      PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, true);
 37      declare(electrons, "electrons");
 38
 39      DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true);
 40      declare(dressedelectrons, "dressedelectrons");
 41
 42      DressedLeptons ewdressedelectrons(all_photons, electrons, 0.1, eta_full, true);
 43      declare(ewdressedelectrons, "ewdressedelectrons");
 44
 45      PromptFinalState muons(Cuts::abspid == PID::MUON, true);
 46      declare(muons, "muons");
 47
 48      DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true);
 49      declare(dressedmuons, "dressedmuons");
 50
 51      DressedLeptons ewdressedmuons(all_photons, muons, 0.1, eta_full, true);
 52      declare(ewdressedmuons, "ewdressedmuons");
 53
 54      const InvisibleFinalState neutrinos(true, true);
 55
 56      VetoedFinalState vfs(fs);
 57      vfs.addVetoOnThisFinalState(dressedelectrons);
 58      vfs.addVetoOnThisFinalState(dressedmuons);
 59      vfs.addVetoOnThisFinalState(neutrinos);
 60      FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::ALL);
 61      declare(jets, "boosted_jets");
 62
 63      VetoedFinalState vfs_res(fs);
 64      vfs_res.addVetoOnThisFinalState(ewdressedelectrons);
 65      vfs_res.addVetoOnThisFinalState(ewdressedmuons);
 66      vfs_res.addVetoOnThisFinalState(neutrinos);
 67      FastJets jets_res(vfs_res, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::ALL);
 68      declare(jets_res, "resolved_jets");
 69
 70      declare(MissingMomentum(), "MissingMomentum");
 71
 72      // Bins for 2D resolved
 73      std::vector<double> ttbar_m_2D_bins = {200,400,550,700,1000,2000};
 74      std::vector<double> top_had_pt_2D_bins = {0,60,120,200,300,1000};
 75      std::vector<double> ttbar_pt_2D_bins = {0,30,80,190,800};
 76      std::vector<double> top_had_abs_y_2D_bins = {0,0.7,1.4,2.5};
 77      std::vector<double> ttbar_abs_y_2D_bins = {0.0,0.4,0.8,1.2,2.5};
 78
 79      std::vector<double> n_jet_bins = {3.5, 4.5, 5.5, 6.5, 7.5};
 80      std::vector<double> n_jet_bins_for_ttbar_m = {3.5, 4.5, 5.5, 6.5};
 81      std::vector<double> n_extrajet_bins = {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5};
 82
 83      //Bins for 2D boosted
 84      std::vector<double> eta_2D_bins = {0,1,2};
 85      std::vector<double> etattbar_2D_bins = {0,60,120,200,300,1000};
 86      std::vector<double> pttbar_2D_bins = {0.0, 40.0, 150.0, 1000.0};
 87      std::vector<double> mtt_2D_bins = {490.0, 1160, 3000.0};
 88      std::vector<double> eta_external_2D_bins = {0.0, 0.65, 1.3, 2.0};
 89      std::vector<double> ptt_external_mtt_2D_bins = {0.0, 40.0, 150.0, 1000.0};
 90      std::vector<double> Htt_external_2D_bins = {350.0, 780.0, 2500.0};
 91      std::vector<double> eta_external_ptt_2D_bins = {0.0, 0.65, 2.0};
 92
 93      std::vector<double> n_jet_pttop_bins = {-0.5,1.5,2.5,3.5};
 94      std::vector<double> n_jet_ptttbar_bins = {-0.5,1.5,3.5};
 95      std::vector<double> n_jet_Pout_bins = {-0.5,1.5,3.5};
 96      std::vector<double> n_jet_mtt_bins = {-0.5,0.5,1.5,2.5};
 97
 98      //Resolved histograms (digits correspond to "Table ID" from HepData)
 99      book2D("ttbar_m_top_had_pt_multi_norm", ttbar_m_2D_bins,54);
100      book2D("ttbar_m_top_had_pt_multi", ttbar_m_2D_bins, 74);
101
102      book2D("ttbar_m_ttbar_pt_multi_norm", ttbar_m_2D_bins, 94);
103      book2D("ttbar_m_ttbar_pt_multi", ttbar_m_2D_bins, 114);
104
105      book2D("top_had_pt_absPout_multi_norm", top_had_pt_2D_bins, 134);
106      book2D("top_had_pt_absPout_multi", top_had_pt_2D_bins,154);
107
108      book2D("top_had_pt_jet_n_multi_norm", n_jet_bins,174);
109      book2D("top_had_pt_jet_n_multi", n_jet_bins, 188 );
110
111      book2D("ttbar_m_jet_n_multi_norm", n_jet_bins_for_ttbar_m,202);
112      book2D("ttbar_m_jet_n_multi", n_jet_bins_for_ttbar_m,211);
113
114      book2D("ttbar_pt_jet_n_multi_norm", n_jet_bins, 220 );
115      book2D("ttbar_pt_jet_n_multi", n_jet_bins, 234 );
116
117      book2D("absPout_jet_n_multi_norm", n_jet_bins, 248 );
118      book2D("absPout_jet_n_multi", n_jet_bins, 262 );
119
120      book2D("deltaPhi_tt_jet_n_multi_norm", n_jet_bins, 276 );
121      book2D("deltaPhi_tt_jet_n_multi", n_jet_bins, 290 );
122
123      book2D("HT_tt_jet_n_multi_norm", n_jet_bins, 304 );
124      book2D("HT_tt_jet_n_multi", n_jet_bins, 318 );
125
126      book2D("top_had_abs_y_jet_n_multi_norm", n_jet_bins, 332 );
127      book2D("top_had_abs_y_jet_n_multi", n_jet_bins, 346 );
128
129      book2D("ttbar_abs_y_jet_n_multi_norm", n_jet_bins, 360 );
130      book2D("ttbar_abs_y_jet_n_multi", n_jet_bins, 374 );
131
132      book2D("chi_tt_jet_n_multi_norm", n_jet_bins, 388 );
133      book2D("chi_tt_jet_n_multi", n_jet_bins, 402 );
134
135      book2D("top_had_abs_y_top_had_pt_multi_norm", top_had_abs_y_2D_bins,416 );
136      book2D("top_had_abs_y_top_had_pt_multi", top_had_abs_y_2D_bins, 425);
137
138      book2D("ttbar_abs_y_ttbar_pt_multi_norm", ttbar_abs_y_2D_bins, 434);
139      book2D("ttbar_abs_y_ttbar_pt_multi", ttbar_abs_y_2D_bins, 448);
140
141      book2D("ttbar_abs_y_ttbar_m_multi_norm", ttbar_abs_y_2D_bins, 462);
142      book2D("ttbar_abs_y_ttbar_m_multi", ttbar_abs_y_2D_bins, 476);
143
144      book2D("ttbar_pt_top_had_pt_multi_norm", ttbar_pt_2D_bins, 490);
145      book2D("ttbar_pt_top_had_pt_multi", ttbar_pt_2D_bins, 504);
146
147      book_hist("top_had_pt",1);
148      book_hist("top_had_abs_y_fine",5);
149      book_hist("leading_top_pt",9);
150      book_hist("subleading_top_pt",13);
151      book_hist("ttbar_m",17);
152      book_hist("ttbar_pt",21);
153      book_hist("absPout",25);
154      book_hist("deltaPhi_tt",29);
155      book_hist("HT_tt",33);
156      book_hist("extrajet_n",37);
157      book_hist("ttbar_abs_y_fine",41);
158      book_hist("abs_y_boost",45);
159      book_hist("chi_tt",49);
160
161      //Boosted histograms (digits correspond to "Table ID" from HepData)
162      book2D("boosted_rc_pttop_etatop_multi", eta_2D_bins, 922);
163      book2D("boosted_rc_pttop_etattbar_multi", eta_2D_bins, 912);
164      book2D("boosted_rc_pttop_ptttbar_multi", pttbar_2D_bins, 898);
165      book2D("boosted_rc_pttop_mttbar_multi", mtt_2D_bins, 932);
166      book2D("boosted_rc_mttbar_etattbar_multi", eta_external_2D_bins, 974);
167      book2D("boosted_rc_mttbar_ptttbar_multi", ptt_external_mtt_2D_bins, 956);
168      book2D("boosted_rc_mttbar_HT_multi", Htt_external_2D_bins, 942);
169      book2D("boosted_rc_pttop_extrajet_multi", n_jet_pttop_bins, 992);
170      book2D("boosted_rc_ptttbar_extrajet_multi", n_jet_ptttbar_bins, 1006);
171      book2D("boosted_rc_mttbar_extrajet_multi", n_jet_mtt_bins, 1020);
172
173      book2D("boosted_rc_pttop_etatop_multi_norm", eta_2D_bins, 917);
174      book2D("boosted_rc_pttop_etattbar_multi_norm", eta_2D_bins, 907);
175      book2D("boosted_rc_pttop_ptttbar_multi_norm", pttbar_2D_bins, 889);
176      book2D("boosted_rc_pttop_mttbar_multi_norm", mtt_2D_bins, 927);
177      book2D("boosted_rc_mttbar_etattbar_multi_norm", eta_external_2D_bins, 965);
178      book2D("boosted_rc_mttbar_ptttbar_multi_norm", ptt_external_mtt_2D_bins, 947);
179      book2D("boosted_rc_mttbar_HT_multi_norm", Htt_external_2D_bins, 937);
180      book2D("boosted_rc_pttop_extrajet_multi_norm", n_jet_pttop_bins, 983);
181      book2D("boosted_rc_ptttbar_extrajet_multi_norm", n_jet_ptttbar_bins, 1001);
182      book2D("boosted_rc_mttbar_extrajet_multi_norm", n_jet_mtt_bins, 1011);
183
184      book_hist("hadTop_boosted_rc_pt",840);
185      book_hist("hadTop_boosted_rc_y",844);
186      book_hist("LeadingTop_boosted_rc_pt",848);
187      book_hist("SubLeadingTop_boosted_rc_pt",852);
188      book_hist("boosted_rc_Pout_lep",872);
189      book_hist("boosted_rc_chi_tt",868);
190      book_hist("boosted_rc_HT",876);
191      book_hist("hadTop_boosted_rc_subjets",884);
192      book_hist("boosted_rc_extrajet",880);
193      book_hist("ttbar_boosted_rc_m",864);
194      book_hist("ttbar_boosted_rc_pt",856);
195      book_hist("ttbar_boosted_rc_Rapidity",860);
196
197    }
198
199
200    void analyze(const Event& event) {
201      if (_doResolved)  Resolved_selection(event);
202      if (_doBoosted)   Boosted_selection(event);
203    }
204
205
206    void Resolved_selection(const Event& event) {
207
208      // Get the selected objects, using the projections.
209      vector<DressedLepton> electrons = apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
210      vector<DressedLepton> muons     = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
211      const Jets& jets = apply<FastJets>(event, "resolved_jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
212      FourMomentum met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
213
214      Jets bjets, lightjets;
215
216      // OVERLAP REMOVAL
217      idiscardIfAnyDeltaRLess(muons, jets, 0.4);
218      idiscardIfAnyDeltaRLess(electrons, jets, 0.4);
219
220      // b-tagging
221      // If there are more than 2 b-tagged jets, the extra b-tagged jets will be treat as light jets
222      for (const Jet& jet : jets) {
223        bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
224        if ( b_tagged && bjets.size() < 2)  bjets +=jet;
225        else lightjets += jet;
226      }
227
228      bool single_electron = electrons.size() == 1 && muons.empty();
229      bool single_muon     = muons.size() == 1 && electrons.empty();
230
231      DressedLepton *lepton = NULL;
232      if (single_electron)   lepton = &electrons[0];
233      else if (single_muon)  lepton = &muons[0];
234
235      if (!single_electron && !single_muon) vetoEvent;
236      bool num_b_tagged_jets = (bjets.size() == 2);
237      if (!num_b_tagged_jets) vetoEvent;
238
239      if (lightjets.size() < 2) vetoEvent;
240
241      FourMomentum pbjet1; //Momentum of bjet1
242      FourMomentum pbjet2; //Momentum of bjet
243      if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) {
244        pbjet1 = bjets[0].momentum();
245        pbjet2 = bjets[1].momentum();
246      } else {
247        pbjet1 = bjets[1].momentum();
248        pbjet2 = bjets[0].momentum();
249      }
250
251      double bestWmass = 1000.0*TeV;
252      double mWPDG = 80.399*GeV;
253      int Wj1index = -1, Wj2index = -1;
254      for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) {
255        for (unsigned int j = i + 1; j < lightjets.size(); ++j) {
256          double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass();
257          if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) {
258            bestWmass = wmass;
259            Wj1index = i;
260            Wj2index = j;
261          }
262        }
263      }
264
265      FourMomentum pjet1 = lightjets[Wj1index].momentum();
266      FourMomentum pjet2 = lightjets[Wj2index].momentum();
267
268      // compute hadronic W boson
269      FourMomentum pWhadron = pjet1 + pjet2;
270      double pz = computeneutrinoz(lepton->momentum(), met);
271      FourMomentum ppseudoneutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz);
272
273      //compute leptonic, hadronic, combined pseudo-top
274      FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1;
275      FourMomentum ppseudotophadron = pbjet2 + pWhadron;
276      FourMomentum pttbar = ppseudotoplepton + ppseudotophadron;
277
278      Vector3 z_versor(0,0,1);
279      Vector3 vpseudotophadron = ppseudotophadron.vector3();
280      Vector3 vpseudotoplepton = ppseudotoplepton.vector3();
281
282      // Variables
283      double ystar = (ppseudotophadron.pt() > ppseudotoplepton.pt()) ? 0.5 * (ppseudotophadron.rap()-ppseudotoplepton.rap()) : 0.5*(ppseudotoplepton.rap()-ppseudotophadron.rap());
284      double chi_ttbar = exp(2 * fabs(ystar));
285      double deltaPhi_ttbar = deltaPhi(ppseudotoplepton,ppseudotophadron);
286      double HT_ttbar = ppseudotophadron.pt() + ppseudotoplepton.pt();
287      double Yboost = 0.5 * (ppseudotophadron.rapidity() + ppseudotoplepton.rapidity());
288      double Pout = vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod()));
289      double absPout = fabs(Pout);
290      double Leading_top_pt = (ppseudotophadron.pt() > ppseudotoplepton.pt()) ? ppseudotophadron.pt() : ppseudotoplepton.pt();
291      double Subleading_top_pt = (ppseudotophadron.pt() > ppseudotoplepton.pt()) ? ppseudotoplepton.pt() : ppseudotophadron.pt();
292      int jet_multiplicity = jets.size();
293      int extrajet_n = jet_multiplicity - 4;
294      int new_jet_multi = TransformJetMultiplicity(jet_multiplicity);
295      int new_jet_multi_for_ttbar_m = TransformJetMultiplicity_for_ttbar_m(jet_multiplicity);
296      int new_extrajet_multi =TransformExtrajetMultiplicity(extrajet_n);
297
298      _h_multi["top_had_pt_Pout_multi"].fill(ppseudotophadron.pt()/GeV, Pout);
299      _h_multi["top_had_pt_absPout_multi"].fill(ppseudotophadron.pt()/GeV, absPout);
300      _h_multi["ttbar_m_top_had_pt_multi"].fill(pttbar.mass()/GeV, ppseudotophadron.pt()/GeV);
301      _h_multi["ttbar_m_ttbar_pt_multi"].fill(pttbar.mass()/GeV, pttbar.pt()/GeV);
302      _h_multi["ttbar_pt_top_had_pt_multi"].fill(pttbar.pt()/GeV, ppseudotophadron.pt()/GeV);
303      _h_multi["ttbar_abs_y_ttbar_pt_multi"].fill(pttbar.absrap(), pttbar.pt()/GeV);
304      _h_multi["ttbar_abs_y_ttbar_m_multi"].fill(pttbar.absrap(), pttbar.mass()/GeV);
305      _h_multi["top_had_abs_y_top_had_pt_multi"].fill(ppseudotophadron.absrap(), ppseudotophadron.pt()/GeV);
306
307      _h_multi["ttbar_pt_jet_n_multi"].fill(new_jet_multi, pttbar.pt()/GeV);
308      _h_multi["ttbar_m_jet_n_multi"].fill(new_jet_multi_for_ttbar_m, pttbar.mass()/GeV);
309      _h_multi["chi_tt_jet_n_multi"].fill(new_jet_multi, chi_ttbar);
310      _h_multi["absPout_jet_n_multi"].fill(new_jet_multi, absPout);
311      _h_multi["deltaPhi_tt_jet_n_multi"].fill(new_jet_multi, deltaPhi_ttbar);
312      _h_multi["HT_tt_jet_n_multi"].fill(new_jet_multi, HT_ttbar/GeV);
313      _h_multi["top_had_pt_jet_n_multi"].fill(new_jet_multi, ppseudotophadron.pt()/GeV);
314      _h_multi["top_had_abs_y_jet_n_multi"].fill(new_jet_multi, ppseudotophadron.absrap());
315      _h_multi["ttbar_abs_y_jet_n_multi"].fill(new_jet_multi, pttbar.absrap());
316
317      _h_multi["top_had_pt_Pout_multi_norm"].fill(ppseudotophadron.pt()/GeV, Pout);
318      _h_multi["top_had_pt_absPout_multi_norm"].fill(ppseudotophadron.pt()/GeV, absPout);
319      _h_multi["ttbar_m_top_had_pt_multi_norm"].fill(pttbar.mass()/GeV, ppseudotophadron.pt()/GeV);
320      _h_multi["ttbar_m_ttbar_pt_multi_norm"].fill(pttbar.mass()/GeV, pttbar.pt()/GeV);
321      _h_multi["ttbar_pt_top_had_pt_multi_norm"].fill(pttbar.pt()/GeV, ppseudotophadron.pt()/GeV);
322      _h_multi["ttbar_abs_y_ttbar_pt_multi_norm"].fill(pttbar.absrap(), pttbar.pt()/GeV);
323      _h_multi["ttbar_abs_y_ttbar_m_multi_norm"].fill(pttbar.absrap(), pttbar.mass()/GeV);
324      _h_multi["top_had_abs_y_top_had_pt_multi_norm"].fill(ppseudotophadron.absrap(), ppseudotophadron.pt()/GeV);
325
326      _h_multi["ttbar_pt_jet_n_multi_norm"].fill(new_jet_multi, pttbar.pt());
327      _h_multi["ttbar_m_jet_n_multi_norm"].fill(new_jet_multi_for_ttbar_m, pttbar.mass());
328      _h_multi["chi_tt_jet_n_multi_norm"].fill(new_jet_multi, chi_ttbar);
329      _h_multi["absPout_jet_n_multi_norm"].fill(new_jet_multi, absPout);
330      _h_multi["deltaPhi_tt_jet_n_multi_norm"].fill(new_jet_multi, deltaPhi_ttbar);
331      _h_multi["HT_tt_jet_n_multi_norm"].fill(new_jet_multi, HT_ttbar/GeV);
332      _h_multi["top_had_pt_jet_n_multi_norm"].fill(new_jet_multi, ppseudotophadron.pt()/GeV);
333      _h_multi["top_had_abs_y_jet_n_multi_norm"].fill(new_jet_multi, ppseudotophadron.absrap());
334      _h_multi["ttbar_abs_y_jet_n_multi_norm"].fill(new_jet_multi, pttbar.absrap());
335
336      _h["chi_tt"]->fill(chi_ttbar);
337      _h["deltaPhi_tt"]->fill(deltaPhi_ttbar);
338      _h["HT_tt"]->fill(HT_ttbar/GeV);
339      _h["absPout"]->fill(absPout);
340      _h["abs_y_boost"]->fill(fabs(Yboost));
341      _h["top_had_pt"]->fill(ppseudotophadron.pt()/GeV);
342      _h["top_had_abs_y_fine"]->fill(ppseudotophadron.absrap());
343      _h["ttbar_pt"]->fill(pttbar.pt()/GeV);
344      _h["ttbar_m"]->fill(pttbar.mass()/GeV);
345      _h["ttbar_abs_y_fine"]->fill(pttbar.absrap());
346      _h["leading_top_pt"]->fill(Leading_top_pt/GeV);
347      _h["subleading_top_pt"]->fill(Subleading_top_pt/GeV);
348      _h["extrajet_n"]->fill(new_extrajet_multi+1);
349
350      _h["chi_tt_norm"]->fill(chi_ttbar);
351      _h["deltaPhi_tt_norm"]->fill(deltaPhi_ttbar);
352      _h["HT_tt_norm"]->fill(HT_ttbar/GeV);
353      _h["absPout_norm"]->fill(absPout);
354      _h["abs_y_boost_norm"]->fill(fabs(Yboost));
355      _h["top_had_pt_norm"]->fill(ppseudotophadron.pt()/GeV);
356      _h["top_had_abs_y_fine_norm"]->fill(ppseudotophadron.absrap());
357      _h["ttbar_pt_norm"]->fill(pttbar.pt()/GeV);
358      _h["ttbar_m_norm"]->fill(pttbar.mass()/GeV);
359      _h["ttbar_abs_y_fine_norm"]->fill(pttbar.absrap());
360      _h["leading_top_pt_norm"]->fill(Leading_top_pt/GeV);
361      _h["subleading_top_pt_norm"]->fill(Subleading_top_pt/GeV);
362      _h["extrajet_n_norm"]->fill(new_extrajet_multi+1);
363    }
364
365    void Boosted_selection(const Event& event) {
366
367      //Projections
368      vector<DressedLepton> electrons = apply<DressedLeptons>(event, "dressedelectrons").dressedLeptons();
369      vector<DressedLepton> muons = apply<DressedLeptons>(event, "dressedmuons").dressedLeptons();
370      const Jets& jets = apply<FastJets>(event, "boosted_jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta <= 2.5);
371      const FourMomentum& met = apply<MissingMomentum>(event, "MissingMomentum").missingMomentum();
372
373      if (jets.size() < 2)  vetoEvent;
374      PseudoJets smallRjets;
375      for (const Jet& jet : jets) {
376        smallRjets += jet.pseudojet();
377        bool b_tagged = jet.bTagged(Cuts::pT > 5*GeV);
378        smallRjets[smallRjets.size()-1].set_user_index(b_tagged); // cheeky, but works
379      }
380
381      idiscardIfAnyDeltaRLess(muons, jets, 0.4);
382      idiscardIfAnyDeltaRLess(electrons, jets, 0.4);
383
384      fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0), fastjet::SelectorPtFractionMin(0.05));
385      fastjet::ClusterSequence antikt_cs(smallRjets, fastjet::JetDefinition(fastjet::antikt_algorithm, 1.0));
386      PseudoJets reclustered_jets = antikt_cs.inclusive_jets();
387
388      // trim the jets
389      Jets TrimmedJets;
390      for (PseudoJet pjet : reclustered_jets) {
391        PseudoJet ptrim = trimmer(pjet);
392        if (ptrim.perp() < 350*GeV)  continue;
393        if (fabs(ptrim.eta()) > 2.0) continue;
394        bool bTagged  = false;
395        Particles constituents;
396        for (const PseudoJet& c : ptrim.constituents()) {
397          // we only care about the number of subjets, so
398          // fine to treat as Particles with dummy PID
399          constituents += Particle(0, momentum(c));
400          bTagged |= c.user_index();
401        }
402        ptrim.set_user_index(bTagged);
403        TrimmedJets += Jet(ptrim, constituents);
404      }
405      Cut trim_selection = Cuts::abseta < 2.0 && Cuts::pT > 200*GeV && Cuts::massIn(120*GeV, 220*GeV);
406      ifilter_select(isortByPt(TrimmedJets), trim_selection);
407      if (TrimmedJets.empty())  vetoEvent;
408
409
410      // SINGLE LEPTON
411      bool single_electron=(electrons.size() == 1) && (muons.empty());
412      bool single_muon=(muons.size() == 1) && (electrons.empty());
413
414      DressedLepton *lepton = NULL;
415      if (single_electron)   lepton = &electrons[0];
416      else if (single_muon)  lepton = &muons[0];
417      if (!single_electron && !single_muon) vetoEvent;
418
419      //MET
420      if (met.pT() < 20*GeV)  vetoEvent;
421
422      //MET+MWT
423      double transmass = TransMass(lepton->pt(), lepton->phi(), met.pt(), met.phi());
424      if ((met.pT() + transmass) < 60*GeV)   vetoEvent;
425
426      size_t subjets = 0;
427      bool btag_hadside=false;
428      bool hasHadTopCandidate = false;
429      FourMomentum HadTopCandidate;
430      for (const Jet& rc_jet : TrimmedJets) {
431        FourMomentum rc_jet_mom = rc_jet.mom();
432        if (rc_jet_mom.pt() < 350*GeV)   continue;
433        double dPhi_lepJet = fabs(deltaPhi(rc_jet_mom.phi(), lepton->phi()));
434        if (dPhi_lepJet < 1.)  continue;
435        if (rc_jet.pseudojet().user_index()) {
436          btag_hadside=true;
437        }
438        HadTopCandidate = momentum(rc_jet);
439        subjets = rc_jet.constituents().size();
440        hasHadTopCandidate = true;
441        break;
442      }
443      if (!hasHadTopCandidate)  vetoEvent;
444
445      Jets LepTopCandidates = filter_discard(jets, [&](const Jet& j) {
446          return deltaR(j, HadTopCandidate) < 1.5 || deltaR(j, *lepton) > 2.0;
447        });
448      if (LepTopCandidates.empty()) vetoEvent;
449
450      FourMomentum ltop;
451      bool btag_lepside=false;
452      for (const Jet& jet : LepTopCandidates) {
453        if (jet.bTagged(Cuts::pT > 5*GeV)) {
454          btag_lepside = true;
455          ltop = jet.mom();
456          break;
457        }
458      }
459      if (!btag_hadside && !btag_lepside)  vetoEvent;
460      if (!btag_lepside)   ltop = LepTopCandidates[0].momentum();
461      double pz = computeneutrinoz(lepton->momentum(), met);
462      FourMomentum neutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz);
463      FourMomentum LeptonicTop = lepton->momentum() + neutrino + ltop;
464      FourMomentum HadronicTop = HadTopCandidate;
465      FourMomentum pttbar = HadronicTop + LeptonicTop;
466
467      Vector3 z_versor(0,0,1);
468      Vector3 vpseudotophadron = HadronicTop.vector3();
469      Vector3 vpseudotoplepton = LeptonicTop.vector3();
470      // Variables
471      double ystar = (HadronicTop.pt() > LeptonicTop.pt()) ? 0.5 * (HadronicTop.rap()-LeptonicTop.rap()) : 0.5*(LeptonicTop.rap()-HadronicTop.rap());
472      double chi_ttbar = exp(2 * fabs(ystar));
473      double pt_leading = (HadronicTop.pt() > LeptonicTop.pt()) ? HadronicTop.pt() : LeptonicTop.pt();
474      double pt_subleading = (HadronicTop.pt() > LeptonicTop.pt()) ? LeptonicTop.pt() : HadronicTop.pt();
475      double HT_ttbar = HadronicTop.pt() + LeptonicTop.pt();
476      double absPout_lep = fabs(vpseudotoplepton.dot((vpseudotophadron.cross(z_versor))/(vpseudotophadron.cross(z_versor).mod())));
477      size_t extrajet = smallRjets.size() - subjets -1;
478
479      size_t new_subjets_multi = TransformExtrajetMultiplicity_boosted(subjets);
480      size_t new_extrajet_multi = TransformExtrajetMultiplicity_boosted(extrajet);
481      size_t new_extrajet_multi_pttop = TransformJetMultiplicity_pttop(extrajet);
482      size_t new_extrajet_multi_ptttbar = TransformJetMultiplicity_ptttbar(extrajet);
483      size_t new_extrajet_multi_mttbar = TransformJetMultiplicity_mttbar(extrajet);
484
485      _h_multi["boosted_rc_pttop_etatop_multi"].fill(HadronicTop.absrap(), HadronicTop.pt()/GeV);
486      _h_multi["boosted_rc_pttop_etattbar_multi"].fill(pttbar.absrap(), HadronicTop.pt()/GeV);
487      _h_multi["boosted_rc_pttop_ptttbar_multi"].fill(pttbar.pt()/GeV, HadronicTop.pt()/GeV);
488      _h_multi["boosted_rc_pttop_mttbar_multi"].fill(pttbar.mass()/GeV, HadronicTop.pt()/GeV);
489      _h_multi["boosted_rc_mttbar_etattbar_multi"].fill(pttbar.absrap(), pttbar.mass()/GeV);
490      _h_multi["boosted_rc_mttbar_ptttbar_multi"].fill(pttbar.pt()/GeV, pttbar.mass()/GeV);
491      _h_multi["boosted_rc_mttbar_HT_multi"].fill(HT_ttbar, pttbar.mass()/GeV);
492
493      _h_multi["boosted_rc_pttop_extrajet_multi"].fill(new_extrajet_multi_pttop, HadronicTop.pt()/GeV);
494      _h_multi["boosted_rc_ptttbar_extrajet_multi"].fill(new_extrajet_multi_ptttbar, pttbar.pt()/GeV);
495      _h_multi["boosted_rc_mttbar_extrajet_multi"].fill(new_extrajet_multi_mttbar, pttbar.mass()/GeV);
496
497      _h_multi["boosted_rc_pttop_etatop_multi_norm"].fill(HadronicTop.absrap(), HadronicTop.pt()/GeV);
498      _h_multi["boosted_rc_pttop_etattbar_multi_norm"].fill(pttbar.absrap(), HadronicTop.pt()/GeV);
499      _h_multi["boosted_rc_pttop_ptttbar_multi_norm"].fill(pttbar.pt()/GeV, HadronicTop.pt()/GeV);
500      _h_multi["boosted_rc_pttop_mttbar_multi_norm"].fill(pttbar.mass()/GeV, HadronicTop.pt()/GeV);
501      _h_multi["boosted_rc_mttbar_etattbar_multi_norm"].fill(pttbar.absrap(), pttbar.mass()/GeV);
502      _h_multi["boosted_rc_mttbar_ptttbar_multi_norm"].fill(pttbar.pt()/GeV, pttbar.mass()/GeV);
503      _h_multi["boosted_rc_mttbar_HT_multi_norm"].fill(HT_ttbar/GeV, pttbar.mass()/GeV);
504      _h_multi["boosted_rc_pttop_extrajet_multi_norm"].fill(new_extrajet_multi_pttop, HadronicTop.pt()/GeV);
505      _h_multi["boosted_rc_ptttbar_extrajet_multi_norm"].fill(new_extrajet_multi_ptttbar, pttbar.pt()/GeV);
506      _h_multi["boosted_rc_mttbar_extrajet_multi_norm"].fill(new_extrajet_multi_mttbar, pttbar.mass()/GeV);
507
508      _h["hadTop_boosted_rc_pt"]->fill(HadronicTop.pt()/GeV);
509      _h["hadTop_boosted_rc_y"]->fill(HadronicTop.absrap());
510      _h["LeadingTop_boosted_rc_pt"]->fill(pt_leading/GeV);
511      _h["SubLeadingTop_boosted_rc_pt"]->fill(pt_subleading/GeV);
512      _h["boosted_rc_Pout_lep"]->fill(absPout_lep);
513      _h["boosted_rc_chi_tt"]->fill(chi_ttbar);
514      _h["boosted_rc_HT"]->fill(HT_ttbar/GeV);
515      _h["hadTop_boosted_rc_subjets"]->fill(new_subjets_multi);
516      _h["boosted_rc_extrajet"]->fill(new_extrajet_multi+1);
517      _h["ttbar_boosted_rc_m"]->fill(pttbar.mass()/GeV);
518      _h["ttbar_boosted_rc_pt"]->fill(pttbar.pt()/GeV);
519      _h["ttbar_boosted_rc_Rapidity"]->fill(pttbar.absrapidity());
520
521      _h["hadTop_boosted_rc_pt_norm"]->fill(HadronicTop.pt()/GeV);
522      _h["hadTop_boosted_rc_y_norm"]->fill(HadronicTop.absrap());
523      _h["LeadingTop_boosted_rc_pt_norm"]->fill(pt_leading/GeV);
524      _h["SubLeadingTop_boosted_rc_pt_norm"]->fill(pt_subleading/GeV);
525      _h["boosted_rc_Pout_lep_norm"]->fill(absPout_lep);
526      _h["boosted_rc_chi_tt_norm"]->fill(chi_ttbar);
527      _h["boosted_rc_HT_norm"]->fill(HT_ttbar/GeV);
528      _h["hadTop_boosted_rc_subjets_norm"]->fill(new_subjets_multi);
529      _h["boosted_rc_extrajet_norm"]->fill(new_extrajet_multi+1);
530      _h["ttbar_boosted_rc_m_norm"]->fill(pttbar.mass()/GeV);
531      _h["ttbar_boosted_rc_pt_norm"]->fill(pttbar.pt()/GeV);
532      _h["ttbar_boosted_rc_Rapidity_norm"]->fill(pttbar.absrap()) ;
533    }
534
535
536    void finalize() {
537      // Normalize to cross-section
538      const double sf = crossSection() / sumOfWeights();
539      for (auto& hit : _h) {
540        scale(hit.second, sf);
541        if (hit.first.find("_norm") != string::npos)  normalize(hit.second, 1.0, false);
542      }
543      for (auto& hit : _h_multi) {
544        if (hit.first.find("_norm") != string::npos) {
545          for (Histo1DPtr& hist : hit.second.histos()) { scale(hist, sf); }
546          const double norm2D = integral2D(hit.second);
547          hit.second.scale(safediv(1.0, norm2D), this);
548        }
549        else {
550          hit.second.scale(sf, this);
551        }
552      }
553    }
554
555
556  private:
557
558    bool _doBoosted, _doResolved;
559
560
561    double TransMass(double ptLep, double phiLep, double met, double phiMet) {
562      return std::sqrt(2.0*ptLep*met*( 1 - std::cos( phiLep-phiMet ) ) );
563    }
564
565
566    double computeneutrinoz(const FourMomentum& lepton, const FourMomentum& met) const {
567      //computing z component of neutrino momentum given lepton and met
568      double pzneutrino;
569      double m_W = 80.399; // in GeV, given in the paper
570      double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py());
571      double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
572      double b = -2*k*lepton.pz();
573      double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k );
574      double discriminant = sqr(b) - 4 * a * c;
575      double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
576      if (discriminant < 0)  pzneutrino = - b / (2 * a); //if the discriminant is negative
577      else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
578        double absquad[2];
579        for (int n=0; n<2; ++n)  absquad[n] = fabs(quad[n]);
580        if (absquad[0] < absquad[1])  pzneutrino = quad[0];
581        else                          pzneutrino = quad[1];
582      }
583      return pzneutrino;
584    }
585
586
587    double integral2D(BinnedHistogram& h_multi) {
588      double total_integral = 0;
589      for  (Histo1DPtr& h : h_multi.histos()) {
590        total_integral += h->integral(false);
591      }
592      return total_integral;
593    }
594
595
596    void book2D(string name, std::vector<double>& doubleDiff_bins, size_t table){
597      for (size_t i = 0; i < doubleDiff_bins.size() - 1; ++i) {
598        string nbin = std::to_string(i);
599        // HepData entry has dummy "Table of Contents", 
600        // so need to offset everything by one unit
601        { Histo1DPtr tmp; _h_multi[name].add(doubleDiff_bins[i], doubleDiff_bins[i+1], book(tmp, table+1+i, 1, 1)); }
602      }
603    }
604
605
606    void book_hist(string name, size_t table) {
607      // HepData entry has dummy "Table of Contents", 
608      // so need to offset everything by one unit
609      book(_h[name], table+3, 1, 1);
610      book(_h[name+"_norm"], table+1, 1, 1);
611    }
612
613    size_t TransformJetMultiplicity(size_t jet_n) { return jet_n > 7 ? 7 : jet_n; }
614
615    size_t TransformExtrajetMultiplicity(size_t jet_n) { return jet_n > 6 ? 6 : jet_n; }
616
617    size_t TransformExtrajetMultiplicity_boosted(size_t jet_n) { return jet_n > 4 ? 4 : jet_n; }
618
619    size_t TransformJetMultiplicity_for_ttbar_m(size_t jet_n) { return jet_n > 6 ? 6 : jet_n; }
620
621    size_t TransformJetMultiplicity_pttop(size_t jet_n) {
622      if (jet_n >= 0 && jet_n < 2)  return 0;
623      if (jet_n == 2)               return 2;
624      if (jet_n > 2)                return 3;
625      return jet_n;
626    }
627
628    size_t TransformJetMultiplicity_ptttbar(size_t jet_n) {
629      if (jet_n >= 0 && jet_n < 2)  return 0;
630      if (jet_n >= 2)               return 2;
631      return jet_n;
632    }
633
634    size_t TransformJetMultiplicity_mttbar(size_t jet_n) {
635      if (jet_n == 0)  return 0;
636      if (jet_n == 1)  return 1;
637      if (jet_n >= 2)  return 2;
638      return jet_n;
639    }
640
641    /// @name Objects that are used by the event selection decisions
642    ///@{
643    map<string, Histo1DPtr> _h;
644    map<string, BinnedHistogram> _h_multi;
645    ///@}
646
647  };
648
649
650  RIVET_DECLARE_PLUGIN(ATLAS_2019_I1750330);
651
652}