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