Rivet analyses

title: ATLAS_2015_I1376945

Colour flow in hadronic top decay at 8 TeV

Experiment: ATLAS (LHC)

Inspire ID: 1376945

Status: VALIDATED

Authors: - Ben Nachman - Christian Gutschow

References: - Expt page: ATLAS-TOPQ-2014-09 - arXiv: 1506.05629

Beams: p+ p+

Beam energies: (4000.0, 4000.0)GeV

Run details: - ttbar production with one W decaying leptonically, the other one hadronically

The distribution and orientation of energy inside jets is predicted to be an experimental handle on colour connections between the hard-scatter quarks and gluons initiating the jets. This is a measurement of the distribution of one such variable, the jet pull angle. The pull angle is measured for jets produced in $t\bar{t}$ events with one $W$ boson decaying leptonically and the other decaying to jets using 20.3\,$\text{fb}^{-1}$ of data recorded with the ATLAS detector at a centre-of-mass energy of $\sqrt{s} = 8$\,TeV at the LHC. The jet pull angle distribution is corrected for detector resolution and acceptance effects.

Source code:ATLAS_2015_I1376945.cc

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

include "Rivet/Analysis.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/IdentifiedFinalState.hh"

include "Rivet/Projections/LeptonFinder.hh"

include "Rivet/Projections/PromptFinalState.hh"

include "Rivet/Projections/VetoedFinalState.hh"

namespace Rivet {

/// @brief Colour flow in hadronic top decay at 8 TeV class ATLAS_2015_I1376945 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2015_I1376945);


/// @name Analysis methods
/// @{

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

  const FinalState fs;

  PromptFinalState promptFs(fs);
  promptFs.acceptTauDecays(true);
  promptFs.acceptMuonDecays(false);

  IdentifiedFinalState neutrino_fs(promptFs);
  neutrino_fs.acceptNeutrinos();
  declare(neutrino_fs, "NEUTRINO_FS");

  IdentifiedFinalState Photon(fs);
  Photon.acceptIdPair(PID::PHOTON);

  IdentifiedFinalState bare_muons_fs(promptFs);
  bare_muons_fs.acceptIdPair(PID::MUON);

  IdentifiedFinalState bare_elecs_fs(promptFs);
  bare_elecs_fs.acceptIdPair(PID::ELECTRON);

  Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 1 * MeV);
  LeptonFinder muons(bare_muons_fs, Photon, 0.1, lep_cuts);
  declare(muons, "MUONS");

  LeptonFinder elecs(bare_elecs_fs, Photon, 0.1, lep_cuts);
  declare(elecs, "ELECS");

  VetoedFinalState vfs;
  vfs.addVetoOnThisFinalState(muons);
  vfs.addVetoOnThisFinalState(elecs);
  vfs.addVetoOnThisFinalState(neutrino_fs);

  FastJets fjets(vfs, JetAlg::ANTIKT, 0.4);
  fjets.useInvisibles();
  declare(fjets, "jets");

  book(h_pull_all, 4, 1, 1);
  book(h_pull_charged, 5, 1, 1);
}


/// Perform the per-event analysis
void analyze(const Event& event) {

  /**************
   *    JETS    *
   **************/
  const Jets& allJets = apply<FastJets>(event, "jets")
                            .jetsByPt(Cuts::pT > 25.0 * GeV && Cuts::absrap < 2.5);
  const DressedLeptons& all_elecs = apply<LeptonFinder>(event, "ELECS").dressedLeptons();
  const DressedLeptons& all_muons = apply<LeptonFinder>(event, "MUONS").dressedLeptons();
  Jets goodJets;
  for (const Jet& j : allJets) {
    bool keep = true;
    for (const DressedLepton& el : all_elecs) keep &= deltaR(j, el) >= 0.2;
    if (keep) goodJets += j;
  }
  if (goodJets.size() < 4) vetoEvent;

  /****************
   *    LEPTONS   *
   ****************/
  DressedLeptons muons, vetoMuons;
  for (const DressedLepton& mu : all_muons) {
    bool keep = true;
    for (const Jet& j : goodJets) keep &= deltaR(j, mu) >= 0.4;
    if (keep && mu.pt() > 15 * GeV) {
      vetoMuons.push_back(mu);
      if (mu.pt() > 25 * GeV) muons.push_back(mu);
    }
  }

  DressedLeptons elecs, vetoElecs;
  for (const DressedLepton& el : all_elecs) {
    bool keep = true;
    for (const Jet& j : goodJets) keep &= deltaR(j, el) >= 0.4;
    if (keep && el.pt() > 15 * GeV) {
      vetoElecs.push_back(el);
      if (el.pt() > 25 * GeV) elecs.push_back(el);
    }
  }

  if (muons.empty() && elecs.empty()) vetoEvent;

  bool muCandidate = !(muons.size() < 1 || vetoMuons.size() > 1 || vetoElecs.size() > 0);
  bool elCandidate = !(elecs.size() < 1 || vetoElecs.size() > 1 || vetoMuons.size() > 0);

  if (!elCandidate && !muCandidate) vetoEvent;

  /******************************
   *    ELECTRON-MUON OVERLAP   *
   ******************************/
  for (const DressedLepton& electron : elecs) {
    for (const DressedLepton& muon : muons) {
      double d_theta = fabs(muon.theta() - electron.theta());
      double d_phi = deltaPhi(muon.phi(), electron.phi());
      if (d_theta < 0.005 && d_phi < 0.005) vetoEvent;
    }
  }

  /****************
   *  NEUTRINOS   *
   ****************/
  const Particles& neutrinos = apply<IdentifiedFinalState>(event, "NEUTRINO_FS").particlesByPt();
  FourMomentum metVector = FourMomentum(0., 0., 0., 0.);
  for (const Particle& n : neutrinos) {
    metVector += n.momentum();
  }
  double met = metVector.pt();
  if (met <= 20 * GeV) vetoEvent;

  if ((_mT(muCandidate ? muons[0] : elecs[0], metVector) + met) <= 60.) vetoEvent;

  /****************
   *    B-JETS    *
   ****************/
  Jets bJets, wJets;
  for (Jet& j : goodJets) {
    bool b_tagged = false;
    Particles bTags = j.bTags();
    for (Particle& b : bTags) {
      b_tagged |= b.pT() > 5 * GeV;
    }
    if (b_tagged) bJets += j;
    if (!b_tagged && j.abseta() < 2.1) wJets += j;
  }

  if (bJets.size() < 2 || wJets.size() < 2) vetoEvent;

  double pull_angle = fabs(CalculatePullAngle(wJets[0], wJets[1], 0));
  h_pull_all->fill(pull_angle / Rivet::PI);

  double pull_angle_charged = fabs(CalculatePullAngle(wJets[0], wJets[1], 1));
  h_pull_charged->fill(pull_angle_charged / Rivet::PI);
}

Vector3 CalculatePull(Jet& jet, bool& isCharged) {
  Vector3 pull(0.0, 0.0, 0.0);
  double PT = jet.pT();
  Particles& constituents = jet.particles();
  Particles charged_constituents;
  if (isCharged) {
    for (Particle& p : constituents) {
      if (p.charge3() != 0) charged_constituents += p;
    }
    constituents = charged_constituents;
  }
  // calculate axis
  FourMomentum axis;
  for (Particle& p : constituents) axis += p.momentum();
  Vector3 J(axis.rap(), axis.phi(MINUSPI_PLUSPI), 0.0);
  // calculate pull
  for (Particle& p : constituents) {
    Vector3 ri = Vector3(p.rap(), p.phi(MINUSPI_PLUSPI), 0.0) - J;
    while (ri.y() > Rivet::PI) ri.setY(ri.y() - Rivet::TWOPI);
    while (ri.y() < -Rivet::PI) ri.setY(ri.y() + Rivet::TWOPI);
    pull.setX(pull.x() + (ri.mod() * ri.x() * p.pT()) / PT);
    pull.setY(pull.y() + (ri.mod() * ri.y() * p.pT()) / PT);
  }
  return pull;
}

double CalculatePullAngle(Jet& jet1, Jet& axisjet, bool isCharged) {
  Vector3 pull_vector = CalculatePull(jet1, isCharged);
  pull_vector = Vector3(1000. * pull_vector.x(), 1000. * pull_vector.y(), 0.);
  double drap = axisjet.rap() - jet1.rap();
  double dphi = axisjet.phi(MINUSPI_PLUSPI) - jet1.phi(MINUSPI_PLUSPI);
  Vector3 j2_vector(drap, dphi, 0.0);
  return mapAngleMPiToPi(deltaPhi(pull_vector, j2_vector));
}

double _mT(const FourMomentum& l, FourMomentum& nu) const {
  return sqrt(2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))));
}

/// Normalise histograms etc., after the run
void finalize() {
  normalize(h_pull_all);
  normalize(h_pull_charged);
}

/// @}

private:

Histo1DPtr h_pull_all;
Histo1DPtr h_pull_charged;

};

RIVET_DECLARE_PLUGIN(ATLAS_2015_I1376945);

} ```