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