diff --git a/Generators/GeneratorFilters/GeneratorFilters/VBFMjjIntervalFilter.h b/Generators/GeneratorFilters/GeneratorFilters/VBFMjjIntervalFilter.h
index 04b4519f27d5e1ba942e445cfd512397f27b377b..9318b38e07d88a439e8f67e49c123e509f5c8b19 100644
--- a/Generators/GeneratorFilters/GeneratorFilters/VBFMjjIntervalFilter.h
+++ b/Generators/GeneratorFilters/GeneratorFilters/VBFMjjIntervalFilter.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#ifndef GENERATORFILTERSVBFMJJINTERVALFILTER_H
@@ -48,9 +48,9 @@ private:
bool m_taujetoverlap;
double m_alpha;
- bool checkOverlap(double, double, std::vector<HepMC::GenParticlePtr>);
- bool checkOverlap(double, double, std::vector<CLHEP::HepLorentzVector*>);
- CLHEP::HepLorentzVector sumDaughterNeutrinos( HepMC::GenParticlePtr );
+ bool checkOverlap(double, double, std::vector<HepMC::ConstGenParticlePtr>);
+ bool checkOverlap(double, double, std::vector<CLHEP::HepLorentzVector>);
+ CLHEP::HepLorentzVector sumDaughterNeutrinos( HepMC::ConstGenParticlePtr );
public:
diff --git a/Generators/GeneratorFilters/GeneratorFilters/VHtoVVDiLepFilter.h b/Generators/GeneratorFilters/GeneratorFilters/VHtoVVDiLepFilter.h
index 2886150e441acc47fc2462a4682ce8b0021acc05..3397add11221b5a1493381de6223651a84c3dd8e 100644
--- a/Generators/GeneratorFilters/GeneratorFilters/VHtoVVDiLepFilter.h
+++ b/Generators/GeneratorFilters/GeneratorFilters/VHtoVVDiLepFilter.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#ifndef GENERATORFILTERS_VHTOVVFILTERDILEP_H
@@ -36,7 +36,7 @@ private:
int m_nVHtoVV;
int m_nGoodVHtoVV;
- void findAncestor(const HepMC::GenVertexPtr searchvertex,
+ void findAncestor(HepMC::ConstGenVertexPtr searchvertex,
int targetPDGID, int& n_okPDGChild);
};
diff --git a/Generators/GeneratorFilters/src/MissingEtFilter.cxx b/Generators/GeneratorFilters/src/MissingEtFilter.cxx
index fa478d14add76066edb0514feafa43e4211cc53f..0415f0b50be2e81b2aeff116f43c934542711080 100644
--- a/Generators/GeneratorFilters/src/MissingEtFilter.cxx
+++ b/Generators/GeneratorFilters/src/MissingEtFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/MissingEtFilter.h"
@@ -24,7 +24,7 @@ StatusCode MissingEtFilter::filterEvent() {
for (auto pitr: *genEvt) {
if (!MC::isGenStable(pitr)) continue;
if (!MC::isNonInteracting(pitr)) continue;
- if(!m_useHadronicNu && MC::isNeutrino(pitr) && !(fromWZ(pitr) || fromTau(pitr)) ) continue;
+ if(!m_useHadronicNu && MC::PID::isNeutrino(pitr->pdg_id()) && !(fromWZ(pitr) || fromTau(pitr)) ) continue;
ATH_MSG_VERBOSE("Found noninteracting particle: ID = " << pitr->pdg_id() << " PX = " << pitr->momentum().px() << " PY = "<< pitr->momentum().py());
sumx += pitr->momentum().px();
sumy += pitr->momentum().py();
@@ -32,7 +32,7 @@ StatusCode MissingEtFilter::filterEvent() {
}
// Now see what the total missing Et is and compare to minimum
- double met = sqrt(sumx*sumx + sumy*sumy);
+ double met = std::sqrt(sumx*sumx + sumy*sumy);
ATH_MSG_DEBUG("Totals for event: EX = " << sumx << ", EY = "<< sumy << ", ET = " << met);
setFilterPassed(met >= m_METmin);
return StatusCode::SUCCESS;
@@ -66,7 +66,7 @@ bool MissingEtFilter::fromWZ( HepMC::ConstGenParticlePtr part ) const
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 ( abs( parent_pdgid ) < 9 ) return true;
+ if ( std::abs( parent_pdgid ) < 9 ) return true;
if ( parent_pdgid == part->pdg_id() ) return fromWZ( *iter );
}
#endif
@@ -89,15 +89,15 @@ bool MissingEtFilter::fromTau( HepMC::ConstGenParticlePtr part ) const
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 (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();
- if ( abs( parent_pdgid ) == 15 && fromWZ(*iter)) return true;
- if (MC::PID::isHadron( parent_pdgid ) || abs( parent_pdgid ) < 9 ) return false;
+ if ( std::abs( parent_pdgid ) == 15 && fromWZ(*iter)) return true;
+ if (MC::PID::isHadron( parent_pdgid ) || std::abs( parent_pdgid ) < 9 ) return false;
if ( parent_pdgid == part->pdg_id() ) return fromTau( *iter );
}
#endif
diff --git a/Generators/GeneratorFilters/src/TTbarMassFilter.cxx b/Generators/GeneratorFilters/src/TTbarMassFilter.cxx
index 64f2a727c82c89245ba2ca761ac2bfd69ec9e1e1..2ab630d8d255c1f8768d033a419c37ee97f75007 100644
--- a/Generators/GeneratorFilters/src/TTbarMassFilter.cxx
+++ b/Generators/GeneratorFilters/src/TTbarMassFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/TTbarMassFilter.h"
@@ -20,11 +20,113 @@ StatusCode TTbarMassFilter::filterEvent() {
bool isLastTop = false;
bool isFirstTop = false;
double topPairInvariantMass = 0.;
- std::vector<HepMC::GenParticlePtr> tops;
- std::vector<HepMC::GenVertexPtr> top_vtxs;
+ std::vector<HepMC::ConstGenParticlePtr> tops;
+ std::vector<HepMC::ConstGenVertexPtr> top_vtxs;
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent* genEvt = (*itr);
+#ifdef HEPMC3
+ for ( auto mcpart: *genEvt) {
+ // Reset the flag. it become true if the top with a final ('last') status is found
+ isLastTop = false;
+ // Reset the flag. becomes true if the top with the initial ('first') status is found
+ isFirstTop = false;
+
+ // Work only with tops
+ if (std::abs(mcpart->pdg_id()) == 6) {
+ // Assume that this is the 'last' top
+ isLastTop=true;
+ auto decayVtx = mcpart->end_vertex();
+ // Unusual case...
+ if (!decayVtx) {
+ ATH_MSG_WARNING("top particle with a status "<<mcpart->status()<<" has no valid decay vertex. ");
+ ATH_MSG_WARNING("It looks like a Pythia history particle if it has a status=3. Skip this particle ");
+ continue;
+ }
+
+ // Find out whether this is the top particle with final statuscode, so, just before its decay.
+ /// @todo How generator-portable is this status code assumption?
+ for (auto child_mcpart: decayVtx->particles_out() ) {
+ if (std::abs(child_mcpart->pdg_id()) == 6) {
+ // This is not a 'last' top: break the loop over the children, and do nothing with this top particle
+ isLastTop = false;
+ break;
+ }
+ }
+
+ // Store the 'last' top
+ if (isLastTop) {
+ ATH_MSG_DEBUG("Top particle with a status " << mcpart->status() << " is found and stored");
+ tops.push_back(mcpart);
+ }
+
+ // Find and store production vertices for the 'last' tops. This will help to easily couple top and anti-top
+ // pairs later if there are 2 pairs of them. It's a very rare case but needs to be properly handled.
+ if (isLastTop) {
+ // Investigate whether this top particle is the 'first' one.
+ isFirstTop = false;
+ // Retrieve the production vertex of the current 'last' top particle
+ auto prodVtx = mcpart->production_vertex();
+ if (!prodVtx) {
+ ATH_MSG_WARNING("Top particle with a status " << mcpart->status() << " has no valid production vertex");
+ //save null pointer for consistency
+ top_vtxs.push_back(nullptr);
+ } else {
+ // Loop until the 'first' top particle production vertex is not reached
+ while (!isFirstTop && prodVtx) {
+ // Loop over the top mother particles
+ for (auto mother_mcpart: prodVtx->particles_in()) {
+ // One of the mother particles is still top quark. needed to go up in the hierarchy
+ if (mother_mcpart->pdg_id() == 6) {
+ isFirstTop = false;
+ prodVtx = mother_mcpart->production_vertex();
+ if (!prodVtx) {
+ ATH_MSG_WARNING("mother particle is still a top with a status " << mcpart->status() << ", but has no valid production vertex");
+ top_vtxs.push_back(nullptr);
+ }
+ break;
+ } else {
+ // The current (in loop over the mother particles) mother particle is not a top.
+ // maybe the found top particle is the initial one? set the flag true. it will be
+ // set false if one of the other mother particles is top quark.
+ isFirstTop = true;
+ }
+ }
+ }
+
+ //Store the production vertex for the given top particle. it's a production vertex
+ // of the initial ('first') top particle for the current final ('last') top.
+ if (isFirstTop) {
+ top_vtxs.push_back(prodVtx);
+ }
+ }
+ }
+
+ // Investigate how many top and anti-tops are found and make a decision
+ // to continue or stop looping over GenParticles in the event
+ if (isLastTop) {
+ // One more 'last' top particle was stored into the vector --> increment the number of top or anti-top particles.
+ if (tops[tops.size()-1]->pdg_id() == 6) top++; else topbar++;
+
+ // One top and one ati-top 'last' particles are found. Stop looping
+ // over the particles
+ if (top==1 && topbar==1) break;
+
+ // This should not happen. But depends on a style how a given
+ // generator records the particles. In very rare cases there are more
+ // than 1 top-pairs created. If occasionally the same signed tops from
+ // the different couples are ordered first in the generator record,
+ // then the situation below is possible
+ if (top > 1 || topbar > 1) ATH_MSG_INFO("More than one top-pair exist in the event.");
+
+ // Both top-pairs are found. don't expect more than that. So, it's
+ // time to break the loop over particles. this fairly assumes that
+ // there are no more top-pairs
+ if (top==2 && topbar==2) break;
+ }
+ }
+ }
+#else
HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin();
for (; pitr!=genEvt->particles_end(); ++pitr) {
// Reset the flag. it become true if the top with a final ('last') status is found
@@ -34,7 +136,7 @@ StatusCode TTbarMassFilter::filterEvent() {
HepMC::GenParticle* mcpart = (*pitr);
// Work only with tops
- if (abs(mcpart->pdg_id()) == 6) {
+ if (std::abs(mcpart->pdg_id()) == 6) {
// Assume that this is the 'last' top
isLastTop=true;
HepMC::GenVertex* decayVtx = mcpart->end_vertex();
@@ -135,6 +237,7 @@ StatusCode TTbarMassFilter::filterEvent() {
}
}
+#endif
}
// Main case: there is one top-pair in the event
@@ -158,7 +261,7 @@ StatusCode TTbarMassFilter::filterEvent() {
double topPairInvMass_2=0.;
// There must be only two different vertecies out of the four
- const HepMC::GenVertexPtr prodVtx = top_vtxs[0];
+ auto prodVtx = top_vtxs[0];
int top_11 = 0;
int top_12 = -1;
int top_21 = -1;
diff --git a/Generators/GeneratorFilters/src/TauFilter.cxx b/Generators/GeneratorFilters/src/TauFilter.cxx
index 66a77781272b8f5cbc453b9a47d28f612cd9f710..253607e6d5709464185fd1e979aebfa9f6a7a118 100644
--- a/Generators/GeneratorFilters/src/TauFilter.cxx
+++ b/Generators/GeneratorFilters/src/TauFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/TauFilter.h"
@@ -58,35 +58,30 @@ CLHEP::HepLorentzVector TauFilter::sumDaughterNeutrinos( HepMC::ConstGenParticle
StatusCode TauFilter::filterEvent() {
- HepMC::GenParticlePtr tau;
CLHEP::HepLorentzVector mom_tauprod; // will contain the momentum of the products of the tau decay
CLHEP::HepLorentzVector tauvis;
CLHEP::HepLorentzVector nutau;
- tau = 0;
int ntau = 0;
McEventCollection::const_iterator itr;
for (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) {
+ for ( auto tau: *genEvt) {
// Look for the first tau with genstat != 3
- if (abs((*pitr)->pdg_id()) == 15 && (*pitr)->status() != 3) {
- tau = (*pitr);
- ATH_MSG_DEBUG("found tau with barcode " << tau->barcode() << " status " << (*pitr)->status());
+ if (std::abs(tau->pdg_id()) != 15 || tau->status() == 3) continue;
+ ATH_MSG_DEBUG("found tau with barcode " << HepMC::barcode(tau) << " status " << tau->status());
ATH_MSG_DEBUG("pT\t\teta\tphi\tid");
ATH_MSG_DEBUG(tau->momentum().perp() << "\t" <<
tau->momentum().eta() << "\t" <<
tau->momentum().phi() << "\t" <<
tau->pdg_id() << "\t");
-
- HepMC::GenVertex::particles_out_const_iterator begin = tau->end_vertex()->particles_out_const_begin();
- HepMC::GenVertex::particles_out_const_iterator end = tau->end_vertex()->particles_out_const_end();
int leptonic = 0;
- for ( ; begin != end; begin++ ) {
- if ( (*begin)->production_vertex() != tau->end_vertex() ) continue;
- if ( abs( (*begin)->pdg_id() ) == 12 ) leptonic = 1;
- if ( abs( (*begin)->pdg_id() ) == 14 ) leptonic = 2;
- if ( abs( (*begin)->pdg_id() ) == 15 ) leptonic = 11;
+ if ( tau->production_vertex()) continue;
+ for (auto beg: *(tau->end_vertex())) {
+ if ( beg->production_vertex() != tau->end_vertex() ) continue;
+ if ( std::abs( beg->pdg_id() ) == 12 ) leptonic = 1;
+ if ( std::abs( beg->pdg_id() ) == 14 ) leptonic = 2;
+ if ( std::abs( beg->pdg_id() ) == 15 ) leptonic = 11;
}
if (leptonic == 11) {
@@ -124,7 +119,6 @@ StatusCode TauFilter::filterEvent() {
ntau++;
m_eventshadacc++;
}
- }
}
}
diff --git a/Generators/GeneratorFilters/src/TopCKMFilter.cxx b/Generators/GeneratorFilters/src/TopCKMFilter.cxx
index 8b749b18453a2e306cd80cfdb94d4259c27eb384..6363fc2b8b3849bb2cc5d55d9f877a89b142818c 100644
--- a/Generators/GeneratorFilters/src/TopCKMFilter.cxx
+++ b/Generators/GeneratorFilters/src/TopCKMFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/TopCKMFilter.h"
@@ -29,6 +29,37 @@ StatusCode TopCKMFilter::filterEvent() {
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent * genEvt = (*itr);
+#ifdef HEPMC3
+ // Get the number of parents on the event (ttbar --> 2 tops!)
+ std::vector<HepMC::ConstGenParticlePtr> Parents;
+ for (auto partItr: *genEvt) if (std::abs(partItr->pdg_id()) == 6) Parents.push_back(partItr);
+ if (Parents.size() != 2) continue;
+ // Check that the first and second parents have the corresponding childs
+ if (!Parents.at(0) || !Parents.at(1)) continue;
+ int Wbosons[2]={0,0};
+ int quarks[2]={0,0};
+ bool isTau[2] ={false,false};
+ for (int i=0;i<2;i++)
+ {
+ for (auto Child: *(Parents.at(i)->end_vertex()))
+ {
+ if ((std::abs(Child->pdg_id()) == m_PDGChild[0] || std::abs(Child->pdg_id()) == m_PDGChild[1])
+ && (Child->momentum().perp() > m_PtMinChild)
+ && (std::abs(Child->momentum().eta()) < m_EtaRangeChild))
+ quarks[i]=Child->pdg_id();
+ if (std::abs(Child->pdg_id()) == 24)
+ {
+ 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;
+ }
+ }
+ }
+ if ( quarks[0]!=0 && quarks[1]!=0 && Wbosons[0]!=0 && Wbosons[1]!=0 && !isTau[0] && !isTau[1] && (std::abs(quarks[0]) != std::abs(quarks[1])) )
+ {
+ return StatusCode::SUCCESS;
+ }
+#else
HepMC::GenEvent::particle_const_iterator part1 = genEvt->particles_begin();
HepMC::GenEvent::particle_const_iterator partItr = part1;
HepMC::GenEvent::particle_const_iterator partE = genEvt->particles_end();
@@ -77,17 +108,17 @@ StatusCode TopCKMFilter::filterEvent() {
isOK1 = 1;
quark1 = (*firstChildItr);
}
- if (abs((*firstChildItr)->pdg_id()) == 24) wBoson1 = (*firstChildItr);
+ if (std::abs((*firstChildItr)->pdg_id()) == 24) wBoson1 = (*firstChildItr);
}
for (secondChildItr = secondChild1; secondChildItr != secondChildE; ++secondChildItr) {
- if ((abs((*secondChildItr)->pdg_id()) == m_PDGChild[0] || abs((*secondChildItr)->pdg_id()) == m_PDGChild[1])
+ if ((std::abs((*secondChildItr)->pdg_id()) == m_PDGChild[0] || std::abs((*secondChildItr)->pdg_id()) == m_PDGChild[1])
&& (*secondChildItr)->momentum().perp() > m_PtMinChild
&& std::abs((*secondChildItr)->momentum().eta()) < m_EtaRangeChild){
isOK2 = 1;
quark2 = (*secondChildItr);
}
- if (abs((*secondChildItr)->pdg_id()) == 24) wBoson2 = (*secondChildItr);
+ if (std::abs((*secondChildItr)->pdg_id()) == 24) wBoson2 = (*secondChildItr);
}
// Now avoid the decay W --> tau nu
@@ -105,16 +136,17 @@ StatusCode TopCKMFilter::filterEvent() {
w2SonE = wVertex2->particles_end(HepMC::children);
for (w1SonItr = w1Son1; w1SonItr != w1SonE; ++w1SonItr) {
- if (abs((*w1SonItr)->pdg_id()) == 15) isTau1 = 1;
+ if (std::abs((*w1SonItr)->pdg_id()) == 15) isTau1 = 1;
}
for (w2SonItr = w2Son1; w2SonItr != w2SonE; ++w2SonItr) {
- if (abs((*w2SonItr)->pdg_id()) == 15) isTau2 = 1;
+ if (std::abs((*w2SonItr)->pdg_id()) == 15) isTau2 = 1;
}
- if (isOK1 && isOK2 && !isTau1 && !isTau2 && abs(quark1->pdg_id()) != abs(quark2->pdg_id())) {
+ if (isOK1 && isOK2 && !isTau1 && !isTau2 && std::abs(quark1->pdg_id()) != std::abs(quark2->pdg_id())) {
return StatusCode::SUCCESS;
}
}
}
+#endif
}
setFilterPassed(false);
return StatusCode::SUCCESS;
diff --git a/Generators/GeneratorFilters/src/TopCKMKinFilter.cxx b/Generators/GeneratorFilters/src/TopCKMKinFilter.cxx
index f91f010f2755d939b685a615e4975779dc9c629b..cd534e9948145ee16e939bc015658c56f29d97f1 100644
--- a/Generators/GeneratorFilters/src/TopCKMKinFilter.cxx
+++ b/Generators/GeneratorFilters/src/TopCKMKinFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/TopCKMKinFilter.h"
@@ -32,6 +32,44 @@ StatusCode TopCKMKinFilter::filterInitialize() {
StatusCode TopCKMKinFilter::filterEvent() {
for (McEventCollection::const_iterator itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent * genEvt = (*itr);
+#ifdef HEPMC3
+ // Get the number of parents on the event (ttbar --> 2 tops!)
+ std::vector<HepMC::ConstGenParticlePtr> Parents;
+ for (auto partItr: *genEvt) if (std::abs(partItr->pdg_id()) == 6) Parents.push_back(partItr);
+ if (Parents.size() != 2) continue;
+ // Check that the first and second parents have the corresponding childs
+ if (!Parents.at(0) || !Parents.at(1)) continue;
+ int Wbosons[2]={0,0};
+ int quarks[2]={0,0};
+ bool isTau[2] ={false,false};
+ bool isKinLep[2] ={false,false};
+ for (int i=0;i<2;i++)
+ {
+ for (auto Child: *(Parents.at(i)->end_vertex()))
+ {
+ if ((std::abs(Child->pdg_id()) == m_PDGChild[0] || std::abs(Child->pdg_id()) == m_PDGChild[1])
+ && (Child->momentum().perp() > m_PtMinChild)
+ && (std::abs(Child->momentum().eta()) < m_EtaRangeChild))
+ quarks[i]=Child->pdg_id();
+ if (std::abs(Child->pdg_id()) == 24)
+ {
+ 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;
+ if (std::abs(gChild->pdg_id()) == 11 || std::abs(gChild->pdg_id()) == 13){
+ if (gChild->momentum().perp() > m_PtMinChild && std::abs(gChild->momentum().eta()) < m_EtaRangeChild) isKinLep[i] = true;
+ }
+ }
+ }
+ }
+ if (isKinLep[0] && isKinLep[1] && quarks[0]!=0 && quarks[1]!=0 && Wbosons[0]!=0 && Wbosons[1]!=0 && !isTau[0] && !isTau[1] && (std::abs(quarks[0]) != std::abs(quarks[1]))){
+ return StatusCode::SUCCESS;
+ }
+
+ }
+#else
HepMC::GenEvent::particle_const_iterator part1 = genEvt->particles_begin();
HepMC::GenEvent::particle_const_iterator partItr = part1;
HepMC::GenEvent::particle_const_iterator partE = genEvt->particles_end();
@@ -133,6 +171,7 @@ StatusCode TopCKMKinFilter::filterEvent() {
}
}
+#endif
}
setFilterPassed(false);
diff --git a/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx b/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx
index b920071b026db39776020431b1d2cac829ad331b..71a088e7b246d899e8cecbc64b1bd89eb2593848 100644
--- a/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx
+++ b/Generators/GeneratorFilters/src/TtHtoVVDecayFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/TtHtoVVDecayFilter.h"
@@ -64,10 +64,42 @@ StatusCode TtHtoVVDecayFilter::filterEvent() {
for (McEventCollection::const_iterator 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) {
+ for (size_t iParent=0; iParent < m_PDGParent.size(); iParent++) {
+ if ( std::abs(pitr->pdg_id()) == m_PDGParent[iParent] && pitr->status() == m_StatusParent[iParent]) { // if find a parent here
+ bool isGrandParentOK = false;
+ for (auto thisMother: pitr->production_vertex()->particles_in()) { // Loop over all grandparent of this founded parent
+ ATH_MSG_DEBUG("Parent " << pitr->pdg_id() << " barcode = " << HepMC::barcode(pitr) << " status = " << pitr->status());
+ ATH_MSG_DEBUG("A parent mother " << thisMother->pdg_id()<< " barcode = " << HepMC::barcode(thisMother));
+ for (size_t iGrandParent = 0; iGrandParent < m_PDGGrandParent.size(); iGrandParent++) { // Loop over all Grandparents we want to check
+ if ( thisMother->pdg_id() == m_PDGGrandParent[iGrandParent] ) isGrandParentOK = true;
+ }
+ }
+ ATH_MSG_DEBUG(" Grandparent is OK? " << isGrandParentOK);
+ if (!isGrandParentOK) continue;
+
+ ++nGoodParent; // good parent from the indicated GrandParents
+ const bool okPDGChild = findAncestor(pitr->end_vertex(), m_PDGParent[iParent]);
+ ATH_MSG_DEBUG("Result " << nGoodParent << " " << okPDGChild );
+ if (!(okPDGChild)) continue;
+
+ ++nGoodDecayedParent; // good parent has the indicated child and child2
+ if (nGoodDecayedParent == 1) {
+ if ( pitr->pdg_id() < 0) charge = -1;
+ else if ( pitr->pdg_id() > 0) charge = 1;
+ } else {
+ if (charge == -1 && pitr->pdg_id() > 0) sameCharge = false;
+ if (charge == 1 && pitr->pdg_id() < 0) sameCharge = false;
+ }
+ }
+ } // End loop over all Parents we want to check
+ } // End loop over all particles in the event
+#else
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++) {
- if ( abs((*pitr)->pdg_id()) == m_PDGParent[iParent] && (*pitr)->status() == m_StatusParent[iParent]) { // if find a parent here
+ if ( std::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);
HepMC::GenVertex::particle_iterator thisMother = firstMother;
@@ -100,6 +132,7 @@ StatusCode TtHtoVVDecayFilter::filterEvent() {
} // End loop over all Parents we want to check
} // End loop over all particles in the event
+#endif
}
ATH_MSG_DEBUG("Selected # good parents = " << nGoodParent << "; # good decayed parents = " << nGoodDecayedParent);
diff --git a/Generators/GeneratorFilters/src/VBFMjjIntervalFilter.cxx b/Generators/GeneratorFilters/src/VBFMjjIntervalFilter.cxx
index 0d1ccfb52ff066b61e7a2591b1ebdb72efd67ae5..f7b9c04f81bc60fef946fac09229c0d22d243626 100644
--- a/Generators/GeneratorFilters/src/VBFMjjIntervalFilter.cxx
+++ b/Generators/GeneratorFilters/src/VBFMjjIntervalFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
// Header for this module
@@ -80,52 +80,48 @@ StatusCode VBFMjjIntervalFilter::filterEvent() {
}
// Find overlap objects
- std::vector<HepMC::GenParticlePtr> MCTruthPhotonList;
- std::vector<HepMC::GenParticlePtr> MCTruthElectronList;
- std::vector<CLHEP::HepLorentzVector*> MCTruthTauList;
+ std::vector<HepMC::ConstGenParticlePtr> MCTruthPhotonList;
+ std::vector<HepMC::ConstGenParticlePtr> MCTruthElectronList;
+ std::vector<CLHEP::HepLorentzVector> MCTruthTauList;
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) {
+ for (auto pitr: *genEvt) {
if (m_photonjetoverlap==true) {
// photon - copied from VBFForwardJetsFilter.cxx
- if ( (*pitr)->pdg_id() == 22 && (*pitr)->status() == 1 &&
- (*pitr)->momentum().perp() >= m_olapPt &&
- std::abs((*pitr)->momentum().pseudoRapidity()) <= m_yMax) {
- MCTruthPhotonList.push_back((*pitr));
+ if ( pitr->pdg_id() == 22 && pitr->status() == 1 &&
+ pitr->momentum().perp() >= m_olapPt &&
+ std::abs(pitr->momentum().pseudoRapidity()) <= m_yMax) {
+ MCTruthPhotonList.push_back(pitr);
}
}
if (m_electronjetoverlap==true) {
// electron
- if (abs((*pitr)->pdg_id()) == 11 && (*pitr)->status() == 1 &&
- (*pitr)->momentum().perp() >= m_olapPt &&
- std::abs((*pitr)->momentum().pseudoRapidity()) <= m_yMax) {
- MCTruthElectronList.push_back((*pitr));
+ if (std::abs(pitr->pdg_id()) == 11 && pitr->status() == 1 &&
+ pitr->momentum().perp() >= m_olapPt &&
+ std::abs(pitr->momentum().pseudoRapidity()) <= m_yMax) {
+ MCTruthElectronList.push_back(pitr);
}
}
if (m_taujetoverlap==true) {
// tau - copied from VBFForwardJetsFilter.cxx
- if ( abs((*pitr)->pdg_id()) == 15 && (*pitr)->status() != 3 ) {
- HepMC::GenParticle *tau = (*pitr);
- HepMC::GenVertex::particles_out_const_iterator begin = tau->end_vertex()->particles_out_const_begin();
- HepMC::GenVertex::particles_out_const_iterator end = tau->end_vertex()->particles_out_const_end();
+ if ( std::abs(pitr->pdg_id()) == 15 && pitr->status() != 3 ) {
+ auto tau = pitr;
int leptonic = 0;
- for ( ; begin != end; begin++ ) {
- if ( (*begin)->production_vertex() != tau->end_vertex() ) continue;
- if ( abs( (*begin)->pdg_id() ) == 12 ) leptonic = 1;
- if ( abs( (*begin)->pdg_id() ) == 14 ) leptonic = 2;
- if ( abs( (*begin)->pdg_id() ) == 15 ) leptonic = 11;
+ for ( auto beg: *(tau->end_vertex()) ) {
+ if ( beg->production_vertex() != tau->end_vertex() ) continue;
+ if ( std::abs( beg->pdg_id() ) == 12 ) leptonic = 1;
+ if ( std::abs( beg->pdg_id() ) == 14 ) leptonic = 2;
+ if ( std::abs( beg->pdg_id() ) == 15 ) leptonic = 11;
}
if (leptonic == 0) {
CLHEP::HepLorentzVector nutau = sumDaughterNeutrinos( tau );
- CLHEP::HepLorentzVector *tauvis = new CLHEP::HepLorentzVector(tau->momentum().px()-nutau.px(),
+ CLHEP::HepLorentzVector tauvis = CLHEP::HepLorentzVector(tau->momentum().px()-nutau.px(),
tau->momentum().py()-nutau.py(),
tau->momentum().pz()-nutau.pz(),
tau->momentum().e()-nutau.e());
- if (tauvis->vect().perp() >= m_olapPt && std::abs(tauvis->vect().pseudoRapidity()) <= m_yMax) {
+ if (tauvis.vect().perp() >= m_olapPt && std::abs(tauvis.vect().pseudoRapidity()) <= m_yMax) {
MCTruthTauList.push_back(tauvis);
- } else {
- delete tauvis;
}
}
}
@@ -193,7 +189,7 @@ StatusCode VBFMjjIntervalFilter::filterEvent() {
}
-bool VBFMjjIntervalFilter::checkOverlap(double eta, double phi, std::vector<HepMC::GenParticlePtr> list) {
+bool VBFMjjIntervalFilter::checkOverlap(double eta, double phi, std::vector<HepMC::ConstGenParticlePtr> list) {
for (size_t i = 0; i < list.size(); ++i) {
double pt = list[i]->momentum().perp();
if (pt > m_olapPt) {
@@ -202,7 +198,7 @@ bool VBFMjjIntervalFilter::checkOverlap(double eta, double phi, std::vector<HepM
double deta = eta-list[i]->momentum().pseudoRapidity();
if (dphi > M_PI) { dphi -= 2.*M_PI; }
if (dphi < -M_PI) { dphi += 2.*M_PI; }
- double dr = sqrt(deta*deta+dphi*dphi);
+ double dr = std::sqrt(deta*deta+dphi*dphi);
if (dr < 0.3) return true;
}
}
@@ -211,16 +207,16 @@ bool VBFMjjIntervalFilter::checkOverlap(double eta, double phi, std::vector<HepM
-bool VBFMjjIntervalFilter::checkOverlap(double eta, double phi, std::vector<CLHEP::HepLorentzVector*> list) {
+bool VBFMjjIntervalFilter::checkOverlap(double eta, double phi, std::vector<CLHEP::HepLorentzVector> list) {
for (size_t i = 0; i < list.size(); ++i) {
- double pt = list[i]->vect().perp();
+ double pt = list[i].vect().perp();
if (pt > m_olapPt) {
/// @todo Provide a helper function for this (and similar)
- double dphi = phi-list[i]->phi();
- double deta = eta-list[i]->vect().pseudoRapidity();
+ double dphi = phi-list[i].phi();
+ double deta = eta-list[i].vect().pseudoRapidity();
if (dphi > M_PI) { dphi -= 2.*M_PI; }
if (dphi < -M_PI) { dphi += 2.*M_PI; }
- double dr = sqrt(deta*deta+dphi*dphi);
+ double dr = std::sqrt(deta*deta+dphi*dphi);
if (dr < 0.3) return true;
}
}
@@ -263,10 +259,10 @@ double VBFMjjIntervalFilter::getEventWeight(xAOD::JetContainer *jets) {
}
-CLHEP::HepLorentzVector VBFMjjIntervalFilter::sumDaughterNeutrinos( HepMC::GenParticlePtr part ) {
+CLHEP::HepLorentzVector VBFMjjIntervalFilter::sumDaughterNeutrinos( HepMC::ConstGenParticlePtr part ) {
CLHEP::HepLorentzVector nu( 0, 0, 0, 0);
- if ( ( abs( part->pdg_id() ) == 12 ) || ( abs( part->pdg_id() ) == 14 ) || ( abs( part->pdg_id() ) == 16 ) ) {
+ if ( ( std::abs( part->pdg_id() ) == 12 ) || ( std::abs( part->pdg_id() ) == 14 ) || ( std::abs( part->pdg_id() ) == 16 ) ) {
nu.setPx(part->momentum().px());
nu.setPy(part->momentum().py());
nu.setPz(part->momentum().pz());
@@ -274,10 +270,8 @@ CLHEP::HepLorentzVector VBFMjjIntervalFilter::sumDaughterNeutrinos( HepMC::GenPa
return nu;
}
- if ( part->end_vertex() == 0 ) return nu;
+ if ( !part->end_vertex() ) return nu;
- HepMC::GenVertex::particles_out_const_iterator begin = part->end_vertex()->particles_out_const_begin();
- HepMC::GenVertex::particles_out_const_iterator end = part->end_vertex()->particles_out_const_end();
- for ( ; begin != end; begin++ ) nu += sumDaughterNeutrinos( *begin );
+ for (auto beg: *(part->end_vertex()) ) nu += sumDaughterNeutrinos( beg );
return nu;
}
diff --git a/Generators/GeneratorFilters/src/VHtoVVDiLepFilter.cxx b/Generators/GeneratorFilters/src/VHtoVVDiLepFilter.cxx
index 0d5372a9979dc085f6a3d7fddae441c6c303689c..0ef5945fc21df91a71ee3d57f4882acd0b671a87 100644
--- a/Generators/GeneratorFilters/src/VHtoVVDiLepFilter.cxx
+++ b/Generators/GeneratorFilters/src/VHtoVVDiLepFilter.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
#include "GeneratorFilters/VHtoVVDiLepFilter.h"
@@ -48,10 +48,42 @@ StatusCode VHtoVVDiLepFilter::filterEvent() {
int nVParent = 0;
for (McEventCollection::const_iterator ievt = events()->begin(); ievt != events()->end(); ++ievt) {
// Loop over all particles in the event
+#ifdef HEPMC3
+ for (auto pitr: **ievt) {
+
+ // Loop over particles from the primary interaction that match the PDG Parent
+ if ( (std::abs(pitr->pdg_id()) == m_PDGParent ||std::abs(pitr->pdg_id()) == m_PDGAssoc) && std::abs(pitr->status()) == m_StatusCode) {
+ // Check if originated from the Higgs or is the associated vector boson
+ bool isGrandParentHiggs = false;
+ bool isGrandParentV = false;
+ for (auto thisMother: pitr->production_vertex()->particles_in()) { // loop over chain of grandparents
+ 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() == m_PDGGrandParent ) isGrandParentHiggs = true;
+ else isGrandParentV = true;
+ }
+ ATH_MSG_DEBUG(" Grand Parent is Higgs? " << isGrandParentHiggs );
+ ATH_MSG_DEBUG(" Grand Parent is V? " << isGrandParentV );
+
+ if (!isGrandParentHiggs && !isGrandParentV) continue;
+
+ if (isGrandParentHiggs) {
+ ++nHiggsParent;
+ if (!pitr->end_vertex()) continue;
+ findAncestor(pitr->end_vertex(), m_PDGParent, n_okPDGHVChildren);
+ } // end of higgs grandparent loop
+
+ if (isGrandParentV) {
+ ++nVParent;
+ findAncestor(pitr->end_vertex(), m_PDGAssoc, n_okPDGAssocVChild);
+ } // end of v grandparent loop
+ } // end good parent loop
+ }
+#else
for (HepMC::GenEvent::particle_const_iterator pitr = (*ievt)->particles_begin(); pitr != (*ievt)->particles_end(); ++pitr) {
// Loop over particles from the primary interaction that match the PDG Parent
- if ( (abs((*pitr)->pdg_id()) == m_PDGParent ||abs((*pitr)->pdg_id()) == m_PDGAssoc) && abs((*pitr)->status()) == m_StatusCode) {
+ if ( (std::abs((*pitr)->pdg_id()) == m_PDGParent ||std::abs((*pitr)->pdg_id()) == m_PDGAssoc) && std::abs((*pitr)->status()) == m_StatusCode) {
HepMC::GenVertex::particle_iterator firstMother = (*pitr)->production_vertex()->particles_begin(HepMC::parents);
HepMC::GenVertex::particle_iterator endMother = (*pitr)->production_vertex()->particles_end(HepMC::parents);
HepMC::GenVertex::particle_iterator thisMother = firstMother;
@@ -82,6 +114,7 @@ StatusCode VHtoVVDiLepFilter::filterEvent() {
} // end good parent loop
}
+#endif
}
ATH_MSG_DEBUG("Result " << nHiggsParent << " " << n_okPDGHVChildren );
@@ -99,24 +132,21 @@ StatusCode VHtoVVDiLepFilter::filterEvent() {
return StatusCode::SUCCESS;
}
-void VHtoVVDiLepFilter::findAncestor(const HepMC::GenVertexPtr searchvertex,
+void VHtoVVDiLepFilter::findAncestor(HepMC::ConstGenVertexPtr searchvertex,
int targetPDGID, int& n_okPDGChild) {
std::vector<int> foundCodes;
if (!searchvertex) return;
- 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();
- HepMC::GenVertex::particles_out_const_iterator thisAncestor = firstAncestor;
- for (; thisAncestor != endAncestor; ++thisAncestor) {
- if (abs((*thisAncestor)->pdg_id()) == targetPDGID) { // same particle as parent
- findAncestor((*thisAncestor)->end_vertex(), targetPDGID, n_okPDGChild);
+ for (auto thisAncestor: *searchvertex) {
+ if (std::abs(thisAncestor->pdg_id()) == targetPDGID) { // same particle as parent
+ findAncestor(thisAncestor->end_vertex(), targetPDGID, n_okPDGChild);
} else {
for (size_t i = 0; i < m_PDGChildren.size(); ++i) {
- int testPdgID = (*thisAncestor)->pdg_id();
- if (abs(testPdgID) == m_PDGChildren[i]) {
+ int testPdgID = thisAncestor->pdg_id();
+ if (std::abs(testPdgID) == m_PDGChildren[i]) {
const bool alreadyFound = (std::find(foundCodes.begin(), foundCodes.end(), testPdgID) != foundCodes.end());
if (!alreadyFound) {
n_okPDGChild++;
- foundCodes.push_back((*thisAncestor)->pdg_id()); // add to list of found particles and check to avoid double counting
+ foundCodes.push_back(thisAncestor->pdg_id()); // add to list of found particles and check to avoid double counting
}
else break;
}