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