70 #include "G4ProductionCutsTable.hh"
71 #include "G4ParticleTable.hh"
72 #include "G4ChordFinder.hh"
73 #include "G4SafetyHelper.hh"
74 #include "G4FieldManagerStore.hh"
75 #include "G4EquationOfMotion.hh"
76 #include "G4MagIntegratorDriver.hh"
77 #include "G4MagIntegratorStepper.hh"
78 #include "G4Version.hh"
80 #include "CLHEP/Units/SystemOfUnits.h"
81 #include "CLHEP/Units/PhysicalConstants.h"
83 class G4VSensitiveDetector;
90 : G4VProcess( G4String(
"mplAtlasTransportation"), fTransportation ),
91 fTransportEndKineticEnergy (0.),
92 fMomentumChanged( false ),
94 fParticleIsLooping( false ),
95 fCurrentTouchableHandle(),
96 fGeometryLimitedStep( false ),
97 fPreviousSftOrigin (0.,0.,0.),
98 fPreviousSafety ( 0.0 ),
99 endpointDistance ( 0.0 ),
100 fThreshold_Warning_Energy( 100 *
CLHEP::
MeV ),
101 fThreshold_Important_Energy( 250 *
CLHEP::
MeV ),
102 fThresholdTrials( 10 ),
105 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
106 fShortStepOptimisation(false),
107 fVerboseLevel( verboseLevel ),
113 G4TransportationManager* transportMgr ;
115 G4cout <<
"!!! G4mplAtlasTransportation constructor " << G4endl;
117 transportMgr = G4TransportationManager::GetTransportationManager() ;
143 G4cout <<
"!!! G4mplAtlasTransportation destructor " << G4endl;
146 G4cout <<
" G4mplAtlasTransportation: Statistics for looping particles " << G4endl;
147 G4cout <<
" Sum of energy of loopers killed: " <<
fSumEnergyKilled << G4endl;
162 G4double currentMinimumStep,
163 G4double& currentSafety,
166 G4double geometryStepLength, newSafety ;
182 const G4DynamicParticle* pParticle =
track.GetDynamicParticle() ;
183 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
184 const G4ThreeVector& startMomentumDir = pParticle->GetMomentumDirection() ;
185 G4ThreeVector startPosition =
track.GetPosition() ;
194 G4double MagSqShift = OriginShift.mag2() ;
197 currentSafety = 0.0 ;
207 G4double particleElCharge = pParticle->GetCharge() ;
222 G4FieldManager* fieldMgr=0;
223 G4bool fieldExertsForce = false ;
224 if( (particleElCharge != 0.0) || (particleMagCharge!=0.0) )
232 if (fieldMgr != oldFieldMgr) {
238 if (fieldMgr !=
nullptr) {
240 fieldMgr->ConfigureForTrack( &
track );
245 if (particleMagCharge!=0.0) fieldMgr->SetFieldChangesEnergy(
true);
248 fieldExertsForce = (fieldMgr->GetDetectorField() !=
nullptr);
250 if( fieldExertsForce)
251 fEquationSetup->SwitchStepperAndChordFinder( (particleMagCharge != 0.0), fieldMgr );
260 if( !fieldExertsForce )
262 G4double linearStepLength ;
267 geometryStepLength = currentMinimumStep ;
286 currentSafety = newSafety ;
292 geometryStepLength = linearStepLength ;
297 geometryStepLength = currentMinimumStep ;
317 G4ThreeVector EndUnitMomentum ;
318 G4double lengthAlongCurve ;
319 G4double restMass = pParticleDef->GetPDGMass() ;
320 #if G4VERSION_NUMBER > 1009
321 G4double momentumMagnitude = pParticle->GetTotalMomentum();
322 G4ChargeState chargeState(particleElCharge,
323 pParticleDef->GetPDGSpin(),
327 G4EquationOfMotion* equationOfMotion = (
fFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper())->GetEquationOfMotion();
328 equationOfMotion->SetChargeMomentumMass(chargeState,
337 G4ThreeVector
spin =
track.GetPolarization() ;
338 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
339 track.GetMomentumDirection(),
341 track.GetKineticEnergy(),
344 track.GetGlobalTime(),
345 track.GetProperTime(),
347 if( currentMinimumStep > 0 )
355 track.GetVolume() ) ;
358 geometryStepLength = lengthAlongCurve ;
360 geometryStepLength = currentMinimumStep ;
365 geometryStepLength = lengthAlongCurve= 0.0 ;
407 G4double startEnergy=
track.GetKineticEnergy();
410 static std::atomic<G4int> no_inexact_steps=0, no_large_ediff;
411 G4double absEdiff = std::fabs(startEnergy- endEnergy);
421 static std::atomic<G4int> no_warnings= 0, warnModulo=1, moduloFactor= 10;
423 if( (no_large_ediff% warnModulo) == 0 )
426 G4cout <<
"WARNING - G4mplAtlasTransportation::AlongStepGetPIL() "
427 <<
" Energy change in Step is above 1^-3 relative value. " << G4endl
428 <<
" Relative change in 'tracking' step = "
429 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
430 <<
" Starting E= " << std::setw(12) << startEnergy /
CLHEP::MeV <<
" MeV " << G4endl
431 <<
" Ending E= " << std::setw(12) << endEnergy /
CLHEP::MeV <<
" MeV " << G4endl;
432 G4cout <<
" Energy has been corrected -- however, review"
433 <<
" field propagation parameters for accuracy." << G4endl;
434 if( (
fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){
435 G4cout <<
" These include EpsilonStepMax(/Min) in G4FieldManager "
436 <<
" which determine fractional error per step for integrated quantities. " << G4endl
437 <<
" Note also the influence of the permitted number of integration steps."
440 G4cerr <<
"ERROR - G4mplAtlasTransportation::AlongStepGetPIL()" << G4endl
441 <<
" Bad 'endpoint'. Energy change detected"
442 <<
" and corrected. "
443 <<
" Has occurred already "
444 << no_large_ediff <<
" times." << G4endl;
445 if( no_large_ediff == warnModulo * moduloFactor )
447 warnModulo = warnModulo * moduloFactor;
463 fieldMgr->SetFieldChangesEnergy(
false);
470 if( currentMinimumStep == 0.0 )
487 currentSafety = endSafety ;
497 #ifdef G4DEBUG_TRANSPORT
498 G4cout.precision(12) ;
499 G4cout <<
"***G4mplAtlasTransportation::AlongStepGPIL ** " << G4endl ;
501 <<
" and it returned safety= " << endSafety << G4endl ;
503 <<
" to obtain pseudo-safety= " << currentSafety << G4endl ;
532 return geometryStepLength ;
541 const G4Step& stepData )
544 static std::atomic<G4int> noCalls=0;
548 static const G4ParticleDefinition*
const fOpticalPhoton =
549 G4ParticleTable::GetParticleTable()->FindParticle(
"opticalphoton");
564 G4double deltaTime = 0.0 ;
580 G4double finalVelocity =
track.GetVelocity() ;
581 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
582 G4double stepLength =
track.GetStepLength() ;
585 const G4DynamicParticle* fpDynamicParticle =
track.GetDynamicParticle();
586 if (fpDynamicParticle->GetDefinition()== fOpticalPhoton)
590 deltaTime = stepLength/finalVelocity ;
592 else if (finalVelocity > 0.0)
594 G4double meanInverseVelocity ;
596 meanInverseVelocity = 0.5
597 * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ;
598 deltaTime = stepLength * meanInverseVelocity ;
600 else if( initialVelocity > 0.0 )
602 deltaTime = stepLength/initialVelocity ;
615 G4double restMass =
track.GetDynamicParticle()->GetMass() ;
616 G4double deltaProperTime = deltaTime*( restMass/
track.GetTotalEnergy() ) ;
650 G4cout <<
" G4mplAtlasTransportation is killing track that is looping or stuck "
653 <<
" MeV energy." << G4endl;
655 <<
" No of calls to AlongStepDoIt = " << noCalls
665 G4cout <<
" G4mplAtlasTransportation::AlongStepDoIt(): Particle looping - "
667 <<
" No of calls to = " << noCalls
696 G4ForceCondition* pForceCond )
698 *pForceCond = Forced ;
708 G4TouchableHandle retCurrentTouchable ;
730 LocateGlobalPointAndUpdateTouchableHandle(
track.GetPosition(),
731 track.GetMomentumDirection(),
756 retCurrentTouchable =
track.GetTouchableHandle() ;
759 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
760 G4Material* pNewMaterial = 0 ;
761 G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
765 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
766 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
773 fParticleChange.SetSensitiveDetectorInTouchable( pNewSensitiveDetector ) ;
775 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
778 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
781 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
785 pNewMaterialCutsCouple =
786 G4ProductionCutsTable::GetProductionCutsTable()
787 ->GetMaterialCutsCouple(pNewMaterial,
788 pNewMaterialCutsCouple->GetProductionCuts());
790 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
809 G4VProcess::StartTracking(aTrack);
836 G4FieldManagerStore::GetInstance()->ClearAllChordFindersState();
848 fEquationSetup->SwitchStepperAndChordFinder( (particleMagCharge != 0.0), fieldMgr );
860 G4VProcess::EndTracking();