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