HepMC3/src/HepMC3InputConverter.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Io/HepMC3/HepMC3InputConverter.hpp"
0010
0011 #include "Acts/Utilities/ScopedTimer.hpp"
0012
0013 #include <HepMC3/GenEvent.h>
0014 #include <HepMC3/GenParticle.h>
0015 #include <HepMC3/GenVertex.h>
0016
0017 // This is a hack to make HepMC3::Print::listing public
0018 // It's pretty evil but should have no side-effects
0019 #if defined(__clang__)
0020 #pragma clang diagnostic push
0021 #pragma clang diagnostic ignored "-Wkeyword-macro"
0022 #endif
0023 #define private public
0024 #include <HepMC3/Print.h>
0025 #undef private
0026 #if defined(__clang__)
0027 #pragma clang diagnostic pop
0028 #endif
0029
0030 using namespace Acts::UnitLiterals;
0031
0032 namespace ActsExamples {
0033
0034 HepMC3InputConverter::HepMC3InputConverter(const Config& config,
0035                                            Acts::Logging::Level level)
0036     : IAlgorithm("HepMC3InputConverter", level), m_cfg(config) {
0037   m_outputParticles.initialize(m_cfg.outputParticles);
0038   m_outputVertices.initialize(m_cfg.outputVertices);
0039   m_inputEvent.initialize(m_cfg.inputEvent);
0040 }
0041
0042 ProcessCode HepMC3InputConverter::execute(const AlgorithmContext& ctx) const {
0043   ACTS_DEBUG("HepMC3InputConverter::execute");
0044   const auto& genEvent = *m_inputEvent(ctx);
0045   ACTS_DEBUG("Have " << genEvent.vertices().size() << " vertices");
0046   ACTS_DEBUG("Have " << genEvent.particles().size() << " particles");
0047   ACTS_DEBUG("Have " << genEvent.event_number() << " event number");
0048
0049   convertHepMC3ToInternalEdm(ctx, genEvent);
0050
0051   return ProcessCode::SUCCESS;
0052 }
0053
0054 namespace {
0055
0056 constexpr int kBeamParticleStatus = 4;
0057 constexpr int kUndecayedParticleStatus = 1;
0058
0059 Acts::Vector4 convertPosition(const HepMC3::FourVector& vec) {
0060   return Acts::Vector4(vec.x() * 1_mm, vec.y() * 1_mm, vec.z() * 1_mm,
0061                        vec.t() * 1_mm);
0062 }
0063
0064 std::string printListing(const auto& vertices, const auto& particles) {
0065   auto findParticle = [&](SimBarcode particleId) {
0066     if (auto it = std::ranges::find_if(particles,
0067                                        [&](const auto& particle) {
0068                                          return particle.particleId() ==
0069                                                 particleId;
0070                                        });
0071         it != particles.end()) {
0072       return *it;
0073     }
0074     throw std::runtime_error("Particle not found");
0075   };
0076   std::stringstream ss;
0077
0078   for (const auto& vertex : vertices) {
0079     ss << "Vtx:    " << vertex.vertexId() << " at "
0080        << vertex.position4.transpose() << "\n";
0081     for (const auto& [idx, particleId] : Acts::enumerate(vertex.incoming)) {
0082       const auto& particle = findParticle(particleId);
0083       if (idx == 0) {
0084         ss << " I:     ";
0085       } else {
0086         ss << "        ";
0087       }
0088       ss << particle;
0089       ss << "\n";
0090     }
0091
0092     for (const auto& [idx, particleId] : Acts::enumerate(vertex.outgoing)) {
0093       const auto& particle = findParticle(particleId);
0094       if (idx == 0) {
0095         ss << " O:     ";
0096       } else {
0097         ss << "        ";
0098       }
0099       ss << particle;
0100       ss << "\n";
0101     }
0102   }
0103
0104   return ss.str();
0105 };
0106 }  // namespace
0107
0108 void HepMC3InputConverter::handleVertex(const HepMC3::GenVertex& genVertex,
0109                                         SimVertex& vertex,
0110                                         std::vector<SimVertex>& vertices,
0111                                         std::vector<SimParticle>& particles,
0112                                         std::size_t& nSecondaryVertices,
0113                                         std::size_t& nParticles,
0114                                         std::vector<bool>& seenVertices) const {
0115   for (const auto& particle : genVertex.particles_out()) {
0116     if (particle->end_vertex() != nullptr) {
0117       // This particle has an end vertex, we need to handle that vertex
0118       const auto& endVertex = *particle->end_vertex();
0119       if (seenVertices.at(std::abs(endVertex.id()) - 1)) {
0120         // We've seen this vertex before, we don't need to create it again
0121         continue;
0122       }
0123       seenVertices.at(std::abs(endVertex.id()) - 1) = true;
0124
0125       const auto endVertexPosition = convertPosition(endVertex.position());
0126
0127       ACTS_VERBOSE("Found secondary vertex at "
0128                    << endVertexPosition.transpose());
0129       double distance = (endVertexPosition.template head<3>() -
0130                          vertex.position4.template head<3>())
0131                             .squaredNorm();
0132
0133       if (distance <=
0134               m_cfg.vertexSpatialThreshold * m_cfg.vertexSpatialThreshold &&
0135           m_cfg.mergeSecondaries) {
0136         handleVertex(endVertex, vertex, vertices, particles, nSecondaryVertices,
0137                      nParticles, seenVertices);
0138       } else {
0139         // Over threshold, make a new vertex
0140         SimVertex secondaryVertex;
0141         nSecondaryVertices += 1;
0142         secondaryVertex.id =
0143             SimVertexBarcode{vertex.id}.withVertexSecondary(nSecondaryVertices);
0144         secondaryVertex.position4 = convertPosition(endVertex.position());
0145
0146         handleVertex(endVertex, secondaryVertex, vertices, particles,
0147                      nSecondaryVertices, nParticles, seenVertices);
0148
0149         // Only keep the secondary vertex if it has outgoing particles
0150         if (!secondaryVertex.outgoing.empty()) {
0151           vertices.push_back(secondaryVertex);
0152         }
0153       }
0154     } else {
0155       if (particle->status() != kUndecayedParticleStatus) {
0156         ACTS_ERROR("Undecayed particle has status "
0157                    << particle->status() << "(and not "
0158                    << kUndecayedParticleStatus << ")");
0159       }
0160       // This particle is a final state particle
0161       nParticles += 1;
0162       SimBarcode particleId =
0163           SimBarcode()
0164               .withVertexPrimary(vertex.vertexId().vertexPrimary())
0165               .withVertexSecondary(vertex.vertexId().vertexSecondary())
0166               .withParticle(nParticles);
0167
0168       Acts::PdgParticle pdg{particle->pdg_id()};
0169       double mass = 0.0;
0170       double charge = 0.0;
0171
0172       if (pdg != Acts::PdgParticle::eInvalid) {
0173         if (auto massOpt = Acts::findMass(pdg); massOpt.has_value()) {
0174           mass = massOpt.value();
0175         } else {
0176           ACTS_ERROR("No mass found for PDG ID " << pdg);
0177           throw std::bad_optional_access{};
0178         }
0179
0180         if (auto chargeOpt = Acts::findCharge(pdg); chargeOpt.has_value()) {
0181           charge = chargeOpt.value();
0182         } else {
0183           ACTS_ERROR("No charge found for PDG ID " << pdg);
0184           throw std::bad_optional_access{};
0185         }
0186       }
0187
0188       if (std::isnan(mass) || std::isnan(charge)) {
0189         // This particle does not have a mass or a charge based on our data
0190         // table:
0191         // => We can't handle it
0192         continue;
0193       }
0194
0195       SimParticle simParticle{particleId, pdg, charge, mass};
0196       simParticle.initialState().setPosition4(vertex.position4);
0197
0198       const HepMC3::FourVector& genMomentum = particle->momentum();
0199       Acts::Vector3 momentum{genMomentum.px() * 1_GeV, genMomentum.py() * 1_GeV,
0200                              genMomentum.pz() * 1_GeV};
0201       double p = momentum.norm();
0202
0203       simParticle.initialState().setDirection(momentum.normalized());
0204       simParticle.initialState().setAbsoluteMomentum(p);
0205
0206       particles.push_back(simParticle);
0207       vertex.outgoing.insert(particleId);
0208     }
0209   }
0210 }
0211
0212 void HepMC3InputConverter::convertHepMC3ToInternalEdm(
0213     const AlgorithmContext& ctx, const HepMC3::GenEvent& genEvent) const {
0214   ACTS_DEBUG("Converting HepMC3 event to internal EDM");
0215
0216   ACTS_VERBOSE("Have " << genEvent.vertices().size() << " vertices");
0217   ACTS_VERBOSE("Have " << genEvent.particles().size() << " particles");
0218
0219   using VertexCluster = std::vector<HepMC3::ConstGenVertexPtr>;
0220   std::vector<VertexCluster> vertexClusters;
0221
0222   ACTS_VERBOSE("Finding primary vertex clusters with threshold "
0223                << m_cfg.primaryVertexSpatialThreshold);
0224
0225   // Find all vertices whose incoming particles are either all beam particles or
0226   // that don't have any incoming particles
0227   {
0228     Acts::ScopedTimer timer("Finding primary vertex clusters", logger(),
0229                             Acts::Logging::DEBUG);
0230     for (const auto& vertex : genEvent.vertices()) {
0231       if (vertex->particles_in().empty() ||
0232           std::ranges::all_of(vertex->particles_in(), [](const auto& particle) {
0233             return particle->status() == kBeamParticleStatus;
0234           })) {
0235         // Check if primary vertex is within tolerance of an existing primary
0236         // vertex
0237
0238         ACTS_VERBOSE("Found primary vertex candidate at:\n"
0239                      << [&]() {
0240                           std::stringstream ss;
0241                           HepMC3::Print::listing(ss, vertex);
0242                           return ss.str();
0243                         }());
0244
0245         auto position = convertPosition(vertex->position());
0246
0247         if (auto it = std::ranges::find_if(
0248                 vertexClusters,
0249                 [&](const auto& cluster) {
0250                   const auto clusterPosition =
0251                       convertPosition(cluster.at(0)->position());
0252                   return (position - clusterPosition)
0253                              .template head<3>()
0254                              .cwiseAbs()
0255                              .maxCoeff() < m_cfg.primaryVertexSpatialThreshold;
0256                 });
0257             it != vertexClusters.end() && m_cfg.mergePrimaries) {
0258           // Add the vertex to the cluster
0259           it->push_back(vertex);
0260         } else {
0261           // Make a new cluster
0262           vertexClusters.push_back({vertex});
0263         }
0264       }
0265     }
0266   }
0267
0268   ACTS_DEBUG("Found " << vertexClusters.size()
0269                       << " primary vertex clusters (~ primary vertices)");
0270
0271   std::vector<SimVertex> verticesUnordered;
0272   std::vector<SimParticle> particlesUnordered;
0273
0274   std::vector<bool> seenVertices;
0275   seenVertices.resize(genEvent.vertices().size());
0276
0277   for (const auto& cluster : vertexClusters) {
0278     for (const auto& vertex : cluster) {
0279       seenVertices.at(std::abs(vertex->id()) - 1) = true;
0280     }
0281   }
0282
0283   std::size_t nPrimaryVertices = 0;
0284   {
0285     Acts::ScopedTimer timer("Converting HepMC3 vertices to internal EDM",
0286                             logger(), Acts::Logging::DEBUG);
0287
0288     std::size_t nUndecayedParticles = 0;
0289     for (const auto& particle : genEvent.particles()) {
0290       if (particle->end_vertex() == nullptr) {
0291         nUndecayedParticles += 1;
0292       }
0293     }
0294     ACTS_DEBUG("Found " << nUndecayedParticles
0295                         << " undecayed particles in the event");
0296     ACTS_DEBUG("Reserving " << genEvent.vertices().size() / 2
0297                             << " vertices and " << nUndecayedParticles
0298                             << " particles");
0299     verticesUnordered.reserve(genEvent.vertices().size() / 2);
0300     particlesUnordered.reserve(nUndecayedParticles);
0301
0302     for (auto& cluster : vertexClusters) {
0303       ACTS_VERBOSE("Primary vertex cluster at "
0304                    << convertPosition(cluster.at(0)->position()).transpose()
0305                    << " containing " << cluster.size() << " vertices:\n"
0306                    <<
0307                    [&]() {
0308                      std::stringstream ss;
0309                      for (auto& vertex : cluster) {
0310                        HepMC3::Print::listing(ss, vertex);
0311                      }
0312                      return ss.str();
0313                    }()
0314                    << "--------------");
0315
0316       nPrimaryVertices += 1;
0317       SimVertex primaryVertex;
0318       primaryVertex.id = SimVertexBarcode().withVertexPrimary(nPrimaryVertices);
0319       primaryVertex.position4 = convertPosition(cluster.at(0)->position());
0320
0321       std::size_t nSecondaryVertices = 0;
0322       std::size_t nParticles = 0;
0323
0324       for (auto& genVertex : cluster) {
0325         handleVertex(*genVertex, primaryVertex, verticesUnordered,
0326                      particlesUnordered, nSecondaryVertices, nParticles,
0327                      seenVertices);
0328       }
0329       verticesUnordered.push_back(primaryVertex);
0330     }
0331   }
0332
0333   ACTS_DEBUG("Converted " << particlesUnordered.size() << " particles and "
0334                           << verticesUnordered.size() << " vertices");
0335
0336   if (m_cfg.printListing) {
0337     ACTS_VERBOSE("Converted event record:\n"
0338                  << printListing(verticesUnordered, particlesUnordered));
0339   }
0340
0341   {
0342     Acts::ScopedTimer timer("Sorting particles and vertices", logger(),
0343                             Acts::Logging::DEBUG);
0344
0345     std::ranges::sort(particlesUnordered, [](const auto& a, const auto& b) {
0346       return a.particleId() < b.particleId();
0347     });
0348     std::ranges::sort(verticesUnordered, [](const auto& a, const auto& b) {
0349       return a.vertexId() < b.vertexId();
0350     });
0351   }
0352
0353   if (m_cfg.checkConsistency) {
0354     ACTS_DEBUG("Checking consistency of particles");
0355     auto equalParticleIds = [](const auto& a, const auto& b) {
0356       return a.particleId() == b.particleId();
0357     };
0358     if (auto it =
0359             std::ranges::adjacent_find(particlesUnordered, equalParticleIds);
0360         it != particlesUnordered.end()) {
0361       ACTS_ERROR("Found duplicate particle id " << it->particleId());
0362     }
0363     auto equalVertexIds = [](const auto& a, const auto& b) {
0364       return a.vertexId() == b.vertexId();
0365     };
0366     if (auto it = std::ranges::adjacent_find(verticesUnordered, equalVertexIds);
0367         it != verticesUnordered.end()) {
0368       ACTS_ERROR("Found duplicate vertex id " << it->vertexId());
0369     }
0370   }
0371
0372   SimParticleContainer particles{particlesUnordered.begin(),
0373                                  particlesUnordered.end()};
0374   SimVertexContainer vertices{verticesUnordered.begin(),
0375                               verticesUnordered.end()};
0376
0377   // move generated event to the store
0378   m_outputParticles(ctx, std::move(particles));
0379   m_outputVertices(ctx, std::move(vertices));
0380 }
0381
0382 }  // namespace ActsExamples