diff --git a/Generators/GeneratorFilters/GeneratorFilters/DstD0K3piFilter.h b/Generators/GeneratorFilters/GeneratorFilters/DstD0K3piFilter.h
index 84adc3da8d67f70a563f24e2c4101f044c2eb33c..a5076dfee0862537f58963a05fce2c8f9833e8da 100644
--- a/Generators/GeneratorFilters/GeneratorFilters/DstD0K3piFilter.h
+++ b/Generators/GeneratorFilters/GeneratorFilters/DstD0K3piFilter.h
@@ -37,8 +37,8 @@ public:
private:
- bool CheckChildLundId(HepMC::GenParticlePtr mcpart, unsigned int nth, int chLundId); // Check lundId of Nth(0 or 1) children from mcpart...
- bool IsCandidate(std::vector<float>& lundIds, std::vector<HepMC::GenParticlePtr>& genParticles); // Is Candidate among 6 intermediate state??
+ bool CheckChildLundId(HepMC::ConstGenParticlePtr mcpart, unsigned int nth, int chLundId); // Check lundId of Nth(0 or 1) children from mcpart...
+ bool IsCandidate(std::vector<float>& lundIds, std::vector<HepMC::ConstGenParticlePtr>& genParticles); // Is Candidate among 6 intermediate state??
double m_Ptmin;
double m_EtaRange;
diff --git a/Generators/GeneratorFilters/src/DstD0K3piFilter.cxx b/Generators/GeneratorFilters/src/DstD0K3piFilter.cxx
index 4c086ac70e359e27a4ff1b55327367772db5cf7f..8a9e571ce06519bba66e1c405985df83f9d4aeb8 100644
--- a/Generators/GeneratorFilters/src/DstD0K3piFilter.cxx
+++ b/Generators/GeneratorFilters/src/DstD0K3piFilter.cxx
@@ -15,7 +15,16 @@ DstD0K3piFilter::DstD0K3piFilter(const std::string& name, ISvcLocator* pSvcLocat
// If mcpart first child lundId equals chLundId, return true
/// @todo No need for this to be a member function... useful in a HepMC utils library?
-bool DstD0K3piFilter::CheckChildLundId(HepMC::GenParticlePtr mcpart, unsigned int nth, int chLundId) {
+bool DstD0K3piFilter::CheckChildLundId(HepMC::ConstGenParticlePtr mcpart, unsigned int nth, int chLundId) {
+#ifdef HEPMC3
+ auto DecayVtx = mcpart->end_vertex();
+ if (!DecayVtx) return false;
+ auto children=DecayVtx->particles_out();
+ if (children.size() < 2) return false;
+ if (children.size() < nth) return false;
+ if (std::abs(children.at(nth)->pdg_id()) == chLundId) return true;
+ return false;
+#else
int nChild = 0;
const HepMC::GenVertex* DecayVtx = mcpart->end_vertex();
if (DecayVtx != 0) nChild = DecayVtx->particles_out_size();
@@ -33,11 +42,59 @@ bool DstD0K3piFilter::CheckChildLundId(HepMC::GenParticlePtr mcpart, unsigned in
if (std::abs(child_mcpart->pdg_id()) == chLundId) return true;
}
return false;
+#endif
}
/// @todo No need for this to be a member function
-bool DstD0K3piFilter::IsCandidate(std::vector<float>& lundIds, std::vector<HepMC::GenParticlePtr>& genParticles) {
+bool DstD0K3piFilter::IsCandidate(std::vector<float>& lundIds, std::vector<HepMC::ConstGenParticlePtr>& genParticles) {
+#ifdef HEPMC3
+ unsigned int nDecay = lundIds.size();
+ if (nDecay == 2) {
+ unsigned int id0 = std::abs( static_cast<int>(lundIds[0]) );
+ unsigned int id1 = std::abs( static_cast<int>(lundIds[1]) );
+ // 10323 211 k_1+ CLHEP::pi-
+ if ( id0 == 10323 && id1 == 211 ) { // // K_1+ CLHEP::pi-
+ if ( CheckChildLundId(genParticles[0], 0, 313) ) { // K*0 pi+
+ int nChild = 0;
+ auto DecayVtx = genParticles[0]->end_vertex();
+ if ( !DecayVtx) return false;
+ auto children=DecayVtx->particles_out();
+ if ( children.size() != 2 ) return false;
+ if ( CheckChildLundId(children[0], 0, 321) ) return true;
+ } else if ( CheckChildLundId(genParticles[0], 0, 321) ) { // K+ rho0 or K+ CLHEP::pi CLHEP::pi
+ if ( CheckChildLundId(genParticles[0], 1, 113) || CheckChildLundId(genParticles[0], 1, 211) ) return true;
+ }
+ } else if (id0 == 313 && id1 == 113) { // // K^*0 rho0
+ if ( CheckChildLundId(genParticles[0], 0, 321) ) // rho0 decays to CLHEP::pi+CLHEP::pi- 100% in decay.dec
+ return true;
+ } else if (id0 == 321 && id1 == 20213) { // // K+ a_1+
+ if ( CheckChildLundId(genParticles[1], 0, 113) ) // a_1+ -> rho0 (113) CLHEP::pi+
+ return true;
+ }
+
+ } else if ( nDecay == 3 ) {
+ unsigned int id0 = std::abs( static_cast<int>(lundIds[0]) );
+ unsigned int id1 = std::abs( static_cast<int>(lundIds[1]) );
+ unsigned int id2 = std::abs( static_cast<int>(lundIds[2]) );
+ if (id0 == 321 && id1 == 211 && id2 == 113) { // K+ CLHEP::pi- rho0
+ return true; // rho0 decays to CLHEP::pi+CLHEP::pi- 100% in decay.dec
+ }
+ if (id0 == 313 && id1 == 211 && id2 == 211 ) { // K^*+ CLHEP::pi- CLHEP::pi+
+ if ( CheckChildLundId(genParticles[0], 0, 321) ) return true;
+ }
+
+ } else if (nDecay ==4) { // K CLHEP::pi CLHEP::pi CLHEP::pi non-resonant
+ unsigned int id0 = std::abs( static_cast<int>(lundIds[0]) );
+ unsigned int id1 = std::abs( static_cast<int>(lundIds[1]) );
+ unsigned int id2 = std::abs( static_cast<int>(lundIds[2]) );
+ unsigned int id3 = std::abs( static_cast<int>(lundIds[3]) );
+ if (id0 == 321 && id1 == 211 && id2 == 211 && id3 == 211) {
+ return true;
+ }
+ }
+ return false;
+#else
unsigned int nDecay = lundIds.size();
if (nDecay == 2) {
unsigned int id0 = std::abs( static_cast<int>(lundIds[0]) );
@@ -85,6 +142,7 @@ bool DstD0K3piFilter::IsCandidate(std::vector<float>& lundIds, std::vector<HepMC
}
}
return false;
+#endif
}
@@ -92,6 +150,36 @@ StatusCode DstD0K3piFilter::filterEvent() {
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr != events()->end(); ++itr) { // Loop to select D*
const HepMC::GenEvent* genEvt = *itr;
+#ifdef HEPMC3
+ for (auto pitr: *genEvt) {
+ if (std::abs(pitr->pdg_id()) != 413) continue; // D*+
+ if (pitr->momentum().perp() < m_Ptmin) continue;
+ if (std::abs(pitr->momentum().pseudoRapidity()) > m_EtaRange) continue;
+ auto DstDecayVtx = pitr->end_vertex();
+ if (!DstDecayVtx) continue; // Check that we got a valid pointer and retrieve the number of daughters
+ auto DstChild = DstDecayVtx->particles_out();
+ ATH_MSG_DEBUG("D*+- meson found with Nchild = " << DstChild.size() << " and PDG ID = " << pitr->pdg_id());
+ if (DstChild.size() != 2) continue;
+ for ( auto mcpartD0: DstChild ) {
+ if (std::abs(mcpartD0->pdg_id()) != 421) continue; // D0
+ auto D0DecayVtx = mcpartD0->end_vertex();
+ if (!D0DecayVtx) continue; // Check that we got a valid pointer and retrieve the number of daughters
+ auto D0Child = D0DecayVtx->particles_out();
+ if (D0Child.size() > 4) continue; // For this analysis we are only interested in D*->D0 pi+
+ ATH_MSG_DEBUG("D0 meson found with Nchild = " << D0Child.size() << " and PDF ID = " << mcpartD0->pdg_id());
+
+ std::vector<float> lundIds;
+ std::vector<HepMC::ConstGenParticlePtr> genParticles;
+ for (auto grandchild: *D0DecayVtx){
+ genParticles.push_back(grandchild);
+ lundIds.push_back(grandchild->pdg_id());
+ }
+ if (IsCandidate(lundIds, genParticles)) {
+ return StatusCode::SUCCESS;
+ }
+ }
+ }
+#else
for (HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr) {
// Work only with D*
if (std::abs((*pitr)->pdg_id()) != 413) continue; // D*+
@@ -123,7 +211,7 @@ StatusCode DstD0K3piFilter::filterEvent() {
ATH_MSG_DEBUG("D0 meson found with Nchild = " << nD0Child << " and PDF ID = " << child_mcpart->pdg_id());
std::vector<float> lundIds;
- std::vector<HepMC::GenParticle*> genParticles;
+ std::vector<const HepMC::GenParticle*> genParticles;
if (D0DecayVtx) {
HepMC::GenVertex::particles_in_const_iterator grandchild_mcpartItr = D0DecayVtx->particles_out_const_begin();
@@ -146,6 +234,7 @@ StatusCode DstD0K3piFilter::filterEvent() {
}
}
}
+#endif
}
setFilterPassed(false);
return StatusCode::SUCCESS;
diff --git a/Generators/GeneratorFilters/src/TTbarBoostCatFilter.cxx b/Generators/GeneratorFilters/src/TTbarBoostCatFilter.cxx
index 7f3b3c9322f4e0783ed0606b01239c1f436dc30f..50e14e42b34970112662eb0c176455417f11b15e 100644
--- a/Generators/GeneratorFilters/src/TTbarBoostCatFilter.cxx
+++ b/Generators/GeneratorFilters/src/TTbarBoostCatFilter.cxx
@@ -61,13 +61,72 @@ StatusCode TTbarBoostCatFilter::filterEvent() {
if(m_LepPtmin*m_LepPtmax <0 && m_LepPtmin < 0 ) m_LepPtmin = 0.;
if(m_LepPtmin*m_LepPtmax <0 && m_LepPtmax < 0 ) m_LepPtmax = 14000000.; // 14 TeV
- std::vector<HepMC::GenParticlePtr> tops;
- std::vector<HepMC::GenParticlePtr> ws; // W from top decay (from tops)
- std::vector<HepMC::GenParticlePtr> leps; // e, mu, tau from W decay (from ws)
- std::vector<HepMC::GenParticlePtr> nus; // nutrino from W decay (from ws)
+ std::vector<HepMC::ConstGenParticlePtr> tops;
+ std::vector<HepMC::ConstGenParticlePtr> ws; // W from top decay (from tops)
+ std::vector<HepMC::ConstGenParticlePtr> leps; // e, mu, tau from W decay (from ws)
+ std::vector<HepMC::ConstGenParticlePtr> nus; // nutrino from W decay (from ws)
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent* genEvt = (*itr);
+#ifdef HEPMC3
+ for (auto pitr: *genEvt) {
+ if (std::abs(pitr->pdg_id()) != 6) continue;
+ if ( pitr->pdg_id() == 6 ) N_quark_t_all++;
+ if ( pitr->pdg_id() == -6 ) N_quark_tbar_all++;
+ auto decayVtx = pitr->end_vertex();
+ // Verify if we got a valid pointer and retrieve the number of daughters
+ if (!decayVtx) continue;
+ // For this analysis we are not interested in t->t MC structures, only in decays
+ for (auto child_mcpart: *decayVtx ) {
+ // Implicitly assume that tops always decay to W X
+ if (std::abs(child_mcpart->pdg_id()) == 24) {
+ if ( pitr->pdg_id() == 6 ){
+ N_quark_t++;
+ tops.push_back(pitr);
+ ws.push_back(child_mcpart);
+ }
+ if ( pitr->pdg_id() == -6 ){
+ N_quark_tbar++;
+ tops.push_back(pitr);
+ ws.push_back(child_mcpart);
+ }
+ bool useNextVertex = false;
+ auto w_decayVtx = child_mcpart->end_vertex();
+
+ while (w_decayVtx) {
+ useNextVertex = false;
+ for (auto grandchild_mcpart: *w_decayVtx) {
+ int grandchild_pid = grandchild_mcpart->pdg_id();
+ ATH_MSG_DEBUG("W (t/tbar) has " << w_decayVtx->particles_out().size() << " children and the pdg_id of the next is " << grandchild_pid);
+ // Check if the W's child is W again. If yes, then move to its next decay vertex in a decay tree
+ if (std::abs(grandchild_pid) == 24) {
+ w_decayVtx = grandchild_mcpart->end_vertex();
+
+ // If something wrong comes from truth...
+ if (!w_decayVtx) {
+ ATH_MSG_ERROR("A stable W is found... ");
+ break;
+ }
+ useNextVertex = true;
+ break;
+ }
+ if (std::abs(grandchild_pid) == 11 || std::abs(grandchild_pid) == 13 || std::abs(grandchild_pid) == 15) {
+ leps.push_back(grandchild_mcpart);
+ if (grandchild_mcpart->momentum().perp() >= m_Ptmin) N_pt_above_cut++;
+ // W decay lepton is found. Break loop over the decay product particles
+ // break;
+ }
+ if (abs(grandchild_pid) == 12 || abs(grandchild_pid) == 14 || abs(grandchild_pid) == 16) {
+ nus.push_back(grandchild_mcpart);
+ }
+ }
+ // If investigation of W's next decay vertex is not required then finish looking for leptons
+ if (!useNextVertex) break;
+ }
+ }
+ }
+ }
+#else
for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) {
if (std::abs((*pitr)->pdg_id()) == 6) {
if ( (*pitr)->pdg_id() == 6 ) N_quark_t_all++;
@@ -154,6 +213,7 @@ StatusCode TTbarBoostCatFilter::filterEvent() {
}
}
}
+#endif
}
if(tops.size()==2){
@@ -231,29 +291,21 @@ StatusCode TTbarBoostCatFilter::filterEvent() {
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
event++;
const HepMC::GenEvent* genEvt = (*itr);
- HepMC::GenEvent::particle_const_iterator mcpartItr = genEvt->particles_begin();
- HepMC::GenEvent::particle_const_iterator mcpartItrE = genEvt->particles_end();
- int part ( 0 );
- for (; mcpartItr != mcpartItrE; ++mcpartItr) {
+ int part=0;
+ for (auto mcpart: *genEvt ) {
part++;
- HepMC::GenParticle * mcpart = (*mcpartItr);
int pid = mcpart->pdg_id();
ATH_MSG_ERROR("In event (from MC collection) " << event << " particle number " << part << " has pdg_id = " << pid);
// retrieve decay vertex
- const HepMC::GenVertex * decayVtx = mcpart->end_vertex();
-
+ auto decayVtx = mcpart->end_vertex();
// verify if we got a valid pointer
- if ( decayVtx != 0 ) {
- 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();
- int part_child ( 0 );
- for (; child_mcpartItr != child_mcpartItrE; ++child_mcpartItr) {
- part_child++;
- HepMC::GenParticle * child_mcpart = (*child_mcpartItr);
- int child_pid = child_mcpart->pdg_id();
- ATH_MSG_ERROR(" child " << part_child << " with pdg_id = " << child_pid);
- }
+ if ( !decayVtx) continue;
+ int part_child =0 ;
+ for ( auto child_mcpart: *decayVtx) {
+ part_child++;
+ int child_pid = child_mcpart->pdg_id();
+ ATH_MSG_ERROR(" child " << part_child << " with pdg_id = " << child_pid);
}
}
}
diff --git a/Generators/GeneratorFilters/src/TTbarPlusHeavyFlavorFilter.cxx b/Generators/GeneratorFilters/src/TTbarPlusHeavyFlavorFilter.cxx
index 5092d2403ce78c454e302ca62946097e09d472a6..11f46f22778159d6c542eec75602166fe498685b 100644
--- a/Generators/GeneratorFilters/src/TTbarPlusHeavyFlavorFilter.cxx
+++ b/Generators/GeneratorFilters/src/TTbarPlusHeavyFlavorFilter.cxx
@@ -3,7 +3,7 @@
*/
#include "GeneratorFilters/TTbarPlusHeavyFlavorFilter.h"
-
+#include "AtlasHepMC/Relatives.h"
#include "GaudiKernel/MsgStream.h"
//--------------------------------------------------------------------------
@@ -226,6 +226,15 @@ bool TTbarPlusHeavyFlavorFilter::isInitialHadron(HepMC::ConstGenParticlePtr part
auto prod = part->production_vertex();
if(prod){
int type = hadronType(part->pdg_id());
+#ifdef HEPMC3
+ for(auto firstParent: prod->particles_in()){
+ if( HepMC::barcode(part) < HepMC::barcode(firstParent) ) continue; /// protection for sherpa
+ int mothertype = hadronType( firstParent->pdg_id() );
+ if( mothertype == type ){
+ return false;
+ }
+ }
+#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){
@@ -235,6 +244,7 @@ bool TTbarPlusHeavyFlavorFilter::isInitialHadron(HepMC::ConstGenParticlePtr part
return false;
}
}
+#endif
}
@@ -247,6 +257,15 @@ bool TTbarPlusHeavyFlavorFilter::isFinalHadron(HepMC::ConstGenParticlePtr part)
auto end = part->end_vertex();
if(end){
int type = hadronType(part->pdg_id());
+#ifdef HEPMC3
+ for(auto firstChild: end->particles_in()){
+ if( HepMC::barcode(part) > HepMC::barcode(firstChild) ) continue; /// protection for sherpa
+ int childtype = hadronType( firstChild->pdg_id() );
+ if( childtype == type ){
+ return false;
+ }
+ }
+#else
HepMC::GenVertex::particle_iterator firstChild = end->particles_begin(HepMC::children);
HepMC::GenVertex::particle_iterator endChild = end->particles_end(HepMC::children);
for(;firstChild!=endChild; ++firstChild){
@@ -256,6 +275,7 @@ bool TTbarPlusHeavyFlavorFilter::isFinalHadron(HepMC::ConstGenParticlePtr part)
return false;
}
}
+#endif
}
@@ -269,6 +289,15 @@ bool TTbarPlusHeavyFlavorFilter::isQuarkFromHadron(HepMC::ConstGenParticlePtr pa
auto prod = part->production_vertex();
if(prod){
+#ifdef HEPMC3
+ for(auto firstParent: HepMC::ancestor_particles(prod)){
+ if( HepMC::barcode(part) < HepMC::barcode(firstParent) ) continue; /// protection for sherpa
+ int mothertype = hadronType( firstParent->pdg_id() );
+ if( 4 == mothertype || 5 == mothertype ){
+ return true;
+ }
+ }
+#else
HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::ancestors);
HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::ancestors);
for(;firstParent!=endParent; ++firstParent){
@@ -278,6 +307,7 @@ bool TTbarPlusHeavyFlavorFilter::isQuarkFromHadron(HepMC::ConstGenParticlePtr pa
return true;
}
}
+#endif
}
@@ -291,6 +321,14 @@ bool TTbarPlusHeavyFlavorFilter::isCHadronFromB(HepMC::ConstGenParticlePtr part)
auto prod = part->production_vertex();
if(prod){
+#ifdef HEPMC3
+ for(auto firstParent:HepMC::ancestor_particles(prod)){
+ if( HepMC::barcode(part) < HepMC::barcode(firstParent) ) continue; /// protection for sherpa
+ if( isBHadron(firstParent) ){
+ return true;
+ }
+ }
+#else
HepMC::GenVertex::particle_iterator firstParent = prod->particles_begin(HepMC::ancestors);
HepMC::GenVertex::particle_iterator endParent = prod->particles_end(HepMC::ancestors);
for(;firstParent!=endParent; ++firstParent){
@@ -299,6 +337,7 @@ bool TTbarPlusHeavyFlavorFilter::isCHadronFromB(HepMC::ConstGenParticlePtr part)
return true;
}
}
+#endif
}
@@ -314,12 +353,19 @@ bool TTbarPlusHeavyFlavorFilter::isLooping(HepMC::ConstGenParticlePtr part, std:
init_part.insert(part);
+#ifdef HEPMC3
+ for(auto firstParent: prod->particles_in()){
+ if( init_part.find(firstParent) != init_part.end() ) return true;
+ if( isLooping(firstParent, init_part) ) 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( init_part.find(*firstParent) != init_part.end() ) return true;
if( isLooping(*firstParent, init_part) ) return true;
}
+#endif
return false;
@@ -333,6 +379,14 @@ HepMC::ConstGenParticlePtr TTbarPlusHeavyFlavorFilter::findInitial(HepMC::Const
if(!prod) return part;
+#ifdef HEPMC3
+ for(auto firstParent: prod->particles_in()){
+ if( HepMC::barcode(part) < HepMC::barcode(firstParent) && looping) continue; /// protection for sherpa
+ if( part->pdg_id() == firstParent->pdg_id() ){
+ return findInitial(firstParent, looping);
+ }
+ }
+#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){
@@ -341,6 +395,7 @@ HepMC::ConstGenParticlePtr TTbarPlusHeavyFlavorFilter::findInitial(HepMC::Const
return findInitial(*firstParent, looping);
}
}
+#endif
return part;
@@ -359,12 +414,19 @@ bool TTbarPlusHeavyFlavorFilter::isDirectlyFromTop(HepMC::ConstGenParticlePtr pa
if(!prod) return false;
+#ifdef HEPMC3
+ for( auto firstParent: prod->particles_in()){
+ if( HepMC::barcode(part) < HepMC::barcode(firstParent) && looping ) continue; /// protection for sherpa
+ if( std::abs( firstParent->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() && looping ) continue; /// protection for sherpa
if( abs( (*firstParent)->pdg_id() ) == 6 ) return true;
}
+#endif
return false;
}
@@ -377,6 +439,14 @@ bool TTbarPlusHeavyFlavorFilter::isDirectlyFromWTop(HepMC::ConstGenParticlePtr p
if(!prod) return false;
+#ifdef HEPMC3
+ for(auto firstParent: prod->particles_in()){
+ if( HepMC::barcode(part) < HepMC::barcode(firstParent) && looping ) continue; /// protection for sherpa
+ if( abs( firstParent->pdg_id() ) == 24 ){
+ if( isFromTop(firstParent, looping) ) 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){
@@ -385,6 +455,7 @@ bool TTbarPlusHeavyFlavorFilter::isDirectlyFromWTop(HepMC::ConstGenParticlePtr p
if( isFromTop(*firstParent, looping) ) return true;
}
}
+#endif
return false;
diff --git a/Generators/GeneratorFilters/src/TTbarWToLeptonFilter.cxx b/Generators/GeneratorFilters/src/TTbarWToLeptonFilter.cxx
index f23134a46c15ba880d70d7efb3a1798fc3a6be93..0276a6e5a5948c314bfa81e73b46c8b31cc0a80f 100644
--- a/Generators/GeneratorFilters/src/TTbarWToLeptonFilter.cxx
+++ b/Generators/GeneratorFilters/src/TTbarWToLeptonFilter.cxx
@@ -21,6 +21,53 @@ StatusCode TTbarWToLeptonFilter::filterEvent() {
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent* genEvt = (*itr);
+#ifdef HEPMC3
+ for (auto pitr: *genEvt) {
+ if (std::abs(pitr->pdg_id()) != 6) continue;
+ if ( pitr->pdg_id() == 6 ) N_quark_t_all++;
+ if ( pitr->pdg_id() == -6 ) N_quark_tbar_all++;
+ auto decayVtx = pitr->end_vertex();
+ // Verify if we got a valid pointer and retrieve the number of daughters
+ if (!decayVtx) continue;
+ // For this analysis we are not interested in t->t MC structures, only in decays
+ if (decayVtx->particles_out().size() < 2) continue;
+ for (auto child_mcpart: decayVtx->particles_out()) {
+ // Implicitly assume that tops always decay to W X
+ if (std::abs(child_mcpart->pdg_id()) != 24) continue;
+ if ( pitr->pdg_id() == 6 ) N_quark_t++;
+ if ( pitr->pdg_id() == -6 ) N_quark_tbar++;
+
+ bool useNextVertex = false;
+ auto w_decayVtx = child_mcpart->end_vertex();
+ while (w_decayVtx) {
+ useNextVertex = false;
+ for (auto grandchild_mcpart: *w_decayVtx) {
+ int grandchild_pid = grandchild_mcpart->pdg_id();
+ ATH_MSG_DEBUG("W (t/tbar) has " << w_decayVtx->particles_out().size() << " children and the pdg_id of the next is " << grandchild_pid);
+ // Check if the W's child is W again. If yes, then move to its next decay vertex in a decay tree
+ if (std::abs(grandchild_pid) == 24) {
+ w_decayVtx = grandchild_mcpart->end_vertex();
+ // If something wrong comes from truth...
+ if (!w_decayVtx) {
+ ATH_MSG_ERROR("A stable W is found... ");
+ break;
+ }
+ useNextVertex = true;
+ break;
+ }
+
+ if (std::abs(grandchild_pid) == 11 || std::abs(grandchild_pid) == 13 || abs(grandchild_pid) == 15) {
+ if (grandchild_mcpart->momentum().perp() >= m_Ptmin) N_pt_above_cut++;
+ // W decay lepton is found. Break loop over the decay product particles
+ break;
+ }
+ }
+ // If investigation of W's next decay vertex is not required then finish looking for leptons
+ if (!useNextVertex) break;
+ }
+ }
+ }
+#else
for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) {
if (std::abs((*pitr)->pdg_id()) == 6) {
if ( (*pitr)->pdg_id() == 6 ) N_quark_t_all++;
@@ -93,6 +140,7 @@ StatusCode TTbarWToLeptonFilter::filterEvent() {
}
}
}
+#endif
}
ATH_MSG_INFO("Found " << N_quark_t_all << " t quarks in event record");
@@ -113,30 +161,21 @@ StatusCode TTbarWToLeptonFilter::filterEvent() {
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
event++;
const HepMC::GenEvent* genEvt = (*itr);
- HepMC::GenEvent::particle_const_iterator mcpartItr = genEvt->particles_begin();
- HepMC::GenEvent::particle_const_iterator mcpartItrE = genEvt->particles_end();
- int part ( 0 );
- for (; mcpartItr != mcpartItrE; ++mcpartItr) {
+ int part=0 ;
+ for (auto mcpart: *genEvt) {
part++;
- HepMC::GenParticle * mcpart = (*mcpartItr);
int pid = mcpart->pdg_id();
ATH_MSG_ERROR("In event (from MC collection) " << event << " particle number " << part << " has pdg_id = " << pid);
-
// retrieve decay vertex
- const HepMC::GenVertex * decayVtx = mcpart->end_vertex();
-
+ auto decayVtx = mcpart->end_vertex();
// verify if we got a valid pointer
- if ( decayVtx != 0 ) {
- 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();
- int part_child ( 0 );
- for (; child_mcpartItr != child_mcpartItrE; ++child_mcpartItr) {
+ if ( !decayVtx ) continue;
+ int part_child=0;
+ for ( auto child_mcpart: *decayVtx) {
part_child++;
- HepMC::GenParticle * child_mcpart = (*child_mcpartItr);
int child_pid = child_mcpart->pdg_id();
ATH_MSG_ERROR(" child " << part_child << " with pdg_id = " << child_pid);
}
- }
}
}
setFilterPassed(false);
diff --git a/Generators/GeneratorFilters/src/TopCKMFilter.cxx b/Generators/GeneratorFilters/src/TopCKMFilter.cxx
index 6363fc2b8b3849bb2cc5d55d9f877a89b142818c..deec55083eb0eec6a8d4f738337fda6164f1156e 100644
--- a/Generators/GeneratorFilters/src/TopCKMFilter.cxx
+++ b/Generators/GeneratorFilters/src/TopCKMFilter.cxx
@@ -51,7 +51,7 @@ StatusCode TopCKMFilter::filterEvent() {
{
Wbosons[i]=Child->pdg_id();
if (!Child->end_vertex()) continue;
- for (auto gChild: *(Child->end_vertex()) if (std::abs(gChild->pdg_id()) == 15) isTau[i] = true;
+ for (auto gChild: *(Child->end_vertex())) if (std::abs(gChild->pdg_id()) == 15) isTau[i] = true;
}
}
}
diff --git a/Generators/GeneratorFilters/src/WZtoLeptonFilter.cxx b/Generators/GeneratorFilters/src/WZtoLeptonFilter.cxx
index c91f1f2002fd2175abe2caf3974cbf41a1df31e4..a93729e6880bd7988604cad6ee6c9ad5b517bef2 100644
--- a/Generators/GeneratorFilters/src/WZtoLeptonFilter.cxx
+++ b/Generators/GeneratorFilters/src/WZtoLeptonFilter.cxx
@@ -8,6 +8,7 @@
#include "McParticleEvent/TruthParticleContainer.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "GaudiKernel/ObjectList.h"
+#include "AtlasHepMC/Relatives.h"
WZtoLeptonFilter::WZtoLeptonFilter(const std::string& name, ISvcLocator* pSvcLocator)
: GenFilter( name,pSvcLocator ),
@@ -203,8 +204,6 @@ StatusCode WZtoLeptonFilter::filterFinalize() {
StatusCode WZtoLeptonFilter::filterEvent() {
- HepMC::GenVertexPtr LePrdVrt;
- HepMC::GenVertexPtr TauPrdVrt;
// Momentum of the products of the tau decay
CLHEP::HepLorentzVector mom_hadrons;
@@ -244,6 +243,103 @@ StatusCode WZtoLeptonFilter::filterEvent() {
if (wgtsC.size() > 0) wght = wgtsC[0];
m_tot_wghts += wght;
+#ifdef HEPMC3
+HepMC::ConstGenVertexPtr LePrdVrt;
+HepMC::ConstGenVertexPtr TauPrdVrt;
+ for ( auto pitr: *genEvt) {
+ LePrdVrt = 0;
+ TauPrdVrt = 0;
+ lepid = pitr->pdg_id();
+ abslepid = std::abs( lepid );
+
+ if ((abslepid != 11 && abslepid != 13) || pitr->status() != 1) continue;
+ double leppt = pitr->momentum().perp();
+ double lepeta = std::abs( pitr->momentum().pseudoRapidity() );
+
+ LePrdVrt = pitr->production_vertex();
+ int anceWZ = 0;
+ int momWZ = 0;
+ int taufromWZ = 0;
+
+ if (LePrdVrt != NULL) {
+ int ancecnt = 0;
+ for (auto lepanc: HepMC::ancestor_particles(LePrdVrt) ) {
+ int ancepid = lepanc->pdg_id();
+ int ancestatus = lepanc->status();
+ if ( m_tesit == 1 ) {
+ ATH_MSG_DEBUG("lepton=" << lepid << " " << ancecnt <<
+ "'th ancestors=" << ancepid << " status =" << ancestatus);
+ }
+ if ( std::abs(ancepid) == 24 || std::abs(ancepid) == 23 ) anceWZ ++;
+ ancecnt ++;
+ } // end of lepton ancestors( mother, grandmother ... ) test
+
+
+ int momcnt = 0;
+ taufromWZ = 0;
+ for (auto lepanc: LePrdVrt->particles_in() ) {
+ int mompid = lepanc->pdg_id();
+ int momstatus = lepanc->status();
+ if ( m_tesit == 1 ) {
+ ATH_MSG_DEBUG(momcnt << "'th mom with pid= " << mompid << " mom status = " << momstatus);
+ }
+ if ( std::abs(mompid) ==15 ) {
+ TauPrdVrt = lepanc->production_vertex();
+ for (auto taumom: TauPrdVrt->particles_in() ) {
+ int wzpdg = taumom->pdg_id();
+ if (abs(wzpdg) == 24 || abs(wzpdg) == 23 || (wzpdg == mompid && anceWZ > 0)) taufromWZ++;
+ if ( m_tesit == 1 ) ATH_MSG_DEBUG("tau mother =" << wzpdg);
+ }
+ }
+
+ if (abs(mompid)==24 || abs(mompid)==23 || (anceWZ > 0 && mompid==lepid) || (abs(mompid) == 15 && taufromWZ > 0)) momWZ++;
+ } // end of lepton mother test loop
+
+ if ( momWZ > 0 && anceWZ > 0 ) iWL++;
+ else iBL++;
+
+ if ( ( abslepid == 11 && leppt >= m_Pt_e && lepeta <= m_Eta_e ) ||( abslepid == 13 && leppt >= m_Pt_mu && lepeta <= m_Eta_mu ) ) {
+ if ( m_tesit == 1 ) ATH_MSG_DEBUG("Phase space OK");
+ if ( momWZ > 0 && anceWZ > 0 ) {
+ iwl++;
+ if (m_signal != 1) iwls++;
+ if ( taufromWZ > 0 ) taulep++;
+ else {
+ if ( abslepid == 11 ) etronCT++; else muonCT++;
+ }
+ } else {
+ ibl++;
+ if ( m_signal != 1 ) {
+ if ( abslepid == 11 ) etronCT ++;
+ else muonCT ++;
+ }
+ }
+ if ( lepid > 0 ) posilep ++;
+ else negalep ++;
+ }
+
+ if (m_tesit == 1) {
+ ATH_MSG_DEBUG("iWL=" << iWL << " iwl=" << iwl << " iBL=" << iBL << " ibl=" << ibl);
+ if ( iWL == 0 && iBL == 0 ) ATH_MSG_WARNING("Check !!! Unexpected filter LEAKAGE !");
+ }
+ } else { // I've prayed for the upstream generators to give less chaos/duplications
+ nullvertex ++;
+ if ( m_tesit == 1 ) ATH_MSG_DEBUG(" NULL production vertex is met !");
+ if ( m_signal != 1 ) {
+ if ( ( abslepid == 11 && leppt >= m_Pt_e && lepeta <= m_Eta_e ) ||( abslepid == 13 && leppt >= m_Pt_mu && lepeta <= m_Eta_mu ) ) {
+ if ( m_tesit == 1 ) ATH_MSG_DEBUG("Phase space OK, NULL prod-vertex");
+ if ( abslepid == 11 ) etronCT ++;
+ else muonCT ++;
+
+ if ( lepid > 0 ) posilep ++;
+ else negalep ++;
+ } // anyway need for kinematical requirements
+ } // only work for conservation for backgrounds
+ } // End if for muon/electron mom vertex test
+ } // end of all mc_particles
+#else
+ HepMC::GenVertex* LePrdVrt;
+ HepMC::GenVertex* TauPrdVrt;
for ( HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) {
LePrdVrt = 0;
TauPrdVrt = 0;
@@ -346,6 +442,7 @@ StatusCode WZtoLeptonFilter::filterEvent() {
} // End if for muon/electron mom vertex test
} // end of leptonic pid
} // end of all mc_particles
+#endif
leps = (m_signal == 1) ? etronCT + muonCT : etronCT + muonCT + taulep;