Skip to content
Snippets Groups Projects
Commit d52465c7 authored by Andrii Verbytskyi's avatar Andrii Verbytskyi Committed by Adam Edward Barton
Browse files

Migration of ISF_HepMC_Tools to HepmC3

parent 841c4ce8
No related branches found
No related tags found
No related merge requests found
Showing
with 182 additions and 1 deletion
......@@ -37,8 +37,13 @@ namespace ISF {
/// Creates the InterfaceID and interfaceID() method
DeclareInterfaceID(IGenParticleFilter, 1, 0);
#ifdef HEPMC3
/** Returns a boolean if the particle has passed or not */
virtual bool pass(HepMC::ConstGenParticlePtr particle) const = 0;
#else
/** Returns a boolean if the particle has passed or not */
virtual bool pass(const HepMC::GenParticle& particle) const = 0;
#endif
};
......
......@@ -39,6 +39,17 @@ StatusCode ISF::GenParticleFinalStateFilter::initialize()
/** returns true if the the particle is considered stable */
#ifdef HEPMC3
bool ISF::GenParticleFinalStateFilter::pass(HepMC::ConstGenParticlePtr particle) const
{
bool passFilter = true;
passFilter &= isFinalState(particle);
passFilter &= (!m_checkGenSimStable) || MC::isSimStable(particle);
passFilter &= (!m_checkGenInteracting) || MC::isSimInteracting(particle);
return passFilter;
}
#else
bool ISF::GenParticleFinalStateFilter::pass(const HepMC::GenParticle& particle) const
{
bool passFilter = true;
......@@ -48,6 +59,7 @@ bool ISF::GenParticleFinalStateFilter::pass(const HepMC::GenParticle& particle)
return passFilter;
}
#endif
StatusCode ISF::GenParticleFinalStateFilter::finalize()
......@@ -57,9 +69,18 @@ StatusCode ISF::GenParticleFinalStateFilter::finalize()
}
/** checks if the particle is in its final state (no end vertex) */
#ifdef HEPMC3
bool ISF::GenParticleFinalStateFilter::isFinalState(HepMC::ConstGenParticlePtr p) const {
// particle is in its final state if both:
// * no end_vertex
// * status==1
return ( !p->end_vertex() && p->status()==1 );
}
#else
bool ISF::GenParticleFinalStateFilter::isFinalState(const HepMC::GenParticle &p) const {
// particle is in its final state if both:
// * no end_vertex
// * status==1
return ( !p.end_vertex() && p.status()==1 );
}
#endif
......@@ -42,12 +42,21 @@ namespace ISF {
StatusCode initialize();
StatusCode finalize();
#ifdef HEPMC3
/** Returns the Particle Stack, should register truth */
bool pass(HepMC::ConstGenParticlePtr particle) const;
private:
/** checks if the particle is in its final state (no end vertex) */
bool isFinalState( HepMC::ConstGenParticlePtr p) const;
#else
/** Returns the Particle Stack, should register truth */
bool pass(const HepMC::GenParticle& particle) const;
private:
/** checks if the particle is in its final state (no end vertex) */
bool isFinalState( const HepMC::GenParticle& p) const;
#endif
bool m_checkGenSimStable; //!< boolean switch to check on sim stable
bool m_checkGenInteracting; //!< boolean switch to check on gen interacting
......
......@@ -81,17 +81,29 @@ StatusCode ISF::GenParticleGenericFilter::finalize()
/** Returns whether the given particle passes all cuts or not */
#ifdef HEPMC3
bool ISF::GenParticleGenericFilter::pass(HepMC::ConstGenParticlePtr particle) const
#else
bool ISF::GenParticleGenericFilter::pass(const HepMC::GenParticle& particle) const
#endif
{
bool pass = true;
#ifdef HEPMC3
HepMC::ConstGenVertexPtr productionVertex = particle->production_vertex();
#else
HepMC::ConstGenVertexPtr productionVertex = particle.production_vertex();
#endif
const auto* position = productionVertex ? &productionVertex->position() : nullptr;
if (!position || position->perp()<=m_maxApplicableRadius) {
pass = check_cuts_passed(particle);
}
#ifdef HEPMC3
const auto& momentum = particle->momentum();
#else
const auto& momentum = particle.momentum();
#endif
ATH_MSG_VERBOSE( "GenParticle '" << particle << "' with "
<< (position ? "pos: r=" + std::to_string(position->perp()) : "")
<< ", mom: eta=" << momentum.eta() << " phi=" << momentum.phi()
......@@ -102,12 +114,18 @@ bool ISF::GenParticleGenericFilter::pass(const HepMC::GenParticle& particle) con
/** Check whether the given particle passes all configure cuts or not */
#ifdef HEPMC3
bool ISF::GenParticleGenericFilter::check_cuts_passed(HepMC::ConstGenParticlePtr particle) const {
const auto momentum = particle->momentum();
int pdg = particle->pdg_id();
#else
bool ISF::GenParticleGenericFilter::check_cuts_passed(const HepMC::GenParticle &particle) const {
const auto& momentum = particle.momentum();
int pdg = particle.pdg_id();
#endif
double mom = std::sqrt(momentum.x()*momentum.x()+momentum.y()*momentum.y()+momentum.z()*momentum.z());
double eta = momentum.eta();
double phi = momentum.phi();
int pdg = particle.pdg_id();
// check the particle pdg code
if( m_pdgs.size() && std::find(std::begin(m_pdgs), std::end(m_pdgs), pdg) == std::end(m_pdgs) ) {
......
......@@ -52,11 +52,19 @@ typedef std::vector<int> PDGCodes;
StatusCode finalize();
/// Interface method that returns whether the given particle passes all cuts or not
#ifdef HEPMC3
bool pass(HepMC::ConstGenParticlePtr particle) const;
#else
bool pass(const HepMC::GenParticle& particle) const;
#endif
private:
/// Check whether the given particle passes all configure cuts or not
#ifdef HEPMC3
bool check_cuts_passed(HepMC::ConstGenParticlePtr particle) const;
#else
bool check_cuts_passed(const HepMC::GenParticle& particle) const;
#endif
/// the cuts defined by the use
double m_minEta; //!< min pseudorapidity cut
......
......@@ -11,6 +11,7 @@
// HepMC includes
#include "AtlasHepMC/GenParticle.h"
#include "AtlasHepMC/Flow.h"
// Helper function
#include "TruthUtils/HepMCHelpers.h"
......@@ -67,6 +68,19 @@ StatusCode ISF::GenParticleInteractingFilter::initialize()
}
/** passes through to the private version of the filter */
#ifdef HEPMC3
bool ISF::GenParticleInteractingFilter::pass(HepMC::ConstGenParticlePtr particle) const
{
const int pdg_id = particle->pdg_id();
const bool isInteracting = find(m_additionalInteractingParticleTypes.begin(),
m_additionalInteractingParticleTypes.end(),
pdg_id) != m_additionalInteractingParticleTypes.end();
const bool isNonInteracting = find(m_additionalNonInteractingParticleTypes.begin(),
m_additionalNonInteractingParticleTypes.end(),
pdg_id) != m_additionalNonInteractingParticleTypes.end();
return !(MC::isNonInteracting( particle ) || isNonInteracting) || isInteracting;
}
#else
bool ISF::GenParticleInteractingFilter::pass(const HepMC::GenParticle& particle) const
{
const int& pdg_id = particle.pdg_id();
......@@ -78,4 +92,5 @@ bool ISF::GenParticleInteractingFilter::pass(const HepMC::GenParticle& particle)
pdg_id) != m_additionalNonInteractingParticleTypes.end();
return !(MC::isNonInteracting( &particle ) || isNonInteracting) || isInteracting;
}
#endif
......@@ -18,6 +18,7 @@
// STL includes
#include <string>
#include "AtlasHepMC/GenParticle.h"
namespace ISF {
class ISFParticle;
......@@ -40,9 +41,14 @@ namespace ISF {
/** Framework methods */
virtual StatusCode initialize() override;
#ifdef HEPMC3
/** passes through to the private version */
virtual bool pass(HepMC::ConstGenParticlePtr particle ) const override;
#else
/** passes through to the private version */
virtual bool pass(const HepMC::GenParticle& particle ) const override;
#endif
/** Additional PDG codes to classify as interacting */
std::vector<int> m_additionalInteractingParticleTypes;
......
......@@ -27,10 +27,18 @@ ISF::GenParticleLifetimeFilter::GenParticleLifetimeFilter( const std::string& t,
/** does the given particle pass the filter? */
#ifdef HEPMC3
bool ISF::GenParticleLifetimeFilter::pass(HepMC::ConstGenParticlePtr particle) const
#else
bool ISF::GenParticleLifetimeFilter::pass(const HepMC::GenParticle& particle) const
#endif
{
// the GenParticle end vertex
#ifdef HEPMC3
auto endVtx = particle->end_vertex();
#else
auto endVtx = particle.end_vertex();
#endif
// no production vertex?
if (!endVtx) {
ATH_MSG_DEBUG("GenParticle does not have an end vertex, this is fine");
......@@ -40,7 +48,11 @@ bool ISF::GenParticleLifetimeFilter::pass(const HepMC::GenParticle& particle) co
const auto& end4Vec = endVtx->position();
// the GenParticle production vertex
#ifdef HEPMC3
auto prodVtx = particle->production_vertex();
#else
auto prodVtx = particle.production_vertex();
#endif
// no production vertex?
if (!prodVtx) {
ATH_MSG_DEBUG("GenParticle does not have a production vertex, filtering it out");
......
......@@ -34,7 +34,11 @@ namespace ISF {
~GenParticleLifetimeFilter(){}
/** does the given particle pass the filter? */
#ifdef HEPMC3
bool pass(HepMC::ConstGenParticlePtr particle) const;
#else
bool pass(const HepMC::GenParticle& particle) const;
#endif
private:
double m_minLifetime{0.000001}; //units of c*ns
......
......@@ -12,6 +12,7 @@
// HepMC includes
#include "AtlasHepMC/GenParticle.h"
#include "AtlasHepMC/GenVertex.h"
#include "AtlasHepMC/SimpleVector.h"
/** Constructor **/
ISF::GenParticlePositionFilter::GenParticlePositionFilter( const std::string& t,
......@@ -48,10 +49,17 @@ StatusCode ISF::GenParticlePositionFilter::initialize()
/** does the given particle pass the filter? */
#ifdef HEPMC3
bool ISF::GenParticlePositionFilter::pass(HepMC::ConstGenParticlePtr particle) const
{
// the GenParticle production vertex
auto vtx = particle->production_vertex();
#else
bool ISF::GenParticlePositionFilter::pass(const HepMC::GenParticle& particle) const
{
// the GenParticle production vertex
HepMC::GenVertexPtr vtx = particle.production_vertex();
#endif
// no production vertex?
if (!vtx) {
......
......@@ -45,7 +45,11 @@ namespace ISF {
StatusCode finalize();
/** does the given particle pass the filter? */
#ifdef HEPMC3
bool pass(HepMC::ConstGenParticlePtr particle) const;
#else
bool pass(const HepMC::GenParticle& particle) const;
#endif
private:
ServiceHandle<IGeoIDSvc> m_geoIDSvc;
......
......@@ -72,6 +72,33 @@ StatusCode ISF::GenParticleSimWhiteList::initialize()
}
/** passes through to the private version of the filter */
#ifdef HEPMC3
bool ISF::GenParticleSimWhiteList::pass(HepMC::ConstGenParticlePtr particle) const
{
ATH_MSG_VERBOSE( "Checking whether " << particle << " passes the filter." );
static std::vector<int> vertices(500);
vertices.clear();
bool so_far_so_good = pass( particle , vertices );
// Test all parent particles
if (so_far_so_good && particle->production_vertex() && m_qs){
for (auto pit: particle->production_vertex()->particles_in()){
// Loop breaker
if ( HepMC::barcode(pit) == HepMC::barcode(particle) ) continue;
// Check this particle
vertices.clear();
bool parent_all_clear = pass( pit , vertices );
ATH_MSG_VERBOSE( "Parent all clear: " << parent_all_clear <<
"\nIf true, will not pass the daughter because it should have been picked up through the parent already (to avoid multi-counting)." );
so_far_so_good = so_far_so_good && !parent_all_clear;
} // Loop over parents
} // particle had parents
return so_far_so_good;
}
#else
bool ISF::GenParticleSimWhiteList::pass(const HepMC::GenParticle& particle) const
{
......@@ -98,8 +125,42 @@ bool ISF::GenParticleSimWhiteList::pass(const HepMC::GenParticle& particle) cons
return so_far_so_good;
}
#endif
/** returns true if the the particle and all daughters are on the white list */
#ifdef HEPMC3
bool ISF::GenParticleSimWhiteList::pass(HepMC::ConstGenParticlePtr particle , std::vector<int> & used_vertices ) const
{
// See if the particle is in the white list
bool passFilter = std::binary_search( m_pdgId.begin() , m_pdgId.end() , particle->pdg_id() ) || MC::PID::isNucleus( particle->pdg_id() );
// Remove documentation particles
passFilter = passFilter && particle->status()<3;
// Test all daughter particles
if (particle->end_vertex() && m_qs && passFilter){
// Break loops
if ( std::find( used_vertices.begin() , used_vertices.end() , HepMC::barcode(particle->end_vertex()) )==used_vertices.end() ){
used_vertices.push_back( HepMC::barcode(particle->end_vertex()) );
for (auto pit: particle->end_vertex()->particles_out()){
passFilter = passFilter && pass( pit , used_vertices );
if (!passFilter) {
ATH_MSG_VERBOSE( "Daughter particle " << pit << " does not pass." );
break;
}
} // Loop over daughters
} // Break loops
} // particle had daughters
else if (!particle->end_vertex() && !passFilter && particle->status()<3) { // no daughters... No end vertex... Check if this isn't trouble
ATH_MSG_ERROR( "Found a particle with no end vertex that does not appear in the white list." );
ATH_MSG_ERROR( "This is VERY likely pointing to a problem with either the configuration you ");
ATH_MSG_ERROR( "are using, or a bug in the generator. Either way it should be fixed. The");
ATH_MSG_ERROR( "particle will come next, and then we will throw.");
ATH_MSG_ERROR( particle );
throw;
}
return passFilter;
}
#else
bool ISF::GenParticleSimWhiteList::pass(const HepMC::GenParticle& particle , std::vector<int> & used_vertices ) const
{
// See if the particle is in the white list
......@@ -132,6 +193,7 @@ bool ISF::GenParticleSimWhiteList::pass(const HepMC::GenParticle& particle , std
return passFilter;
}
#endif
StatusCode ISF::GenParticleSimWhiteList::finalize()
{
......
......@@ -19,6 +19,7 @@
#include <string>
#include <vector>
#include "AtlasHepMC/GenParticle.h"
namespace ISF {
class ISFParticle;
......@@ -44,11 +45,19 @@ namespace ISF {
StatusCode finalize();
/** passes through to the private version */
#ifdef HEPMC3
bool pass(HepMC::ConstGenParticlePtr particle ) const;
#else
bool pass(const HepMC::GenParticle& particle ) const;
#endif
private:
/** returns true if the the particle and all daughters are on the white list */
#ifdef HEPMC3
bool pass(HepMC::ConstGenParticlePtr particle , std::vector<int> & used_vertices ) const;
#else
bool pass(const HepMC::GenParticle& particle , std::vector<int> & used_vertices ) const;
#endif
std::vector<std::string> m_whiteLists; //!< The location of the white lists
std::vector<long int> m_pdgId; //!< Allowed PDG IDs
bool m_qs; //!< Switch for quasi-stable particle simulation
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment