From aeb4e15866fab450a163c1cf6ec976bcd23529f0 Mon Sep 17 00:00:00 2001
From: Andrii Verbytskyi <averbyts@cern.ch>
Date: Thu, 23 Jul 2020 20:25:11 +0200
Subject: [PATCH 1/4] Migrate EvgenProdTools to HepMC3
---
.../EvgenProdTools/EvgenProdTools/TestHepMC.h | 1 +
Generators/EvgenProdTools/src/FixHepMC.cxx | 77 ++++++++++++++++---
Generators/EvgenProdTools/src/TestHepMC.cxx | 13 ++--
3 files changed, 71 insertions(+), 20 deletions(-)
diff --git a/Generators/EvgenProdTools/EvgenProdTools/TestHepMC.h b/Generators/EvgenProdTools/EvgenProdTools/TestHepMC.h
index fc661c193f70..ae8306628f92 100644
--- a/Generators/EvgenProdTools/EvgenProdTools/TestHepMC.h
+++ b/Generators/EvgenProdTools/EvgenProdTools/TestHepMC.h
@@ -13,6 +13,7 @@
#include "TFile.h"
#include "TH1.h"
#include "AtlasHepMC/GenEvent.h"
+#include "AtlasHepMC/Relatives.h"
#include<cmath>
#include<fstream>
diff --git a/Generators/EvgenProdTools/src/FixHepMC.cxx b/Generators/EvgenProdTools/src/FixHepMC.cxx
index cee6992cfa65..c16e7ee7a6ac 100644
--- a/Generators/EvgenProdTools/src/FixHepMC.cxx
+++ b/Generators/EvgenProdTools/src/FixHepMC.cxx
@@ -27,6 +27,68 @@ StatusCode FixHepMC::execute() {
for (McEventCollection::const_iterator ievt = events()->begin(); ievt != events()->end(); ++ievt) {
// FIXME: const_cast
HepMC::GenEvent* evt = const_cast<HepMC::GenEvent*>(*ievt);
+#ifdef HEPMC3
+ // Add a unit entry to the event weight vector if it's currently empty
+ if (evt->weights().empty()) {
+ ATH_MSG_DEBUG("Adding a unit weight to empty event weight vector");
+ evt->weights().push_back(1);
+ }
+ // Set a (0,0,0) vertex to be the signal vertex if not already set
+ if (!HepMC::signal_process_vertex(evt)) {
+ const HepMC::FourVector nullpos;
+ for (auto iv: evt->vertices()) {
+ if (iv->position() == nullpos) {
+ ATH_MSG_DEBUG("Setting representative event position vertex");
+ HepMC::set_signal_process_vertex(evt,iv);
+ break;
+ }
+ }
+ }
+ // Event particle content cleaning -- remove "bad" structures
+ std::vector<HepMC::GenParticlePtr> toremove;
+ long seenThisEvent = 0;
+ /// @todo Use nicer particles accessor from TruthUtils / HepMC3 when it exists
+ for (auto ip: evt->particles()) {
+ // Skip this particle if (somehow) its pointer is null
+ if (!ip) continue;
+ m_totalSeen += 1;
+ seenThisEvent += 1;
+ // Flag to declare if a particle should be removed
+ bool bad_particle = false;
+ // Check for loops
+ if ( m_killLoops && isLoop(ip) ) {
+ bad_particle = true;
+ m_loopKilled += 1;
+ ATH_MSG_DEBUG( "Found a looper : " );
+ if ( msgLvl( MSG::DEBUG ) ) HepMC::Print::line(ip);
+ }
+ // Check on PDG ID 0
+ if ( m_killPDG0 && isPID0(ip) ) {
+ bad_particle = true;
+ m_pdg0Killed += 1;
+ ATH_MSG_DEBUG( "Found PDG ID 0 : " );
+ if ( msgLvl( MSG::DEBUG ) )HepMC::Print::line(ip);
+ }
+ // Clean decays
+ if ( m_cleanDecays && isNonTransportableInDecayChain(ip) ) {
+ bad_particle = true;
+ m_decayCleaned += 1;
+ ATH_MSG_DEBUG( "Found a bad particle in a decay chain : " );
+ if ( msgLvl( MSG::DEBUG ) ) HepMC::Print::line(ip);
+ }
+ // Only add to the toremove vector once, even if multiple tests match
+ if (bad_particle) toremove.push_back(ip);
+ }
+ // Escape here if there's nothing more to do, otherwise do the cleaning
+ if (toremove.empty()) continue;
+ ATH_MSG_DEBUG("Cleaning event record of " << toremove.size() << " bad particles");
+ // Properties before cleaning
+ const int num_particles_orig = evt->particles().size();
+ for (auto part: toremove) evt->remove_particle(part);
+ const int num_particles_filt = evt->particles().size();
+ // Write out the change in the number of particles
+ ATH_MSG_INFO("Particles filtered: " << num_particles_orig << " -> " << num_particles_filt);
+ #else
// Add a unit entry to the event weight vector if it's currently empty
if (evt->weights().empty()) {
@@ -120,6 +182,7 @@ StatusCode FixHepMC::execute() {
if (num_nochild_vtxs_filt != num_nochild_vtxs_orig)
ATH_MSG_WARNING("Change in no-parent vertices: " << num_nochild_vtxs_orig << " -> " << num_nochild_vtxs_filt);
+#endif
}
return StatusCode::SUCCESS;
}
@@ -143,7 +206,7 @@ bool FixHepMC::isPID0(const HepMC::GenParticlePtr p) {
// Identify non-transportable stuff _after_ hadronisation
bool FixHepMC::isNonTransportableInDecayChain(const HepMC::GenParticlePtr p) {
- return !MC::isTransportable(p) && MC::fromDecay(p);
+ return !MC::isTransportable(p->pdg_id()) && MC::fromDecay(p);
}
// Identify internal "loop" particles
@@ -151,22 +214,12 @@ bool FixHepMC::isLoop(const HepMC::GenParticlePtr p) {
if (p->production_vertex() == p->end_vertex() && p->end_vertex() != NULL) return true;
if (m_loopByBC && p->production_vertex()) {
/// @todo Use new particle MC::parents(...) tool
-#ifdef HEPMC3
- for (auto itrParent: p->production_vertex()->particles_out()) {
+ for (auto itrParent: *(p->production_vertex())) {
if ( HepMC::barcode(itrParent) > HepMC::barcode(p) ) {
ATH_MSG_VERBOSE("Found a loop (a la Sherpa sample) via barcode.");
return true; // Cannot vectorize, but this is a pretty short loop
} // Check on barcodes
} // Loop over parent particles
-#else
- for (HepMC::GenVertex::particle_iterator itrParent = p->production_vertex()->particles_begin(HepMC::parents);
- itrParent != p->production_vertex()->particles_end(HepMC::parents); ++itrParent) {
- if ( (*itrParent)->barcode() > p->barcode() ) {
- ATH_MSG_VERBOSE("Found a loop (a la Sherpa sample) via barcode.");
- return true; // Cannot vectorize, but this is a pretty short loop
- } // Check on barcodes
- } // Loop over parent particles
-#endif
} // Has a production vertex
return false;
}
diff --git a/Generators/EvgenProdTools/src/TestHepMC.cxx b/Generators/EvgenProdTools/src/TestHepMC.cxx
index c8dd06ff725e..f1ca2e1b0e09 100644
--- a/Generators/EvgenProdTools/src/TestHepMC.cxx
+++ b/Generators/EvgenProdTools/src/TestHepMC.cxx
@@ -541,8 +541,12 @@ StatusCode TestHepMC::execute() {
auto vtx = pitr->end_vertex();
if (vtx) {
double p_energy = 0;
+#ifdef HEPMC3
+ for (auto desc: HepMC::descendant_particles(vtx)) {
+#else
for (auto desc_it = vtx->particles_begin(HepMC::descendants); desc_it != vtx->particles_end(HepMC::descendants); ++desc_it) {
auto desc=(*desc_it);
+#endif
if (std::abs(desc->pdg_id()) == m_pdg) tau_child = 1;
if (desc->status() == 1) p_energy += desc->momentum().e();
}
@@ -572,14 +576,7 @@ StatusCode TestHepMC::execute() {
const HepMC::FourVector decaypos = decayvtx->position();
const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
if (displacement > 1e-6) {
-#ifdef HEPMC3
- for (auto ip:decayvtx->particles_out()) {
-#else
- HepMC::GenVertex::particle_iterator ipbegin = decayvtx->particles_begin(HepMC::children);
- HepMC::GenVertex::particle_iterator ipend = decayvtx->particles_end(HepMC::children);
- for (HepMC::GenVertex::particle_iterator ip_it = ipbegin; ip_it != ipend; ++ip_it) {
- auto ip=(*ip_it);
-#endif
+ for (auto ip: *decayvtx
const HepMC::FourVector pos2 = ip->production_vertex()->position();
const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
if (displacement2 < 1e-6) {
--
GitLab
From 9d7dd8e47ce2f3af04add65ff4988329cf933996 Mon Sep 17 00:00:00 2001
From: Andrii Verbytskyi <averbyts@cern.ch>
Date: Thu, 23 Jul 2020 20:27:05 +0200
Subject: [PATCH 2/4] Migrate EvgenProdTools to HepMC3
---
Generators/AtlasHepMC/AtlasHepMC/Relatives.h | 24 ++++++++++++++++++++
1 file changed, 24 insertions(+)
create mode 100644 Generators/AtlasHepMC/AtlasHepMC/Relatives.h
diff --git a/Generators/AtlasHepMC/AtlasHepMC/Relatives.h b/Generators/AtlasHepMC/AtlasHepMC/Relatives.h
new file mode 100644
index 000000000000..261c4f48690a
--- /dev/null
+++ b/Generators/AtlasHepMC/AtlasHepMC/Relatives.h
@@ -0,0 +1,24 @@
+/* Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ Author: Andrii Verbytskyi andrii.verbytskyi@mpp.mpg.de
+*/
+#ifndef ATLASHEPMC_RELATIVES_H
+#define ATLASHEPMC_RELATIVES_H
+#ifdef HEPMC3
+#include "HepMC3/Relatives.h"
+namespace HepMC {
+typedef HepMC3::Relatives Relatives;
+using HepMC3::children_particles;
+using HepMC3::children_vertices;
+using HepMC3::grandchildren_particles;
+using HepMC3::grandchildren_vertices;
+using HepMC3::parent_particles;
+using HepMC3::parent_vertices;
+using HepMC3::grandparent_particles;
+using HepMC3::grandparent_vertices;
+using HepMC3::descendant_particles;
+using HepMC3::descendant_vertices;
+using HepMC3::ancestor_particles;
+using HepMC3::ancestor_vertices;
+}
+#endif
+#endif
\ No newline at end of file
--
GitLab
From e5d82b9588adee66fbfd8fa9e07154578304d54e Mon Sep 17 00:00:00 2001
From: Andrii Verbytskyi <averbyts@cern.ch>
Date: Fri, 24 Jul 2020 01:09:13 +0200
Subject: [PATCH 3/4] Fix
---
Generators/EvgenProdTools/src/TestHepMC.cxx | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Generators/EvgenProdTools/src/TestHepMC.cxx b/Generators/EvgenProdTools/src/TestHepMC.cxx
index f1ca2e1b0e09..1aa5f8906c44 100644
--- a/Generators/EvgenProdTools/src/TestHepMC.cxx
+++ b/Generators/EvgenProdTools/src/TestHepMC.cxx
@@ -576,7 +576,7 @@ StatusCode TestHepMC::execute() {
const HepMC::FourVector decaypos = decayvtx->position();
const double displacement = decaypos.x()*decaypos.x() + decaypos.y()*decaypos.y() + decaypos.z()*decaypos.z();
if (displacement > 1e-6) {
- for (auto ip: *decayvtx
+ for (auto ip: *decayvtx) {
const HepMC::FourVector pos2 = ip->production_vertex()->position();
const double displacement2 = pos2.x()*pos2.x() + pos2.y()*pos2.y() + pos2.z()*pos2.z();
if (displacement2 < 1e-6) {
--
GitLab
From 6a6fc95a36381f5be010557b19dc7030f9a75494 Mon Sep 17 00:00:00 2001
From: Andrii Verbytskyi <averbyts@cern.ch>
Date: Fri, 24 Jul 2020 01:11:01 +0200
Subject: [PATCH 4/4] Fix
---
Generators/EvgenProdTools/src/FixHepMC.cxx | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/Generators/EvgenProdTools/src/FixHepMC.cxx b/Generators/EvgenProdTools/src/FixHepMC.cxx
index c16e7ee7a6ac..a6e4c42216c0 100644
--- a/Generators/EvgenProdTools/src/FixHepMC.cxx
+++ b/Generators/EvgenProdTools/src/FixHepMC.cxx
@@ -206,7 +206,7 @@ bool FixHepMC::isPID0(const HepMC::GenParticlePtr p) {
// Identify non-transportable stuff _after_ hadronisation
bool FixHepMC::isNonTransportableInDecayChain(const HepMC::GenParticlePtr p) {
- return !MC::isTransportable(p->pdg_id()) && MC::fromDecay(p);
+ return !MC::PID::isTransportable(p->pdg_id()) && MC::fromDecay(p);
}
// Identify internal "loop" particles
--
GitLab