diff --git a/Generators/AtlasHepMC/AtlasHepMC/GenEvent.h b/Generators/AtlasHepMC/AtlasHepMC/GenEvent.h index 729e7e5637a65498725b53fe141cade954e0d990..9fcc46a700fd4ae363dbe7a0e5fa626d13bead44 100644 --- a/Generators/AtlasHepMC/AtlasHepMC/GenEvent.h +++ b/Generators/AtlasHepMC/AtlasHepMC/GenEvent.h @@ -11,6 +11,7 @@ #include "HepMC3/PrintStreams.h" #include "AtlasHepMC/GenVertex.h" #include "AtlasHepMC/GenParticle.h" +#include "AtlasHepMC/SimpleVector.h" namespace HepMC3 { inline std::vector<HepMC3::GenParticlePtr>::const_iterator begin(HepMC3::GenEvent& e){ return e.particles().begin(); } diff --git a/Generators/AtlasHepMC/AtlasHepMC/GenVertex.h b/Generators/AtlasHepMC/AtlasHepMC/GenVertex.h index d3810d3d002a5075b7321c683b792e715a03c67a..d1c4dcc375bbdb00da439a7f38f11034cc37773c 100644 --- a/Generators/AtlasHepMC/AtlasHepMC/GenVertex.h +++ b/Generators/AtlasHepMC/AtlasHepMC/GenVertex.h @@ -6,6 +6,11 @@ #ifdef HEPMC3 #include "HepMC3/GenVertex.h" #include "HepMC3/PrintStreams.h" +namespace HepMC3 +{ +inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator begin(const HepMC3::GenVertex& v){ return v.particles_out().begin(); } +inline std::vector<HepMC3::ConstGenParticlePtr>::const_iterator end(const HepMC3::GenVertex& v){ return v.particles_out().end(); } +} namespace HepMC { typedef HepMC3::GenVertexPtr GenVertexPtr; typedef HepMC3::ConstGenVertexPtr ConstGenVertexPtr; @@ -36,12 +41,16 @@ inline std::vector<HepMC3::ConstGenVertexPtr> DESCENDANTS(HepMC3::ConstGenVertex return std::vector<HepMC3::ConstGenVertexPtr>(); } inline void* raw_pointer(GenVertexPtr p){ return p.get();} +inline const void* raw_pointer(ConstGenVertexPtr p){ return p.get();} } #else #include "HepMC/GenVertex.h" namespace HepMC { typedef HepMC::GenVertex* GenVertexPtr; typedef const HepMC::GenVertex* ConstGenVertexPtr; +inline GenVertex::particles_out_const_iterator begin(const HepMC::GenVertex& v){ return v.particles_out_const_begin(); } +inline GenVertex::particles_out_const_iterator end(const HepMC::GenVertex& v){ return v.particles_out_const_end(); } + inline GenVertexPtr newGenVertexPtr(const HepMC::FourVector &pos = HepMC::FourVector(0.0,0.0,0.0,0.0), const int i=0) { return new HepMC::GenVertex(pos,i); } diff --git a/Generators/AtlasHepMC/AtlasHepMC/Polarization.h b/Generators/AtlasHepMC/AtlasHepMC/Polarization.h index c637be312b2f081984b71f51e47ed4106728057e..59ecdd8d743220b6dec821f0cccceed4cddd12f8 100644 --- a/Generators/AtlasHepMC/AtlasHepMC/Polarization.h +++ b/Generators/AtlasHepMC/AtlasHepMC/Polarization.h @@ -26,6 +26,13 @@ inline Polarization polarization(HepMC3::GenParticlePtr a){ double theta=(theta_A?theta_A->value():0.0); return Polarization(theta,phi); } +inline Polarization polarization(HepMC3::ConstGenParticlePtr a){ + std::shared_ptr<HepMC3::DoubleAttribute> phi_A =a->attribute<HepMC3::DoubleAttribute>("phi"); + std::shared_ptr<HepMC3::DoubleAttribute> theta_A=a->attribute<HepMC3::DoubleAttribute>("theta"); + double phi=(phi_A?phi_A->value():0.0); + double theta=(theta_A?theta_A->value():0.0); + return Polarization(theta,phi); +} template<class T> void set_polarization( T a, Polarization b) { a->add_attribute("phi",std::make_shared<HepMC3::DoubleAttribute>(b.phi())); diff --git a/Generators/GeneratorFilters/GeneratorFilters/HTFilter.h b/Generators/GeneratorFilters/GeneratorFilters/HTFilter.h index e89acf18ad75c2151e8264e5b89a1f7e86bd6180..9bbab097e62a78c39dbe0e6bb7e9d10f3d2093b7 100644 --- a/Generators/GeneratorFilters/GeneratorFilters/HTFilter.h +++ b/Generators/GeneratorFilters/GeneratorFilters/HTFilter.h @@ -22,8 +22,8 @@ public: virtual StatusCode filterFinalize(); virtual StatusCode filterEvent(); - bool fromTau( const HepMC::GenParticlePtr tp ) const; - bool fromWZ( const HepMC::GenParticlePtr tp ) const; + bool fromTau( HepMC::ConstGenParticlePtr tp ) const; + bool fromWZ( HepMC::ConstGenParticlePtr tp ) const; private: diff --git a/Generators/GeneratorFilters/GeneratorFilters/TtHtoVVDecayFilter.h b/Generators/GeneratorFilters/GeneratorFilters/TtHtoVVDecayFilter.h index 55012480b5bb302b93e0691ddb9ce97e01bdaf5f..e135e0f14df684c8da83e94823b5b1ea30291485 100644 --- a/Generators/GeneratorFilters/GeneratorFilters/TtHtoVVDecayFilter.h +++ b/Generators/GeneratorFilters/GeneratorFilters/TtHtoVVDecayFilter.h @@ -40,7 +40,7 @@ private: int m_nGoodHtoVV; int m_nGoodHtoVVSameCharge; - bool findAncestor(const HepMC::GenVertexPtr searchvertex, + bool findAncestor(HepMC::ConstGenVertexPtr searchvertex, int targetPDGID); }; diff --git a/Generators/GeneratorFilters/src/BoostedHadTopAndTopPair.cxx b/Generators/GeneratorFilters/src/BoostedHadTopAndTopPair.cxx index 4b799ed9f476934b7d32fc59cc9d408a65f824cc..14ceaf9977d38c3e2425f8c9b917866d3de755b4 100644 --- a/Generators/GeneratorFilters/src/BoostedHadTopAndTopPair.cxx +++ b/Generators/GeneratorFilters/src/BoostedHadTopAndTopPair.cxx @@ -103,12 +103,16 @@ bool BoostedHadTopAndTopPair::isFromTop(HepMC::ConstGenParticlePtr part) const{ if(!prod) return false; +#ifdef HEPMC3 + for (auto p: prod->particles_in()) if (std::abs(p->pdg_id()) == 6) return true; +#else HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::parents); HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::parents); for(;firstParent!=endParent; ++firstParent){ //if( part->barcode() < (*firstParent)->barcode() ) continue; /// protection for sherpa if( abs( (*firstParent)->pdg_id() ) == 6 ) return true; } +#endif return false; } @@ -119,12 +123,16 @@ HepMC::ConstGenParticlePtr BoostedHadTopAndTopPair::findInitial(HepMC::ConstGen if(!prod) return part; +#ifdef HEPMC3 + for (auto p: prod->particles_in()) if (part->pdg_id() == p->pdg_id()) return findInitial(part); +#else HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::parents); HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::parents); for(;firstParent!=endParent; ++firstParent){ //if( part->barcode() < (*firstParent)->barcode() ) continue; /// protection for sherpa if( part->pdg_id() == (*firstParent)->pdg_id() ) return findInitial(*firstParent); } +#endif return part; } @@ -132,12 +140,10 @@ HepMC::ConstGenParticlePtr BoostedHadTopAndTopPair::findInitial(HepMC::ConstGen bool BoostedHadTopAndTopPair::isHadronic(HepMC::ConstGenParticlePtr part) const{ auto end = part->end_vertex(); - - HepMC::GenVertex::particle_iterator firstChild = end->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator endChild = end->particles_end (HepMC::children); - for(;firstChild!=endChild; ++firstChild){ - //if( part->barcode() > (*firstChild)->barcode() ) continue; /// protection for sherpa - if( abs((*firstChild)->pdg_id()) <= 5 ) return true; + if (end) { + for(auto firstChild: *end){ + if( std::abs(firstChild->pdg_id()) <= 5 ) return true; + } } return false; } @@ -148,12 +154,8 @@ bool BoostedHadTopAndTopPair::isFinalParticle(HepMC::ConstGenParticlePtr part) c auto end = part->end_vertex(); if(end){ int type = part->pdg_id(); - HepMC::GenVertex::particle_iterator firstChild = end->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator endChild = end->particles_end(HepMC::children); - for(;firstChild!=endChild; ++firstChild){ - //if( part->barcode() > (*firstChild)->barcode() ) continue; /// protection for sherpa - int childtype = (*firstChild)->pdg_id(); - if( childtype == type ) return false; + for(auto firstChild: *end){ + if( firstChild->pdg_id() == type ) return false; } } return true; @@ -166,14 +168,12 @@ HepMC::FourVector BoostedHadTopAndTopPair::momentumBofW(HepMC::ConstGenParticleP auto prod = initpart->production_vertex(); HepMC::FourVector b(0,0,0,0); - - HepMC::GenVertex::particle_iterator firstChild = prod->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator endChild = prod->particles_end(HepMC::children); - for(;firstChild!=endChild; ++firstChild){ - //if( part->barcode() > (*firstChild)->barcode() ) continue; /// protection for sherpa - if( abs( (*firstChild)->pdg_id() ) == 5 ){ - b.set((*firstChild)->momentum().x(), (*firstChild)->momentum().y(), (*firstChild)->momentum().z(), (*firstChild)->momentum().t()); +if (prod) { + for( auto firstChild: *prod){ + if( std::abs( firstChild->pdg_id() ) == 5 ){ + b.set(firstChild->momentum().x(), firstChild->momentum().y(), firstChild->momentum().z(), firstChild->momentum().t()); } } +} return b; } diff --git a/Generators/GeneratorFilters/src/DecayModeFilter.cxx b/Generators/GeneratorFilters/src/DecayModeFilter.cxx index d8624c3c88f6ba809cd23907dc9220e8e2cdce13..3ac09b477f8cdd478fce67ad851d5886922cb7fc 100644 --- a/Generators/GeneratorFilters/src/DecayModeFilter.cxx +++ b/Generators/GeneratorFilters/src/DecayModeFilter.cxx @@ -116,13 +116,24 @@ StatusCode DecayModeFilter::filterEvent() { McEventCollection::const_iterator itr; for (itr = events()->begin(); itr != events()->end(); ++itr) { const HepMC::GenEvent* theGenEvent = *(events()->begin()); +#ifdef HEPMC3 + auto vtx_begin = ((HepMC::GenEvent*)theGenEvent)->vertices().begin(); + auto vtx_end = ((HepMC::GenEvent*)theGenEvent)->vertices().end(); +#else HepMC::GenEvent::vertex_const_iterator vtx_begin = theGenEvent->vertices_begin(); HepMC::GenEvent::vertex_const_iterator vtx_end = theGenEvent->vertices_end(); +#endif for (auto vtx_iter=vtx_begin; vtx_iter!=vtx_end; ++vtx_iter) { // Look for initial vertex +#ifdef HEPMC3 + if ((*vtx_iter)->particles_in().size() != 2) continue; + if ((*vtx_iter)->particles_out().size() < 2) continue; + auto outParticle = (*vtx_iter)->particles_out().begin(); +#else if ((*vtx_iter)->particles_in_size() != 2) continue; if ((*vtx_iter)->particles_out_size() < 2) continue; auto outParticle = (*vtx_iter)->particles_out_const_begin(); +#endif auto parent1 = *outParticle; auto parent2 = *(++outParticle); ATH_MSG_DEBUG("Two in, two out: " << parent1->pdg_id() << " " << parent2->pdg_id()); @@ -185,15 +196,15 @@ StatusCode DecayModeFilter::filterEvent() { string DecayModeFilter::printChain(HepMC::GenParticlePtr parent) const { + if (!parent) return std::string(""); std::stringstream ss; - ss << " " << abs(parent->pdg_id()) << " -> "; + ss << " " << std::abs(parent->pdg_id()) << " -> "; HepMC::GenParticlePtr foundChild(NULL); int SMchild_PDG(0); - HepMC::GenVertex::particle_iterator child = parent->end_vertex()->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator end_child = parent->end_vertex()->particles_end(HepMC::children); - for (; child != end_child; ++child) { - if (abs((*child)->pdg_id()) < 1000) SMchild_PDG = abs((*child)->pdg_id()); - else foundChild = *child; + if (!parent->end_vertex()) return ss.str(); + for ( auto child: *(parent->end_vertex())) { + if (std::abs(child->pdg_id()) < 1000) SMchild_PDG = std::abs(child->pdg_id()); + else foundChild = child; } ss << (foundChild ? abs(foundChild->pdg_id()) : -9999) << " (" << SMchild_PDG << ") "; if (foundChild && abs(foundChild->pdg_id()) != 1000022) ss << printChain(foundChild); @@ -242,14 +253,16 @@ void DecayModeFilter::countChain(HepMC::GenParticlePtr parent, int& length, int& Nchi2, int& NW,int& NZ,int& NH, int& Nse, int& Nsmu, int& Nstau, int& nChargedLeptons, int& nSMParticles) const { HepMC::GenParticlePtr foundChild=nullptr; int SMchild_PDG(0); - HepMC::GenVertex::particle_iterator child = parent->end_vertex()->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator end_child = parent->end_vertex()->particles_end(HepMC::children); - for (; child != end_child; ++child) { - if (abs((*child)->pdg_id()) < 1000){ - SMchild_PDG = abs((*child)->pdg_id()); + if (!parent) return; + if (parent->end_vertex()) + { + for (auto child: *(parent->end_vertex())) { + if (std::abs(child->pdg_id()) < 1000){ + SMchild_PDG = std::abs(child->pdg_id()); nSMParticles++; } - if (abs((*child)->pdg_id()) > 1000000) foundChild = *child; + if (std::abs(child->pdg_id()) > 1000000) foundChild = child; + } } if (!foundChild) { length = 0; diff --git a/Generators/GeneratorFilters/src/DecayPositionFilter.cxx b/Generators/GeneratorFilters/src/DecayPositionFilter.cxx index dec2072405d6743a54dee291d10005201776afd8..c1005e3f6f7b73215d51fcc748fedd3b2719b356 100644 --- a/Generators/GeneratorFilters/src/DecayPositionFilter.cxx +++ b/Generators/GeneratorFilters/src/DecayPositionFilter.cxx @@ -35,32 +35,25 @@ StatusCode DecayPositionFilter::filterEvent() { int nPass = 0; for (McEventCollection::const_iterator itr = events()->begin(); itr != events()->end(); ++itr) { const HepMC::GenEvent* genEvt = (*itr); - for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) { - if (abs((*pitr)->pdg_id()) == m_PDGID) { - if ( fabs((*pitr)->momentum().pseudoRapidity()) ) { - + for (auto pitr: *genEvt) { + if (std::abs(pitr->pdg_id()) != m_PDGID) continue; + if ( !std::abs( pitr->momentum().pseudoRapidity()) ) continue;// AV:WHAT IS THIS? // Count only particles not decaying to themselves - bool notSelfDecay = true; - if ((*pitr)->end_vertex()) { - HepMC::GenVertex::particle_iterator child = (*pitr)->end_vertex()->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator childE = (*pitr)->end_vertex()->particles_end(HepMC::children); - for (; child != childE; ++child) { - if ( (*child)->pdg_id() == (*pitr)->pdg_id() && (*child)->barcode()!=(*pitr)->barcode() && (*child)->barcode() < 100000) { + bool notSelfDecay = true; + if (!pitr->end_vertex()) continue; + for ( auto child: *(pitr->end_vertex())) { + if ( child->pdg_id() == pitr->pdg_id() && HepMC::barcode(child)!=HepMC::barcode(pitr) && HepMC::barcode(child) < 100000) { notSelfDecay = false; break; } - } - if(notSelfDecay){ - - HepMC::GenVertexPtr vtx = (*pitr)->end_vertex(); - float x = vtx->position().x(); - float y = vtx->position().y(); - float Rdecay = sqrt(x*x + y*y); - if(Rdecay < m_RCut ) nPass++; - } - } - } - } + } + if(notSelfDecay){ + auto vtx = pitr->end_vertex(); + double x = vtx->position().x(); + double y = vtx->position().y(); + double Rdecay = std::sqrt(x*x + y*y); + if(Rdecay < m_RCut ) nPass++; + } } } setFilterPassed(nPass == m_MinPass); diff --git a/Generators/GeneratorFilters/src/DecaysFinalStateFilter.cxx b/Generators/GeneratorFilters/src/DecaysFinalStateFilter.cxx index f68e50cee483de745e6a9fd1b17d8df1498996e6..667dd2cc81f11934684f5b206de849276e42d463 100644 --- a/Generators/GeneratorFilters/src/DecaysFinalStateFilter.cxx +++ b/Generators/GeneratorFilters/src/DecaysFinalStateFilter.cxx @@ -56,23 +56,19 @@ StatusCode DecaysFinalStateFilter::filterEvent() { for (itr = events()->begin(); itr != events()->end(); ++itr) { const HepMC::GenEvent* genEvt = *itr; - // Loop over all particles in event - HepMC::GenEvent::particle_const_iterator pitr; - for (pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr ) { + for (auto part: *genEvt){ // look only at the allowed parents (e.g. W, Z) bool allowedParent = false; for (size_t i=0; i<m_PDGAllowedParents.size(); ++i) { - if ( (*pitr)->pdg_id() == m_PDGAllowedParents[i]) allowedParent = true; + if ( part->pdg_id() == m_PDGAllowedParents[i]) allowedParent = true; } if (!allowedParent) continue; - auto vtx = (*pitr)->end_vertex(); - if (!vtx) continue; + if (!part->end_vertex()) continue; - for (HepMC::GenVertex::particles_out_const_iterator opitr = vtx->particles_out_const_begin(); - opitr != vtx->particles_out_const_end(); ++opitr) { - int apid = abs((*opitr)->pdg_id()); + for (auto opitr: *(part->end_vertex())) { + int apid = std::abs(opitr->pdg_id()); if (apid == 1 || apid == 2 || apid == 3 || apid == 4 || apid ==5) nQuarks++; if (apid == 11 || apid == 13 || apid == 15) nChargedLeptons++; if (apid == 12 || apid == 14 || apid == 16) nNeutrinos++; diff --git a/Generators/GeneratorFilters/src/HTFilter.cxx b/Generators/GeneratorFilters/src/HTFilter.cxx index 2ba7c8e4ebeadcefb677a905d8a5e49b8758e558..8f161fd3db6ae97e409cdc807241ee464c4e698c 100644 --- a/Generators/GeneratorFilters/src/HTFilter.cxx +++ b/Generators/GeneratorFilters/src/HTFilter.cxx @@ -110,25 +110,25 @@ StatusCode HTFilter::filterEvent() { std::vector<HepMC::GenParticlePtr> WZleptons; WZleptons.reserve(10); - for (HepMC::GenEvent::particle_const_iterator iter=(*mecc)[0]->particles_begin(); iter!=(*mecc)[0]->particles_end();++iter){ - if ( !(*iter) ) continue; - int pdgid = (*iter)->pdg_id(); - if (m_UseNu && MC::PID::isNeutrino(pdgid) && MC::isGenStable(*iter)) { - if( fromWZ(*iter) || fromTau(*iter) ) { - HT += (*iter)->momentum().perp(); + for (auto iter: *((*mecc)[0])){ + if ( !iter ) continue; + int pdgid = iter->pdg_id(); + if (m_UseNu && MC::PID::isNeutrino(pdgid) && MC::isGenStable(iter)) { + if( fromWZ(iter) || fromTau(iter) ) { + HT += iter->momentum().perp(); } } // pick muons and electrons specifically -- isLepton selects both charged leptons and neutrinos - if (m_UseLep && (abs(pdgid)==11 || abs(pdgid)==13) && MC::isGenStable(*iter) - && (*iter)->momentum().perp()>m_MinLepPt*CLHEP::GeV && fabs((*iter)->momentum().eta())<m_MaxLepEta) { - bool isFromWZ = fromWZ(*iter); - if( isFromWZ || fromTau(*iter) ) { - ATH_MSG_VERBOSE("Adding W/Z/tau lepton with pt " << (*iter)->momentum().perp() - << ", eta " << (*iter)->momentum().eta() - << ", phi " << (*iter)->momentum().phi() - << ", status " << (*iter)->status() + if (m_UseLep && (std::abs(pdgid)==11 || std::abs(pdgid)==13) && MC::isGenStable(iter) + && (iter)->momentum().perp()>m_MinLepPt*CLHEP::GeV && std::abs(iter->momentum().eta())<m_MaxLepEta) { + bool isFromWZ = fromWZ(iter); + if( isFromWZ || fromTau(iter) ) { + ATH_MSG_VERBOSE("Adding W/Z/tau lepton with pt " << iter->momentum().perp() + << ", eta " << iter->momentum().eta() + << ", phi " << iter->momentum().phi() + << ", status " << iter->status() << ", pdgId " << pdgid); - HT += (*iter)->momentum().perp(); + HT += iter->momentum().perp(); } } } @@ -149,7 +149,7 @@ StatusCode HTFilter::filterEvent() { return StatusCode::SUCCESS; } -bool HTFilter::fromWZ( const HepMC::GenParticlePtr part ) const +bool HTFilter::fromWZ(HepMC::ConstGenParticlePtr part ) const { // !!! IMPORTANT !!! This is a TEMPORARY function // it's used in place of code in MCTruthClassifier as long as this package is not dual-use @@ -163,6 +163,15 @@ bool HTFilter::fromWZ( const HepMC::GenParticlePtr part ) const // generators do not include the W or Z in the truth record (like Sherpa) // This code, like the code before it, really assumes one incoming particle per vertex... if (!part->production_vertex()) return false; +#ifdef HEPMC3 + for (auto iter: part->production_vertex()->particles_in()){ + int parent_pdgid = iter->pdg_id(); + if (MC::PID::isW(parent_pdgid) || MC::PID::isZ(parent_pdgid)) return true; + if (MC::PID::isHadron( parent_pdgid ) ) return false; + if ( std::abs( parent_pdgid ) < 9 ) return true; + if ( parent_pdgid == part->pdg_id() ) return fromWZ( iter ); + } +#else for (HepMC::GenVertex::particles_in_const_iterator iter=part->production_vertex()->particles_in_const_begin(); iter!=part->production_vertex()->particles_in_const_end();++iter){ int parent_pdgid = (*iter)->pdg_id(); @@ -171,10 +180,11 @@ bool HTFilter::fromWZ( const HepMC::GenParticlePtr part ) const if ( abs( parent_pdgid ) < 9 ) return true; if ( parent_pdgid == part->pdg_id() ) return fromWZ( *iter ); } +#endif return false; } -bool HTFilter::fromTau( const HepMC::GenParticlePtr part ) const +bool HTFilter::fromTau(HepMC::ConstGenParticlePtr part ) const { // !!! IMPORTANT !!! This is a TEMPORARY function // it's used in place of code in MCTruthClassifier as long as this package is not dual-use @@ -186,6 +196,14 @@ bool HTFilter::fromTau( const HepMC::GenParticlePtr part ) const // Find a hadron or parton -> return false // This code, like the code before it, really assumes one incoming particle per vertex... if (!part->production_vertex()) return false; +#ifdef HEPMC3 + for (auto iter: part->production_vertex()->particles_in()){ + int parent_pdgid = iter->pdg_id(); + if ( std::abs( parent_pdgid ) == 15 ) return true; + if (MC::PID::isHadron( parent_pdgid ) || std::abs( parent_pdgid ) < 9 ) return false; + if ( parent_pdgid == part->pdg_id() ) return fromTau( iter ); + } +#else for (HepMC::GenVertex::particles_in_const_iterator iter=part->production_vertex()->particles_in_const_begin(); iter!=part->production_vertex()->particles_in_const_end();++iter){ int parent_pdgid = (*iter)->pdg_id(); @@ -193,5 +211,6 @@ bool HTFilter::fromTau( const HepMC::GenParticlePtr part ) const if (MC::PID::isHadron( parent_pdgid ) || abs( parent_pdgid ) < 9 ) return false; if ( parent_pdgid == part->pdg_id() ) return fromTau( *iter ); } +#endif return false; } diff --git a/Generators/GeneratorFilters/src/HiggsFilter.cxx b/Generators/GeneratorFilters/src/HiggsFilter.cxx index bf11e5b3a072fd0e91b95d16b942380e7df25993..be99a51b7f67ecbcfe86e02ea97311279ebc3b7f 100644 --- a/Generators/GeneratorFilters/src/HiggsFilter.cxx +++ b/Generators/GeneratorFilters/src/HiggsFilter.cxx @@ -58,49 +58,27 @@ StatusCode HiggsFilter::filterEvent() { // Loop over all particles in the event const HepMC::GenEvent* genEvt = (*itr); - HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); - for(; pitr!=genEvt->particles_end(); ++pitr ){ - - // Work only with tops - if( std::abs((*pitr)->pdg_id()) == 25 ){ - - N_Higgs_all++; - - int n_daughters = 0; - - HepMC::GenParticle * mcpart = (*pitr); - const HepMC::GenVertex * decayVtx = mcpart->end_vertex(); - + for(auto pitr: *genEvt){ + if( std::abs(pitr->pdg_id()) != 25 ) continue; + N_Higgs_all++; + auto decayVtx = pitr->end_vertex(); // verify if we got a valid pointer and retrieve the number of daughters - if ( decayVtx != 0 ) { - n_daughters = decayVtx->particles_out_size(); - } - - // For this analysis we are not interested in t->t MC structures, only in decays - if( n_daughters >= 2 ){ - - HepMC::GenVertex::particles_in_const_iterator child_mcpartItr = decayVtx->particles_out_const_begin(); - HepMC::GenVertex::particles_in_const_iterator child_mcpartItrE = decayVtx->particles_out_const_end(); - - for (; child_mcpartItr != child_mcpartItrE; ++child_mcpartItr) { - - HepMC::GenParticle * child_mcpart = (*child_mcpartItr); - - // Implicitly, I assume that tops always decay to W X - if ( std::abs(child_mcpart->pdg_id()) == 5 ){ - - if ( (*pitr)->pdg_id() == 25 ) { + if (!decayVtx ) continue; +#ifdef HEPMC3 + int n_daughters = decayVtx->particles_out().size(); +#else + int n_daughters = decayVtx->particles_out_size(); +#endif + if( n_daughters < 2 ) continue; + for ( auto child_mcpart: *decayVtx) { + if ( std::abs(child_mcpart->pdg_id()) != 5 ) continue; + if ( pitr->pdg_id() == 25 ) { N_Higgs++; } - - if( ((*pitr)->momentum().perp() >=m_Ptmin && (*pitr)->momentum().perp() <m_Ptmax ) ){ + if( (pitr->momentum().perp() >=m_Ptmin && pitr->momentum().perp() <m_Ptmax ) ){ N_pt_above_cut++; } - - }//b bbar decay - }//child loop - }// only decaying Higgs - }//Higgs + }//child loop } } diff --git a/Generators/GeneratorFilters/src/MissingEtFilter.cxx b/Generators/GeneratorFilters/src/MissingEtFilter.cxx index c15a979a9fb6b4bc72e6c453c0148828b66f67df..fa478d14add76066edb0514feafa43e4211cc53f 100644 --- a/Generators/GeneratorFilters/src/MissingEtFilter.cxx +++ b/Generators/GeneratorFilters/src/MissingEtFilter.cxx @@ -52,6 +52,15 @@ bool MissingEtFilter::fromWZ( HepMC::ConstGenParticlePtr part ) const // generators do not include the W or Z in the truth record (like Sherpa) // This code, like the code before it, really assumes one incoming particle per vertex... if (!part->production_vertex()) return false; +#ifdef HEPMC3 + for (auto iter: part->production_vertex()->particles_in()){ + int parent_pdgid = iter->pdg_id(); + if (MC::PID::isW(parent_pdgid) || MC::PID::isZ(parent_pdgid)) return true; + if (MC::PID::isHadron( parent_pdgid ) ) return false; + if ( std::abs( parent_pdgid ) < 9 ) return true; + if ( parent_pdgid == part->pdg_id() ) return fromWZ( iter ); + } +#else for (HepMC::GenVertex::particles_in_const_iterator iter=part->production_vertex()->particles_in_const_begin(); iter!=part->production_vertex()->particles_in_const_end();++iter){ int parent_pdgid = (*iter)->pdg_id(); @@ -60,6 +69,7 @@ bool MissingEtFilter::fromWZ( HepMC::ConstGenParticlePtr part ) const if ( abs( parent_pdgid ) < 9 ) return true; if ( parent_pdgid == part->pdg_id() ) return fromWZ( *iter ); } +#endif return false; } @@ -75,6 +85,14 @@ bool MissingEtFilter::fromTau( HepMC::ConstGenParticlePtr part ) const // Find a hadron or parton -> return false // This code, like the code before it, really assumes one incoming particle per vertex... if (!part->production_vertex()) return false; +#ifdef HEPMC3 + for (auto iter: part->production_vertex()->particles_in()){ + int parent_pdgid = iter->pdg_id(); + if ( std::abs( parent_pdgid ) == 15 && fromWZ(iter)) return true; + if (MC::PID::isHadron( parent_pdgid ) || abs( parent_pdgid ) < 9 ) return false; + if ( parent_pdgid == part->pdg_id() ) return fromTau( iter ); + } +#else for (HepMC::GenVertex::particles_in_const_iterator iter=part->production_vertex()->particles_in_const_begin(); iter!=part->production_vertex()->particles_in_const_end();++iter){ int parent_pdgid = (*iter)->pdg_id(); @@ -82,5 +100,6 @@ bool MissingEtFilter::fromTau( HepMC::ConstGenParticlePtr part ) const if (MC::PID::isHadron( parent_pdgid ) || abs( parent_pdgid ) < 9 ) return false; if ( parent_pdgid == part->pdg_id() ) return fromTau( *iter ); } +#endif return false; } diff --git a/Generators/GeneratorFilters/src/MultiHiggsFilter.cxx b/Generators/GeneratorFilters/src/MultiHiggsFilter.cxx index 7b22b594f757f6e829f60e37e1d7e7f635f9607e..13adf90914596ab15c40a86996ec99ab11e5b2ba 100644 --- a/Generators/GeneratorFilters/src/MultiHiggsFilter.cxx +++ b/Generators/GeneratorFilters/src/MultiHiggsFilter.cxx @@ -75,6 +75,20 @@ StatusCode MultiHiggsFilter::filterEvent() { for (itr = events()->begin(); itr!=events()->end(); ++itr) { // Loop over all particles in the event const HepMC::GenEvent* genEvt = (*itr); +#ifdef HEPMC3 + for(auto pitr: *genEvt) { + if (!( std::abs(pitr->pdg_id()) == 25 && ( pitr->status()==m_Status || !m_UseStatus))) continue; + bool hasHiggsParent = false; + for(auto thisMother: pitr->production_vertex()->particles_in()){ + ATH_MSG_DEBUG(" Parent " << pitr->pdg_id() << " barcode = " << HepMC::barcode(pitr) << " status = " << pitr->status()); + ATH_MSG_DEBUG(" a Parent mother " << thisMother->pdg_id()<< " barc = " << HepMC::barcode(thisMother)); + if( thisMother->pdg_id() == 25 ) hasHiggsParent = true; + } + ATH_MSG_DEBUG(" has Higgs parent? " << hasHiggsParent); + if (hasHiggsParent) continue; + ++nGoodParent; + } +#else for(HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr ) { if( abs((*pitr)->pdg_id()) == 25 && ( (*pitr)->status()==m_Status || !m_UseStatus)){ @@ -94,6 +108,7 @@ StatusCode MultiHiggsFilter::filterEvent() { ++nGoodParent; } } +#endif } ATH_MSG_INFO("Result " << nGoodParent << " Higgs "); diff --git a/Generators/GeneratorFilters/src/ParentTwoChildrenFilter.cxx b/Generators/GeneratorFilters/src/ParentTwoChildrenFilter.cxx index 555a27ce2698b0803c91ca0d5524b4770e4505db..0e760cb1546ad0c1b9faf70078ce1922389c5a04 100644 --- a/Generators/GeneratorFilters/src/ParentTwoChildrenFilter.cxx +++ b/Generators/GeneratorFilters/src/ParentTwoChildrenFilter.cxx @@ -40,35 +40,33 @@ StatusCode ParentTwoChildrenFilter::filterEvent() { } for (McEventCollection::const_iterator itr = events()->begin(); itr != events()->end(); ++itr) { const HepMC::GenEvent* genEvt = (*itr); - for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) { - int id = (*pitr)->pdg_id(); - if (abs(id) != m_PDGParent[0]) continue; - if ((*pitr)->momentum().perp() < m_PtMinParent) continue; + for (auto pitr: *genEvt) { + int id = pitr->pdg_id(); + if (std::abs(id) != m_PDGParent[0]) continue; + if (pitr->momentum().perp() < m_PtMinParent) continue; n_parents++; - - int n_daughters = 0; - HepMC::GenParticle* mcpart = (*pitr); - const HepMC::GenVertex* decayVtx = mcpart->end_vertex(); + HepMC::ConstGenVertexPtr decayVtx = pitr->end_vertex(); // Verify if we got a valid pointer and retrieve the number of daughters - if (decayVtx != 0) n_daughters = decayVtx->particles_out_size(); - if (n_daughters >= 2) { - HepMC::GenVertex::particle_iterator firstChild = (*pitr)->end_vertex()->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator endChild = (*pitr)->end_vertex()->particles_end(HepMC::children); - HepMC::GenVertex::particle_iterator thisChild = firstChild; - for(; thisChild != endChild; ++thisChild) { - ATH_MSG_DEBUG(" ParentTwoChildrenFilter: parent ==> " <<(*pitr)->pdg_id() << " child ===> " <<(*thisChild)->pdg_id()); + if (!decayVtx) continue; +#ifdef HEPMC3 + int n_daughters = decayVtx->particles_out().size(); +#else + int n_daughters = decayVtx->particles_out_size(); +#endif + if (n_daughters < 2) continue; + for( auto thisChild: *decayVtx ) { + ATH_MSG_DEBUG(" ParentTwoChildrenFilter: parent ==> " <<pitr->pdg_id() << " child ===> " <<thisChild->pdg_id()); for (int i = 0; i < 2; i++) { - if ( abs((*thisChild)->pdg_id()) == m_PDGChild[i] ){ - if ( (*thisChild)->pdg_id() == m_PDGChild[i] ){ - if( ((*thisChild)->momentum().perp() >= m_PtMinChild) )N_Child[i][0]++; + if ( std::abs(thisChild->pdg_id()) == m_PDGChild[i] ){ + if ( thisChild->pdg_id() == m_PDGChild[i] ){ + if (thisChild->momentum().perp() >= m_PtMinChild) N_Child[i][0]++; } - if ( (*thisChild)->pdg_id() == -m_PDGChild[i] ){ - if( ((*thisChild)->momentum().perp() >= m_PtMinChild) )N_Child[i][1]++; + if ( thisChild->pdg_id() == -m_PDGChild[i] ){ + if (thisChild->momentum().perp() >= m_PtMinChild) N_Child[i][1]++; } } } } - } } } setFilterPassed(N_Child[0][0] >= 1 && N_Child[0][1] >= 1 && N_Child[1][0] >= 1 && N_Child[1][1] >= 1); diff --git a/Generators/GeneratorFilters/src/ParticleFilter.cxx b/Generators/GeneratorFilters/src/ParticleFilter.cxx index 7964aa90e54b7fb8cf9f5f62dd4b8a7f45a0e615..72727d9a290ee4ba0c81c723ea652c187c4afa97 100644 --- a/Generators/GeneratorFilters/src/ParticleFilter.cxx +++ b/Generators/GeneratorFilters/src/ParticleFilter.cxx @@ -35,21 +35,16 @@ StatusCode ParticleFilter::filterEvent() { int nParts = 0; for (McEventCollection::const_iterator itr = events()->begin(); itr != events()->end(); ++itr) { const HepMC::GenEvent* genEvt = (*itr); - for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) { - if (abs((*pitr)->pdg_id()) == m_PDGID && (m_StatusReq == -1 || (*pitr)->status() == m_StatusReq)) { - if ((*pitr)->momentum().perp() >= m_Ptmin && - fabs((*pitr)->momentum().pseudoRapidity()) <= m_EtaRange && - (*pitr)->momentum().e() <= m_EnergyRange) { - if ((!m_Exclusive)&&(m_MinParts == 1)) // Found at least one particle and we have an inclusive requirement - return StatusCode::SUCCESS; - else { + for (auto pitr: *genEvt) { + if (std::abs(pitr->pdg_id()) != m_PDGID || !(m_StatusReq == -1 || pitr->status() == m_StatusReq)) 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()) { - HepMC::GenVertex::particle_iterator child = (*pitr)->end_vertex()->particles_begin(HepMC::children); - HepMC::GenVertex::particle_iterator childE = (*pitr)->end_vertex()->particles_end(HepMC::children); - for (; child != childE; ++child) { - if ( (*child)->pdg_id() == (*pitr)->pdg_id() && (*child)->barcode()!=(*pitr)->barcode() && (*child)->barcode() < 100000) { + if (pitr->end_vertex()) { + for (auto child: *(pitr->end_vertex())) { + if ( child->pdg_id() == pitr->pdg_id() && HepMC::barcode(child)!=HepMC::barcode(pitr) && HepMC::barcode(child) < 100000) { notSelfDecay = false; break; } @@ -57,8 +52,6 @@ StatusCode ParticleFilter::filterEvent() { } if (notSelfDecay) nParts++; } - } - } } } if (m_Exclusive) diff --git a/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx b/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx index 822c2ac92eebd80bf557db6685a82cf2dfc30ae0..b920071b026db39776020431b1d2cac829ad331b 100644 --- a/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx +++ b/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx @@ -67,9 +67,6 @@ StatusCode TtHtoVVDecayFilter::filterEvent() { for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) { for (size_t iParent=0; iParent < m_PDGParent.size(); iParent++) { - //ATH_MSG_INFO(" Parent to match? pdg = " << m_PDGParent[iParent] << " status = " << m_StatusParent[iParent] ); - //if ( abs((*pitr)->pdg_id()) == m_PDGParent[iParent] ) ATH_MSG_INFO(" particle to check pdg = " << abs((*pitr)->pdg_id()) << " status = " << (*pitr)->status() ); - if ( abs((*pitr)->pdg_id()) == m_PDGParent[iParent] && (*pitr)->status() == m_StatusParent[iParent]) { // if find a parent here HepMC::GenVertex::particle_iterator firstMother = (*pitr)->production_vertex()->particles_begin(HepMC::parents); HepMC::GenVertex::particle_iterator endMother = (*pitr)->production_vertex()->particles_end(HepMC::parents); @@ -126,18 +123,16 @@ StatusCode TtHtoVVDecayFilter::filterEvent() { } -bool TtHtoVVDecayFilter::findAncestor(const HepMC::GenVertexPtr searchvertex, +bool TtHtoVVDecayFilter::findAncestor(HepMC::ConstGenVertexPtr searchvertex, int targetPDGID) { if (!searchvertex) return false; - const HepMC::GenVertex::particles_out_const_iterator firstAncestor = searchvertex->particles_out_const_begin(); - const HepMC::GenVertex::particles_out_const_iterator endAncestor = searchvertex->particles_out_const_end(); - for (auto anc = firstAncestor; anc != endAncestor; ++anc) { - if (abs((*anc)->pdg_id()) == targetPDGID) { // same particle as parent - return findAncestor((*anc)->end_vertex(), + for (HepMC::ConstGenParticlePtr anc: *searchvertex) { + if (std::abs(anc->pdg_id()) == targetPDGID) { // same particle as parent + return findAncestor(anc->end_vertex(), targetPDGID); } else { for (size_t i = 0; i < m_PDGChild.size(); ++i) { - if (abs((*anc)->pdg_id()) == m_PDGChild[i]) return true; + if (std::abs(anc->pdg_id()) == m_PDGChild[i]) return true; } } } diff --git a/Generators/GeneratorModules/GeneratorModules/GenBase.h b/Generators/GeneratorModules/GeneratorModules/GenBase.h index f10e0aa5ae8be98d71e63d39b163090a18c058f5..854f2301dfc3f7c1573c896f1106ae9e1cd30fe6 100644 --- a/Generators/GeneratorModules/GeneratorModules/GenBase.h +++ b/Generators/GeneratorModules/GeneratorModules/GenBase.h @@ -15,6 +15,7 @@ #include "GeneratorObjects/McEventCollection.h" #include "AtlasHepMC/GenEvent.h" +#include "AtlasHepMC/GenVertex.h" #include "HepPDT/ParticleData.hh" #include "HepPDT/ParticleDataTable.hh"