CellFinder.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
 
#include "DiTauRec/CellFinder.h"
#include "fastjet/PseudoJet.hh"
 
CellFinder::CellFinder(const std::string& type,
               const std::string& name,
               const IInterface * parent) :
  DiTauToolBase(type, name, parent)
{
  declareInterface<DiTauToolBase > (this);
}
 
 
CellFinder::~CellFinder() = default;
 
 
StatusCode CellFinder::initialize() {
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode CellFinder::execute(DiTauCandidateData * data,
                               const EventContext& /*ctx*/) const {
 
  ATH_MSG_DEBUG("execute CellFinder...");
 
  // get ditau and its seed jet
 
  xAOD::DiTauJet* pDiTau = data->xAODDiTau;
  if (!pDiTau) {
    ATH_MSG_ERROR("no di-tau candidate given");
    return StatusCode::FAILURE;
  }
 
  const xAOD::Jet* pSeed = data->seed;
  if (!pSeed) {
    ATH_MSG_WARNING("No jet seed given.");
    return StatusCode::FAILURE;
  }
 
  std::vector<fastjet::PseudoJet> vSubjets = data->subjets;
  if (vSubjets.empty()) {
    ATH_MSG_WARNING("No subjets given. Continue without cell information.");
    return StatusCode::SUCCESS;
  }
 
  // get clusters linked to the seed jet. Loop over clusters to get linked cells
  // Local variable cellSeen uses 25000 bytes of stack space,
  //coverity[STACK_USE]
  std::bitset<200000> cellSeen;
  std::vector<const CaloCell*> subjetCells;
 
  // loop over seed jet constituents
  for (const auto *const seedConst: pSeed->getConstituents()) {
    // cast jet constituent to cluster object
    const xAOD::CaloCluster* cluster = dynamic_cast<const xAOD::CaloCluster*>( seedConst->rawConstituent() );
 
    // loop over cells which are linked to the cluster
    for (const auto *const cc : *(cluster->getCellLinks())) {
      // skip if pt<0 or cell already encountered
      if (cc->pt() < 0) continue;
      if (cellSeen.test(cc->caloDDE()->calo_hash())) continue;
      // register cell hash as already seen
      cellSeen.set(cc->caloDDE()->calo_hash());            
 
      TLorentzVector temp_cc_p4;
      temp_cc_p4.SetPtEtaPhiM(cc->pt(), cc->eta(), cc->phi(), cc->m());
 
      // check if cell is in one of the subjets cones
      for (const auto& subjet : vSubjets) {
        TLorentzVector temp_sub_p4;
        temp_sub_p4.SetPtEtaPhiM(subjet.pt(), subjet.eta(), subjet.phi_std(), subjet.m());        
    if (temp_cc_p4.DeltaR(temp_sub_p4) < m_Rsubjet) {
      subjetCells.push_back(cc);
    }
      }
    }
  }
 
  ATH_MSG_DEBUG("subjetCells.size()=" << subjetCells.size());
  data->subjetCells = std::move(subjetCells);
 
  // write f_core
  float f_core{};
  for (unsigned int i = 0; i < vSubjets.size(); i++) {
    const fastjet::PseudoJet& subjet = vSubjets.at(i);
    float ptAll = 0.;
    float ptCore = 0.;
 
    TLorentzVector temp_sub_p4;
    temp_sub_p4.SetPtEtaPhiM(subjet.pt(), subjet.eta(), subjet.phi_std(), subjet.m());
 
    for (const auto& cc : data->subjetCells) {
     
      TLorentzVector temp_cc_p4;
      temp_cc_p4.SetPtEtaPhiM(cc->pt(), cc->eta(), cc->phi(), cc->m()); 
 
      if (temp_cc_p4.DeltaR(temp_sub_p4) < data->Rsubjet) {
    ptAll += cc->pt();
      }
 
      if (temp_cc_p4.DeltaR(temp_sub_p4) < data->Rcore) {
    ptCore += cc->pt();
      }
    }
 
    if (ptAll != 0.)
      f_core = ptCore/ptAll;
    else 
      f_core = -999.;
 
    ATH_MSG_DEBUG("subjet "<< i << ": f_core=" << f_core);
    pDiTau->setfCore(i, f_core);
  }
 
  return StatusCode::SUCCESS;
}