HepMC/src/HepMCProcessExtractor.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/HepMC/HepMCProcessExtractor.hpp"
0010
0011 #include "ActsExamples/Framework/WhiteBoard.hpp"
0012 #include "ActsExamples/Io/HepMC3/HepMC3Particle.hpp"
0013
0014 #include <algorithm>
0015 #include <stdexcept>
0016
0017 #include <HepMC3/GenEvent.h>
0018 #include <HepMC3/GenParticle.h>
0019 #include <HepMC3/GenVertex.h>
0020
0021 namespace ActsExamples {
0022
0023 namespace {
0024
0025 /// @brief This method searches for an outgoing particle from a vertex
0026 ///
0027 /// @param [in] vertex The vertex
0028 /// @param [in] id The track ID of the particle
0029 ///
0030 /// @return The particle pointer if found, else nullptr
0031 HepMC3::ConstGenParticlePtr searchProcessParticleById(
0032     const HepMC3::ConstGenVertexPtr& vertex, const int id) {
0033   // Loop over all outgoing particles
0034   for (const auto& particle : vertex->particles_out()) {
0035     const int trackid =
0036         particle->attribute<HepMC3::IntAttribute>("TrackID")->value();
0037     // Compare ID
0038     if (trackid == id) {
0039       return particle;
0040     }
0041   }
0042   return nullptr;
0043 }
0044
0045 /// @brief This method collects the material in X_0 and L_0 a particle has
0046 /// passed from its creation up to a certain vertex.
0047 ///
0048 /// @param [in] vertex The end vertex of the collection
0049 /// @param [in] id The track ID
0050 /// @param [in, out] particle The particle that get the passed material attached
0051 void setPassedMaterial(const HepMC3::ConstGenVertexPtr& vertex, const int id,
0052                        SimParticle& particle) {
0053   double x0 = 0.;
0054   double l0 = 0.;
0055   HepMC3::ConstGenParticlePtr currentParticle = nullptr;
0056   HepMC3::ConstGenVertexPtr currentVertex = vertex;
0057   // Loop backwards and test whether the track still exists
0058   while (currentVertex && !currentVertex->particles_in().empty() &&
0059          currentVertex->particles_in()[0]->attribute<HepMC3::IntAttribute>(
0060              "TrackID") &&
0061          currentVertex->particles_in()[0]
0062                  ->attribute<HepMC3::IntAttribute>("TrackID")
0063                  ->value() == id) {
0064     // Get the step length
0065     currentParticle = currentVertex->particles_in()[0];
0066     const double stepLength =
0067         currentParticle->attribute<HepMC3::DoubleAttribute>("StepLength")
0068             ->value();
0069     // Add the passed material
0070     x0 +=
0071         stepLength /
0072         currentParticle->attribute<HepMC3::DoubleAttribute>("NextX0")->value();
0073     l0 +=
0074         stepLength /
0075         currentParticle->attribute<HepMC3::DoubleAttribute>("NextL0")->value();
0076     currentVertex = currentParticle->production_vertex();
0077   }
0078   // Assign the passed material to the particle
0079   particle.final().setMaterialPassed(x0, l0);
0080 }
0081
0082 /// @brief This function collects outgoing particles from a vertex while keeping
0083 /// track of the future of the ingoing particle.
0084 ///
0085 /// @param [in] vertex The vertex
0086 /// @param [in] trackID The track ID of the ingoing particle
0087 ///
0088 /// @return Vector containing the outgoing particles from a vertex
0089 std::vector<SimParticle> selectOutgoingParticles(
0090     const HepMC3::ConstGenVertexPtr& vertex, const int trackID) {
0091   std::vector<SimParticle> finalStateParticles;
0092
0093   // Identify the ingoing particle in the outgoing particles
0094   HepMC3::ConstGenParticlePtr procPart =
0095       searchProcessParticleById(vertex, trackID);
0096
0097   // Test whether this particle survives or dies
0098   HepMC3::ConstGenVertexPtr endVertex = procPart->end_vertex();
0099   if (endVertex
0100           ->attribute<HepMC3::StringAttribute>("NextProcessOf-" +
0101                                                std::to_string(trackID))
0102           ->value() != "Death") {
0103     // Store the particle if it survives
0104     finalStateParticles.push_back(HepMC3Particle::particle(procPart));
0105   } else {
0106     // Store the leftovers if it dies
0107     for (const HepMC3::ConstGenParticlePtr& procPartOut :
0108          endVertex->particles_out()) {
0109       if (procPartOut->attribute<HepMC3::IntAttribute>("TrackID")->value() ==
0110               trackID &&
0111           procPartOut->end_vertex()) {
0112         for (const HepMC3::ConstGenParticlePtr& dyingPartOut :
0113              procPartOut->end_vertex()->particles_out()) {
0114           finalStateParticles.push_back(HepMC3Particle::particle(dyingPartOut));
0115         }
0116       }
0117     }
0118   }
0119
0120   // Record the particles produced in this process
0121   const std::vector<std::string> attributes = endVertex->attribute_names();
0122   for (const auto& att : attributes) {
0123     // Search for initial parameters
0124     if (att.find("InitialParametersOf") != std::string::npos) {
0125       const std::vector<double> mom4 =
0126           endVertex->attribute<HepMC3::VectorDoubleAttribute>(att)->value();
0127       const HepMC3::FourVector& pos4 = endVertex->position();
0128       const int id = stoi(att.substr(att.find("-") + 1));
0129       HepMC3::ConstGenParticlePtr genParticle =
0130           searchProcessParticleById(endVertex, id);
0131       ActsFatras::Barcode barcode = ActsFatras::Barcode().setParticle(id);
0132       auto pid = static_cast<Acts::PdgParticle>(genParticle->pid());
0133
0134       // Build an Acts particle out of the data
0135       SimParticleState simParticle(barcode, pid);
0136       simParticle.setPosition4(pos4.x(), pos4.y(), pos4.z(), pos4.t());
0137       Acts::Vector3 mom3(mom4[0], mom4[1], mom4[2]);
0138       simParticle.setDirection(mom3.normalized());
0139       simParticle.setAbsoluteMomentum(mom3.norm());
0140
0141       // Store the particle
0142       finalStateParticles.push_back(SimParticle(simParticle, simParticle));
0143     }
0144   }
0145
0146   return finalStateParticles;
0147 }
0148
0149 /// @brief This method filters and sorts the recorded interactions.
0150 ///
0151 /// @param [in] cfg Configuration of the filtering
0152 /// @param [in, out] interactions The recorded interactions
0153 void filterAndSort(const HepMCProcessExtractor::Config& cfg,
0154                    ExtractedSimulationProcessContainer& interactions) {
0155   for (auto& interaction : interactions) {
0156     for (auto cit = interaction.after.cbegin();
0157          cit != interaction.after.cend();) {
0158       // Test whether a particle fulfills the conditions
0159       if (cit->pdg() < cfg.absPdgMin || cit->pdg() > cfg.absPdgMax ||
0160           cit->absoluteMomentum() < cfg.pMin) {
0161         interaction.after.erase(cit);
0162       } else {
0163         cit++;
0164       }
0165     }
0166   }
0167
0168   // Sort the particles based on their momentum
0169   for (auto& interaction : interactions) {
0170     std::ranges::sort(interaction.after, std::greater{},
0171                       [](const auto& a) { return a.absoluteMomentum(); });
0172   }
0173 }
0174
0175 }  // namespace
0176
0177 HepMCProcessExtractor::~HepMCProcessExtractor() = default;
0178
0179 HepMCProcessExtractor::HepMCProcessExtractor(
0180     HepMCProcessExtractor::Config config, Acts::Logging::Level level)
0181     : IAlgorithm("HepMCProcessExtractor", level), m_cfg(std::move(config)) {
0182   if (m_cfg.inputEvents.empty()) {
0183     throw std::invalid_argument("Missing input event collection");
0184   }
0185   if (m_cfg.outputSimulationProcesses.empty()) {
0186     throw std::invalid_argument("Missing output collection");
0187   }
0188   if (m_cfg.extractionProcess.empty()) {
0189     throw std::invalid_argument("Missing extraction process");
0190   }
0191
0192   m_inputEvents.initialize(m_cfg.inputEvents);
0193   m_outputSimulationProcesses.initialize(m_cfg.outputSimulationProcesses);
0194 }
0195
0196 ProcessCode HepMCProcessExtractor::execute(
0197     const AlgorithmContext& context) const {
0198   // Retrieve the initial particles
0199   const auto events = m_inputEvents(context);
0200
0201   ExtractedSimulationProcessContainer fractions;
0202   for (const HepMC3::GenEvent& event : events) {
0203     // Fast exit
0204     if (event.particles().empty() || event.vertices().empty()) {
0205       break;
0206     }
0207
0208     // Get the initial particle
0209     HepMC3::ConstGenParticlePtr initialParticle = event.particles()[0];
0210     SimParticle simParticle = HepMC3Particle::particle(initialParticle);
0211
0212     // Get the final state particles
0213     SimParticle particleToInteraction;
0214     std::vector<SimParticle> finalStateParticles;
0215     // Search the process vertex
0216     bool vertexFound = false;
0217     for (const auto& vertex : event.vertices()) {
0218       const std::vector<std::string> attributes = vertex->attribute_names();
0219       for (const auto& attribute : attributes) {
0220         if (vertex->attribute_as_string(attribute).find(
0221                 m_cfg.extractionProcess) != std::string::npos) {
0222           const int procID = stoi(attribute.substr(attribute.find("-") + 1));
0223           // Get the particle before the interaction
0224           particleToInteraction =
0225               HepMC3Particle::particle(vertex->particles_in()[0]);
0226           // Attach passed material to the particle
0227           setPassedMaterial(vertex, procID, particleToInteraction);
0228           // Record the final state particles
0229           finalStateParticles = selectOutgoingParticles(vertex, procID);
0230           vertexFound = true;
0231           break;
0232         }
0233       }
0234       if (vertexFound) {
0235         break;
0236       }
0237     }
0238     fractions.push_back(ExtractedSimulationProcess{
0239         simParticle, particleToInteraction, finalStateParticles});
0240   }
0241
0242   // Filter and sort the record
0243   filterAndSort(m_cfg, fractions);
0244
0245   ACTS_INFO(events.size() << " processed");
0246
0247   // Write the recorded material to the event store
0248   m_outputSimulationProcesses(context, std::move(fractions));
0249
0250   return ProcessCode::SUCCESS;
0251 }
0252
0253 }  // namespace ActsExamples