67 fParticleIsLooping( false ),
68 fPreviousSftOrigin (0.,0.,0.),
69 fPreviousSafety ( 0.0 ),
70 fThreshold_Warning_Energy( 100 *
MeV ),
71 fThreshold_Important_Energy( 250 *
MeV ),
72 fThresholdTrials( 10 ),
75 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
76 fShortStepOptimisation(false),
105 fCurrentTouchableHandle = nullTouchableHandle;
107 fEndGlobalTimeComputed =
false;
108 fCandidateEndGlobalTime = 0;
116 G4cout <<
" G4MonopoleTransportation: Statistics for looping particles "
118 G4cout <<
" Sum of energy of loopers killed: " << fSumEnergyKilled <<
G4endl;
119 G4cout <<
" Max energy of loopers killed: " << fMaxEnergyKilled <<
G4endl;
140 G4double geometryStepLength, newSafety ;
141 fParticleIsLooping = false ;
166 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
168 if( MagSqShift >=
sqr(fPreviousSafety) )
170 currentSafety = 0.0 ;
174 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
182 fGeometryLimitedStep = false ;
192 G4bool fieldExertsForce = false ;
194 if( (particleMagneticCharge != 0.0) )
214 if( !fieldExertsForce )
217 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
221 geometryStepLength = currentMinimumStep ;
222 fGeometryLimitedStep = false ;
228 linearStepLength = fLinearNavigator->
ComputeStep( startPosition,
234 fPreviousSftOrigin = startPosition ;
235 fPreviousSafety = newSafety ;
240 currentSafety = newSafety ;
242 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
243 if( fGeometryLimitedStep )
246 geometryStepLength = linearStepLength ;
251 geometryStepLength = currentMinimumStep ;
254 endpointDistance = geometryStepLength ;
258 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
262 fTransportEndMomentumDir = startMomentumDir ;
265 fParticleIsLooping = false ;
266 fMomentumChanged = false ;
267 fEndGlobalTimeComputed = false ;
280 particleMagneticCharge );
284 ->GetEquationOfMotion();
301 if( currentMinimumStep > 0 )
305 lengthAlongCurve = fFieldPropagator->
ComputeStep( aFieldTrack,
309 fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep;
310 if( fGeometryLimitedStep ) {
311 geometryStepLength = lengthAlongCurve ;
313 geometryStepLength = currentMinimumStep ;
318 geometryStepLength = lengthAlongCurve= 0.0 ;
319 fGeometryLimitedStep = false ;
324 fPreviousSftOrigin = startPosition ;
325 fPreviousSafety = currentSafety ;
330 fTransportEndPosition = aFieldTrack.
GetPosition() ;
334 fMomentumChanged = true ;
340 fEndGlobalTimeComputed =
true;
342 fTransportEndSpin = aFieldTrack.
GetSpin();
344 endpointDistance = (fTransportEndPosition - startPosition).mag() ;
350 if( currentMinimumStep == 0.0 )
352 if( currentSafety == 0.0 ) fGeometryLimitedStep = true ;
358 if( currentSafety < endpointDistance )
363 if( particleMagneticCharge != 0.0 ) {
367 currentSafety = endSafety ;
368 fPreviousSftOrigin = fTransportEndPosition ;
369 fPreviousSafety = currentSafety ;
375 currentSafety += endpointDistance ;
377 #ifdef G4DEBUG_TRANSPORT
379 G4cout <<
"***G4MonopoleTransportation::AlongStepGPIL ** " <<
G4endl ;
380 G4cout <<
" Called Navigator->ComputeSafety at " << fTransportEndPosition
381 <<
" and it returned safety= " << endSafety <<
G4endl ;
382 G4cout <<
" Adding endpoint distance " << endpointDistance
383 <<
" to obtain pseudo-safety= " << currentSafety <<
G4endl ;
393 return geometryStepLength ;
404 static G4int noCalls=0;
429 if (!fEndGlobalTimeComputed)
443 deltaTime = stepLength/finalVelocity ;
445 else if (finalVelocity > 0.0)
449 meanInverseVelocity = 0.5
450 * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ;
451 deltaTime = stepLength * meanInverseVelocity ;
453 else if( initialVelocity > 0.0 )
455 deltaTime = stepLength/initialVelocity ;
457 fCandidateEndGlobalTime = startTime + deltaTime ;
461 deltaTime = fCandidateEndGlobalTime - startTime ;
478 if ( fParticleIsLooping )
480 G4double endEnergy= fTransportEndKineticEnergy;
482 if( (endEnergy < fThreshold_Important_Energy)
483 || (fNoLooperTrials >= fThresholdTrials ) ){
489 fSumEnergyKilled += endEnergy;
490 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
494 ( endEnergy > fThreshold_Warning_Energy ) ) {
495 G4cout <<
" G4MonopoleTransportation is killing track that is looping or stuck "
498 <<
" MeV energy." <<
G4endl;
499 G4cout <<
" Number of trials = " << fNoLooperTrials
500 <<
" No of calls to AlongStepDoIt = " << noCalls
510 G4cout <<
" G4MonopoleTransportation::AlongStepDoIt(): Particle looping - "
511 <<
" Number of trials = " << fNoLooperTrials
512 <<
" No of calls to = " << noCalls
529 return &fParticleChange ;
563 if(fGeometryLimitedStep)
571 LocateGlobalPointAndUpdateTouchableHandle( track.
GetPosition(),
573 fCurrentTouchableHandle,
578 if( fCurrentTouchableHandle->
GetVolume() == 0 )
582 retCurrentTouchable = fCurrentTouchableHandle ;
622 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->
GetMaterial()!=pNewMaterial )
626 pNewMaterialCutsCouple =
641 return &fParticleChange ;
660 fPreviousSafety = 0.0 ;
void SetMaterialInTouchable(G4Material *fMaterial)
void SetTouchableHandle(const G4TouchableHandle &fTouchable)
const G4ThreeVector & GetPolarization() const
G4SafetyHelper * GetSafetyHelper() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double ¤tSafety, G4GPILSelection *selection)
void SetMaterialCutsCoupleInTouchable(const G4MaterialCutsCouple *fMaterialCutsCouple)
virtual void SetChargeMomentumMass(G4ChargeState particleCharge, G4double MomentumXc, G4double MassXc2)=0
G4double GetProperTime() const
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetVelocity() const
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
const G4MagIntegratorStepper * GetStepper() const
void ClearAllChordFindersState()
G4Material * GetMaterial() const
const G4ThreeVector & GetPosition() const
const G4ThreeVector & GetMomentumDir() const
G4TrackStatus GetTrackStatus() const
G4ThreeVector GetSpin() const
virtual G4double ComputeStep(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double &pNewSafety)
G4Navigator * GetNavigatorForTracking() const
virtual void StartTracking(G4Track *aTrack)
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &stepData)
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
~G4MonopoleTransportation()
G4ParticleDefinition * GetDefinition() const
G4double GetVelocity() const
Definition of the G4MonopoleTransportation class.
void ProposePosition(G4double x, G4double y, G4double z)
G4double GetTotalMomentum() const
G4StepPoint * GetPreStepPoint() const
virtual void Initialize(const G4Track &)
virtual void StartTracking(G4Track *)
virtual void ConfigureForTrack(const G4Track *)
Definition of the G4Monopole class.
void SetGeometricallyLimitedStep()
G4double GetKineticEnergy() const
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
G4double MagneticCharge() const
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &, G4double previousStepSize, G4ForceCondition *pForceCond)
const G4ThreeVector & GetMomentumDirection() const
G4double GetCharge() const
void ProposeTrueStepLength(G4double truePathLength)
void SetProcessSubType(G4int)
static G4MonopoleFieldSetup * GetMonopoleFieldSetup()
G4double GetGlobalTime() const
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
const G4TouchableHandle & GetTouchableHandle() const
G4bool DoesGlobalFieldExist()
static G4TransportationManager * GetTransportationManager()
G4MonopoleTransportation(const G4Monopole *p, G4int verbosityLevel=1)
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ThreeVector & GetMomentumDirection() const
void SetStepperAndChordFinder(G4int val)
G4LogicalVolume * GetLogicalVolume() const
void ProposeProperTime(G4double finalProperTime)
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4FieldManagerStore * GetInstance()
void SetSensitiveDetectorInTouchable(G4VSensitiveDetector *fSensitiveDetector)
static G4ParticleTable * GetParticleTable()
G4bool IsParticleLooping() const
G4double GetLabTimeOfFlight() const
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
void SetPointerToVectorOfAuxiliaryPoints(std::vector< G4ThreeVector > *theNewVectorPointer)
void ProposeGlobalTime(G4double t)
G4ChordFinder * GetChordFinder()
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
void ProposeEnergy(G4double finalEnergy)
G4VPhysicalVolume * GetVolume() const
G4double GetTotalEnergy() const
G4double GetPDGSpin() const
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
void ClearPropagatorState()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
void ProposeTrackStatus(G4TrackStatus status)
const G4Field * GetDetectorField() const
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
G4MagInt_Driver * GetIntegrationDriver()
void SetMomentumChanged(G4bool b)
G4PropagatorInField * GetPropagatorInField() const
G4ProductionCuts * GetProductionCuts() const
G4VSensitiveDetector * GetSensitiveDetector() const
G4double GetStepLength() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
const G4Material * GetMaterial() const