MuonSensitiveDetector.cxx

Go to the documentation of this file.

 
/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#include "MuonSensitiveDetector.h"
 
#include "MuonSensitiveDetectorsR4/Utils.h"
#include "GeoPrimitives/CLHEPtoEigenConverter.h"
#include "GeoModelKernel/throwExcept.h"
#include "xAODMuonSimHit/MuonSimHitAuxContainer.h"
#include "StoreGate/ReadHandle.h"
 
#include "G4Geantino.hh"
#include "G4ChargedGeantino.hh"
 
#include "MCTruth/TrackHelper.h"
#include "MCTruth/TrackInformation.h"
#include "MCTruth/AtlasG4EventUserInfo.h"
#include "TruthUtils/HepMCHelpers.h"
 
using namespace ActsTrk;
 
namespace {
    static const SG::Decorator<int> dec_G4TrkId{"MuonSim_G4TrkId"};
}
 
namespace MuonG4R4 {
    MuonSensitiveDetector::MuonSensitiveDetector(const std::string& name, 
                                                 const std::string& output_key,
                                                 const std::string& trfStore_key,
                                                 const MuonGMR4::MuonDetectorManager* detMgr):
        G4VSensitiveDetector{name},
        AthMessaging{name},
        m_writeKey{output_key},
        m_trfCacheKey{trfStore_key},
        m_detMgr{detMgr} {
        m_trfCacheKey.initialize().ignore();
    }
    void MuonSensitiveDetector::Initialize(G4HCofThisEvent*) {
        m_g4UserEventInfo = AtlasG4EventUserInfo::GetEventUserInfo();
        if (m_g4UserEventInfo) {
            auto* hitVec = m_g4UserEventInfo->GetHitCollectionMap()->Find<MuonSimHitsVec>(m_writeKey);
            if (!hitVec) {
                THROW_EXCEPTION("The event does not contain a MuonSimHit container called '"<<m_writeKey<<"'");
            }
            m_outContainer = hitVec->container.get();
            ATH_MSG_DEBUG(__func__<<"() "<<__LINE__<<" - Retrieved '"<<m_writeKey<<"'.");
        } else {
            m_outContainer = nullptr;
        }
    }
    AtlasG4EventUserInfo* MuonSensitiveDetector::eventInfo() const {
        return m_g4UserEventInfo ? m_g4UserEventInfo : AtlasG4EventUserInfo::GetEventUserInfo();
    }
    const EventContext& MuonSensitiveDetector::eventContext() const {
        AtlasG4EventUserInfo* g4EventInfo = eventInfo();
        if (!g4EventInfo) {
            THROW_EXCEPTION("No AtlasG4EventUserInfo available");
        }
        return g4EventInfo->GetEventContext();
    }
    ActsTrk::GeometryContext MuonSensitiveDetector::getGeoContext() const {
        ActsTrk::GeometryContext gctx{};
        const ActsTrk::DetectorAlignStore* alignment{nullptr};
        if (!SG::get(alignment, m_trfCacheKey, eventContext()).isSuccess()) {
            THROW_EXCEPTION("Failed to retrieve "<<m_trfCacheKey.fullKey()<<".");
        }
        gctx.setStore(std::make_unique<DetectorAlignStore>(*alignment));
        return gctx;
    }
    bool MuonSensitiveDetector::processStep(const G4Step* aStep) const {
        const G4Track* currentTrack = aStep->GetTrack();
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Check whether step pdgId: "<<(*currentTrack)<<" will be processed.");
        constexpr double velCutOff = 10.*Gaudi::Units::micrometer / Gaudi::Units::second;
        if (aStep->GetStepLength() < std::numeric_limits<float>::epsilon() || currentTrack->GetVelocity() < velCutOff) {
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Step length is too short ");
            return false;
        }
        if (currentTrack->GetDefinition()->GetPDGCharge() == 0.0) {
            ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Particle is neutral");
            return currentTrack->GetDefinition() == G4Geantino::GeantinoDefinition() ||
                   currentTrack->GetDefinition() == G4ChargedGeantino::ChargedGeantinoDefinition();
        }
        return true;
    }
    HepMcParticleLink MuonSensitiveDetector::genParticleLink(const G4Track* track) const {
        TrackHelper trHelper{track};
        AtlasG4EventUserInfo* g4EventInfo = eventInfo();
        return trHelper.GenerateParticleLink(g4EventInfo ? g4EventInfo->GetEventStore() : nullptr);
    }
    xAOD::MuonSimHit* MuonSensitiveDetector::propagateAndSaveStrip(const Identifier& hitID,
                                                                   const Acts::Transform3& toGasGap,
                                                                   const G4Step* aStep) {
        
        const G4Track* currentTrack = aStep->GetTrack();
        const Amg::Vector3D locPostStep{toGasGap*Amg::Hep3VectorToEigen(aStep->GetPostStepPoint()->GetPosition())};
        Amg::Vector3D locPreStep{toGasGap*Amg::Hep3VectorToEigen(aStep->GetPreStepPoint()->GetPosition())};

        Amg::Vector3D locHitDir = toGasGap.linear() * Amg::Hep3VectorToEigen(currentTrack->GetMomentumDirection());

        xAOD::MuonSimHit* prevHit = lastSnapShot(hitID, aStep);
        if (MC::isElectron(currentTrack->GetParticleDefinition()->GetPDGEncoding()) && prevHit) {
            locPreStep = xAOD::toEigen(prevHit->localPosition()) 
                       - 0.5* prevHit->stepLength() * xAOD::toEigen(prevHit->localDirection());
 
            locHitDir = (locPostStep - locPreStep).unit();
        }
        const Amg::Vector3D locHitPos = 0.5* (locPreStep + locPostStep);
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - "<<m_detMgr->idHelperSvc()->toStringGasGap(hitID)
            <<" - track: "<<(*currentTrack) <<", deposit: "<<aStep->GetTotalEnergyDeposit()
            <<", -- local coords -- prestep: "<<Amg::toString(locPreStep)
            <<",  post step: "<<Amg::toString(locPostStep) 
            <<" mid point: "<< Amg::toString(locHitPos)<<", direction: "<<Amg::toString(locHitDir)
            <<", deposit: "<<aStep->GetTotalEnergyDeposit());
 
        const double globalTime = currentTrack->GetGlobalTime() + locHitDir.dot(locPostStep - locHitPos) / currentTrack->GetVelocity();
  
        xAOD::MuonSimHit* newHit = saveHit(hitID, locHitPos, locHitDir, globalTime, aStep);
        if (newHit && prevHit) {
            newHit->setStepLength((locPostStep - locPreStep).mag());
        }
        return newHit;
    }
    xAOD::MuonSimHit* MuonSensitiveDetector::saveHit(const Identifier& hitId,
                                                     const Amg::Vector3D& hitPos,
                                                     const Amg::Vector3D& hitDir,
                                                     const double globTime,
                                                     const G4Step* aStep) {
        const G4Track* currentTrack = aStep->GetTrack();
        const HepMcParticleLink particleLink{genParticleLink(currentTrack)};
        // If Geant4 propagates the same track through the same volume just update the last hit and don't write a new one
        xAOD::MuonSimHit* hit = lastSnapShot(hitId, aStep, particleLink);
        bool newHit{false};
        if (!hit) {
            if (!m_outContainer) {
                THROW_EXCEPTION("Cannot find container '"<<m_writeKey<<"'.");
            }
            hit = m_outContainer->push_back(std::make_unique<xAOD::MuonSimHit>());
            newHit = true;
        }
        dec_G4TrkId(*hit) = currentTrack->GetTrackID();
        hit->setIdentifier(hitId); 
        hit->setLocalPosition(xAOD::toStorage(hitPos));  
        hit->setLocalDirection(xAOD::toStorage(hitDir));
        hit->setMass(currentTrack->GetDefinition()->GetPDGMass());
        hit->setGlobalTime(globTime);
        int pdgId = currentTrack->GetDefinition()->GetPDGEncoding();
        if (pdgId ==0) {
            if (currentTrack->GetDefinition() == G4ChargedGeantino::ChargedGeantinoDefinition()){
              pdgId = std::copysign(MC::GEANTINOPLUS, currentTrack->GetDynamicParticle()->GetCharge());
            } else if (currentTrack->GetDefinition() == G4Geantino::GeantinoDefinition()) {
                pdgId = MC::GEANTINO0;
            } else {
                ATH_MSG_WARNING(__func__<<"() "<<__LINE__<<" - The track has an undefined pdgId "<<aStep);
                m_outContainer->pop_back();
                return nullptr;
            }
        }
        hit->setPdgId(pdgId);
        hit->setEnergyDeposit(aStep->GetTotalEnergyDeposit() + (newHit ? 0. : hit->energyDeposit()));
        hit->setKineticEnergy(currentTrack->GetKineticEnergy());
        hit->setGenParticleLink(particleLink);
        hit->setStepLength(aStep->GetStepLength());
        
        ATH_MSG_VERBOSE(__func__<<"() "<<__LINE__<<" - Save new hit "<<m_detMgr->idHelperSvc()->toString(hitId)
                        <<", pdgId: "<<hit->pdgId()<<", "<<particleLink
                        <<", trackId: "<<currentTrack->GetTrackID()
                        <<", "<<particleLink.cptr()<<std::endl
                        <<"pos: "<<Amg::toString(hitPos)<<", dir: "<<Amg::toString(hitDir)<<", time: "<<globTime
                        <<", energy: "<<hit->kineticEnergy()<<", stepLength: "<<hit->stepLength()<<", "
                        <<", deposit energy: "<<hit->energyDeposit());
        return hit;
    }
    xAOD::MuonSimHit* MuonSensitiveDetector::lastSnapShot(const Identifier& hitId,
                                                          const G4Step* hitStep) {
        const HepMcParticleLink particleLink{genParticleLink(hitStep->GetTrack())};
        return lastSnapShot(hitId, hitStep, particleLink);
    }
    xAOD::MuonSimHit* MuonSensitiveDetector::lastSnapShot(const Identifier& hitId,
                                                          const G4Step* hitStep,
                                                          const HepMcParticleLink& particleLink) {
        if (m_outContainer->empty() || 
            m_outContainer->back()->identify() != hitId ||
            particleLink != m_outContainer->back()->genParticleLink() ||
            dec_G4TrkId(*m_outContainer->back()) != hitStep->GetTrack()->GetTrackID()) {
            return nullptr;
        }
        return m_outContainer->back();
    }
}