rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

ATLAS_2018_I1646686

All-hadronic boosted ttbar at 13 TeV
Experiment: ATLAS (LHC)
Inspire ID: 1646686
Status: VALIDATED
Authors:
  • Kyle Cormier
  • Christian Gutschow
References: Beams: p+ p+
Beam energies: (6500.0, 6500.0) GeV
Run details:
  • p + p -> ttbar (all-hadronic, boosted)

Measurements are made of differential cross-sections of highly boosted pair-produced top quarks as a function of top-quark and $t\bar{t}$ system kinematic observables using proton--proton collisions at a center-of-mass energy of $\sqrt{s}=13$ TeV. The data set corresponds to an integrated luminosity of 36.1 fb${}^{-1}$, recorded in 2015 and 2016 with the ATLAS detector at the CERN Large Hadron Collider. Events with two large-radius jets in the final state, one with transverse momentum pT>500 GeV and a second with $p_\text{T} > 350$ GeV, are used for the measurement. The top-quark candidates are separated from the multijet background using jet substructure information and association with a $b$-tagged jet. The measured spectra are corrected for detector effects to a particle-level fiducial phase space and a parton-level limited phase space, and are compared to several Monte Carlo simulations by means of calculated $\chi^2$ values. The cross-section for $t\bar{t}$ production in the fiducial phase-space region is 292$\pm$7 (stat)$\pm$76(syst) fb, to be compared to the theoretical prediction of 384$\pm$36 fb.

Source code: ATLAS_2018_I1646686.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/VetoedFinalState.hh"
  5#include "Rivet/Projections/IdentifiedFinalState.hh"
  6#include "Rivet/Projections/PromptFinalState.hh"
  7#include "Rivet/Projections/DressedLeptons.hh"
  8#include "Rivet/Projections/FastJets.hh"
  9#include "Rivet/Projections/PartonicTops.hh"
 10#include "Rivet/Math/LorentzTrans.hh"
 11
 12namespace Rivet {
 13
 14
 15  /// @brief All-hadronic ttbar at 13 TeV
 16  class ATLAS_2018_I1646686 : public Analysis {
 17  public:
 18
 19      /// Constructor
 20      RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1646686);
 21
 22      /// Book histograms and initialise projections before the run
 23      void init() {
 24
 25        // Get options particle-level only.
 26        _mode = 0;
 27        if ( getOption("TMODE") == "PARTICLE" ) _mode = 0;
 28        if ( getOption("TMODE") == "BOTH" ) _mode = 1;
 29
 30        //histogram booking
 31        book(_h["inclusive"],1,1,1);
 32        bookHistograms("t_pt", 0, true);
 33        bookHistograms("t_y",  1, true);
 34        bookHistograms("t1_pt",         2);
 35        bookHistograms("t1_y",          3);
 36        bookHistograms("t2_pt",         4);
 37        bookHistograms("t2_y",          5);
 38        bookHistograms("tt_m",          6);
 39        bookHistograms("tt_pt",         7);
 40        bookHistograms("tt_y",          8);
 41        bookHistograms("tt_chi",        9);
 42        bookHistograms("tt_yboost",    10);
 43        bookHistograms("tt_pout",      11);
 44        bookHistograms("tt_dPhi",      12);
 45        bookHistograms("tt_Ht",        13);
 46        bookHistograms("tt_cosThStar", 14);
 47
 48        // Projections
 49        Cut dressed_lep = (Cuts::abseta < 2.5) && (Cuts::pT >= 25*GeV);
 50        Cut eta_full = (Cuts::abseta < 5.0);
 51
 52        // All final state particles
 53        FinalState fs(eta_full);
 54
 55        // Get photons to dress leptons
 56        IdentifiedFinalState photons(fs);
 57        photons.acceptIdPair(PID::PHOTON);
 58
 59        // Projection to find the electrons
 60        PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, true);
 61        DressedLeptons dressedelectrons(photons, electrons, 0.1, dressed_lep);
 62        declare(dressedelectrons, "elecs");
 63        DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full);
 64
 65        // Projection to find the muons
 66        PromptFinalState muons(Cuts::abspid == PID::MUON, true);
 67        DressedLeptons dressedmuons(photons, muons, 0.1, dressed_lep);
 68        declare(dressedmuons, "muons");
 69        DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full);
 70
 71        // Jet clustering.
 72        VetoedFinalState vfs;
 73        vfs.addVetoOnThisFinalState(ewdressedelectrons);
 74        vfs.addVetoOnThisFinalState(ewdressedmuons);
 75
 76        FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::DECAY, JetAlg::Invisibles::DECAY);
 77        declare(jets, "jets");
 78
 79        FastJets ljets(fs, FastJets::ANTIKT, 1.0, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE);
 80        declare(ljets, "ljets" );
 81
 82        if (_mode != 0 ){
 83          PartonicTops partonTops;
 84          declare(partonTops, "partonicTops");
 85        }
 86      }
 87
 88
 89      void analyze(const Event& event) {
 90
 91	if (_mode != 0){
 92
 93          // Parton-level top quarks
 94          const Particles partonicTops = apply<PartonicTops>( event, "partonicTops").particlesByPt();
 95          FourMomentum top, tbar;
 96          bool foundT = false, foundTBar = false;
 97          for (const Particle& ptop : partonicTops) {
 98            const int pid = ptop.pid();
 99            if (pid == PID::TQUARK) {
100              top = ptop.momentum();
101              foundT = true;
102            } else if (pid == -PID::TQUARK) {
103              tbar = ptop.momentum();
104              foundTBar = true;
105            }
106          }
107
108	  FourMomentum t1_parton, t2_parton, ttbar_parton;
109	  if ( foundT && foundTBar ) {
110	    t1_parton = top.pT2() > tbar.pT2() ? top : tbar;
111	    t2_parton = top.pT2() > tbar.pT2() ? tbar : top;
112	    ttbar_parton = t1_parton + t2_parton;
113	    
114	    if ( t1_parton.pT() > 500*GeV && t2_parton.pT() > 350*GeV) {
115	      
116	      const double chi_parton = calcChi(t1_parton, t2_parton);
117	      const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton));
118	      if (cosThetaStar_parton == -99) {
119		MSG_DEBUG("ttbar going faster than light! Vetoing event. Try turning of partonic tops?");
120		vetoEvent;
121	      }
122	      const double pout_parton = abs(calcPout(t1_parton, t2_parton));
123	      const double dPhi_parton = deltaPhi(t1_parton, t2_parton);
124	      
125	      const int randomChoice = rand() % 2;
126	      const FourMomentum& randomTopParton = (randomChoice == 0) ? t1_parton : t2_parton;
127	      
128	      fillParton("t_pt", randomTopParton.pT()/GeV);
129	      fillParton("t_y",  randomTopParton.absrap());
130	      
131	      fillParton("t1_pt", t1_parton.pT()/GeV);
132	      fillParton("t1_y",  t1_parton.absrap());
133	      fillParton("t2_pt", t2_parton.pT()/GeV);
134	      fillParton("t2_y",  t2_parton.absrap());
135	      
136	      fillParton("tt_m",  ttbar_parton.mass()/TeV);
137	      fillParton("tt_pt", ttbar_parton.pT()/GeV);
138	      fillParton("tt_Ht", (t1_parton.pT() + t2_parton.pT())/GeV);
139	      fillParton("tt_y",  ttbar_parton.absrap());
140
141	      fillParton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity()));
142	      fillParton("tt_chi", chi_parton);
143	      fillParton("tt_cosThStar", cosThetaStar_parton);
144	      fillParton("tt_pout", pout_parton/GeV);
145	      fillParton("tt_dPhi", dPhi_parton);
146	    }
147	  }
148        }
149
150        // Get and veto on dressed leptons
151        const vector<DressedLepton> dressedElectrons = apply<DressedLeptons>(event, "elecs").dressedLeptons();
152        const vector<DressedLepton> dressedMuons     = apply<DressedLeptons>(event, "muons").dressedLeptons();
153        if (!dressedElectrons.empty()) vetoEvent;
154        if (!dressedMuons.empty()) vetoEvent;
155
156        // Get jets
157        const Jets& all_jets  = apply<FastJets>( event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5);
158        const FastJets& ljets_fj = apply<FastJets>( event, "ljets");
159        const Jets all_ljets = ljets_fj.jetsByPt();
160
161        // Trim the large-R jets
162        Jets trimmedJets;
163        fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05));
164        for (const Jet& jet : all_ljets)
165          trimmedJets += ljets_fj.trimJet(jet, trimmer);
166        trimmedJets = sortByPt(trimmedJets);
167
168        // Check large-R jets
169        Jets ljets;
170        vector<bool> b_tagged;
171        for (const Jet& jet : trimmedJets) {
172
173          if (jet.pT() < 250 * GeV)  continue;
174          if (jet.pT() > 3000 * GeV) continue;
175          if (jet.mass() > jet.pT()) continue;
176          if (jet.abseta() > 2.0 )   continue;
177
178          ljets += jet;
179          b_tagged += jet.bTagged();
180        }
181
182        if (all_jets.size() < 2)  vetoEvent;
183        if (ljets.size() < 2)     vetoEvent;
184
185        // Identify top and anti top, compute some event variables
186        const FourMomentum ttbar = ljets[0].momentum() + ljets[1].momentum();
187        const FourMomentum t1 = ljets[0].momentum();
188        const FourMomentum t2 = ljets[1].momentum();
189
190        const double chi = calcChi(t1, t2);
191        const double cosThetaStar = abs(calcCosThetaStar(t1, t2));
192	if (cosThetaStar == -99) {
193	  MSG_DEBUG("real ttbar going faster than light! This should not happen. Vetoing event.");
194	  vetoEvent;
195	}
196        const double pout = abs(calcPout(t1, t2));
197        const double dPhi = deltaPhi(t1, t2);
198
199        if ( t2.pT() < 350*GeV)  vetoEvent;
200        if ( t1.pT() < 500*GeV)  vetoEvent;
201
202        // b-tagging for particle done on large-R jets
203        if (!(b_tagged[0] && b_tagged[1]))  vetoEvent;
204
205        // Continues with signal region cuts
206        if ( abs(t1.mass() - 172.5 * GeV) > 50*GeV )  vetoEvent;
207        if ( abs(t2.mass() - 172.5 * GeV) > 50*GeV )  vetoEvent;
208
209        _h["inclusive"]->fill(0);
210
211        fillHistograms("t1_pt", t1.pT()/GeV);
212        fillHistograms("t1_y",  t1.absrap());
213        fillHistograms("t2_pt", t2.pT()/GeV);
214        fillHistograms("t2_y",  t2.absrap());
215
216        fillHistograms("tt_m",  ttbar.mass()/TeV);
217        fillHistograms("tt_pt", ttbar.pT()/GeV);
218        fillHistograms("tt_Ht", (t1.pT() + t2.pT())/GeV);
219        fillHistograms("tt_y",  ttbar.absrap());
220
221        fillHistograms("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity()));
222        fillHistograms("tt_chi", chi);
223        fillHistograms("tt_cosThStar", cosThetaStar);
224        fillHistograms("tt_pout", pout/GeV);
225        fillHistograms("tt_dPhi", dPhi);
226
227      }
228
229
230      void finalize() {
231        // Normalize histograms
232        const double sf = crossSection() / sumOfWeights();
233        for (auto &hist : _h) {
234          scale(hist.second, sf);
235          if ((hist.first.find("_norm") != string::npos) && hist.second->integral(false)>0) hist.second->normalize(1.0, false);
236        }
237      }
238
239
240      double calcChi(const FourMomentum& t1, const FourMomentum& t2) {
241        double ystar = 0.5 * (t1.rapidity()-t2.rapidity());
242        double chi = exp( 2 * abs(ystar));
243        return chi;
244      }
245
246      double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) {
247        FourMomentum ttbar = t1 + t2;
248        LorentzTransform centreOfMassTrans;
249        ttbar.setX(0);
250        ttbar.setY(0);
251        if (ttbar.betaVec().mod2() > 1) return -99;
252        centreOfMassTrans.setBetaVec( -ttbar.betaVec() );
253        FourMomentum t1_star = centreOfMassTrans.transform(t1);
254        double cosThetaStar;
255        if (t1_star.p3().mod2() >= 0){
256          cosThetaStar = t1_star.pz()/t1_star.p3().mod();
257        } else {
258          return -99;
259        }
260        return cosThetaStar;
261      }
262
263      double calcPout(const FourMomentum& t1, const FourMomentum& t2) {
264        Vector3 t1V = t1.p3();
265        Vector3 t2V = t2.p3();
266        Vector3 zUnit = Vector3(0., 0., 1.);
267        Vector3 vPerp = zUnit.cross(t1V);
268
269        double pout = vPerp.dot(t2V)/vPerp.mod();
270        return pout;
271      }
272
273
274    protected:
275
276      size_t _mode;
277
278    private:
279
280      map<string, Histo1DPtr> _h;
281
282      //some functions for booking, filling and scaling the histograms
283      void fillHistograms(std::string name, double value) {
284        _h[name]->fill(value);
285        _h[name + "_norm"]->fill(value);
286      }
287
288      void fillParton(std::string name, double value) {
289        _h[name + "_parton"]->fill(value);
290        _h[name + "_parton_norm"]->fill(value);
291      }
292
293      void bookHistograms(const std::string name, unsigned int index, bool onlyParton = false) {
294        if (!onlyParton) {
295          book(_h[name], index, 1, 1 );
296          book(_h[name + "_norm"], index + 13, 1, 1 );
297        }
298        if (_mode != 0) {
299          book(_h[name + "_parton"], index + 82, 1, 1 );
300          book(_h[name + "_parton_norm"], index + 97, 1, 1 );
301        }
302      }
303
304  };
305
306
307  RIVET_DECLARE_PLUGIN(ATLAS_2018_I1646686);
308
309}