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