Geant4HepMC/src/EventAction.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "EventAction.hpp"
0010
0011 #include <stdexcept>
0012
0013 #include <G4Event.hh>
0014 #include <G4RunManager.hh>
0015
0016 #include "SteppingAction.hpp"
0017
0018 namespace {
0019
0020 /// @brief This function tests whether a process is available in the record
0021 ///
0022 /// @param [in] vertex The vertex that will be tested
0023 /// @param [in] processFilter List of processes that will be filtered
0024 ///
0025 /// @return True if the process was found, false if not
0026 bool findAttribute(const HepMC3::ConstGenVertexPtr& vertex,
0027                    const std::vector<std::string>& processFilter) {
0028   // Consider only 1->1 vertices to keep a correct history
0029   if ((vertex->particles_in().size() == 1) &&
0030       (vertex->particles_out().size() == 1)) {
0031     // Test for all attributes if one matches the filter pattern
0032     const std::vector<std::string> vertexAttributes = vertex->attribute_names();
0033     for (const auto& att : vertexAttributes) {
0034       const std::string process = vertex->attribute_as_string(att);
0035       if (std::find(processFilter.begin(), processFilter.end(), process) !=
0036           processFilter.end()) {
0037         return true;
0038       }
0039     }
0040   }
0041   return false;
0042 }
0043
0044 /// @brief This function reduces multiple vertices that should be filtered and
0045 /// combines them in a single dummy vertex
0046 ///
0047 /// @param [in, out] event The underlying event
0048 /// @param [in, out] vertex The vertex that will be reduced
0049 /// @param [in] processFilter List of processes that will be filtered
0050 void reduceVertex(HepMC3::GenEvent& event, HepMC3::GenVertexPtr vertex,
0051                   const std::vector<std::string>& processFilter) {
0052   // Store the particles associated to the vertex
0053   HepMC3::GenParticlePtr particleIn = vertex->particles_in()[0];
0054   HepMC3::GenParticlePtr particleOut = vertex->particles_out()[0];
0055
0056   // Walk backwards to find all vertices in a row that match the pattern
0057   while (findAttribute(particleIn->production_vertex(), processFilter)) {
0058     // Take the particles before the vertex and remove particles & vertices in
0059     // between
0060     HepMC3::GenVertexPtr currentVertex = particleOut->production_vertex();
0061     HepMC3::GenParticlePtr nextParticle = currentVertex->particles_in()[0];
0062     // Cut connections from particle and remove it
0063     if (particleIn->end_vertex()) {
0064       particleIn->end_vertex()->remove_particle_in(particleIn);
0065     }
0066     currentVertex->remove_particle_out(particleIn);
0067     event.remove_particle(particleIn);
0068     // Cut connections from vertext and remove it
0069     particleIn = nextParticle;
0070     currentVertex->remove_particle_in(particleIn);
0071     event.remove_vertex(currentVertex);
0072   }
0073   // Walk forwards to find all vertices in a row that match the pattern
0074   while (findAttribute(particleOut->end_vertex(), processFilter)) {
0075     // Take the particles after the vertex and remove particles & vertices in
0076     // between
0077     HepMC3::GenVertexPtr currentVertex = particleOut->end_vertex();
0078     HepMC3::GenParticlePtr nextParticle = currentVertex->particles_out()[0];
0079     // Cut connections from particle and remove it
0080     if (particleOut->production_vertex()) {
0081       particleOut->production_vertex()->remove_particle_out(particleOut);
0082     }
0083     currentVertex->remove_particle_in(particleOut);
0084     event.remove_particle(particleOut);
0085     // Cut connections from vertext and remove it
0086     particleOut = nextParticle;
0087     currentVertex->remove_particle_out(particleOut);
0088     event.remove_vertex(currentVertex);
0089   }
0090
0091   // Build the dummy vertex
0092   auto reducedVertex = std::make_shared<HepMC3::GenVertex>();
0093   event.add_vertex(reducedVertex);
0094
0095   reducedVertex->add_particle_in(particleIn);
0096   reducedVertex->add_particle_out(particleOut);
0097
0098   // Remove and replace the old vertex
0099   event.remove_vertex(vertex);
0100   vertex = reducedVertex;
0101 }
0102
0103 /// @brief This method walks over all vertices and tests whether it should be
0104 /// filtered
0105 ///
0106 /// @param [in, out] event The underlying event
0107 /// @param [in, out] The current vertex under investigation
0108 /// @param [in] processFilter List of processes that will be filtered
0109 void followOutgoingParticles(HepMC3::GenEvent& event,
0110                              const HepMC3::GenVertexPtr& vertex,
0111                              const std::vector<std::string>& processFilter) {
0112   // Replace and reduce vertex if it should be filtered
0113   if (findAttribute(vertex, processFilter)) {
0114     reduceVertex(event, vertex, processFilter);
0115   }
0116   // Move forward to the next vertices
0117   for (const auto& particle : vertex->particles_out()) {
0118     followOutgoingParticles(event, particle->end_vertex(), processFilter);
0119   }
0120 }
0121 }  // namespace
0122
0123 namespace ActsExamples::Geant4::HepMC3 {
0124
0125 EventAction* EventAction::s_instance = nullptr;
0126
0127 EventAction* EventAction::instance() {
0128   // Static access function via G4RunManager
0129   return s_instance;
0130 }
0131
0132 EventAction::EventAction(std::vector<std::string> processFilter)
0133     : G4UserEventAction(), m_processFilter(std::move(processFilter)) {
0134   if (s_instance != nullptr) {
0135     throw std::logic_error("Attempted to duplicate a singleton");
0136   } else {
0137     s_instance = this;
0138   }
0139 }
0140
0141 EventAction::~EventAction() {
0142   s_instance = nullptr;
0143 }
0144
0145 void EventAction::BeginOfEventAction(const G4Event* /*event*/) {
0146   SteppingAction::instance()->clear();
0147   m_event = ::HepMC3::GenEvent(::HepMC3::Units::GEV, ::HepMC3::Units::MM);
0148   m_event.add_beam_particle(std::make_shared<::HepMC3::GenParticle>());
0149 }
0150
0151 void EventAction::EndOfEventAction(const G4Event* /*event*/) {
0152   // Fast exit if the event is empty
0153   if (m_event.vertices().empty()) {
0154     return;
0155   }
0156   // Filter irrelevant processes
0157   auto currentVertex = m_event.vertices()[0];
0158   for (auto& bp : m_event.beams()) {
0159     if (!bp->end_vertex()) {
0160       currentVertex->add_particle_in(bp);
0161     }
0162   }
0163   followOutgoingParticles(m_event, currentVertex, m_processFilter);
0164   // Remove vertices w/o outgoing particles and particles w/o production
0165   // vertices
0166   while (true) {
0167     bool sane = true;
0168     for (const auto& v : m_event.vertices()) {
0169       if (!v) {
0170         continue;
0171       }
0172       if (v->particles_out().empty()) {
0173         m_event.remove_vertex(v);
0174         sane = false;
0175       }
0176     }
0177     for (const auto& p : m_event.particles()) {
0178       if (!p) {
0179         continue;
0180       }
0181       if (!p->production_vertex()) {
0182         m_event.remove_particle(p);
0183         sane = false;
0184       }
0185     }
0186     if (sane) {
0187       break;
0188     }
0189   }
0190 }
0191
0192 void EventAction::clear() {
0193   SteppingAction::instance()->clear();
0194 }
0195
0196 ::HepMC3::GenEvent& EventAction::event() {
0197   return m_event;
0198 }
0199 }  // namespace ActsExamples::Geant4::HepMC3