Geant4/src/SensitiveSteppingAction.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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 "ActsExamples/Geant4/SensitiveSteppingAction.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Utilities/MultiIndex.hpp"
0015 #include "ActsExamples/EventData/SimHit.hpp"
0016 #include "ActsExamples/Geant4/EventStore.hpp"
0017 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0018 #include "ActsFatras/EventData/Barcode.hpp"
0019
0020 #include <cstddef>
0021 #include <string>
0022 #include <unordered_map>
0023 #include <utility>
0024
0025 #include <G4RunManager.hh>
0026 #include <G4Step.hh>
0027 #include <G4StepPoint.hh>
0028 #include <G4Track.hh>
0029 #include <G4UnitsTable.hh>
0030 #include <G4VPhysicalVolume.hh>
0031
0032 class G4PrimaryParticle;
0033
0034 #if BOOST_VERSION >= 107800
0035 #include <boost/describe.hpp>
0036
0037 BOOST_DESCRIBE_ENUM(G4StepStatus, fWorldBoundary, fGeomBoundary,
0038                     fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc,
0039                     fUserDefinedLimit, fExclusivelyForcedProc, fUndefined);
0040
0041 BOOST_DESCRIBE_ENUM(G4ProcessType, fNotDefined, fTransportation,
0042                     fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
0043                     fDecay, fGeneral, fParameterisation, fUserDefined,
0044                     fParallel, fPhonon, fUCN);
0045
0046 BOOST_DESCRIBE_ENUM(G4TrackStatus, fAlive, fStopButAlive, fStopAndKill,
0047                     fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent);
0048 #endif
0049
0050 namespace {
0051
0052 ActsFatras::Hit hitFromStep(const G4StepPoint* preStepPoint,
0053                             const G4StepPoint* postStepPoint,
0054                             ActsFatras::Barcode particleId,
0055                             Acts::GeometryIdentifier geoId, int32_t index) {
0056   static constexpr double convertTime = Acts::UnitConstants::s / CLHEP::s;
0057   static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
0058   static constexpr double convertEnergy = Acts::UnitConstants::GeV / CLHEP::GeV;
0059
0060   G4ThreeVector preStepPosition = convertLength * preStepPoint->GetPosition();
0061   G4double preStepTime = convertTime * preStepPoint->GetGlobalTime();
0062   G4ThreeVector postStepPosition = convertLength * postStepPoint->GetPosition();
0063   G4double postStepTime = convertTime * postStepPoint->GetGlobalTime();
0064
0065   G4ThreeVector preStepMomentum = convertEnergy * preStepPoint->GetMomentum();
0066   G4double preStepEnergy = convertEnergy * preStepPoint->GetTotalEnergy();
0067   G4ThreeVector postStepMomentum = convertEnergy * postStepPoint->GetMomentum();
0068   G4double postStepEnergy = convertEnergy * postStepPoint->GetTotalEnergy();
0069
0070   Acts::ActsScalar hX = 0.5 * (preStepPosition[0] + postStepPosition[0]);
0071   Acts::ActsScalar hY = 0.5 * (preStepPosition[1] + postStepPosition[1]);
0072   Acts::ActsScalar hZ = 0.5 * (preStepPosition[2] + postStepPosition[2]);
0073   Acts::ActsScalar hT = 0.5 * (preStepTime + postStepTime);
0074
0075   Acts::ActsScalar mXpre = preStepMomentum[0];
0076   Acts::ActsScalar mYpre = preStepMomentum[1];
0077   Acts::ActsScalar mZpre = preStepMomentum[2];
0078   Acts::ActsScalar mEpre = preStepEnergy;
0079   Acts::ActsScalar mXpost = postStepMomentum[0];
0080   Acts::ActsScalar mYpost = postStepMomentum[1];
0081   Acts::ActsScalar mZpost = postStepMomentum[2];
0082   Acts::ActsScalar mEpost = postStepEnergy;
0083
0084   Acts::Vector4 particlePosition(hX, hY, hZ, hT);
0085   Acts::Vector4 beforeMomentum(mXpre, mYpre, mZpre, mEpre);
0086   Acts::Vector4 afterMomentum(mXpost, mYpost, mZpost, mEpost);
0087
0088   return ActsFatras::Hit(geoId, particleId, particlePosition, beforeMomentum,
0089                          afterMomentum, index);
0090 }
0091 }  // namespace
0092
0093 ActsExamples::SensitiveSteppingAction::SensitiveSteppingAction(
0094     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0095     : G4UserSteppingAction(), m_cfg(cfg), m_logger(std::move(logger)) {}
0096
0097 void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
0098     const G4Step* step) {
0099   // Unit conversions G4->::ACTS
0100
0101   // The particle after the step
0102   G4Track* track = step->GetTrack();
0103   G4PrimaryParticle* primaryParticle =
0104       track->GetDynamicParticle()->GetPrimaryParticle();
0105
0106   // Bail out if charged & configured to do so
0107   G4double absCharge = std::abs(track->GetParticleDefinition()->GetPDGCharge());
0108   if (!m_cfg.charged && absCharge > 0.) {
0109     return;
0110   }
0111
0112   // Bail out if neutral & configured to do so
0113   if (!m_cfg.neutral && absCharge == 0.) {
0114     return;
0115   }
0116
0117   // Bail out if it is a primary & configured to be ignored
0118   if (!m_cfg.primary && primaryParticle != nullptr) {
0119     return;
0120   }
0121
0122   // Bail out if it is a secondary & configured to be ignored
0123   if (!m_cfg.secondary && primaryParticle == nullptr) {
0124     return;
0125   }
0126
0127   // Get the physical volume & check if it has the sensitive string name
0128   std::string_view volumeName = track->GetVolume()->GetName();
0129
0130   if (volumeName.find(SensitiveSurfaceMapper::mappingPrefix) ==
0131       std::string_view::npos) {
0132     return;
0133   }
0134
0135   // Cast out the GeometryIdentifier
0136   volumeName = volumeName.substr(SensitiveSurfaceMapper::mappingPrefix.size(),
0137                                  std::string_view::npos);
0138   char* end = nullptr;
0139   const Acts::GeometryIdentifier geoId(
0140       std::strtoul(volumeName.data(), &end, 10));
0141
0142   // This is not the case if we have a particle-ID collision
0143   if (eventStore().trackIdMapping.find(track->GetTrackID()) ==
0144       eventStore().trackIdMapping.end()) {
0145     return;
0146   }
0147
0148   const auto particleId = eventStore().trackIdMapping.at(track->GetTrackID());
0149
0150   ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
0151
0152   // Get PreStepPoint and PostStepPoint
0153   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0154   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0155
0156   // Set particle hit count to zero, so we have this entry in the map later
0157   if (eventStore().particleHitCount.find(particleId) ==
0158       eventStore().particleHitCount.end()) {
0159     eventStore().particleHitCount[particleId] = 0;
0160   }
0161
0162   // Extract if we are at volume boundaries
0163   const bool preOnBoundary = preStepPoint->GetStepStatus() == fGeomBoundary;
0164   const bool postOnBoundary = postStepPoint->GetStepStatus() == fGeomBoundary ||
0165                               postStepPoint->GetStepStatus() == fWorldBoundary;
0166   const bool particleStopped = (postStepPoint->GetKineticEnergy() == 0.0);
0167   const bool particleDecayed =
0168       (postStepPoint->GetProcessDefinedStep()->GetProcessType() == fDecay);
0169
0170   auto print = [](auto s) {
0171 #if BOOST_VERSION >= 107800
0172     return boost::describe::enum_to_string(s, "unmatched");
0173 #else
0174     return s;
0175 #endif
0176   };
0177   ACTS_VERBOSE("status: pre="
0178                << print(preStepPoint->GetStepStatus())
0179                << ", post=" << print(postStepPoint->GetStepStatus())
0180                << ", post E_kin=" << std::boolalpha
0181                << postStepPoint->GetKineticEnergy() << ", process_type="
0182                << print(
0183                       postStepPoint->GetProcessDefinedStep()->GetProcessType())
0184                << ", particle="
0185                << track->GetParticleDefinition()->GetParticleName()
0186                << ", process_name="
0187                << postStepPoint->GetProcessDefinedStep()->GetProcessName()
0188                << ", track status=" << print(track->GetTrackStatus()));
0189
0190   // Case A: The step starts at the entry of the volume and ends at the exit.
0191   // Add hit to collection.
0192   if (preOnBoundary && postOnBoundary) {
0193     ACTS_VERBOSE("-> merge single step to hit");
0194     ++eventStore().particleHitCount[particleId];
0195     eventStore().hits.push_back(
0196         hitFromStep(preStepPoint, postStepPoint, particleId, geoId,
0197                     eventStore().particleHitCount.at(particleId) - 1));
0198
0199     eventStore().numberGeantSteps += 1ul;
0200     eventStore().maxStepsForHit = std::max(eventStore().maxStepsForHit, 1ul);
0201     return;
0202   }
0203
0204   // Case B: The step ends at the exit of the volume. Add hit to hit buffer, and
0205   // then combine hit buffer
0206   if (postOnBoundary || particleStopped || particleDecayed) {
0207     ACTS_VERBOSE("-> merge buffer to hit");
0208     auto& buffer = eventStore().hitBuffer;
0209     buffer.push_back(
0210         hitFromStep(preStepPoint, postStepPoint, particleId, geoId, -1));
0211
0212     const auto pos4 =
0213         0.5 * (buffer.front().fourPosition() + buffer.back().fourPosition());
0214
0215     ++eventStore().particleHitCount[particleId];
0216     eventStore().hits.emplace_back(
0217         geoId, particleId, pos4, buffer.front().momentum4Before(),
0218         buffer.back().momentum4After(),
0219         eventStore().particleHitCount.at(particleId) - 1);
0220
0221     assert(std::all_of(buffer.begin(), buffer.end(),
0222                        [&](const auto& h) { return h.geometryId() == geoId; }));
0223     assert(std::all_of(buffer.begin(), buffer.end(), [&](const auto& h) {
0224       return h.particleId() == particleId;
0225     }));
0226
0227     eventStore().numberGeantSteps += buffer.size();
0228     eventStore().maxStepsForHit =
0229         std::max(eventStore().maxStepsForHit, buffer.size());
0230
0231     buffer.clear();
0232     return;
0233   }
0234
0235   // Case C: The step doesn't end at the exit of the volume. Add the hit to the
0236   // hit buffer.
0237   if (!postOnBoundary) {
0238     // ACTS_VERBOSE("-> add hit to buffer");
0239     eventStore().hitBuffer.push_back(
0240         hitFromStep(preStepPoint, postStepPoint, particleId, geoId, -1));
0241     return;
0242   }
0243
0244   assert(false && "should never reach this");
0245 }