HepMC3/src/HepMC3OutputConverter.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/HepMC3OutputConverter.hpp"
0010
0011 #include "ActsExamples/Framework/IAlgorithm.hpp"
0012
0013 #include <HepMC3/GenEvent.h>
0014 #include <HepMC3/GenParticle.h>
0015 #include <HepMC3/GenVertex.h>
0016
0017 using namespace Acts;
0018 using namespace Acts::UnitLiterals;
0019
0020 namespace ActsExamples {
0021
0022 HepMC3OutputConverter::HepMC3OutputConverter(const Config& config,
0023                                              Acts::Logging::Level level)
0024     : IAlgorithm{"HepMC3OutputConverter", level}, m_cfg{config} {
0025   if (m_cfg.inputParticles.empty()) {
0026     throw std::invalid_argument("Missing input particles collection");
0027   }
0028
0029   if (m_cfg.inputVertices.empty()) {
0030     throw std::invalid_argument("Missing input vertices collection");
0031   }
0032
0033   if (m_cfg.outputEvent.empty()) {
0034     throw std::invalid_argument("Missing output event collection");
0035   }
0036
0037   m_inputParticles.initialize(m_cfg.inputParticles);
0038   m_inputVertices.initialize(m_cfg.inputVertices);
0039   m_outputEvent.initialize(m_cfg.outputEvent);
0040 }
0041
0042 ProcessCode HepMC3OutputConverter::execute(const AlgorithmContext& ctx) const {
0043   auto genEvent =
0044       std::make_shared<HepMC3::GenEvent>(HepMC3::Units::GEV, HepMC3::Units::MM);
0045
0046   const auto& inputParticles = m_inputParticles(ctx);
0047   const auto& inputVertices = m_inputVertices(ctx);
0048   ACTS_DEBUG("Converting " << inputParticles.size() << " and "
0049                            << inputVertices.size() << " vertices to HepMC3");
0050
0051   std::unordered_map<SimBarcode, std::shared_ptr<HepMC3::GenParticle>>
0052       barcodeMap;
0053   for (const auto& inParticle : inputParticles) {
0054     const Vector4 momentum =
0055         inParticle.fourMomentum() / 1_GeV;  // Let's ensure we're using GeV
0056     const HepMC3::FourVector vec(momentum[0], momentum[1], momentum[2],
0057                                  momentum[3]);
0058     auto hepmc3Particle =
0059         std::make_shared<HepMC3::GenParticle>(vec, inParticle.pdg());
0060     hepmc3Particle->set_generated_mass(inParticle.mass());
0061     hepmc3Particle->set_status(1);
0062     genEvent->add_particle(hepmc3Particle);
0063     ACTS_VERBOSE("Adding particle with barcode " << inParticle.particleId());
0064     barcodeMap.try_emplace(inParticle.particleId(), hepmc3Particle);
0065   }
0066
0067   for (const auto& inVertex : inputVertices) {
0068     Vector4 position = inVertex.position4;
0069     // Let's make sure we use correct units
0070     position.head<3>() /= 1_mm;
0071     position[3] /= 1_ns;
0072     const HepMC3::FourVector vec(position[0], position[1], position[2],
0073                                  position[3]);
0074     auto hepmc3Vertex = std::make_shared<HepMC3::GenVertex>(vec);
0075
0076     if (inVertex.incoming.empty()) {
0077       // HepMC3 does not like vertices without incoming particles
0078       // Add our dummy incoming particle
0079       ACTS_VERBOSE("Adding dummy incoming particle to vertex");
0080       auto dummyParticle = std::make_shared<HepMC3::GenParticle>(
0081           HepMC3::FourVector(0, 0, 0, 0), 0, 666);
0082       dummyParticle->set_generated_mass(123);
0083       dummyParticle->set_status(12);
0084       genEvent->add_particle(dummyParticle);
0085       hepmc3Vertex->add_particle_in(dummyParticle);
0086     }
0087
0088     hepmc3Vertex->set_status(1);
0089     genEvent->add_vertex(hepmc3Vertex);
0090
0091     for (const auto& particleId : inVertex.incoming) {
0092       auto it = barcodeMap.find(particleId);
0093       if (it == barcodeMap.end()) {
0094         ACTS_ERROR("Particle with barcode " << particleId << " not found");
0095         continue;
0096       }
0097       hepmc3Vertex->add_particle_in(it->second);
0098     }
0099
0100     for (const auto& particleId : inVertex.outgoing) {
0101       auto it = barcodeMap.find(particleId);
0102       if (it == barcodeMap.end()) {
0103         ACTS_ERROR("Particle with barcode " << particleId << " not found");
0104         continue;
0105       }
0106       hepmc3Vertex->add_particle_out(it->second);
0107     }
0108   }
0109
0110   m_outputEvent(ctx, std::move(genEvent));
0111
0112   return ProcessCode::SUCCESS;
0113 }
0114
0115 }  // namespace ActsExamples