PhysicsLists/src/TrackingManagerHelper.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // TrackingManagerHelper
0027 //
0028 // Class description:
0029 //
0030 // Helper class for reducing the effort required to implement a custom tracking
0031 // manager. It implements a stepping loop that calls user actions as the generic
0032 // tracking and stepping managers do, and it implements navigation for charged
0033 // particles in energy-preserving fields and for neutral particles.
0034 //
0035 // Original author: Jonas Hahnfeld, 2021
0036
0037 #include "G4EventManager.hh"
0038 #include "G4Step.hh"
0039 #include "G4StepPoint.hh"
0040 #include "G4Track.hh"
0041 #include "G4TrackVector.hh"
0042 #include "G4UserSteppingAction.hh"
0043 #include "G4UserTrackingAction.hh"
0044 #include "G4VSensitiveDetector.hh"
0045
0046 #include "G4Field.hh"
0047 #include "G4FieldManager.hh"
0048 #include "G4FieldManagerStore.hh"
0049 #include "G4GeometryTolerance.hh"
0050 #include "G4LogicalVolume.hh"
0051 #include "G4Navigator.hh"
0052 #include "G4PropagatorInField.hh"
0053 #include "G4Region.hh"
0054 #include "G4SafetyHelper.hh"
0055 #include "G4TouchableHandle.hh"
0056 #include "G4TouchableHistory.hh"
0057 #include "G4TransportationManager.hh"
0058 #include "G4VPhysicalVolume.hh"
0059
0060 template <typename PhysicsImpl, typename NavigationImpl>
0061 void TrackingManagerHelper::TrackParticle(G4Track* aTrack, PhysicsImpl& physics, NavigationImpl& navigation) {
0062   // Prepare for calling the user action.
0063   auto* evtMgr = G4EventManager::GetEventManager();
0064   auto* userTrackingAction = evtMgr->GetUserTrackingAction();
0065   auto* userSteppingAction = evtMgr->GetUserSteppingAction();
0066
0067   // Locate the track in geometry.
0068   {
0069     auto* transMgr = G4TransportationManager::GetTransportationManager();
0070     auto* linearNavigator = transMgr->GetNavigatorForTracking();
0071
0072     const G4ThreeVector& pos = aTrack->GetPosition();
0073     const G4ThreeVector& dir = aTrack->GetMomentumDirection();
0074
0075     // Do not assign directly, doesn't work if the handle is empty.
0076     G4TouchableHandle touchableHandle;
0077     if (aTrack->GetTouchableHandle()) {
0078       touchableHandle = aTrack->GetTouchableHandle();
0079       // FIXME: This assumes we only ever have G4TouchableHistorys!
0080       auto* touchableHistory = (G4TouchableHistory*)touchableHandle();
0081       G4VPhysicalVolume* oldTopVolume = touchableHandle->GetVolume();
0082       G4VPhysicalVolume* newTopVolume = linearNavigator->ResetHierarchyAndLocate(pos, dir, *touchableHistory);
0083       // TODO: WHY?!
0084       if (newTopVolume != oldTopVolume || oldTopVolume->GetRegularStructureId() == 1) {
0085         touchableHandle = linearNavigator->CreateTouchableHistory();
0086         aTrack->SetTouchableHandle(touchableHandle);
0087       }
0088     } else {
0089       linearNavigator->LocateGlobalPointAndSetup(pos, &dir, false, false);
0090       touchableHandle = linearNavigator->CreateTouchableHistory();
0091       aTrack->SetTouchableHandle(touchableHandle);
0092     }
0093     aTrack->SetNextTouchableHandle(touchableHandle);
0094   }
0095
0096   // Prepare data structures used while tracking.
0097   G4Step step;
0098   step.NewSecondaryVector();
0099   G4StepPoint& preStepPoint = *step.GetPreStepPoint();
0100   step.InitializeStep(aTrack);
0101   aTrack->SetStep(&step);
0102   G4TrackVector secondaries;
0103
0104   // Start of tracking: Inform user and processes.
0105   if (userTrackingAction) {
0106     userTrackingAction->PreUserTrackingAction(aTrack);
0107   }
0108
0109   physics.StartTracking(aTrack);
0110
0111   while (aTrack->GetTrackStatus() == fAlive) {
0112     // Beginning of this step: Prepare data structures.
0113     aTrack->IncrementCurrentStepNumber();
0114
0115     step.CopyPostToPreStepPoint();
0116     step.ResetTotalEnergyDeposit();
0117     aTrack->SetTouchableHandle(aTrack->GetNextTouchableHandle());
0118
0119     auto* lvol = aTrack->GetTouchable()->GetVolume()->GetLogicalVolume();
0120     preStepPoint.SetMaterial(lvol->GetMaterial());
0121     preStepPoint.SetMaterialCutsCouple(lvol->GetMaterialCutsCouple());
0122
0123     // Query step lengths from pyhsics and geometry, decide on limit.
0124     G4double physicalStep = physics.GetPhysicalInteractionLength(*aTrack);
0125     G4double geometryStep = navigation.MakeStep(*aTrack, step, physicalStep);
0126
0127     bool geometryLimitedStep = geometryStep < physicalStep;
0128     G4double finalStep = geometryLimitedStep ? geometryStep : physicalStep;
0129
0130     step.SetStepLength(finalStep);
0131     aTrack->SetStepLength(finalStep);
0132
0133     // Call AlongStepDoIt in every step.
0134     physics.AlongStepDoIt(*aTrack, step, secondaries);
0135     step.UpdateTrack();
0136
0137     if (aTrack->GetTrackStatus() == fAlive && aTrack->GetKineticEnergy() < DBL_MIN) {
0138       if (physics.HasAtRestProcesses()) {
0139         aTrack->SetTrackStatus(fStopButAlive);
0140       } else {
0141         aTrack->SetTrackStatus(fStopAndKill);
0142       }
0143     }
0144
0145     navigation.FinishStep(*aTrack, step);
0146
0147     // Check if the track left the world.
0148     if (aTrack->GetNextVolume() == nullptr) {
0149       aTrack->SetTrackStatus(fStopAndKill);
0150     }
0151
0152     // The check should rather check for == fAlive and avoid calling
0153     // PostStepDoIt for fStopButAlive, but the generic stepping loop
0154     // does it like this...
0155     if (aTrack->GetTrackStatus() != fStopAndKill) {
0156       physics.PostStepDoIt(*aTrack, step, secondaries);
0157     }
0158
0159     // Need to get the true step length, not the geometry step length!
0160     aTrack->AddTrackLength(step.GetStepLength());
0161
0162     // End of this step: Call sensitive detector and stepping actions.
0163     if (step.GetControlFlag() != AvoidHitInvocation) {
0164       auto* sensitive = lvol->GetSensitiveDetector();
0165       if (sensitive) {
0166         sensitive->Hit(&step);
0167       }
0168     }
0169
0170     if (userSteppingAction) {
0171       userSteppingAction->UserSteppingAction(&step);
0172     }
0173
0174     auto* regionalAction = lvol->GetRegion()->GetRegionalSteppingAction();
0175     if (regionalAction) {
0176       regionalAction->UserSteppingAction(&step);
0177     }
0178   }
0179
0180   if (aTrack->GetTrackStatus() == fStopButAlive && aTrack->GetNextVolume() != nullptr) {
0181     // Do one final step.
0182     aTrack->IncrementCurrentStepNumber();
0183
0184     step.CopyPostToPreStepPoint();
0185     step.ResetTotalEnergyDeposit();
0186
0187     physics.AtRestDoIt(*aTrack, step, secondaries);
0188
0189     // End of this step: Call sensitive detector and stepping actions.
0190     auto* lvol = aTrack->GetTouchable()->GetVolume()->GetLogicalVolume();
0191     if (step.GetControlFlag() != AvoidHitInvocation) {
0192       auto sensitive = lvol->GetSensitiveDetector();
0193       if (sensitive) {
0194         sensitive->Hit(&step);
0195       }
0196     }
0197
0198     if (userSteppingAction) {
0199       userSteppingAction->UserSteppingAction(&step);
0200     }
0201
0202     auto* regionalAction = lvol->GetRegion()->GetRegionalSteppingAction();
0203     if (regionalAction) {
0204       regionalAction->UserSteppingAction(&step);
0205     }
0206   }
0207
0208   // End of tracking: Inform processes and user.
0209   physics.EndTracking();
0210
0211   if (userTrackingAction) {
0212     userTrackingAction->PostUserTrackingAction(aTrack);
0213   }
0214
0215   evtMgr->StackTracks(&secondaries);
0216
0217   step.DeleteSecondaryVector();
0218 }
0219
0220 template <typename PhysicsImpl>
0221 void TrackingManagerHelper::TrackChargedParticle(G4Track* aTrack, PhysicsImpl& physics) {
0222   class ChargedNavigation final : public Navigation {
0223   public:
0224     ChargedNavigation() {
0225       auto* transMgr = G4TransportationManager::GetTransportationManager();
0226       fLinearNavigator = transMgr->GetNavigatorForTracking();
0227       fFieldPropagator = transMgr->GetPropagatorInField();
0228       fSafetyHelper = transMgr->GetSafetyHelper();
0229       kCarTolerance = 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
0230
0231       // Reset sstate of field propagator and all chord finders.
0232       fFieldPropagator->ClearPropagatorState();
0233
0234       auto* fieldMgrStore = G4FieldManagerStore::GetInstance();
0235       fieldMgrStore->ClearAllChordFindersState();
0236     }
0237
0238     G4double MakeStep(G4Track& track, G4Step& step, G4double physicalStep) override {
0239       G4ThreeVector pos = track.GetPosition();
0240       G4ThreeVector dir = track.GetMomentumDirection();
0241       G4StepPoint& postStepPoint = *step.GetPostStepPoint();
0242
0243       bool fieldExertsForce = false;
0244       if (auto* fieldMgr = fFieldPropagator->FindAndSetFieldManager(track.GetVolume())) {
0245         fieldMgr->ConfigureForTrack(&track);
0246         if (const G4Field* ptrField = fieldMgr->GetDetectorField()) {
0247           fieldExertsForce = true;
0248         }
0249       }
0250
0251       G4double endpointDistance;
0252       G4double safety = 0.0;
0253       // Setting a fallback value for safety is required in case of where very
0254       // short steps where the field propagator returns immediately without
0255       // calling geometry.
0256       const G4double shiftSquare = (pos - fSafetyOrigin).mag2();
0257       if (shiftSquare < sqr(fSafety)) {
0258         safety = fSafety - std::sqrt(shiftSquare);
0259       }
0260
0261       if (fieldExertsForce) {
0262         const G4DynamicParticle* pParticle = track.GetDynamicParticle();
0263         const G4double particleCharge = pParticle->GetCharge();
0264         const G4double particleMass = pParticle->GetMass();
0265         const G4double magneticMoment = pParticle->GetMagneticMoment();
0266         const G4ThreeVector particleSpin = pParticle->GetPolarization();
0267         const G4double kineticEnergy = pParticle->GetKineticEnergy();
0268         const auto pParticleDef = pParticle->GetDefinition();
0269         const auto particlePDGSpin = pParticleDef->GetPDGSpin();
0270         const auto particlePDGMagM = pParticleDef->GetPDGMagneticMoment();
0271
0272         auto equationOfMotion = fFieldPropagator->GetCurrentEquationOfMotion();
0273         equationOfMotion->SetChargeMomentumMass(G4ChargeState(particleCharge, magneticMoment, particlePDGSpin),
0274                                                 pParticle->GetTotalMomentum(),
0275                                                 particleMass);
0276
0277         const G4ThreeVector startPosition = pos;
0278         const G4ThreeVector startDirection = dir;
0279         G4FieldTrack aFieldTrack(startPosition,
0280                                  track.GetGlobalTime(),  // Lab.
0281                                  dir,
0282                                  kineticEnergy,
0283                                  particleMass,
0284                                  particleCharge,
0285                                  particleSpin,
0286                                  particlePDGMagM,
0287                                  0.0,  // Length along track
0288                                  particlePDGSpin);
0289
0290         // Do the Transport in the field (non recti-linear)
0291         //
0292         fGeometryLimitedStep = false;
0293         const G4double lengthAlongCurve =
0294             fFieldPropagator->ComputeStep(aFieldTrack, physicalStep, safety, track.GetVolume(), kineticEnergy < 250.0);
0295         if (lengthAlongCurve < physicalStep) {
0296           physicalStep = lengthAlongCurve;
0297           fGeometryLimitedStep = true;
0298         }
0299         fSafetyHelper->SetCurrentSafety(safety, pos);
0300         fSafetyOrigin = pos;
0301         fSafety = safety;
0302
0303         if (fFieldPropagator->IsParticleLooping()) {
0304           track.SetTrackStatus(fStopAndKill);
0305         }
0306
0307         pos = aFieldTrack.GetPosition();
0308         dir = aFieldTrack.GetMomentumDir();
0309
0310         postStepPoint.SetPosition(pos);
0311         postStepPoint.SetMomentumDirection(dir);
0312
0313         endpointDistance = (startPosition - pos).mag();
0314       } else {
0315         fGeometryLimitedStep = false;
0316         G4double linearStepLength = fLinearNavigator->ComputeStep(pos, dir, physicalStep, safety);
0317         if (linearStepLength < physicalStep) {
0318           physicalStep = linearStepLength;
0319           fGeometryLimitedStep = true;
0320         }
0321         fSafetyHelper->SetCurrentSafety(safety, pos);
0322         fSafetyOrigin = pos;
0323         fSafety = safety;
0324
0325         // Update the position.
0326         pos += physicalStep * dir;
0327         postStepPoint.SetPosition(pos);
0328
0329         endpointDistance = physicalStep;
0330       }
0331
0332       // Update global, local, and proper time.
0333       double velocity = track.GetVelocity();
0334       double deltaTime = 0;
0335       if (velocity > 0) {
0336         deltaTime = physicalStep / velocity;
0337       }
0338
0339       postStepPoint.AddGlobalTime(deltaTime);
0340       postStepPoint.AddLocalTime(deltaTime);
0341
0342       double restMass = track.GetDynamicParticle()->GetMass();
0343       double deltaProperTime = deltaTime * (restMass / track.GetTotalEnergy());
0344       postStepPoint.AddProperTime(deltaProperTime);
0345
0346       // Compute safety, including the call to safetyHelper, but don't set the
0347       // safety in the post-step point to mimick the generic stepping loop.
0348       if (safety > physicalStep) {
0349         safety -= physicalStep;
0350       } else if (safety < endpointDistance) {
0351         safety = fLinearNavigator->ComputeSafety(pos);
0352         fSafetyHelper->SetCurrentSafety(safety, pos);
0353         fSafetyOrigin = pos;
0354         fSafety = safety;
0355       } else {
0356         safety = 0;
0357       }
0358       if (safety < kCarTolerance) {
0359         fPostStepSafety = kCarTolerance;
0360       } else {
0361         fPostStepSafety = safety;
0362       }
0363
0364       return physicalStep;
0365     }
0366
0367     void FinishStep(G4Track& track, G4Step& step) override {
0368       // Now set the safety that was computed in MakeStep.
0369       G4StepPoint& postStepPoint = *step.GetPostStepPoint();
0370       postStepPoint.SetSafety(fPostStepSafety);
0371
0372       G4TouchableHandle touchableHandle = track.GetTouchableHandle();
0373       const G4ThreeVector& pos = track.GetPosition();
0374       if (fGeometryLimitedStep) {
0375         // Relocate the particle.
0376         fLinearNavigator->SetGeometricallyLimitedStep();
0377         fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0378             pos, track.GetMomentumDirection(), touchableHandle, true);
0379         const G4VPhysicalVolume* newVolume = touchableHandle->GetVolume();
0380         if (newVolume == nullptr) {
0381           postStepPoint.SetStepStatus(fWorldBoundary);
0382         } else {
0383           postStepPoint.SetStepStatus(fGeomBoundary);
0384         }
0385       } else {
0386         // Move the Navigator's location.
0387         fLinearNavigator->LocateGlobalPointWithinVolume(pos);
0388       }
0389
0390       postStepPoint.SetTouchableHandle(touchableHandle);
0391       track.SetNextTouchableHandle(touchableHandle);
0392     }
0393
0394   private:
0395     G4Navigator* fLinearNavigator;
0396     G4PropagatorInField* fFieldPropagator;
0397     G4SafetyHelper* fSafetyHelper;
0398     G4ThreeVector fSafetyOrigin;
0399     G4double fSafety = 0;
0400     G4double fPostStepSafety = 0;
0401     G4double kCarTolerance;
0402     G4bool fGeometryLimitedStep;
0403   };
0404
0405   ChargedNavigation navigation;
0406   TrackParticle(aTrack, physics, navigation);
0407 }
0408
0409 template <typename PhysicsImpl>
0410 void TrackingManagerHelper::TrackNeutralParticle(G4Track* aTrack, PhysicsImpl& physics) {
0411   class NeutralNavigation final : public Navigation {
0412   public:
0413     NeutralNavigation() {
0414       auto* transMgr = G4TransportationManager::GetTransportationManager();
0415       fLinearNavigator = transMgr->GetNavigatorForTracking();
0416       fSafetyHelper = transMgr->GetSafetyHelper();
0417       kCarTolerance = 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
0418     }
0419
0420     G4double MakeStep(G4Track& track, G4Step& step, G4double physicalStep) override {
0421       G4ThreeVector pos = track.GetPosition();
0422       const G4ThreeVector& dir = track.GetMomentumDirection();
0423       G4StepPoint& postStepPoint = *step.GetPostStepPoint();
0424
0425       G4double safety = 0.0;
0426       const G4double shiftSquare = (pos - fSafetyOrigin).mag2();
0427       if (shiftSquare < sqr(fSafety)) {
0428         safety = fSafety - std::sqrt(shiftSquare);
0429       }
0430
0431       fGeometryLimitedStep = false;
0432       G4double linearStepLength = fLinearNavigator->ComputeStep(pos, dir, physicalStep, safety);
0433       if (linearStepLength < physicalStep) {
0434         physicalStep = linearStepLength;
0435         fGeometryLimitedStep = true;
0436       }
0437       fSafetyHelper->SetCurrentSafety(safety, pos);
0438       fSafetyOrigin = pos;
0439       fSafety = safety;
0440
0441       // Update the position.
0442       pos += physicalStep * dir;
0443       postStepPoint.SetPosition(pos);
0444
0445       // Update global, local, and proper time.
0446       double velocity = track.GetVelocity();
0447       double deltaTime = 0;
0448       if (velocity > 0) {
0449         deltaTime = physicalStep / velocity;
0450       }
0451       postStepPoint.AddGlobalTime(deltaTime);
0452       postStepPoint.AddLocalTime(deltaTime);
0453
0454       double restMass = track.GetDynamicParticle()->GetMass();
0455       double deltaProperTime = deltaTime * (restMass / track.GetTotalEnergy());
0456       postStepPoint.AddProperTime(deltaProperTime);
0457
0458       // Compute safety, but don't set the safety in the post-step point to
0459       // mimick the generic stepping loop.
0460       if (safety > physicalStep) {
0461         safety -= physicalStep;
0462       } else {
0463         safety = 0;
0464       }
0465       if (safety < kCarTolerance) {
0466         fPostStepSafety = kCarTolerance;
0467       } else {
0468         fPostStepSafety = safety;
0469       }
0470
0471       return physicalStep;
0472     }
0473
0474     void FinishStep(G4Track& track, G4Step& step) override {
0475       // Now set the safety that was computed in MakeStep.
0476       G4StepPoint& postStepPoint = *step.GetPostStepPoint();
0477       postStepPoint.SetSafety(fPostStepSafety);
0478
0479       G4TouchableHandle touchableHandle = track.GetTouchableHandle();
0480       const G4ThreeVector& pos = track.GetPosition();
0481       if (fGeometryLimitedStep) {
0482         // Relocate the particle.
0483         fLinearNavigator->SetGeometricallyLimitedStep();
0484         fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0485             pos, track.GetMomentumDirection(), touchableHandle, true);
0486         const G4VPhysicalVolume* newVolume = touchableHandle->GetVolume();
0487         if (newVolume == nullptr) {
0488           postStepPoint.SetStepStatus(fWorldBoundary);
0489         } else {
0490           postStepPoint.SetStepStatus(fGeomBoundary);
0491         }
0492       } else {
0493         // Move the Navigator's location.
0494         fLinearNavigator->LocateGlobalPointWithinVolume(pos);
0495       }
0496
0497       postStepPoint.SetTouchableHandle(touchableHandle);
0498       track.SetNextTouchableHandle(touchableHandle);
0499     }
0500
0501   private:
0502     G4Navigator* fLinearNavigator;
0503     G4SafetyHelper* fSafetyHelper;
0504     G4ThreeVector fSafetyOrigin;
0505     G4double fSafety = 0;
0506     G4double fPostStepSafety = 0;
0507     G4double kCarTolerance;
0508     G4bool fGeometryLimitedStep;
0509   };
0510
0511   NeutralNavigation navigation;
0512   TrackParticle(aTrack, physics, navigation);
0513 }