egammaTruthAlg.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "egammaTruthAlg.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "xAODTruth/TruthParticleAuxContainer.h"
#include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
#include "TrkCaloExtension/CaloExtension.h"
#include "TrkCaloExtension/CaloExtensionHelpers.h"
#include "StoreGate/ReadHandle.h"
#include "StoreGate/WriteHandle.h"
#include "GaudiKernel/SystemOfUnits.h"
 
#include "TruthUtils/MagicNumbers.h"
#include "TruthUtils/HepMCHelpers.h"
 
using Gaudi::Units::GeV;
using Gaudi::Units::MeV;
 
 
 
 
 
 
namespace D3PD {
 
 
egammaTruthAlg::egammaTruthAlg (const std::string& name,
                                ISvcLocator* svcloc)
  : AthReentrantAlgorithm (name, svcloc)
{
  declareProperty ("AuxPrefix", m_auxPrefix,
                   "Prefix to add to aux data items.");
 
  declareProperty ("ElectronPtMin", m_electronPtMin = 2*GeV,
                   "Minimum pt for electrons.");
  declareProperty ("PhotonPtMin",   m_photonPtMin = 2*GeV,
                   "Minimum pt for photons.");
  declareProperty ("EtaMax",    m_etaMax = 2.5,
                   "Maximum eta.");
  declareProperty ("IsoCone",    m_isoCone = 0.2,
                   "Isolation cone width.");
  declareProperty ("PhotonEtIsoMax", m_photonEtIsoMax = 2*MeV,
                   "Maximum isolation cone energy allowed to keep a photon.");
}
 
 
StatusCode egammaTruthAlg::initialize()
{
  CHECK( AthReentrantAlgorithm::initialize() );
  CHECK( m_exten.retrieve() );
  CHECK( m_inputKey.initialize() );
  CHECK( m_outputKey.initialize() );
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode egammaTruthAlg::execute (const EventContext& ctx) const
{
  SG::ReadHandle<xAOD::TruthParticleContainer> pin (m_inputKey, ctx);
 
  auto pout = std::make_unique<xAOD::TruthParticleContainer>();
  auto pout_aux = std::make_unique<xAOD::TruthParticleAuxContainer>();
  pout->setStore (pout_aux.get());
 
#define DECOR(TYPE,N) xAOD::TruthParticle::Decorator<TYPE> N (m_auxPrefix + #N)
  DECOR(float,        etaCalo);
  DECOR(float,        phiCalo);
  DECOR(float,        depthCalo);
  DECOR(float,        Etcone20);
#undef DECOR
 
  for (const xAOD::TruthParticle* tp : *pin) {
    float iso = -999;
    if (isAccepted (*tp, *pin, iso)) {
      pout->push_back (std::make_unique<xAOD::TruthParticle>());
      *pout->back() = *tp;
 
      CHECK( findImpact (*tp,
                         etaCalo(*pout->back()),
                         phiCalo(*pout->back()),
                         depthCalo(*pout->back())) );
      Etcone20(*pout->back()) = iso;
    }
  }
 
  SG::WriteHandle<xAOD::TruthParticleContainer> output (m_outputKey, ctx);
  CHECK( output.record (std::move(pout), std::move(pout_aux)) );
 
  return StatusCode::SUCCESS;
}
 
 
bool egammaTruthAlg::isAccepted (const xAOD::TruthParticle& tp,
                                 const xAOD::TruthParticleContainer& cont,
                                 float& iso) const
{
  iso = -999;
 
  int id = tp.pdgId();
  int aid = abs(id);
  int uniqueID = HepMC::uniqueID(tp);
 
  if (aid != abs(MC::ELECTRON) && !MC::isPhoton(aid)) return false;
  if (aid == abs(MC::ELECTRON) && tp.pt() < m_electronPtMin) return false;
  if (MC::isPhoton(aid) && tp.pt() < m_photonPtMin) return false;
 
 
  if (fabs(tp.eta()) > m_etaMax) return false;
  if ( (!((HepMC::is_simulation_particle(&tp) || MC::isZeroEnergyPhoton(&tp)) || ( MC::isStable(&tp) && MC::isSpecialNonInteracting(&tp)))) && MC::isStableOrSimDecayed(&tp)) return false;
 
  // Remove electrons/gammas decaying into themselves
  if( tp.hasDecayVtx() ) {
    const xAOD::TruthVertex* v = tp.decayVtx();
    size_t sz = v->nOutgoingParticles();
    for (size_t i = 0; i < sz; i++) {
      const xAOD::TruthParticle* child = v->outgoingParticle(i);
      if( child && child->pdgId()==id && HepMC::uniqueID(child) != uniqueID
          && (HepMC::generations(child) == 0))
      {
        return false;
      }
    }
  } // end decays into themselves
 
  // Isolation selection for photons.
  iso = computeIso (tp, cont);
  if (aid == abs(MC::PHOTON)) {
    if (iso > m_photonEtIsoMax)
      return false;
  }
 
  return true;
}
 
 
float egammaTruthAlg::computeIso (const xAOD::TruthParticle& tp,
                                  const xAOD::TruthParticleContainer& cont) const
{
  TLorentzVector sum;
  for (const xAOD::TruthParticle* p : cont) {
    if (p == &tp || HepMC::uniqueID(p) == HepMC::uniqueID(tp)) continue;
    if (!((HepMC::is_simulation_particle(p) || MC::isZeroEnergyPhoton(p) || ( MC::isStable(p) && MC::isSpecialNonInteracting(p))) && MC::isStableOrSimDecayed(p))) continue;
    if (tp.p4().DeltaR (p->p4()) < m_isoCone)
      sum += p->p4();
  }
 
  return sum.Pt();
}
 
 
StatusCode egammaTruthAlg::findImpact (const xAOD::TruthParticle& tp,
                                       float& etaCalo,
                                       float& phiCalo,
                                       float& depthCalo) const
{
  etaCalo = -999;
  phiCalo = -999;
  depthCalo = -999;
 
  std::unique_ptr<Trk::CaloExtension> extension =
    m_exten->caloExtension(Gaudi::Hive::currentContext(), tp);
  if (!extension) {
    REPORT_MESSAGE (MSG::ERROR) <<  "Extension to calorimeter failed";
    return StatusCode::FAILURE;
  }
 
  CaloExtensionHelpers::EntryExitPerLayerVector lvec;
  CaloExtensionHelpers::entryExitPerLayerVector (*extension, lvec);
  for (const auto& [sampling, entry, exit] : lvec) {
    if (sampling == CaloSampling::EMB2) {
      etaCalo = entry.eta();
      phiCalo = entry.phi();
      depthCalo = entry.perp();
      break;
    }
    else if (sampling == CaloSampling::EME2) {
      etaCalo = entry.eta();
      phiCalo = entry.phi();
      depthCalo = std::abs(entry.z());
    }
  }
 
  return StatusCode::SUCCESS;
}
 
 
} // namespace D3PD