DiTauTrackFinder.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#include "DiTauRec/DiTauTrackFinder.h"
#include "tauRecTools/TrackSort.h"
#include "StoreGate/ReadHandle.h"
 
 
 
DiTauTrackFinder::DiTauTrackFinder(const std::string& type,
                   const std::string& name,
                   const IInterface * parent) :
  DiTauToolBase(type, name, parent)
{
  declareInterface<DiTauToolBase > (this);
}
 
 
DiTauTrackFinder::~DiTauTrackFinder() = default;
 
 
StatusCode DiTauTrackFinder::initialize() {
 
  ATH_CHECK( m_TrackSelectorTool.retrieve() );
  ATH_CHECK( m_TrackParticleContainerName.initialize() );
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode DiTauTrackFinder::execute(DiTauCandidateData * data,
                                     const EventContext& ctx) const
{
  ATH_MSG_DEBUG("execute DiTauTrackFinder...");
 
  xAOD::DiTauJet *pDiTau = data->xAODDiTau;
 
  if (pDiTau == nullptr) {
    ATH_MSG_ERROR("no di-tau candidate given");
    return StatusCode::FAILURE;
  }
 
  // retrieve track container
  SG::ReadHandle<xAOD::TrackParticleContainer> pTrackParticleCont (m_TrackParticleContainerName, ctx);
 
  std::vector<const xAOD::TrackParticle*> tauTracks;    // good tracks in subjets
  std::vector<const xAOD::TrackParticle*> isoTracks;    // good tracks in isolation region
  std::vector<const xAOD::TrackParticle*> otherTracks;  // tracks failing the track criteria
 
  // get associated primary vertex
  const xAOD::Vertex* pVertex = nullptr;
  if (pDiTau->vertexLink().isValid()) {
    pVertex = *pDiTau->vertexLink();
  }
  else {
    ATH_MSG_WARNING("could not retieve VertexLink in ditau candidate");
  }
 
  // get tracks
  getTracksFromPV(data, pTrackParticleCont.get(), pVertex, tauTracks, isoTracks, otherTracks);
 
  // clear track links before association
  pDiTau->clearTrackLinks();
  pDiTau->clearIsoTrackLinks();
  pDiTau->clearOtherTrackLinks();
 
  // drop subjets without good tracks
  std::vector<fastjet::PseudoJet> vSubjets = data->subjets;
  int nTracks;
  for (auto subjet_itr = vSubjets.begin(); subjet_itr != vSubjets.end(); ) {
    nTracks = 0;
 
    TLorentzVector temp_p4;
    temp_p4.SetPtEtaPhiM(subjet_itr->pt(), subjet_itr->eta(), subjet_itr->phi_std(), subjet_itr->m()); 
 
    for (const auto& track : tauTracks) {
      if ( temp_p4.DeltaR(track->p4()) < m_MaxDrSubjet) nTracks++;
    }
 
    ATH_MSG_DEBUG("number of tracks in subjet: "<< nTracks);
 
    if (nTracks == 0) {
      ATH_MSG_DEBUG("number of tracks is zero. Drop subjet");
      subjet_itr = vSubjets.erase(subjet_itr); //point subjet_itr to the next element/end of the vector
    } else if( m_MaxNTracksSubjet != -1 && nTracks > m_MaxNTracksSubjet){
      ATH_MSG_DEBUG("number of tracks greater than MaxNTracksSubjet threshold. Drop subjet");
      subjet_itr = vSubjets.erase(subjet_itr); //point subjet_itr to the next element/end of the vector
    }       
    else {
      ++subjet_itr;
    }
  }
    
  // check if ditau candidate has still at least 2 subjets
  if (vSubjets.size()<=1) {
    ATH_MSG_DEBUG("Found less than 2 subjets. Reject ditau candidate");
    return StatusCode::FAILURE;
  }   
 
  data->subjets = vSubjets;
  ATH_MSG_DEBUG("number of subjets  after track association: " << data->subjets.size());
  // set subjet p4 in xAODDiTau
  for (unsigned int i = 0; i < vSubjets.size(); i++) {
    const fastjet::PseudoJet& subjet = vSubjets.at(i);
    pDiTau->setSubjetPtEtaPhiE(i, subjet.pt(), subjet.eta(), subjet.phi_std(), subjet.e());
    ATH_MSG_DEBUG("subjet " << i << " pt: " << subjet.pt() << " eta: " << subjet.eta() << " phi: " << subjet.phi_std() << " e: " << subjet.e());
  }
  vSubjets.clear();
 
 
  // associate tau tracks
  for (const auto& track : tauTracks ) {
    ATH_MSG_DEBUG("adding subjet track. eta: " << track->eta() << " phi: " << track->phi());
    pDiTau->addTrack(pTrackParticleCont.get(), track);
  }
 
  // associate isolation tracks
  for (const auto& track : isoTracks ) {
    ATH_MSG_DEBUG("adding iso track. eta: " << track->eta() << " phi: " << track->phi());
    pDiTau->addIsoTrack(pTrackParticleCont.get(), track);
  }
 
  // associate other tracks
  for (const auto& track : otherTracks ) {
    ATH_MSG_DEBUG("adding other track. eta: " << track->eta() << " phi: " << track->phi());
    pDiTau->addOtherTrack(pTrackParticleCont.get(), track);
  }
 
  return StatusCode::SUCCESS;
}
 
 
void DiTauTrackFinder::getTracksFromPV( const DiTauCandidateData* data,
                    const xAOD::TrackParticleContainer* pTrackParticleCont,
                    const xAOD::Vertex* pVertex,
                    std::vector<const xAOD::TrackParticle*> &tauTracks,
                    std::vector<const xAOD::TrackParticle*> &isoTracks,
                    std::vector<const xAOD::TrackParticle*> &otherTracks) const {
    
  for (const auto *const track : *pTrackParticleCont ) {
    DiTauTrackType type = diTauTrackType(data, track, pVertex);
 
    if (type == DiTauSubjetTrack)
      tauTracks.push_back(track);
    else if (type == DiTauIsoTrack)
      isoTracks.push_back(track);
    else if (type == DiTauOtherTrack)
      otherTracks.push_back(track);
  }
 
  std::sort(tauTracks.begin(), tauTracks.end(), TrackSort());
  std::sort(isoTracks.begin(), isoTracks.end(), TrackSort());
  std::sort(otherTracks.begin(), otherTracks.end(), TrackSort());
}
 
 
DiTauTrackFinder::DiTauTrackType DiTauTrackFinder::diTauTrackType(const DiTauCandidateData* data,
                                  const xAOD::TrackParticle* track,
                                  const xAOD::Vertex* pVertex) const {
 
  xAOD::DiTauJet *pDiTau = data->xAODDiTau;
 
  // check if track is outside the jet ditau cone
  if ( pDiTau->p4().DeltaR(track->p4()) > m_MaxDrJet) return OutsideTrack;
    
  // check quality criteria
  bool goodTrack = m_TrackSelectorTool->decision(*track, pVertex);
  if (!goodTrack) return DiTauOtherTrack;
 
  // check if track is inside a subjet
  std::vector<fastjet::PseudoJet> vSubjets = data->subjets;
  for (const auto &subjet : vSubjets) {
 
    TLorentzVector temp_p4;
    temp_p4.SetPtEtaPhiM(subjet.pt(), subjet.eta(), subjet.phi_std(), subjet.m());
 
    if (temp_p4.DeltaR(track->p4()) < m_MaxDrSubjet)
      return DiTauSubjetTrack;
  }
 
  // track is in isolation region
  return DiTauIsoTrack;
}