File indexing completed on 2025-01-31 23:36:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
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
0063 auto* evtMgr = G4EventManager::GetEventManager();
0064 auto* userTrackingAction = evtMgr->GetUserTrackingAction();
0065 auto* userSteppingAction = evtMgr->GetUserSteppingAction();
0066
0067
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
0076 G4TouchableHandle touchableHandle;
0077 if (aTrack->GetTouchableHandle()) {
0078 touchableHandle = aTrack->GetTouchableHandle();
0079
0080 auto* touchableHistory = (G4TouchableHistory*)touchableHandle();
0081 G4VPhysicalVolume* oldTopVolume = touchableHandle->GetVolume();
0082 G4VPhysicalVolume* newTopVolume = linearNavigator->ResetHierarchyAndLocate(pos, dir, *touchableHistory);
0083
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
0097 G4Step step;
0098 step.NewSecondaryVector();
0099 G4StepPoint& preStepPoint = *step.GetPreStepPoint();
0100 step.InitializeStep(aTrack);
0101 aTrack->SetStep(&step);
0102 G4TrackVector secondaries;
0103
0104
0105 if (userTrackingAction) {
0106 userTrackingAction->PreUserTrackingAction(aTrack);
0107 }
0108
0109 physics.StartTracking(aTrack);
0110
0111 while (aTrack->GetTrackStatus() == fAlive) {
0112
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
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
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
0148 if (aTrack->GetNextVolume() == nullptr) {
0149 aTrack->SetTrackStatus(fStopAndKill);
0150 }
0151
0152
0153
0154
0155 if (aTrack->GetTrackStatus() != fStopAndKill) {
0156 physics.PostStepDoIt(*aTrack, step, secondaries);
0157 }
0158
0159
0160 aTrack->AddTrackLength(step.GetStepLength());
0161
0162
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
0182 aTrack->IncrementCurrentStepNumber();
0183
0184 step.CopyPostToPreStepPoint();
0185 step.ResetTotalEnergyDeposit();
0186
0187 physics.AtRestDoIt(*aTrack, step, secondaries);
0188
0189
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
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
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 (fieldMgr->GetDetectorField()) {
0247 fieldExertsForce = true;
0248 }
0249 }
0250
0251 G4double endpointDistance;
0252 G4double safety = 0.0;
0253
0254
0255
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(),
0281 dir,
0282 kineticEnergy,
0283 particleMass,
0284 particleCharge,
0285 particleSpin,
0286 particlePDGMagM,
0287 0.0,
0288 particlePDGSpin);
0289
0290
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
0326 pos += physicalStep * dir;
0327 postStepPoint.SetPosition(pos);
0328
0329 endpointDistance = physicalStep;
0330 }
0331
0332
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
0347
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
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
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
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
0442 pos += physicalStep * dir;
0443 postStepPoint.SetPosition(pos);
0444
0445
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
0459
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
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
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
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 }