TtGenEvent.cc

Go to the documentation of this file.

//
//

#include "FWCore/Utilities/interface/EDMException.h"
#include "CommonTools/CandUtils/interface/pdgIdUtils.h"
#include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"

TtGenEvent::TtGenEvent(reco::GenParticleRefProd& decaySubset, reco::GenParticleRefProd& initSubset) {
  parts_ = decaySubset;
  initPartons_ = initSubset;
  if (top() && topBar())
    topPair_ = math::XYZTLorentzVector(top()->p4() + topBar()->p4());
}

bool TtGenEvent::fromGluonFusion() const {
  const reco::GenParticleCollection& initPartsColl = *initPartons_;
  if (initPartsColl.size() == 2)
    if (initPartsColl[0].pdgId() == 21 && initPartsColl[1].pdgId() == 21)
      return true;
  return false;
}

bool TtGenEvent::fromQuarkAnnihilation() const {
  const reco::GenParticleCollection& initPartsColl = *initPartons_;
  if (initPartsColl.size() == 2)
    if (std::abs(initPartsColl[0].pdgId()) < TopDecayID::tID && initPartsColl[0].pdgId() == -initPartsColl[1].pdgId())
      return true;
  return false;
}

WDecay::LeptonType TtGenEvent::semiLeptonicChannel() const {
  WDecay::LeptonType type = WDecay::kNone;
  if (isSemiLeptonic() && singleLepton()) {
    if (std::abs(singleLepton()->pdgId()) == TopDecayID::elecID)
      type = WDecay::kElec;
    if (std::abs(singleLepton()->pdgId()) == TopDecayID::muonID)
      type = WDecay::kMuon;
    if (std::abs(singleLepton()->pdgId()) == TopDecayID::tauID)
      type = WDecay::kTau;
  }
  return type;
}

std::pair<WDecay::LeptonType, WDecay::LeptonType> TtGenEvent::fullLeptonicChannel() const {
  WDecay::LeptonType typeA = WDecay::kNone, typeB = WDecay::kNone;
  if (isFullLeptonic()) {
    if (lepton()) {
      if (std::abs(lepton()->pdgId()) == TopDecayID::elecID)
        typeA = WDecay::kElec;
      if (std::abs(lepton()->pdgId()) == TopDecayID::muonID)
        typeA = WDecay::kMuon;
      if (std::abs(lepton()->pdgId()) == TopDecayID::tauID)
        typeA = WDecay::kTau;
    }
    if (leptonBar()) {
      if (std::abs(leptonBar()->pdgId()) == TopDecayID::elecID)
        typeB = WDecay::kElec;
      if (std::abs(leptonBar()->pdgId()) == TopDecayID::muonID)
        typeB = WDecay::kMuon;
      if (std::abs(leptonBar()->pdgId()) == TopDecayID::tauID)
        typeB = WDecay::kTau;
    }
  }
  return (std::pair<WDecay::LeptonType, WDecay::LeptonType>(typeA, typeB));
}

const reco::GenParticle* TtGenEvent::lepton(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  const reco::GenParticleCollection& partsColl = *parts_;
  for (unsigned int i = 0; i < partsColl.size(); ++i) {
    if (reco::isLepton(partsColl[i]) && partsColl[i].mother() &&
        std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
      if (reco::flavour(partsColl[i]) > 0) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::leptonBar(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  const reco::GenParticleCollection& partsColl = *parts_;
  for (unsigned int i = 0; i < partsColl.size(); ++i) {
    if (reco::isLepton(partsColl[i]) && partsColl[i].mother() &&
        std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
      if (reco::flavour(partsColl[i]) < 0) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::singleLepton(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (isSemiLeptonic(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (reco::isLepton(partsColl[i]) && partsColl[i].mother() &&
          std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::neutrino(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  const reco::GenParticleCollection& partsColl = *parts_;
  for (unsigned int i = 0; i < partsColl.size(); ++i) {
    if (reco::isNeutrino(partsColl[i]) && partsColl[i].mother() &&
        std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
      if (reco::flavour(partsColl[i]) > 0) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::neutrinoBar(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  const reco::GenParticleCollection& partsColl = *parts_;
  for (unsigned int i = 0; i < partsColl.size(); ++i) {
    if (reco::isNeutrino(partsColl[i]) && partsColl[i].mother() &&
        std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
      if (reco::flavour(partsColl[i]) < 0) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::singleNeutrino(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (isSemiLeptonic(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (reco::isNeutrino(partsColl[i]) && partsColl[i].mother() &&
          std::abs(partsColl[i].mother()->pdgId()) == TopDecayID::WID) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::hadronicDecayQuark(bool invertFlavor) const {
  const reco::GenParticle* cand = nullptr;
  // catch W boson and check its daughters for a quark;
  // make sure the decay is semi-leptonic first; this
  // only makes sense if taus are not excluded from the
  // decision
  if (singleLepton(false)) {
    for (reco::GenParticleCollection::const_iterator w = parts_->begin(); w != parts_->end(); ++w) {
      if (std::abs(w->pdgId()) == TopDecayID::WID) {
        // make sure that the particle is a W daughter
        for (reco::GenParticle::const_iterator wd = w->begin(); wd != w->end(); ++wd) {
          // make sure that the parton is a quark
          if (std::abs(wd->pdgId()) < TopDecayID::tID) {
            if (invertFlavor ? reco::flavour(*wd) < 0 : reco::flavour(*wd) > 0) {
              cand = dynamic_cast<const reco::GenParticle*>(&(*wd));
              if (cand == nullptr) {
                throw edm::Exception(edm::errors::InvalidReference, "Not a GenParticle");
              }
              break;
            }
          }
        }
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::hadronicDecayB(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (singleLepton(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (std::abs(partsColl[i].pdgId()) == TopDecayID::bID &&
          reco::flavour(singleLep) == reco::flavour(partsColl[i])) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::hadronicDecayW(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (singleLepton(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (std::abs(partsColl[i].pdgId()) == TopDecayID::WID &&
          reco::flavour(singleLep) != -reco::flavour(partsColl[i])) {
        // PDG Id:13=mu- 24=W+ (+24)->(-13) (-24)->(+13) opposite sign
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::hadronicDecayTop(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (singleLepton(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (std::abs(partsColl[i].pdgId()) == TopDecayID::tID &&
          reco::flavour(singleLep) == reco::flavour(partsColl[i])) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::leptonicDecayB(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (singleLepton(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (std::abs(partsColl[i].pdgId()) == TopDecayID::bID &&
          reco::flavour(singleLep) != reco::flavour(partsColl[i])) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::leptonicDecayW(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (singleLepton(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (std::abs(partsColl[i].pdgId()) == TopDecayID::WID &&
          reco::flavour(singleLep) == -reco::flavour(partsColl[i])) {
        // PDG Id:13=mu- 24=W+ (+24)->(-13) (-24)->(+13) opposite sign
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

const reco::GenParticle* TtGenEvent::leptonicDecayTop(bool excludeTauLeptons) const {
  const reco::GenParticle* cand = nullptr;
  if (singleLepton(excludeTauLeptons)) {
    const reco::GenParticleCollection& partsColl = *parts_;
    const reco::GenParticle& singleLep = *singleLepton(excludeTauLeptons);
    for (unsigned int i = 0; i < partsColl.size(); ++i) {
      if (std::abs(partsColl[i].pdgId()) == TopDecayID::tID &&
          reco::flavour(singleLep) != reco::flavour(partsColl[i])) {
        cand = &partsColl[i];
      }
    }
  }
  return cand;
}

std::vector<const reco::GenParticle*> TtGenEvent::leptonicDecayTopRadiation(bool excludeTauLeptons) const {
  if (leptonicDecayTop(excludeTauLeptons)) {
    return (leptonicDecayTop(excludeTauLeptons)->pdgId() > 0 ? radiatedGluons(TopDecayID::tID)
                                                             : radiatedGluons(-TopDecayID::tID));
  }
  std::vector<const reco::GenParticle*> rad;
  return (rad);
}

std::vector<const reco::GenParticle*> TtGenEvent::hadronicDecayTopRadiation(bool excludeTauLeptons) const {
  if (hadronicDecayTop(excludeTauLeptons)) {
    return (hadronicDecayTop(excludeTauLeptons)->pdgId() > 0 ? radiatedGluons(TopDecayID::tID)
                                                             : radiatedGluons(-TopDecayID::tID));
  }
  std::vector<const reco::GenParticle*> rad;
  return (rad);
}