Rivet analyses

title: ATLAS_2016_I1449082

Charge asymmetry in top quark pair production in dilepton channel

Experiment: ATLAS (LHC)

Inspire ID: 1449082

Status: VALIDATED

Authors: - Roman Lysak

References: - Expt page: ATLAS-TOPQ-2015-08 - Phys.Rev. D94 (2016) no.3, 032006 - DOI: 10.1103/PhysRevD.94.032006 - arXiv: 1604.05538

Beams: p+ p+

Beam energies: (4000.0, 4000.0)GeV

Run details: - pp --> top + antitop events at 8 TeV

Measurements of the top-antitop quark pair production charge asymmetry in the dilepton channel are presented using data corresponding to an integrated luminosity of $20.3~\text{fb}^{-1}$ from $pp$ collisions at a center-of-mass energy of $\sqrt{s} = 8$ TeV collected with the ATLAS detector at the Large Hadron Collider at CERN. Inclusive and differential measurements as a function of the invariant mass, transverse momentum, and longitudinal boost of the $t\bar{t}$ system are performed both in the full phase space and in a fiducial phase space closely matching the detector acceptance (at least 2 leptons and 2 jets with $p_\text{T} > 25$ GeV and $|\eta| < 2.5$). Two observables are studied: $A_C^{\ell\ell}$ based on the selected leptons and $A_C^{t\bar{t}}$ based on the reconstructed $t\bar{t}$ final state. The unfolded distributions of $\Delta|\eta| = |\eta|{\ell^+} - |\eta|{\ell^-}$ and $\Delta|y| = |y|\text{top} - |y|\text{antitop}$ are provided. USERS SHOULD NOTE THAT EXPLICIT RECONSTRUCTION OF INDIVIDUAL STANDARD MODEL NEUTRINOS IS USED IN THIS ANALYSIS ROUTINE TO MATCH THE MONTE-CARLO-BASED CORRECTION TO THE FIDUCIAL REGION APPLIED IN THE PAPER.

Source code:ATLAS_2016_I1449082.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 {

/// Charge asymmetry in top quark pair production in dilepton channel class ATLAS_2016_I1449082 : public Analysis { public:

const double MW = 80.300 * GeV;
const double MTOP = 172.5 * GeV;

enum MeasureType { kInclMeas, kmttMeas, kbetaMeas, kptMeas, kNmeas };
const size_t kNbins = 2;
const double bins[kNmeas][3];
const string measStr[kNmeas] = {"incl", "mtt", "beta", "ptt"};
const string rangeStr[kNmeas][2];


/// Constructor
//RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1449082);
ATLAS_2016_I1449082()
    : Analysis("ATLAS_2016_I1449082"),
      // inclusive (dummy), mtt [GeV], beta, pTtt
      bins{{0., 1., 2.}, {0., 500., 2000.}, {0., 0.6, 1.0}, {0., 30., 1000.}},
      rangeStr{{"0_1", "1_2"}, {"0_500", "500_2000"}, {"0_0.6", "0.6_1.0"}, {"0_30", "30_1000"}} { }


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

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

  // Cuts
  const Cut eta_full = Cuts::abseta < 5.0;
  const Cut lep_cuts = Cuts::abseta < 2.5 && Cuts::pT > 25 * GeV;
  // All final state particles
  FinalState fs(eta_full);
  // Get photons to dress leptons
  IdentifiedFinalState photons(fs, PID::PHOTON);


  // Electron projections
  // ---------------------
  // Electron/muons are defined from electron/muon and photons within a cone of DR = 0.1.
  // No isolation condition is imposed. The parent of the electron/muon is required not to be a hadron or quark.
  IdentifiedFinalState el_id(fs, {PID::ELECTRON, -PID::ELECTRON});
  PromptFinalState electrons(el_id);
  electrons.acceptTauDecays(true);
  // Electron dressing
  LeptonFinder dressedelectrons(electrons, photons, 0.1, lep_cuts);
  declare(dressedelectrons, "dressedelectrons");
  LeptonFinder dressedelectrons_full(electrons, photons, 0.1, eta_full);

  // Muon projections
  // ---------------------
  IdentifiedFinalState mu_id(fs, {PID::MUON, -PID::MUON});
  PromptFinalState muons(mu_id);
  muons.acceptTauDecays(true);
  // Muon dressing
  LeptonFinder dressedmuons(muons, photons, 0.1, lep_cuts);
  declare(dressedmuons, "dressedmuons");
  LeptonFinder dressedmuons_full(muons, photons, 0.1, eta_full);

  // Neutrino projections
  // ---------------------
  // Missing ET is calculated as the 4–vector sum of neutrinos from W/Z-boson decays. Tau decays are
  // included. A neutrino is treated as a detectable particle and is selected for consideration in the same
  // way as electrons or muons, i.e. the parent is required not to be a hadron or quark (u − b).
  IdentifiedFinalState nu_id;
  nu_id.acceptNeutrinos();
  PromptFinalState neutrinos(nu_id);
  neutrinos.acceptTauDecays(true);
  declare(neutrinos, "neutrinos");

  // Jets projections
  // ---------------------
  // Jets are defined with the anti-kt algorithm, clustering all stable particles excluding the electrons,
  // muons, neutrinos, and photons used in the definition of the selected leptons.
  VetoedFinalState vfs(fs);
  vfs.addVetoOnThisFinalState(dressedelectrons_full);
  vfs.addVetoOnThisFinalState(dressedmuons_full);
  vfs.addVetoOnThisFinalState(neutrinos);
  declare(FastJets(vfs, JetAlg::ANTIKT, 0.4), "Jets");

  // Book histograms
  book(_h_dEta, 1, 1, 1);
  book(_h_dY, 2, 1, 1);
  for (size_t iM = 0; iM < kNmeas; ++iM) {
    book(_h_Acll[iM], 3 + iM, 1, 1);
    book(_h_Actt[iM], 7 + iM, 1, 1);
  }
  for (size_t iM = 0; iM < kNmeas; ++iM) {
    for (size_t iB = 0; iB < kNbins; ++iB) {
      book(_h_dEta_asym[iM][iB], "_dEta_asym_" + measStr[iM] + "_bin" + rangeStr[iM][iB], 2, -10., 10.);
      book(_h_dY_asym[iM][iB], "_dY_asym_" + measStr[iM] + "_bin" + rangeStr[iM][iB], 2, -10., 10.);
    }
  }
}


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

  // Get the electrons and muons
  const DressedLeptons dressedelectrons = apply<LeptonFinder>(event, "dressedelectrons").dressedLeptons();
  const DressedLeptons dressedmuons = apply<LeptonFinder>(event, "dressedmuons").dressedLeptons();
  const DressedLeptons leptons = dressedelectrons + dressedmuons;
  // Require at least 2 leptons in the event
  if (leptons.size() < 2) vetoEvent;

  // Get the neutrinos
  const Particles neutrinos = apply<PromptFinalState>(event, "neutrinos").particlesByPt();
  // Require at least 2 neutrinos in the event (ick)
  if (neutrinos.size() < 2) vetoEvent;

  // Get jets and apply selection
  const Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 25 * GeV && Cuts::abseta < 2.5);
  // Require at least 2 jets in the event
  if (jets.size() < 2) vetoEvent;

  // Remaining selections
  // Events where leptons and jets overlap, within dR = 0.4, are rejected.
  for (const DressedLepton& lepton : leptons) {
    if (any(jets, deltaRLess(lepton, 0.4))) vetoEvent;
  }

  // Construct pseudo-tops
  // Exactly 2 opposite-sign leptons are required (e/mu)
  if (leptons.size() != 2) vetoEvent;
  if ((leptons[0].charge() * leptons[1].charge()) > 0.) vetoEvent;
  const FourMomentum lep_p = (leptons[0].charge3() > 0) ? leptons[0] : leptons[1];
  const FourMomentum lep_n = (leptons[0].charge3() > 0) ? leptons[1] : leptons[0];

  // Only the 2 leading pT selected neutrinos are considered
  const FourMomentum nu1 = neutrinos[0].momentum();
  const FourMomentum nu2 = neutrinos[1].momentum();

  // Two jets correspond to the two leading jets in the event.
  // If there is any b-tagged jet in the event, then the b-tagged jets
  // are preferentially selected over the non-tagged jets without taking into account its pT.
  // A jet is a b–jet if any B–hadron is included in the jet.
  // Only B-hadrons with an initial pT > 5 GeV are considered.
  Jets bjets, lightjets;
  for (const Jet& jet : jets) {
    (jet.bTagged(Cuts::pT > 5 * GeV) ? bjets : lightjets) += jet;
  }
  // Already sorted by construction, since jets is sorted by decreasing pT
  // std::sort(bjets.begin()    , bjets.end()    , cmpMomByPt);
  // std::sort(lightjets.begin(), lightjets.end(), cmpMomByPt);

  // Initially take 2 highest pT jets
  FourMomentum bjet1 = jets[0];
  FourMomentum bjet2 = jets[1];
  if (!bjets.empty()) {
    bjet1 = bjets[0];
    bjet2 = (bjets.size() > 1) ? bjets[1]
                               : lightjets[0]; //< We should have a light jet because >=2 jets requirement
  }
  else {
    // No btagged jets --> should have >= 2 light jets
    bjet1 = lightjets[0];
    bjet2 = lightjets[1];
  }

  // Construct pseudo-W bosons from lepton-neutrino combinations
  // Minimize the difference between the mass computed from each lepton-neutrino combination and the W boson mass
  const double massDiffW1 = fabs((nu1 + lep_p).mass() - MW) + fabs((nu2 + lep_n).mass() - MW);
  const double massDiffW2 = fabs((nu1 + lep_n).mass() - MW) + fabs((nu2 + lep_p).mass() - MW);
  const FourMomentum Wp = (massDiffW1 < massDiffW2) ? nu1 + lep_p : nu2 + lep_p;
  const FourMomentum Wn = (massDiffW1 < massDiffW2) ? nu2 + lep_n : nu1 + lep_n;

  // Construct pseudo-tops from jets and pseudo-W bosons
  // Minimize the difference between the mass computed from each W-boson and b-jet combination and the top mass
  const double massDiffT1 = fabs((Wp + bjet1).mass() * GeV - MTOP)
      + fabs((Wn + bjet2).mass() * GeV - MTOP);
  const double massDiffT2 = fabs((Wp + bjet2).mass() * GeV - MTOP)
      + fabs((Wn + bjet1).mass() * GeV - MTOP);
  const FourMomentum top_p = (massDiffT1 < massDiffT2) ? Wp + bjet1 : Wp + bjet2;
  const FourMomentum top_n = (massDiffT1 < massDiffT2) ? Wn + bjet2 : Wn + bjet1;

  // Calculate d|eta|, d|y|, etc.
  double dEta = lep_p.abseta() - lep_n.abseta();
  double dY = top_p.absrapidity() - top_n.absrapidity();
  double mtt = (top_p + top_n).mass() * GeV;
  double beta = fabs((top_p + top_n).pz()) / (top_p + top_n).E();
  double pttt = (top_p + top_n).pt() * GeV;

  // Fill histos, counters
  _h_dEta->fill(dEta);
  _h_dY->fill(dY);
  // Histos for inclusive and differential asymmetries
  int mttBinID = getBinID(kmttMeas, mtt);
  int betaBinID = getBinID(kbetaMeas, beta);
  int ptttBinID = getBinID(kptMeas, pttt);
  for (int iM = 0; iM < kNmeas; ++iM) {
    int binID = -1;
    switch (iM) {
      case kInclMeas: binID = 0; break;
      case kmttMeas: binID = mttBinID; break;
      case kbetaMeas: binID = betaBinID; break;
      case kptMeas: binID = ptttBinID; break;
      default: binID = -1; break;
    }
    if (binID >= 0) {
      _h_dY_asym[iM][binID]->fill(dY);
      _h_dEta_asym[iM][binID]->fill(dEta);
    }
  }
}


/// Normalise histograms etc., after the run
void finalize() {

  // Calculate charge asymmetries and fill them to hists
  // Just for cross-check, calculate asymmetries from original dEta/dY histos
  double asym = 0, err = 0;
  calcAsymAndError(_h_dEta, asym, err);
  MSG_INFO("Lepton inclusive asymmetry from histo:  = " << asym << " +- " << err);
  calcAsymAndError(_h_dY, asym, err);
  MSG_INFO("ttbar inclusive asymmetry from histo:  = " << asym << " +- " << err);

  // dEta/dY distributions: normalize to unity
  normalize(_h_dEta);
  normalize(_h_dY);

  // Build asymm scatters
  for (size_t iM = 0; iM < kNmeas; ++iM) {
    for (size_t iB = 0; iB < kNbins; ++iB) {
      // Only one bin for inclusive measurement
      if ((iM == kInclMeas) && (iB != 0)) continue;

      calcAsymAndError(_h_dEta_asym[iM][iB], asym, err);
      _h_Acll[iM]->bin(iB + 1).set(asym, err);

      calcAsymAndError(_h_dY_asym[iM][iB], asym, err);
      _h_Actt[iM]->bin(iB + 1).set(asym, err);
    }
  }
}

/// @}

private:

void calcAsymAndError(Histo1DPtr hist, double& asym, double& err) {
  int nBins = hist->numBins();
  if (nBins % 2 != 0) {
    asym = -999;
    err = -999.;
    return;
  }

  double Nneg = 0.;
  double Npos = 0.;
  double dNneg = 0.;
  double dNpos = 0.;
  for (int iB = 0; iB < nBins; ++iB) {
    if (iB < nBins / 2) {
      Nneg += hist->bin(iB + 1).sumW();
      dNneg += hist->bin(iB + 1).sumW2();
    }
    else {
      Npos += hist->bin(iB + 1).sumW();
      dNpos += hist->bin(iB + 1).sumW2();
    }
  }
  dNneg = sqrt(dNneg);
  dNpos = sqrt(dNpos);
  asym = (Npos + Nneg) != 0.0 ? (Npos - Nneg) / (Npos + Nneg) : -999.;

  const double Ntot = Npos + Nneg;
  const double Ntot2 = Ntot * Ntot;
  const double dNpos2 = dNpos * dNpos;
  const double dNneg2 = dNneg * dNneg;
  err = Ntot2 != 0. ? 2. * sqrt((dNneg2 * Npos * Npos + dNpos2 * Nneg * Nneg) / (Ntot2 * Ntot2)) : -999.;
}


int getBinID(MeasureType type, double value) {
  /// @todo Use Rivet index() function
  for (size_t iBin = 0; iBin < kNbins; ++iBin) {
    if (value <= bins[type][iBin + 1]) return iBin;
  }
  return -1;
}


/// @name Histograms
/// @{
Histo1DPtr _h_dEta;
Histo1DPtr _h_dY;
Estimate1DPtr _h_Actt[kNmeas];
Estimate1DPtr _h_Acll[kNmeas];
// Histograms to calculate the asymmetries from
/// @todo Use /TMP histos?
Histo1DPtr _h_dEta_asym[kNmeas][2];
Histo1DPtr _h_dY_asym[kNmeas][2];
/// @}

// Not-scaled histos
Histo1DPtr _h_dEta_notscaled, _h_dY_notscaled;

};

RIVET_DECLARE_PLUGIN(ATLAS_2016_I1449082); } ```