From f8d7b87a1956e0203399dad7e5ab42ceb8fb8c7b Mon Sep 17 00:00:00 2001
From: Andrii Verbytskyi <andrii.verbytskyi@cern.ch>
Date: Tue, 23 Jun 2020 18:51:29 +0000
Subject: [PATCH] Migrate xaodtruthcnv to HepMC3
---
.../xAODTruthCnv/src/HepMCTruthReader.cxx | 31 +++++++-
.../xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.cxx | 71 ++++++++++++++++---
Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.h | 6 ++
3 files changed, 97 insertions(+), 11 deletions(-)
diff --git a/Event/xAOD/xAODTruthCnv/src/HepMCTruthReader.cxx b/Event/xAOD/xAODTruthCnv/src/HepMCTruthReader.cxx
index bb5b253bb65..f0ddf24c69f 100644
--- a/Event/xAOD/xAODTruthCnv/src/HepMCTruthReader.cxx
+++ b/Event/xAOD/xAODTruthCnv/src/HepMCTruthReader.cxx
@@ -74,14 +74,22 @@ StatusCode HepMCTruthReader::execute() {
void HepMCTruthReader::printEvent(const HepMC::GenEvent* event) {
cout << "--------------------------------------------------------------------------------\n";
cout << "GenEvent: #" << "NNN" << "\n";
+#ifdef HEPMC3
+ cout << " Entries this event: " << event->vertices().size() << " vertices, " << event->particles().size() << " particles.\n";
+#else
cout << " Entries this event: " << event->vertices_size() << " vertices, " << event->particles_size() << " particles.\n";
// Particles and vertices
+#endif
cout << " GenParticle Legend\n";
cout << " Barcode PDG ID ( Px, Py, Pz, E ) Stat DecayVtx\n";
cout << "--------------------------------------------------------------------------------\n";
+#ifdef HEPMC3
+ for (auto iv: ((HepMC::GenEvent*)event)->vertices()) { printVertex(iv); }
+#else
for (HepMC::GenEvent::vertex_const_iterator iv = event->vertices_begin(); iv != event->vertices_end(); ++iv) {
printVertex(*iv);
}
+#endif
cout << "--------------------------------------------------------------------------------\n";
}
@@ -130,7 +138,7 @@ void HepMCTruthReader::printVertex(const HepMC::GenVertexPtr vertex) {
// print out gives us a unique tag for the particle.
if (vertex->position().x() != 0.0 && vertex->position().y() != 0.0 && vertex->position().z() != 0.0) {
cout.width(9);
- cout << (void*)vertex;
+ cout << HepMC::raw_pointer(vertex);
cout << " ID:";
cout.width(5);
cout << vertex->id();
@@ -154,7 +162,7 @@ void HepMCTruthReader::printVertex(const HepMC::GenVertexPtr vertex) {
cout << endl;
} else {
cout.width(9);
- cout << (void*)vertex;
+ cout << HepMC::raw_pointer(vertex);
cout << " ID:";
cout.width(5);
cout << vertex->id();
@@ -170,6 +178,24 @@ void HepMCTruthReader::printVertex(const HepMC::GenVertexPtr vertex) {
// cout << endl;
// }
// Print out all the incoming, then outgoing particles
+#ifdef HEPMC3
+ for (auto iPIn: vertex->particles_in()) {
+ if ( iPIn == vertex->particles_in().front() ) {
+ cout << " I: ";
+ cout.width(2);
+ cout << vertex->particles_in().size();
+ } else cout << " ";
+ printParticle(iPIn);
+ }
+ for (auto iPOut: vertex->particles_out()) {
+ if ( iPOut == vertex->particles_out().front()) {
+ cout << " O: ";
+ cout.width(2);
+ cout << vertex->particles_out().size();
+ } else cout << " ";
+ printParticle(iPOut);
+ }
+#else
for (HepMC::GenVertex::particles_in_const_iterator iPIn = vertex->particles_in_const_begin();
iPIn != vertex->particles_in_const_end(); ++iPIn) {
if ( iPIn == vertex->particles_in_const_begin() ) {
@@ -189,6 +215,7 @@ void HepMCTruthReader::printVertex(const HepMC::GenVertexPtr vertex) {
printParticle(*iPOut);
}
+#endif
cout.flags(f);
}
diff --git a/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.cxx b/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.cxx
index 6a3698116eb..a651259d159 100644
--- a/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.cxx
+++ b/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.cxx
@@ -186,8 +186,13 @@ namespace xAODMaker {
xTruthEventContainer->push_back( xTruthEvent );
// Cross-section
auto crossSection = genEvt->cross_section();
+#ifdef HEPMC3
+ xTruthEvent->setCrossSection(crossSection ? (float)crossSection->xsec() : -1);
+ xTruthEvent->setCrossSectionError(crossSection ? (float)crossSection->xsec_err() : -1);
+#else
xTruthEvent->setCrossSection(crossSection ? (float)crossSection->cross_section() : -1);
xTruthEvent->setCrossSectionError(crossSection ? (float)crossSection->cross_section_error() : -1);
+#endif
if (m_writeMetaData) {
//The mcChannelNumber is used as a unique identifier for which truth meta data belongs to
@@ -211,6 +216,22 @@ namespace xAODMaker {
// Heavy ion info
auto const hiInfo = genEvt->heavy_ion();
if (hiInfo) {
+#ifdef HEPMC3
+ /* Please note HepMC3 as well as more recent HePMC2 versions have more Hi parameters */
+ xTruthEvent->setHeavyIonParameter(hiInfo->Ncoll_hard, xAOD::TruthEvent::NCOLLHARD);
+ xTruthEvent->setHeavyIonParameter(hiInfo->Npart_proj, xAOD::TruthEvent::NPARTPROJ);
+ xTruthEvent->setHeavyIonParameter(hiInfo->Npart_targ, xAOD::TruthEvent::NPARTTARG);
+ xTruthEvent->setHeavyIonParameter(hiInfo->Ncoll, xAOD::TruthEvent::NCOLL);
+ xTruthEvent->setHeavyIonParameter(hiInfo->spectator_neutrons, xAOD::TruthEvent::SPECTATORNEUTRONS);
+ xTruthEvent->setHeavyIonParameter(hiInfo->spectator_protons, xAOD::TruthEvent::SPECTATORPROTONS);
+ xTruthEvent->setHeavyIonParameter(hiInfo->N_Nwounded_collisions, xAOD::TruthEvent::NNWOUNDEDCOLLISIONS);
+ xTruthEvent->setHeavyIonParameter(hiInfo->Nwounded_N_collisions, xAOD::TruthEvent::NWOUNDEDNCOLLISIONS);
+ xTruthEvent->setHeavyIonParameter(hiInfo->Nwounded_Nwounded_collisions, xAOD::TruthEvent::NWOUNDEDNWOUNDEDCOLLISIONS);
+ xTruthEvent->setHeavyIonParameter((float)hiInfo->impact_parameter, xAOD::TruthEvent::IMPACTPARAMETER);
+ xTruthEvent->setHeavyIonParameter((float)hiInfo->event_plane_angle, xAOD::TruthEvent::EVENTPLANEANGLE);
+ xTruthEvent->setHeavyIonParameter((float)hiInfo->eccentricity, xAOD::TruthEvent::ECCENTRICITY);
+ xTruthEvent->setHeavyIonParameter((float)hiInfo->sigma_inel_NN, xAOD::TruthEvent::SIGMAINELNN);
+#else
xTruthEvent->setHeavyIonParameter(hiInfo->Ncoll_hard(), xAOD::TruthEvent::NCOLLHARD);
xTruthEvent->setHeavyIonParameter(hiInfo->Npart_proj(), xAOD::TruthEvent::NPARTPROJ);
xTruthEvent->setHeavyIonParameter(hiInfo->Npart_targ(), xAOD::TruthEvent::NPARTTARG);
@@ -224,6 +245,7 @@ namespace xAODMaker {
xTruthEvent->setHeavyIonParameter(hiInfo->event_plane_angle(), xAOD::TruthEvent::EVENTPLANEANGLE);
xTruthEvent->setHeavyIonParameter(hiInfo->eccentricity(), xAOD::TruthEvent::ECCENTRICITY);
xTruthEvent->setHeavyIonParameter(hiInfo->sigma_inel_NN(), xAOD::TruthEvent::SIGMAINELNN);
+#endif
// This doesn't yet exist in our version of HepMC
// xTruthEvent->setHeavyIonParameter(hiInfo->centrality(),xAOD::TruthEvent::CENTRALITY);
}
@@ -232,6 +254,18 @@ namespace xAODMaker {
// This will exist 99% of the time, except for e.g. cosmic or particle gun simulation
auto const pdfInfo = genEvt->pdf_info();
if (pdfInfo) {
+#ifdef HEPMC3
+ xTruthEvent->setPdfInfoParameter(pdfInfo->parton_id[0], xAOD::TruthEvent::PDGID1);
+ xTruthEvent->setPdfInfoParameter(pdfInfo->parton_id[1], xAOD::TruthEvent::PDGID2);
+ xTruthEvent->setPdfInfoParameter(pdfInfo->pdf_id[1], xAOD::TruthEvent::PDFID1);
+ xTruthEvent->setPdfInfoParameter(pdfInfo->pdf_id[1], xAOD::TruthEvent::PDFID2);
+
+ xTruthEvent->setPdfInfoParameter((float)pdfInfo->x[0], xAOD::TruthEvent::X1);
+ xTruthEvent->setPdfInfoParameter((float)pdfInfo->x[1], xAOD::TruthEvent::X2);
+ xTruthEvent->setPdfInfoParameter((float)pdfInfo->scale, xAOD::TruthEvent::Q);
+ xTruthEvent->setPdfInfoParameter((float)pdfInfo->xf[0], xAOD::TruthEvent::XF1);
+ xTruthEvent->setPdfInfoParameter((float)pdfInfo->xf[1], xAOD::TruthEvent::XF2);
+#else
xTruthEvent->setPdfInfoParameter(pdfInfo->id1(), xAOD::TruthEvent::PDGID1);
xTruthEvent->setPdfInfoParameter(pdfInfo->id2(), xAOD::TruthEvent::PDGID2);
xTruthEvent->setPdfInfoParameter(pdfInfo->pdf_id1(), xAOD::TruthEvent::PDFID1);
@@ -242,6 +276,7 @@ namespace xAODMaker {
xTruthEvent->setPdfInfoParameter((float)pdfInfo->scalePDF(), xAOD::TruthEvent::Q);
xTruthEvent->setPdfInfoParameter((float)pdfInfo->pdf1(), xAOD::TruthEvent::XF1);
xTruthEvent->setPdfInfoParameter((float)pdfInfo->pdf2(), xAOD::TruthEvent::XF2);
+#endif
}
}
if (!isSignalProcess) xTruthPileupEventContainer->push_back( xTruthPileupEvent );
@@ -258,8 +293,13 @@ namespace xAODMaker {
// If this is a disconnected vertex, add it manually or won't be added from the loop over particles below.
auto disconnectedSignalProcessVtx = HepMC::signal_process_vertex((HepMC::GenEvent*)genEvt); // Get the signal process vertex
if (disconnectedSignalProcessVtx) {
+#ifdef HEPMC3
+ if (disconnectedSignalProcessVtx->particles_in().size() == 0 &&
+ disconnectedSignalProcessVtx->particles_out().size() == 0 ) {
+#else
if (disconnectedSignalProcessVtx->particles_in_size() == 0 &&
disconnectedSignalProcessVtx->particles_out_size() == 0 ) {
+#endif
//This is a disconnected vertex, add it manually
vertices.push_back (disconnectedSignalProcessVtx);
}
@@ -270,31 +310,37 @@ namespace xAODMaker {
// Get the beam particles
pair<HepMC::GenParticlePtr,HepMC::GenParticlePtr> beamParticles;
bool genEvt_valid_beam_particles=false;
+#ifdef HEPMC3
+ auto beamParticles_vec = ((HepMC::GenEvent*)genEvt)->beams();
+ genEvt_valid_beam_particles=(beamParticles_vec.size()>1);
+ if (genEvt_valid_beam_particles){beamParticles.first=beamParticles_vec[0]; beamParticles.second=beamParticles_vec[1]; }
+#else
genEvt_valid_beam_particles=genEvt->valid_beam_particles();
if ( genEvt_valid_beam_particles ) beamParticles = genEvt->beam_particles();
- for (HepMC::GenEvent::particle_const_iterator pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr) {
+#endif
+ for (auto part: *((HepMC::GenEvent*)genEvt)) {
// (a) create TruthParticle
xAOD::TruthParticle* xTruthParticle = new xAOD::TruthParticle();
xTruthParticleContainer->push_back( xTruthParticle );
- fillParticle(xTruthParticle, *pitr); // (b) Copy HepMC info into the new particle
+ fillParticle(xTruthParticle, part); // (b) Copy HepMC info into the new particle
// (c) Put particle into container; Build Event<->Particle element link
const ElementLink<xAOD::TruthParticleContainer> eltp(*xTruthParticleContainer, xTruthParticleContainer->size()-1);
if (isSignalProcess) xTruthEvent->addTruthParticleLink(eltp);
if (!isSignalProcess) xTruthPileupEvent->addTruthParticleLink(eltp);
// Create link between HepMC and xAOD truth
- if (isSignalProcess) truthLinkVec->push_back(new xAODTruthParticleLink(HepMcParticleLink((*pitr),0,EBC_MAINEVCOLL,HepMcParticleLink::IS_POSITION), eltp));
- if (!isSignalProcess) truthLinkVec->push_back(new xAODTruthParticleLink(HepMcParticleLink((*pitr),genEvt->event_number()), eltp));
+ if (isSignalProcess) truthLinkVec->push_back(new xAODTruthParticleLink(HepMcParticleLink(part,0,EBC_MAINEVCOLL,HepMcParticleLink::IS_POSITION), eltp));
+ if (!isSignalProcess) truthLinkVec->push_back(new xAODTruthParticleLink(HepMcParticleLink(part,genEvt->event_number()), eltp));
// Is this one of the beam particles?
if (genEvt_valid_beam_particles) {
if (isSignalProcess) {
- if (*pitr == beamParticles.first) xTruthEvent->setBeamParticle1Link(eltp);
- if (*pitr == beamParticles.second) xTruthEvent->setBeamParticle2Link(eltp);
+ if (part == beamParticles.first) xTruthEvent->setBeamParticle1Link(eltp);
+ if (part == beamParticles.second) xTruthEvent->setBeamParticle2Link(eltp);
}
}
// (d) Particle's production vertex
- auto productionVertex = (*pitr)->production_vertex();
+ auto productionVertex = part->production_vertex();
if (productionVertex) {
VertexParticles& parts = vertexMap[productionVertex];
if (parts.incoming.empty() && parts.outgoing.empty())
@@ -306,7 +352,7 @@ namespace xAODMaker {
// else maybe want to keep track that this is the production vertex
//
// (e) Particle's decay vertex
- auto decayVertex = (*pitr)->end_vertex();
+ auto decayVertex = part->end_vertex();
if (decayVertex) {
VertexParticles& parts = vertexMap[decayVertex];
if (parts.incoming.empty() && parts.outgoing.empty())
@@ -378,7 +424,7 @@ namespace xAODMaker {
tp->setBarcode(HepMC::barcode(gp));
tp->setStatus(gp->status());
- const HepMC::Polarization& pol = gp->polarization();
+ auto pol = HepMC::polarization(gp);
if (pol.is_defined()) {
tp->setPolarizationParameter(pol.theta(), xAOD::TruthParticle::polarizationTheta);
tp->setPolarizationParameter(pol.phi(), xAOD::TruthParticle::polarizationPhi);
@@ -424,6 +470,12 @@ namespace xAODMaker {
m_tmd->push_back (std::make_unique <xAOD::TruthMetaData>());
xAOD::TruthMetaData* md = m_tmd->back();
+#ifdef HEPMC3
+ ///Here comes the fix. Note that HepMC2.06.11 also contains the fix
+ md->setMcChannelNumber(mcChannelNumber);
+ std::vector<std::string> orderedWeightNameVec=genEvt.weight_names();
+ md->setWeightNames(orderedWeightNameVec);
+#else
// FIXME: class member protection violation here.
// This appears to be because WeightContainer has no public methods
// to get information about the weight names.
@@ -441,6 +493,7 @@ namespace xAODMaker {
md->setMcChannelNumber(mcChannelNumber);
md->setWeightNames( std::move(orderedWeightNameVec) );
+#endif
}
return StatusCode::SUCCESS;
diff --git a/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.h b/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.h
index 79513e702cc..7dbd7ab4166 100644
--- a/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.h
+++ b/Event/xAOD/xAODTruthCnv/src/xAODTruthCnvAlg.h
@@ -7,6 +7,9 @@
#include "AthenaBaseComps/AthReentrantAlgorithm.h"
+#ifdef HEPMC3
+// This form of ifdef is kept for convenience
+#else
// The lines below I don't like. We should fix them when we update the
// the metadata to handles (ATLASRECTS-4162).
// Needs changes in HepMC to resolve.
@@ -20,6 +23,7 @@
#ifdef __clang__
#pragma clang diagnostic pop
#endif
+#endif
#include "GeneratorObjects/xAODTruthParticleLink.h"
#include "GeneratorObjects/McEventCollection.h"
@@ -39,6 +43,8 @@
#include <unordered_set>
+#include "AtlasHepMC/WeightContainer.h"
+#include "AtlasHepMC/Polarization.h"
#include "AtlasHepMC/GenVertex_fwd.h"
#include "AtlasHepMC/GenParticle_fwd.h"
--
GitLab