Rivet analyses

title: ATLAS_2018_I1646686

All-hadronic boosted ttbar at 13 TeV

Experiment: ATLAS (LHC)

Inspire ID: 1646686

Status: VALIDATED

Authors: - Kyle Cormier - Christian Gutschow

References: - Expt page: ATLAS-TOPQ-2016-09 - DOI: 10.17182/hepdata.81709.v1 - arXiv: 1801.02052

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

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Math/LorentzTrans.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/IdentifiedFinalState.hh"

include "Rivet/Projections/LeptonFinder.hh"

include "Rivet/Projections/PartonicTops.hh"

include "Rivet/Projections/PromptFinalState.hh"

include "Rivet/Projections/VetoedFinalState.hh"

include "Rivet/Tools/Random.hh"

namespace Rivet {

/// @brief All-hadronic ttbar at 13 TeV class ATLAS_2018_I1646686 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2018_I1646686);

/// Book histograms and initialise projections before the run
void init() {

  // Get options particle-level only.
  _mode = 0;
  if (getOption("TMODE") == "PARTICLE") _mode = 0;
  if (getOption("TMODE") == "BOTH") _mode = 1;

  //histogram booking
  book(_h["inclusive"], 1, 1, 1);
  bookHistograms("t_pt", 0, true);
  bookHistograms("t_y", 1, true);
  bookHistograms("t1_pt", 2);
  bookHistograms("t1_y", 3);
  bookHistograms("t2_pt", 4);
  bookHistograms("t2_y", 5);
  bookHistograms("tt_m", 6);
  bookHistograms("tt_pt", 7);
  bookHistograms("tt_y", 8);
  bookHistograms("tt_chi", 9);
  bookHistograms("tt_yboost", 10);
  bookHistograms("tt_pout", 11);
  bookHistograms("tt_dPhi", 12);
  bookHistograms("tt_Ht", 13);
  bookHistograms("tt_cosThStar", 14);

  // Projections
  Cut dressed_lep = (Cuts::abseta < 2.5) && (Cuts::pT >= 25 * GeV);
  Cut eta_full = (Cuts::abseta < 5.0);

  // All final state particles
  FinalState fs(eta_full);

  // Get photons to dress leptons
  IdentifiedFinalState photons(fs);
  photons.acceptIdPair(PID::PHOTON);

  // Projection to find the electrons
  PromptFinalState electrons(Cuts::abspid == PID::ELECTRON, TauDecaysAs::PROMPT);
  LeptonFinder dressedelectrons(electrons, photons, 0.1, dressed_lep);
  declare(dressedelectrons, "elecs");
  LeptonFinder ewdressedelectrons(electrons, photons, 0.1, eta_full);

  // Projection to find the muons
  PromptFinalState muons(Cuts::abspid == PID::MUON, TauDecaysAs::PROMPT);
  LeptonFinder dressedmuons(muons, photons, 0.1, dressed_lep);
  declare(dressedmuons, "muons");
  LeptonFinder ewdressedmuons(muons, photons, 0.1, eta_full);

  // Jet clustering.
  VetoedFinalState vfs;
  vfs.addVetoOnThisFinalState(ewdressedelectrons);
  vfs.addVetoOnThisFinalState(ewdressedmuons);

  FastJets jets(vfs, JetAlg::ANTIKT, 0.4, JetMuons::DECAY, JetInvisibles::DECAY);
  declare(jets, "jets");

  FastJets ljets(fs, JetAlg::ANTIKT, 1.0, JetMuons::NONE, JetInvisibles::NONE);
  ljets.addTrf(new fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2),
                                   fastjet::SelectorPtFractionMin(0.05)));
  declare(ljets, "ljets");

  if (_mode != 0) {
    PartonicTops partonTops;
    declare(partonTops, "partonicTops");
  }
}


void analyze(const Event& event) {

  if (_mode != 0) {

    // Parton-level top quarks
    const Particles partonicTops = apply<PartonicTops>(event, "partonicTops").particlesByPt();
    FourMomentum top, tbar;
    bool foundT = false, foundTBar = false;
    for (const Particle& ptop : partonicTops) {
      const int pid = ptop.pid();
      if (pid == PID::TQUARK) {
        top = ptop.momentum();
        foundT = true;
      }
      else if (pid == -PID::TQUARK) {
        tbar = ptop.momentum();
        foundTBar = true;
      }
    }

    FourMomentum t1_parton, t2_parton, ttbar_parton;
    if (foundT && foundTBar) {
      t1_parton = top.pT2() > tbar.pT2() ? top : tbar;
      t2_parton = top.pT2() > tbar.pT2() ? tbar : top;
      ttbar_parton = t1_parton + t2_parton;

      if (t1_parton.pT() > 500 * GeV && t2_parton.pT() > 350 * GeV) {

        const double chi_parton = calcChi(t1_parton, t2_parton);
        const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton));
        if (cosThetaStar_parton == -99) {
          MSG_DEBUG("ttbar going faster than light! Vetoing event. Try turning of partonic tops?");
          vetoEvent;
        }
        const double pout_parton = abs(calcPout(t1_parton, t2_parton));
        const double dPhi_parton = deltaPhi(t1_parton, t2_parton);

        const int randomChoice = int(rand01() < 0.5);
        const FourMomentum& randomTopParton = (randomChoice == 0) ? t1_parton : t2_parton;

        fillParton("t_pt", randomTopParton.pT() / GeV);
        fillParton("t_y", randomTopParton.absrap());

        fillParton("t1_pt", t1_parton.pT() / GeV);
        fillParton("t1_y", t1_parton.absrap());
        fillParton("t2_pt", t2_parton.pT() / GeV);
        fillParton("t2_y", t2_parton.absrap());

        fillParton("tt_m", ttbar_parton.mass() / TeV);
        fillParton("tt_pt", ttbar_parton.pT() / GeV);
        fillParton("tt_Ht", (t1_parton.pT() + t2_parton.pT()) / GeV);
        fillParton("tt_y", ttbar_parton.absrap());

        fillParton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity()));
        fillParton("tt_chi", chi_parton);
        fillParton("tt_cosThStar", cosThetaStar_parton);
        fillParton("tt_pout", pout_parton / GeV);
        fillParton("tt_dPhi", dPhi_parton);
      }
    }
  }

  // Get and veto on dressed leptons
  const DressedLeptons dressedElectrons = apply<LeptonFinder>(event, "elecs").dressedLeptons();
  const DressedLeptons dressedMuons = apply<LeptonFinder>(event, "muons").dressedLeptons();
  if (!dressedElectrons.empty()) vetoEvent;
  if (!dressedMuons.empty()) vetoEvent;

  // Get jets
  const Jets& all_jets = apply<FastJets>(event, "jets")
                             .jetsByPt(Cuts::pT > 25 * GeV && Cuts::abseta < 2.5);
  const FastJets& ljets_fj = apply<FastJets>(event, "ljets");
  const Jets trimmedJets = ljets_fj.jetsByPt();

  // Check large-R jets
  Jets ljets;
  vector<bool> b_tagged;
  for (const Jet& jet : trimmedJets) {

    if (jet.pT() < 250 * GeV) continue;
    if (jet.pT() > 3000 * GeV) continue;
    if (jet.mass() > jet.pT()) continue;
    if (jet.abseta() > 2.0) continue;

    ljets += jet;
    b_tagged += jet.bTagged();
  }

  if (all_jets.size() < 2) vetoEvent;
  if (ljets.size() < 2) vetoEvent;

  // Identify top and anti top, compute some event variables
  const FourMomentum ttbar = ljets[0].momentum() + ljets[1].momentum();
  const FourMomentum t1 = ljets[0].momentum();
  const FourMomentum t2 = ljets[1].momentum();

  const double chi = calcChi(t1, t2);
  const double cosThetaStar = abs(calcCosThetaStar(t1, t2));
  if (cosThetaStar == -99) {
    MSG_DEBUG("real ttbar going faster than light! This should not happen. Vetoing event.");
    vetoEvent;
  }
  const double pout = abs(calcPout(t1, t2));
  const double dPhi = deltaPhi(t1, t2);

  if (t2.pT() < 350 * GeV) vetoEvent;
  if (t1.pT() < 500 * GeV) vetoEvent;

  // b-tagging for particle done on large-R jets
  if (!(b_tagged[0] && b_tagged[1])) vetoEvent;

  // Continues with signal region cuts
  if (abs(t1.mass() - 172.5 * GeV) > 50 * GeV) vetoEvent;
  if (abs(t2.mass() - 172.5 * GeV) > 50 * GeV) vetoEvent;

  _h["inclusive"]->fill(0);

  fillHistograms("t1_pt", t1.pT() / GeV);
  fillHistograms("t1_y", t1.absrap());
  fillHistograms("t2_pt", t2.pT() / GeV);
  fillHistograms("t2_y", t2.absrap());

  fillHistograms("tt_m", ttbar.mass() / TeV);
  fillHistograms("tt_pt", ttbar.pT() / GeV);
  fillHistograms("tt_Ht", (t1.pT() + t2.pT()) / GeV);
  fillHistograms("tt_y", ttbar.absrap());

  fillHistograms("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity()));
  fillHistograms("tt_chi", chi);
  fillHistograms("tt_cosThStar", cosThetaStar);
  fillHistograms("tt_pout", pout / GeV);
  fillHistograms("tt_dPhi", dPhi);
}


void finalize() {
  // Normalize histograms
  const double sf = crossSection() / picobarn / sumOfWeights();
  for (auto& hist : _h) {
    scale(hist.second, sf);
    if ((hist.first.find("_norm") != string::npos) && hist.second->integral(false) > 0)
      hist.second->normalize(1.0, false);
  }
}


double calcChi(const FourMomentum& t1, const FourMomentum& t2) {
  double ystar = 0.5 * (t1.rapidity() - t2.rapidity());
  double chi = exp(2 * abs(ystar));
  return chi;
}

double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) {
  FourMomentum ttbar = t1 + t2;
  LorentzTransform centreOfMassTrans;
  ttbar.setX(0);
  ttbar.setY(0);
  if (ttbar.betaVec().mod2() > 1) return -99;
  centreOfMassTrans.setBetaVec(-ttbar.betaVec());
  FourMomentum t1_star = centreOfMassTrans.transform(t1);
  double cosThetaStar;
  if (t1_star.p3().mod2() >= 0) {
    cosThetaStar = t1_star.pz() / t1_star.p3().mod();
  }
  else {
    return -99;
  }
  return cosThetaStar;
}

double calcPout(const FourMomentum& t1, const FourMomentum& t2) {
  Vector3 t1V = t1.p3();
  Vector3 t2V = t2.p3();
  Vector3 zUnit = Vector3(0., 0., 1.);
  Vector3 vPerp = zUnit.cross(t1V);

  double pout = vPerp.dot(t2V) / vPerp.mod();
  return pout;
}

protected:

size_t _mode;

private:

map<string, Histo1DPtr> _h;

//some functions for booking, filling and scaling the histograms
void fillHistograms(std::string name, double value) {
  _h[name]->fill(value);
  _h[name + "_norm"]->fill(value);
}

void fillParton(std::string name, double value) {
  _h[name + "_parton"]->fill(value);
  _h[name + "_parton_norm"]->fill(value);
}

void bookHistograms(const std::string name, unsigned int index, bool onlyParton = false) {
  if (!onlyParton) {
    book(_h[name], index, 1, 1);
    book(_h[name + "_norm"], index + 13, 1, 1);
  }
  if (_mode != 0) {
    book(_h[name + "_parton"], index + 82, 1, 1);
    book(_h[name + "_parton_norm"], index + 97, 1, 1);
  }
}

};

RIVET_DECLARE_PLUGIN(ATLAS_2018_I1646686);

} ```