File indexing completed on 2026-05-02 07:42:28
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
0038
0039
0040
0041 #include "G4MonopoleTransportation.hh"
0042
0043 #include "DetectorConstruction.hh"
0044
0045 #include "G4ChordFinder.hh"
0046 #include "G4FieldManagerStore.hh"
0047 #include "G4Monopole.hh"
0048 #include "G4ParticleTable.hh"
0049 #include "G4ProductionCutsTable.hh"
0050 #include "G4RunManager.hh"
0051 #include "G4SafetyHelper.hh"
0052 #include "G4SystemOfUnits.hh"
0053 #include "G4TransportationProcessType.hh"
0054
0055 class G4VSensitiveDetector;
0056
0057
0058
0059 G4MonopoleTransportation::G4MonopoleTransportation(const G4Monopole* mpl, G4int verb)
0060 : G4VProcess(G4String("MonopoleTransportation"), fTransportation),
0061 fParticleDef(mpl),
0062 fMagSetup(0),
0063 fLinearNavigator(0),
0064 fFieldPropagator(0),
0065 fParticleIsLooping(false),
0066 fPreviousSftOrigin(0., 0., 0.),
0067 fPreviousSafety(0.0),
0068 fThreshold_Warning_Energy(100 * MeV),
0069 fThreshold_Important_Energy(250 * MeV),
0070 fThresholdTrials(10),
0071
0072 fNoLooperTrials(0),
0073 fSumEnergyKilled(0.0),
0074 fMaxEnergyKilled(0.0),
0075 fShortStepOptimisation(false),
0076 fpSafetyHelper(0),
0077 noCalls(0)
0078 {
0079 verboseLevel = verb;
0080
0081
0082 SetProcessSubType(TRANSPORTATION);
0083
0084
0085 if (G4Threading::IsMasterThread() && G4Threading::IsMultithreadedApplication()) {
0086 return;
0087 }
0088
0089 const DetectorConstruction* detector = static_cast<const DetectorConstruction*>(
0090 G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0091
0092 fMagSetup = detector->GetMonopoleFieldSetup();
0093
0094 G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
0095
0096 fLinearNavigator = transportMgr->GetNavigatorForTracking();
0097
0098 fFieldPropagator = transportMgr->GetPropagatorInField();
0099 fpSafetyHelper = transportMgr->GetSafetyHelper();
0100
0101
0102
0103
0104
0105
0106
0107
0108 fCurrentTouchableHandle = nullptr;
0109
0110 fEndGlobalTimeComputed = false;
0111 fCandidateEndGlobalTime = 0;
0112 }
0113
0114
0115
0116 G4MonopoleTransportation::~G4MonopoleTransportation()
0117 {
0118 if ((verboseLevel > 0) && (fSumEnergyKilled > 0.0)) {
0119 G4cout << " G4MonopoleTransportation: Statistics for looping particles " << G4endl;
0120 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
0121 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
0122 }
0123 }
0124
0125
0126
0127
0128
0129
0130
0131
0132 G4double G4MonopoleTransportation::AlongStepGetPhysicalInteractionLength(
0133 const G4Track& track,
0134 G4double,
0135 G4double currentMinimumStep, G4double& currentSafety, G4GPILSelection* selection)
0136 {
0137 fMagSetup->SetStepperAndChordFinder(1);
0138
0139
0140 G4double geometryStepLength, newSafety;
0141 fParticleIsLooping = false;
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152 *selection = CandidateForSelection;
0153
0154
0155
0156 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
0157 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection();
0158 G4ThreeVector startPosition = track.GetPosition();
0159
0160
0161
0162
0163
0164
0165
0166 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin;
0167 G4double MagSqShift = OriginShift.mag2();
0168 if (MagSqShift >= sqr(fPreviousSafety)) {
0169 currentSafety = 0.0;
0170 }
0171 else {
0172 currentSafety = fPreviousSafety - std::sqrt(MagSqShift);
0173 }
0174
0175
0176
0177 G4double particleMagneticCharge = fParticleDef->MagneticCharge();
0178 G4double particleElectricCharge = pParticle->GetCharge();
0179
0180 fGeometryLimitedStep = false;
0181
0182
0183
0184
0185
0186
0187
0188
0189 G4FieldManager* fieldMgr = 0;
0190 G4bool fieldExertsForce = false;
0191
0192 if ((particleMagneticCharge != 0.0)) {
0193 fieldMgr = fFieldPropagator->FindAndSetFieldManager(track.GetVolume());
0194 if (fieldMgr != 0) {
0195
0196 fieldMgr->ConfigureForTrack(&track);
0197
0198
0199
0200
0201
0202 fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
0203 }
0204 }
0205
0206
0207
0208
0209
0210
0211 if (!fieldExertsForce) {
0212 G4double linearStepLength;
0213 if (fShortStepOptimisation && (currentMinimumStep <= currentSafety)) {
0214
0215
0216 geometryStepLength = currentMinimumStep;
0217 fGeometryLimitedStep = false;
0218 }
0219 else {
0220
0221
0222 linearStepLength = fLinearNavigator->ComputeStep(startPosition, startMomentumDir,
0223 currentMinimumStep, newSafety);
0224
0225
0226 fPreviousSftOrigin = startPosition;
0227 fPreviousSafety = newSafety;
0228
0229
0230
0231
0232 currentSafety = newSafety;
0233
0234 fGeometryLimitedStep = (linearStepLength <= currentMinimumStep);
0235 if (fGeometryLimitedStep) {
0236
0237 geometryStepLength = linearStepLength;
0238 }
0239 else {
0240
0241 geometryStepLength = currentMinimumStep;
0242 }
0243 }
0244 endpointDistance = geometryStepLength;
0245
0246
0247
0248 fTransportEndPosition = startPosition + geometryStepLength * startMomentumDir;
0249
0250
0251
0252 fTransportEndMomentumDir = startMomentumDir;
0253 fTransportEndKineticEnergy = track.GetKineticEnergy();
0254 fTransportEndSpin = track.GetPolarization();
0255 fParticleIsLooping = false;
0256 fMomentumChanged = false;
0257 fEndGlobalTimeComputed = false;
0258 }
0259 else
0260 {
0261 G4double momentumMagnitude = pParticle->GetTotalMomentum();
0262 G4ThreeVector EndUnitMomentum;
0263 G4double lengthAlongCurve;
0264 G4double restMass = fParticleDef->GetPDGMass();
0265
0266 G4ChargeState chargeState(particleElectricCharge,
0267 fParticleDef->GetPDGSpin(),
0268 0,
0269 0,
0270 particleMagneticCharge);
0271
0272 G4EquationOfMotion* equationOfMotion =
0273 fFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetEquationOfMotion();
0274
0275 equationOfMotion->SetChargeMomentumMass(chargeState, momentumMagnitude, restMass);
0276
0277
0278
0279 G4ThreeVector spin = track.GetPolarization();
0280 G4FieldTrack aFieldTrack = G4FieldTrack(startPosition, track.GetMomentumDirection(), 0.0,
0281 track.GetKineticEnergy(), restMass, track.GetVelocity(),
0282 track.GetGlobalTime(),
0283 track.GetProperTime(),
0284 &spin);
0285 if (currentMinimumStep > 0) {
0286
0287
0288 lengthAlongCurve = fFieldPropagator->ComputeStep(aFieldTrack, currentMinimumStep,
0289 currentSafety, track.GetVolume());
0290 fGeometryLimitedStep = lengthAlongCurve < currentMinimumStep;
0291 if (fGeometryLimitedStep) {
0292 geometryStepLength = lengthAlongCurve;
0293 }
0294 else {
0295 geometryStepLength = currentMinimumStep;
0296 }
0297 }
0298 else {
0299 geometryStepLength = lengthAlongCurve = 0.0;
0300 fGeometryLimitedStep = false;
0301 }
0302
0303
0304
0305 fPreviousSftOrigin = startPosition;
0306 fPreviousSafety = currentSafety;
0307
0308
0309
0310
0311 fTransportEndPosition = aFieldTrack.GetPosition();
0312
0313
0314
0315 fMomentumChanged = true;
0316 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir();
0317
0318 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy();
0319
0320 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
0321 fEndGlobalTimeComputed = true;
0322
0323 fTransportEndSpin = aFieldTrack.GetSpin();
0324 fParticleIsLooping = fFieldPropagator->IsParticleLooping();
0325 endpointDistance = (fTransportEndPosition - startPosition).mag();
0326 }
0327
0328
0329
0330
0331 if (currentMinimumStep == 0.0) {
0332 if (currentSafety == 0.0) fGeometryLimitedStep = true;
0333 }
0334
0335
0336
0337
0338 if (currentSafety < endpointDistance) {
0339
0340
0341
0342 if (particleMagneticCharge != 0.0) {
0343 G4double endSafety = fLinearNavigator->ComputeSafety(fTransportEndPosition);
0344 currentSafety = endSafety;
0345 fPreviousSftOrigin = fTransportEndPosition;
0346 fPreviousSafety = currentSafety;
0347 fpSafetyHelper->SetCurrentSafety(currentSafety, fTransportEndPosition);
0348
0349
0350
0351
0352 currentSafety += endpointDistance;
0353
0354 #ifdef G4DEBUG_TRANSPORT
0355 G4cout.precision(12);
0356 G4cout << "***G4MonopoleTransportation::AlongStepGPIL ** " << G4endl;
0357 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition
0358 << " and it returned safety= " << endSafety << G4endl;
0359 G4cout << " Adding endpoint distance " << endpointDistance
0360 << " to obtain pseudo-safety= " << currentSafety << G4endl;
0361 #endif
0362 }
0363 }
0364
0365 fParticleChange.ProposeTrueStepLength(geometryStepLength);
0366
0367 fMagSetup->SetStepperAndChordFinder(0);
0368
0369
0370 return geometryStepLength;
0371 }
0372
0373
0374
0375
0376
0377
0378 G4VParticleChange* G4MonopoleTransportation::AlongStepDoIt(const G4Track& track,
0379 const G4Step& stepData)
0380 {
0381 ++noCalls;
0382
0383 if (fGeometryLimitedStep) {
0384 stepData.GetPostStepPoint()->SetStepStatus(fGeomBoundary);
0385 }
0386
0387 fParticleChange.Initialize(track);
0388
0389
0390
0391 fParticleChange.ProposePosition(fTransportEndPosition);
0392 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir);
0393 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy);
0394 fParticleChange.SetMomentumChanged(fMomentumChanged);
0395
0396 fParticleChange.ProposePolarization(fTransportEndSpin);
0397
0398 G4double deltaTime = 0.0;
0399
0400
0401
0402 G4double startTime = track.GetGlobalTime();
0403
0404 if (!fEndGlobalTimeComputed) {
0405
0406
0407 G4double finalVelocity = track.GetVelocity();
0408 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
0409 G4double stepLength = track.GetStepLength();
0410
0411 deltaTime = 0.0;
0412 if (finalVelocity > 0.0) {
0413 G4double meanInverseVelocity;
0414
0415 meanInverseVelocity = 0.5 * (1.0 / initialVelocity + 1.0 / finalVelocity);
0416 deltaTime = stepLength * meanInverseVelocity;
0417 }
0418 else if (initialVelocity > 0.0) {
0419 deltaTime = stepLength / initialVelocity;
0420 }
0421 fCandidateEndGlobalTime = startTime + deltaTime;
0422 }
0423 else {
0424 deltaTime = fCandidateEndGlobalTime - startTime;
0425 }
0426
0427 fParticleChange.ProposeGlobalTime(fCandidateEndGlobalTime);
0428
0429
0430
0431 G4double restMass = track.GetDynamicParticle()->GetMass();
0432 G4double deltaProperTime = deltaTime * (restMass / track.GetTotalEnergy());
0433
0434 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime);
0435
0436
0437
0438
0439
0440 if (fParticleIsLooping) {
0441 G4double endEnergy = fTransportEndKineticEnergy;
0442
0443 if ((endEnergy < fThreshold_Important_Energy) || (fNoLooperTrials >= fThresholdTrials)) {
0444
0445
0446 fParticleChange.ProposeTrackStatus(fStopAndKill);
0447
0448
0449 fSumEnergyKilled += endEnergy;
0450 if (endEnergy > fMaxEnergyKilled) {
0451 fMaxEnergyKilled = endEnergy;
0452 }
0453
0454 #ifdef G4VERBOSE
0455 if ((verboseLevel > 1) || (endEnergy > fThreshold_Warning_Energy)) {
0456 G4cout << " G4MonopoleTransportation is killing track "
0457 << "that is looping or stuck " << G4endl << " This track has "
0458 << track.GetKineticEnergy() / MeV << " MeV energy." << G4endl;
0459 G4cout << " Number of trials = " << fNoLooperTrials
0460 << " No of calls to AlongStepDoIt = " << noCalls << G4endl;
0461 }
0462 #endif
0463 fNoLooperTrials = 0;
0464 }
0465 else {
0466 fNoLooperTrials++;
0467 #ifdef G4VERBOSE
0468 if ((verboseLevel > 2)) {
0469 G4cout << " G4MonopoleTransportation::AlongStepDoIt(): "
0470 << "Particle looping - "
0471 << " Number of trials = " << fNoLooperTrials << " No of calls to = " << noCalls
0472 << G4endl;
0473 }
0474 #endif
0475 }
0476 }
0477 else {
0478 fNoLooperTrials = 0;
0479 }
0480
0481
0482
0483
0484
0485
0486 fParticleChange.SetPointerToVectorOfAuxiliaryPoints(
0487 fFieldPropagator->GimmeTrajectoryVectorAndForgetIt());
0488
0489 return &fParticleChange;
0490 }
0491
0492
0493
0494
0495
0496
0497
0498 G4double
0499 G4MonopoleTransportation::PostStepGetPhysicalInteractionLength(const G4Track&,
0500 G4double,
0501 G4ForceCondition* pForceCond)
0502 {
0503 *pForceCond = Forced;
0504 return DBL_MAX;
0505 }
0506
0507
0508
0509 G4VParticleChange* G4MonopoleTransportation::PostStepDoIt(const G4Track& track, const G4Step&)
0510 {
0511 G4TouchableHandle retCurrentTouchable;
0512
0513
0514
0515
0516
0517 fParticleChange.ProposeTrackStatus(track.GetTrackStatus());
0518
0519
0520
0521
0522 if (fGeometryLimitedStep) {
0523
0524
0525
0526
0527 fLinearNavigator->SetGeometricallyLimitedStep();
0528 fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0529 track.GetPosition(), track.GetMomentumDirection(), fCurrentTouchableHandle, true);
0530
0531
0532
0533 if (fCurrentTouchableHandle->GetVolume() == 0) {
0534 fParticleChange.ProposeTrackStatus(fStopAndKill);
0535 }
0536 retCurrentTouchable = fCurrentTouchableHandle;
0537 fParticleChange.SetTouchableHandle(fCurrentTouchableHandle);
0538 }
0539 else
0540 {
0541
0542
0543 fLinearNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
0544
0545
0546
0547
0548
0549
0550 fParticleChange.SetTouchableHandle(track.GetTouchableHandle());
0551 retCurrentTouchable = track.GetTouchableHandle();
0552 }
0553
0554 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume();
0555 const G4Material* pNewMaterial = 0;
0556 G4VSensitiveDetector* pNewSensitiveDetector = 0;
0557 if (pNewVol != 0) {
0558 pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial();
0559 pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector();
0560 }
0561
0562
0563
0564 fParticleChange.SetMaterialInTouchable((G4Material*)pNewMaterial);
0565 fParticleChange.SetSensitiveDetectorInTouchable(pNewSensitiveDetector);
0566
0567 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
0568 if (pNewVol != 0) {
0569 pNewMaterialCutsCouple = pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
0570 }
0571
0572 if (pNewVol != 0 && pNewMaterialCutsCouple != 0
0573 && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial)
0574 {
0575
0576
0577 pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()->GetMaterialCutsCouple(
0578 pNewMaterial, pNewMaterialCutsCouple->GetProductionCuts());
0579 }
0580 fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple);
0581
0582
0583
0584
0585
0586
0587
0588 fParticleChange.SetTouchableHandle(retCurrentTouchable);
0589
0590 return &fParticleChange;
0591 }
0592
0593
0594
0595
0596
0597
0598
0599 void G4MonopoleTransportation::StartTracking(G4Track* aTrack)
0600 {
0601 G4VProcess::StartTracking(aTrack);
0602
0603
0604
0605
0606
0607
0608 fPreviousSafety = 0.0;
0609 fPreviousSftOrigin = G4ThreeVector(0., 0., 0.);
0610
0611
0612 fNoLooperTrials = 0;
0613
0614
0615
0616
0617
0618
0619 if (DoesGlobalFieldExist()) {
0620 fFieldPropagator->ClearPropagatorState();
0621
0622
0623
0624
0625
0626
0627 }
0628
0629
0630 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
0631 fieldMgrStore->ClearAllChordFindersState();
0632
0633
0634
0635 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
0636 }
0637
0638