ParticleFilter.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
*/
 
// Allows the user to search for particles with specified kinematics.
// It will pass if there is an particle with pT and eta or E in the specified range
#include "GeneratorFilters/ParticleFilter.h"
 
ParticleFilter::ParticleFilter(const std::string& name, ISvcLocator* pSvcLocator)
  : GenFilter(name, pSvcLocator)
{
}
 
 
StatusCode ParticleFilter::filterInitialize() {
  ATH_MSG_INFO("Ptcut=" << m_Ptmin);
  ATH_MSG_INFO("Etacut=" << m_EtaRange);
  ATH_MSG_INFO("Energycut=" << m_EnergyRange);
  ATH_MSG_INFO("PDG=" << m_PDGID);
  ATH_MSG_INFO("MinParts=" << m_MinParts);
  ATH_MSG_INFO("Exclusive=" << m_Exclusive);
  return StatusCode::SUCCESS;
}
 
 
StatusCode ParticleFilter::filterEvent() {
  int nParts = 0;
  for (McEventCollection::const_iterator itr = events()->begin(); itr != events()->end(); ++itr) {
    const HepMC::GenEvent* genEvt = (*itr);
    for (const auto& pitr: *genEvt) {
      if (std::abs(pitr->pdg_id()) != m_PDGID ) continue;
      if (pitr->momentum().perp() >= m_Ptmin && std::abs(pitr->momentum().pseudoRapidity()) <= m_EtaRange && pitr->momentum().e() <= m_EnergyRange) {
      if ((!m_Exclusive)&&(m_MinParts == 1))  return StatusCode::SUCCESS; // Found at least one particle and we have an inclusive requirement
              
            // Count only particles not decaying to themselves
            bool notSelfDecay = true;
            if (pitr->end_vertex()) {
              for (const auto&  child: *(pitr->end_vertex())) {
                if ( child->pdg_id() != pitr->pdg_id()) continue;
                if ( child == pitr) continue;
                if ( HepMC::is_simulation_particle(child)) continue;
                  notSelfDecay = false;
                  break;
              }
            }
            if (notSelfDecay) nParts++;
          }
    }
  }
   if (m_Exclusive)
     {setFilterPassed(nParts == m_MinParts);}
   else
     {setFilterPassed(nParts >= m_MinParts);}
  return StatusCode::SUCCESS;
}