Geant4HepEmTrackingPhysics.cpp

Go to the documentation of this file.

//==========================================================================
//  AIDA Detector description implementation
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
//==========================================================================
 
#ifndef DDG4_GEANT4HEPEMTRACKINGPHYSICS_H
#define DDG4_GEANT4HEPEMTRACKINGPHYSICS_H 1
 
#include <DDG4/Geant4PhysicsList.h>
 
#include <G4Electron.hh>
#include <G4EventManager.hh>
#include <G4Gamma.hh>
#include <G4Positron.hh>
#include <G4RegionStore.hh>
#include <G4TrackStatus.hh>
#include <G4VUserPhysicsList.hh>
 
#include <G4HepEmConfig.hh>
#include <G4HepEmTrackingManager.hh>
 
#include <memory>
#include <string>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace dd4hep {
 
  namespace sim {
 
    namespace {
 
      class Geant4RegionSelectiveHepEmTrackingManager : public G4VTrackingManager {
      public:
        Geant4RegionSelectiveHepEmTrackingManager(bool trackInAllRegions,
                                                  std::unordered_set<G4int> hepEmRegionIDs,
                                                  G4int verbose)
          : m_trackInAllRegions(trackInAllRegions)
          , m_hepEmRegionIDs(std::move(hepEmRegionIDs))
          , m_hepEmTrackingManager(new G4HepEmTrackingManager(verbose)) {}
 
        virtual ~Geant4RegionSelectiveHepEmTrackingManager() = default;
 
        G4HepEmConfig* GetConfig() { return m_hepEmTrackingManager->GetConfig(); }
 
        void BuildPhysicsTable(const G4ParticleDefinition& part) override {
          m_hepEmTrackingManager->BuildPhysicsTable(part);
        }
 
        void PreparePhysicsTable(const G4ParticleDefinition& part) override {
          m_hepEmTrackingManager->PreparePhysicsTable(part);
        }
 
        void FlushEvent() override {
          m_hepEmTrackingManager->FlushEvent();
        }
 
        void HandOverOneTrack(G4Track* aTrack) override {
          if (useHepEm(aTrack)) {
            m_hepEmTrackingManager->HandOverOneTrack(aTrack);
            return;
          }
          trackWithGeant4(aTrack);
        }
 
      private:
        bool useHepEm(const G4Track* aTrack) const {
          if (m_trackInAllRegions) {
            return true;
          }
          const G4VPhysicalVolume* volume = aTrack->GetVolume();
          if (volume == nullptr && aTrack->GetTouchableHandle()) {
            volume = aTrack->GetTouchableHandle()->GetVolume();
          }
          if (volume == nullptr) {
            return false;
          }
          const G4LogicalVolume* logical = volume->GetLogicalVolume();
          if (logical == nullptr) {
            return false;
          }
          const G4Region* region = logical->GetRegion();
          if (region == nullptr) {
            return false;
          }
          return m_hepEmRegionIDs.find(region->GetInstanceID()) != m_hepEmRegionIDs.end();
        }
 
        void trackWithGeant4(G4Track* track) {
          // We are inside a custom G4VTrackingManager::HandOverOneTrack() path here,
          // but with track hand-off to Geant4 after all (since not in a HepEm region).
          // In this path, G4EventManager does not execute its normal post-ProcessOneTrack()
          // switch (see G4EventManager::DoProcessing, switch(istop) in the non-custom branch),
          // so we must mirror that track-status/secondary/trajectory handling here.
          auto* eventManager = G4EventManager::GetEventManager();
          auto* trackManager = eventManager->GetTrackingManager();
          auto* stackManager = eventManager->GetStackManager();
 
          trackManager->ProcessOneTrack(track);
 
          G4TrackStatus track_status = track->GetTrackStatus();
          G4VTrajectory* trajectory =
              trackManager->GetStoreTrajectory() == 0 ? nullptr : trackManager->GimmeTrajectory();
          G4TrackVector* secondaries = trackManager->GimmeSecondaries();
 
          switch (track_status) {
          case fStopButAlive:
          case fSuspend:
          case fSuspendAndWait:
            stackManager->PushOneTrack(track, trajectory);
            trajectory = nullptr;
            eventManager->StackTracks(secondaries);
            break;
          case fPostponeToNextEvent:
            stackManager->PushOneTrack(track);
            eventManager->StackTracks(secondaries);
            delete trajectory;
            break;
          case fStopAndKill:
            eventManager->StackTracks(secondaries);
            delete trajectory;
            delete track;
            break;
          case fKillTrackAndSecondaries:
            if (secondaries != nullptr) {
              for (auto& secondary : *secondaries) {
                delete secondary;
              }
              secondaries->clear();
            }
            delete trajectory;
            delete track;
            break;
          default:
            G4cerr << "+++ Geant4RegionSelectiveHepEmTrackingManager: unexpected track status from "
                      "G4TrackingManager, killing track."
                   << G4endl;
            eventManager->StackTracks(secondaries);
            delete trajectory;
            delete track;
            break;
          }
        }
 
        bool m_trackInAllRegions;
        std::unordered_set<G4int> m_hepEmRegionIDs;
        std::unique_ptr<G4HepEmTrackingManager> m_hepEmTrackingManager;
      };
 
    }  // namespace
 
 
    class Geant4HepEmTrackingPhysics : public Geant4PhysicsList {
    protected:
      DDG4_DEFINE_ACTION_CONSTRUCTORS(Geant4HepEmTrackingPhysics);
 
      std::vector<std::string> m_woodcockRegions;
      std::vector<std::string> m_hepEmRegions;
      G4int m_verbosity;
 
    public:
      Geant4HepEmTrackingPhysics(Geant4Context* context, const std::string& nam)
        : Geant4PhysicsList(context, nam) {
        declareProperty("WoodcockRegions", m_woodcockRegions);
        declareProperty("HepEmRegions", m_hepEmRegions);
        declareProperty("VerboseLevel", m_verbosity = 0);
      }
 
      virtual ~Geant4HepEmTrackingPhysics() = default;
 
 
      virtual void constructProcesses(G4VUserPhysicsList* /* physics_list */) override {
        const bool trackInAllRegions = m_hepEmRegions.empty();
        std::unordered_set<G4int> hepEmRegionIDs;
        if (trackInAllRegions) {
          info("+++ HepEm region selection: all regions");
        } else {
          for (const auto& regionName : m_hepEmRegions) {
            G4Region* region = G4RegionStore::GetInstance()->GetRegion(regionName, false);
            if (region == nullptr) {
              warning("+++ HepEm region not found in G4RegionStore: %s", regionName.c_str());
              continue;
            }
            hepEmRegionIDs.insert(region->GetInstanceID());
            info("+++ HepEm enabled in G4Region: %s", regionName.c_str());
          }
          if (hepEmRegionIDs.empty()) {
            warning("+++ No valid HepEmRegions were resolved; HepEm tracking will be disabled.");
          }
        }
        const std::unordered_set<G4int> resolvedHepEmRegionIDs = hepEmRegionIDs;
 
        auto* tm =
            new Geant4RegionSelectiveHepEmTrackingManager(trackInAllRegions, std::move(hepEmRegionIDs), m_verbosity);
 
        G4int configuredWoodcockRegions = 0;
        for (const auto& region : m_woodcockRegions) {
          G4Region* woodcockRegion = G4RegionStore::GetInstance()->GetRegion(region, false);
          if (woodcockRegion == nullptr) {
            warning("+++ Woodcock region not found in G4RegionStore: %s", region.c_str());
            continue;
          }
          if (!trackInAllRegions
              && resolvedHepEmRegionIDs.find(woodcockRegion->GetInstanceID()) == resolvedHepEmRegionIDs.end()) {
            warning("+++ Ignoring Woodcock region '%s': region is outside HepEmRegions selection", region.c_str());
            continue;
          }
          info("+++ Enabling Woodcock photon tracking in G4Region: %s", region.c_str());
          tm->GetConfig()->SetWoodcockTrackingRegion(region);
          ++configuredWoodcockRegions;
        }
        if (!trackInAllRegions && !m_woodcockRegions.empty() && configuredWoodcockRegions == 0) {
          warning("+++ No WoodcockRegions remain after applying HepEmRegions intersection.");
        }
 
        // Warn if another custom tracking manager is already installed —
        // SetTrackingManager() is exclusive and will silently replace it.
        auto warnIfConflict = [&](G4ParticleDefinition* particle) {
          if (particle->GetTrackingManager() != nullptr) {
            warning("+++ Replacing existing tracking manager for %s with G4HepEmTrackingManager",
                    particle->GetParticleName().c_str());
          }
        };
        warnIfConflict(G4Electron::Definition());
        warnIfConflict(G4Positron::Definition());
        warnIfConflict(G4Gamma::Definition());
 
        G4Electron::Definition()->SetTrackingManager(tm);
        G4Positron::Definition()->SetTrackingManager(tm);
        G4Gamma::Definition()   ->SetTrackingManager(tm);
 
        info("+++ Installed G4HepEmTrackingManager for e-/e+/gamma");
      }
    };
  }  // namespace sim
}  // namespace dd4hep
 
#endif  // DDG4_GEANT4HEPEMTRACKINGPHYSICS_H
 
#include <DDG4/Factories.h>
using namespace dd4hep::sim;
DECLARE_GEANT4ACTION(Geant4HepEmTrackingPhysics)