Geant4HepMC/src/SteppingAction.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 "SteppingAction.hpp"
0010
0011 #include "Acts/Utilities/Helpers.hpp"
0012
0013 #include <stdexcept>
0014
0015 #include <G4RunManager.hh>
0016 #include <G4Step.hh>
0017 #include <G4VProcess.hh>
0018 #include <HepMC3/Attribute.h>
0019 #include <HepMC3/Units.h>
0020
0021 #include "EventAction.hpp"
0022
0023 namespace ActsExamples::Geant4::HepMC3 {
0024
0025 SteppingAction* SteppingAction::s_instance = nullptr;
0026
0027 SteppingAction* SteppingAction::instance() {
0028   // Static access function via G4RunManager
0029   return s_instance;
0030 }
0031
0032 SteppingAction::SteppingAction(std::vector<std::string> eventRejectionProcess)
0033     : G4UserSteppingAction(),
0034       m_eventRejectionProcess(std::move(eventRejectionProcess)) {
0035   if (s_instance != nullptr) {
0036     throw std::logic_error("Attempted to duplicate a singleton");
0037   } else {
0038     s_instance = this;
0039   }
0040 }
0041
0042 SteppingAction::~SteppingAction() {
0043   s_instance = nullptr;
0044 }
0045
0046 void SteppingAction::UserSteppingAction(const G4Step* step) {
0047   // Test if the event should be aborted
0048   if (Acts::rangeContainsValue(m_eventRejectionProcess,
0049                                step->GetPostStepPoint()
0050                                    ->GetProcessDefinedStep()
0051                                    ->GetProcessName())) {
0052     m_eventAborted = true;
0053     G4RunManager::GetRunManager()->AbortEvent();
0054     return;
0055   }
0056
0057   /// Store the step such that a vertex knows the position and upcoming process
0058   /// for a particle. The particle properties are stored as ingoing before and
0059   /// as outgoing after the step with the process was performed. Therefore the
0060   /// vertex defines the starting point of a process while the next vertex
0061   /// describes the state after the process, including all particles produced
0062   /// along the step. In total the entire event is represented by vertices which
0063   /// describe each step in Geant4.
0064
0065   ::HepMC3::GenEvent& event = EventAction::instance()->event();
0066
0067   // Unit conversions G4->::HepMC3
0068   constexpr double convertLength = 1. / CLHEP::mm;
0069   constexpr double convertEnergy = 1. / CLHEP::GeV;
0070   constexpr double convertTime = 1. / CLHEP::s;
0071
0072   // The particle after the step
0073   auto* track = step->GetTrack();
0074   const std::string trackId = std::to_string(track->GetTrackID());
0075   auto postStepMomentum = track->GetMomentum() * convertEnergy;
0076   auto postStepEnergy = track->GetTotalEnergy() * convertEnergy;
0077   ::HepMC3::FourVector mom4{postStepMomentum[0], postStepMomentum[1],
0078                             postStepMomentum[2], postStepEnergy};
0079   auto postParticle = std::make_shared<::HepMC3::GenParticle>(
0080       mom4, track->GetDynamicParticle()->GetPDGcode());
0081
0082   // The process that led to the current state
0083   auto process = std::make_shared<::HepMC3::StringAttribute>(
0084       step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName());
0085
0086   // Assign particle to a production vertex
0087   if (!m_previousVertex) {
0088     // Get the production position of the particle
0089     auto* preStep = step->GetPreStepPoint();
0090     auto prePosition = preStep->GetPosition() * convertLength;
0091     auto preTime = preStep->GetGlobalTime() * convertTime;
0092     ::HepMC3::FourVector prePos{prePosition[0], prePosition[1], prePosition[2],
0093                                 preTime};
0094
0095     // Handle the first step: No vertices exist
0096     if (event.vertices().empty()) {
0097       auto vertex = std::make_shared<::HepMC3::GenVertex>(prePos);
0098       vertex->add_particle_out(postParticle);
0099       event.add_vertex(vertex);
0100       vertex->set_status(1);
0101       vertex->add_attribute("NextProcessOf" + trackId, process);
0102     } else {
0103       // Search for an existing vertex
0104       for (const auto& vertex : event.vertices()) {
0105         if (vertex->position() == prePos) {
0106           // Add particle to existing vertex
0107           vertex->add_particle_out(postParticle);
0108           vertex->add_attribute("NextProcessOf-" + trackId, process);
0109           auto preStepMomentum =
0110               step->GetPreStepPoint()->GetMomentum() * convertEnergy;
0111           auto preStepEnergy =
0112               step->GetPreStepPoint()->GetTotalEnergy() * convertEnergy;
0113           auto preMom4 = std::make_shared<::HepMC3::VectorDoubleAttribute>(
0114               std::vector<double>{preStepMomentum[0], preStepMomentum[1],
0115                                   preStepMomentum[2], preStepEnergy});
0116           vertex->add_attribute("InitialParametersOf-" + trackId, preMom4);
0117         }
0118       }
0119     }
0120     if (track->GetCreatorProcess() != nullptr) {
0121       postParticle->add_attribute(
0122           "CreatorProcessOf-" + trackId,
0123           std::make_shared<::HepMC3::StringAttribute>(
0124               track->GetCreatorProcess()->GetProcessName()));
0125     }
0126   } else {
0127     // Add particle from same track to vertex
0128     m_previousVertex->add_particle_out(postParticle);
0129     m_previousVertex->add_attribute("NextProcessOf-" + trackId, process);
0130   }
0131
0132   // Build the end vertex
0133   auto* postStep = step->GetPostStepPoint();
0134   auto postPosition = postStep->GetPosition() * convertLength;
0135   auto postTime = postStep->GetGlobalTime() * convertTime;
0136   ::HepMC3::FourVector postPos{postPosition[0], postPosition[1],
0137                                postPosition[2], postTime};
0138   m_previousVertex = std::make_shared<::HepMC3::GenVertex>(postPos);
0139
0140   // Add particle to the vertex
0141   m_previousVertex->add_particle_in(postParticle);
0142
0143   // Store the vertex
0144   event.add_vertex(m_previousVertex);
0145
0146   // Store additional data in the particle
0147   postParticle->add_attribute(
0148       "TrackID", std::make_shared<::HepMC3::IntAttribute>(track->GetTrackID()));
0149   postParticle->add_attribute(
0150       "ParentID",
0151       std::make_shared<::HepMC3::IntAttribute>(track->GetParentID()));
0152   const double X0 = track->GetMaterial()->GetRadlen() * convertLength;
0153   const double L0 =
0154       track->GetMaterial()->GetNuclearInterLength() * convertLength;
0155   const double stepLength = track->GetStepLength();
0156   postParticle->add_attribute("NextX0",
0157                               std::make_shared<::HepMC3::DoubleAttribute>(X0));
0158   postParticle->add_attribute("NextL0",
0159                               std::make_shared<::HepMC3::DoubleAttribute>(L0));
0160   postParticle->add_attribute(
0161       "StepLength", std::make_shared<::HepMC3::DoubleAttribute>(stepLength));
0162   postParticle->set_status(1);
0163
0164   // Stop tracking the vertex if the particle dies
0165   if (track->GetTrackStatus() != fAlive) {
0166     process = std::make_shared<::HepMC3::StringAttribute>("Death");
0167     m_previousVertex->add_attribute("NextProcessOf-" + trackId, process);
0168     m_previousVertex = nullptr;
0169   }
0170 }
0171
0172 void SteppingAction::clear() {
0173   m_previousVertex = nullptr;
0174   m_eventAborted = false;
0175 }
0176
0177 }  // namespace ActsExamples::Geant4::HepMC3