src/DDG4/ActsMaterialRecorder.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 "Acts/Definitions/Units.hpp"
0010 #include "Acts/Material/MaterialSlab.hpp"
0011 #include "ActsPlugins/Geant4/Geant4MaterialStepConverter.hpp"
0012
0013 #include <algorithm>
0014 #include <cmath>
0015 #include <limits>
0016 #include <optional>
0017 #include <string>
0018
0019 #include <CLHEP/Units/SystemOfUnits.h>
0020 #include <DDG4/Factories.h>
0021 #include <DDG4/Geant4Context.h>
0022 #include <DDG4/Geant4EventAction.h>
0023 #include <DDG4/Geant4RunAction.h>
0024 #include <DDG4/Geant4SteppingAction.h>
0025 #include <G4Event.hh>
0026 #include <G4Run.hh>
0027 #include <G4Step.hh>
0028 #include <G4StepPoint.hh>
0029 #include <G4Track.hh>
0030
0031 #include "detail/RecordedMaterialTracks.hpp"
0032
0033 namespace {
0034
0035 constexpr double infinity = std::numeric_limits<double>::infinity();
0036
0037 /// Cylinder around the beam axis that limits the recording, in Geant4 units.
0038 /// It does not need to be symmetric in z.
0039 struct RecordingRegion {
0040   double maxR = infinity;
0041   double minZ = -infinity;
0042   double maxZ = infinity;
0043
0044   bool isLimited() const {
0045     return std::isfinite(maxR) || std::isfinite(minZ) || std::isfinite(maxZ);
0046   }
0047
0048   bool contains(const G4ThreeVector& position) const {
0049     return position.perp() < maxR && position.z() > minZ && position.z() < maxZ;
0050   }
0051
0052   /// Path length from @p position, inside the region, along the unit vector
0053   /// @p direction to the boundary
0054   double pathToBoundary(const G4ThreeVector& position,
0055                         const G4ThreeVector& direction) const {
0056     double pathR = infinity;
0057     if (const double a = direction.perp2(); a > 0. && std::isfinite(maxR)) {
0058       const double b =
0059           position.x() * direction.x() + position.y() * direction.y();
0060       const double c = position.perp2() - maxR * maxR;
0061       pathR = (-b + std::sqrt(b * b - a * c)) / a;
0062     }
0063     double pathZ = infinity;
0064     if (const double zEnd = direction.z() > 0. ? maxZ : minZ;
0065         direction.z() != 0. && std::isfinite(zEnd)) {
0066       pathZ = (zEnd - position.z()) / direction.z();
0067     }
0068     return std::min(pathR, pathZ);
0069   }
0070 };
0071
0072 /// Shorten the last step of @p track to @p length, and its material with it
0073 void clipLastStep(Acts::RecordedMaterialTrack& track, double length) {
0074   Acts::RecordedMaterial& recorded = track.second;
0075   Acts::MaterialInteraction& interaction = recorded.materialInteractions.back();
0076   const Acts::MaterialSlab original = interaction.materialSlab;
0077   interaction.materialSlab =
0078       Acts::MaterialSlab(original.material(), static_cast<float>(length));
0079   interaction.pathCorrection = length;
0080   recorded.materialInX0 +=
0081       interaction.materialSlab.thicknessInX0() - original.thicknessInX0();
0082   recorded.materialInL0 +=
0083       interaction.materialSlab.thicknessInL0() - original.thicknessInL0();
0084 }
0085
0086 }  // namespace
0087
0088 namespace ActsPlugins {
0089
0090 /// @brief DDG4 stepping action that records the material seen by tracks
0091 ///
0092 /// This allows running a material recording with geantinos directly in
0093 /// `ddsim`, on the Geant4 geometry DDG4 builds from the DD4hep detector.
0094 /// Each track is converted with @ref ActsPlugins::Geant4MaterialStepConverter
0095 /// into an @ref Acts::RecordedMaterialTrack, and the tracks of an event are
0096 /// attached to the DDG4 event as a `RecordedMaterialTracks` extension.
0097 /// Tracks that cross no recorded material are kept without steps, so there is
0098 /// exactly one material track per simulated track.
0099 ///
0100 /// The tracks are written by the `ActsMaterialTrackWriter` event action. In
0101 /// multi-threaded running there is one recorder per worker thread.
0102 ///
0103 /// The recording can be limited to a cylinder around the beam axis, e.g. to
0104 /// stop at the calorimeter face: the material mapping assigns every recorded
0105 /// step to a surface, so material beyond the outermost surfaces would end up
0106 /// on them. A track is stopped where it leaves the cylinder, and the step that
0107 /// crosses the boundary is cut there.
0108 ///
0109 /// Properties:
0110 /// - `ExcludeMaterials`: names of materials whose steps are not recorded,
0111 ///   default `["Air", "Vacuum"]`
0112 /// - `RecordElementFractions`: whether to record the element composition of
0113 ///   the stepped materials, default `false`
0114 /// - `MaxR`, `MinZ`, `MaxZ`: radius and z range of the recording region in
0115 ///   Geant4 units, by default infinite, i.e. no limit. The z range does not
0116 ///   need to be symmetric.
0117 class ActsMaterialRecorder final : public dd4hep::sim::Geant4SteppingAction {
0118  public:
0119   /// Standard constructor of a DDG4 action
0120   /// @param context the DDG4 context
0121   /// @param name the name of the action instance
0122   ActsMaterialRecorder(dd4hep::sim::Geant4Context* context,
0123                        const std::string& name);
0124
0125   /// Record the material of one step
0126   /// @param step the Geant4 step
0127   /// @param mgr the Geant4 stepping manager (unused)
0128   void operator()(const G4Step* step, G4SteppingManager* mgr) override;
0129
0130  private:
0131   void beginRun(const G4Run* run);
0132   void beginEvent(const G4Event* event);
0133
0134   /// Configuration of the step conversion, set through the properties
0135   ActsPlugins::Geant4MaterialStepConverter::Config m_stepConverterCfg{
0136       .excludeMaterials = {"Air", "Vacuum"}, .recordElementFractions = false};
0137
0138   /// Created at the start of the run, once the properties are set
0139   std::optional<ActsPlugins::Geant4MaterialStepConverter> m_stepConverter;
0140
0141   /// The recording region, set through the properties
0142   RecordingRegion m_regionCfg;
0143   /// The recording region, if limited
0144   std::optional<RecordingRegion> m_region;
0145
0146   /// Tracks of the current event, owned by the DDG4 event
0147   detail::RecordedMaterialTracks* m_recorded = nullptr;
0148 };
0149
0150 ActsMaterialRecorder::ActsMaterialRecorder(dd4hep::sim::Geant4Context* context,
0151                                            const std::string& name)
0152     : dd4hep::sim::Geant4SteppingAction(context, name) {
0153   declareProperty("ExcludeMaterials", m_stepConverterCfg.excludeMaterials);
0154   declareProperty("RecordElementFractions",
0155                   m_stepConverterCfg.recordElementFractions);
0156   declareProperty("MaxR", m_regionCfg.maxR);
0157   declareProperty("MinZ", m_regionCfg.minZ);
0158   declareProperty("MaxZ", m_regionCfg.maxZ);
0159
0160   runAction().callAtBegin(this, &ActsMaterialRecorder::beginRun);
0161   eventAction().callAtBegin(this, &ActsMaterialRecorder::beginEvent);
0162 }
0163
0164 void ActsMaterialRecorder::beginRun(const G4Run* /*run*/) {
0165   // Properties are only set after construction, so configure here
0166   m_stepConverter.emplace(m_stepConverterCfg);
0167
0168   m_region.reset();
0169   if (m_regionCfg.isLimited()) {
0170     if (!(m_regionCfg.maxR > 0.) || !(m_regionCfg.minZ < m_regionCfg.maxZ)) {
0171       except(
0172           "Invalid recording region: MaxR must be positive and MinZ smaller "
0173           "than MaxZ.");
0174     }
0175     m_region = m_regionCfg;
0176     info("Recording within r < %.1f mm and %.1f mm < z < %.1f mm.",
0177          m_region->maxR / CLHEP::mm, m_region->minZ / CLHEP::mm,
0178          m_region->maxZ / CLHEP::mm);
0179   }
0180 }
0181
0182 void ActsMaterialRecorder::beginEvent(const G4Event* /*event*/) {
0183   // The DDG4 event takes ownership and deletes it after the event
0184   m_recorded =
0185       context()->event().addExtension(new detail::RecordedMaterialTracks());
0186 }
0187
0188 void ActsMaterialRecorder::operator()(const G4Step* step,
0189                                       G4SteppingManager* /*mgr*/) {
0190   // Every track gets an entry, also if it never crosses recorded material
0191   const G4Track& track = *step->GetTrack();
0192   auto& tracks = m_recorded->tracks;
0193   auto it = tracks.find(track.GetTrackID());
0194   if (it == tracks.end()) {
0195     it =
0196         tracks
0197             .try_emplace(track.GetTrackID(), m_stepConverter->startTrack(track))
0198             .first;
0199   }
0200
0201   if (!m_region.has_value()) {
0202     m_stepConverter->addStep(*step, it->second);
0203     return;
0204   }
0205
0206   const G4StepPoint& pre = *step->GetPreStepPoint();
0207   if (!m_region->contains(pre.GetPosition())) {
0208     step->GetTrack()->SetTrackStatus(fStopAndKill);
0209     return;
0210   }
0211   const double inside =
0212       m_region->pathToBoundary(pre.GetPosition(), pre.GetMomentumDirection());
0213   const bool leaves = step->GetStepLength() >= inside;
0214   if (m_stepConverter->addStep(*step, it->second) && leaves) {
0215     clipLastStep(it->second, inside * Acts::UnitConstants::mm / CLHEP::mm);
0216   }
0217   if (leaves) {
0218     // Nothing beyond the region is recorded, so there is no need to go on
0219     step->GetTrack()->SetTrackStatus(fStopAndKill);
0220   }
0221 }
0222
0223 }  // namespace ActsPlugins
0224
0225 DECLARE_GEANT4ACTION_NS(ActsPlugins, ActsMaterialRecorder)