00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 #include "G4Transportation.hh"
00060 #include "G4TransportationProcessType.hh"
00061
00062 #include "G4PhysicalConstants.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "G4ProductionCutsTable.hh"
00065 #include "G4ParticleTable.hh"
00066 #include "G4ChordFinder.hh"
00067 #include "G4SafetyHelper.hh"
00068 #include "G4FieldManagerStore.hh"
00069
00070 class G4VSensitiveDetector;
00071
00073
00074
00075
00076 G4Transportation::G4Transportation( G4int verbosity )
00077 : G4VProcess( G4String("Transportation"), fTransportation ),
00078 fTransportEndPosition( 0.0, 0.0, 0.0 ),
00079 fTransportEndMomentumDir( 0.0, 0.0, 0.0 ),
00080 fTransportEndKineticEnergy( 0.0 ),
00081 fTransportEndSpin( 0.0, 0.0, 0.0 ),
00082 fMomentumChanged(true),
00083 fEndGlobalTimeComputed(false),
00084 fCandidateEndGlobalTime(0.0),
00085 fParticleIsLooping( false ),
00086 fGeometryLimitedStep(true),
00087 fPreviousSftOrigin( 0.,0.,0. ),
00088 fPreviousSafety( 0.0 ),
00089
00090 fEndPointDistance( -1.0 ),
00091 fThreshold_Warning_Energy( 100 * MeV ),
00092 fThreshold_Important_Energy( 250 * MeV ),
00093 fThresholdTrials( 10 ),
00094 fNoLooperTrials( 0 ),
00095 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
00096 fShortStepOptimisation( false ),
00097 fUseMagneticMoment( false ),
00098 fVerboseLevel( verbosity )
00099 {
00100
00101 SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
00102
00103 G4TransportationManager* transportMgr ;
00104
00105 transportMgr = G4TransportationManager::GetTransportationManager() ;
00106
00107 fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
00108
00109 fFieldPropagator = transportMgr->GetPropagatorInField() ;
00110
00111 fpSafetyHelper = transportMgr->GetSafetyHelper();
00112
00113
00114
00115
00116
00117
00118 static G4TouchableHandle nullTouchableHandle;
00119 fCurrentTouchableHandle = nullTouchableHandle;
00120
00121 #ifdef G4VERBOSE
00122 if( fVerboseLevel > 0)
00123 {
00124 G4cout << " G4Transportation constructor> set fShortStepOptimisation to ";
00125 if ( fShortStepOptimisation ) G4cout << "true" << G4endl;
00126 else G4cout << "false" << G4endl;
00127 }
00128 #endif
00129 }
00130
00132
00133 G4Transportation::~G4Transportation()
00134 {
00135 if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) )
00136 {
00137 G4cout << " G4Transportation: Statistics for looping particles " << G4endl;
00138 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
00139 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
00140 }
00141 }
00142
00144
00145
00146
00147
00148
00149
00150 G4double G4Transportation::
00151 AlongStepGetPhysicalInteractionLength( const G4Track& track,
00152 G4double,
00153 G4double currentMinimumStep,
00154 G4double& currentSafety,
00155 G4GPILSelection* selection )
00156 {
00157 G4double geometryStepLength= -1.0, newSafety= -1.0;
00158 fParticleIsLooping = false ;
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 *selection = CandidateForSelection ;
00170
00171
00172
00173 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
00174 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
00175 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
00176 G4ThreeVector startPosition = track.GetPosition() ;
00177
00178
00179
00180
00181
00182
00183
00184 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
00185 G4double MagSqShift = OriginShift.mag2() ;
00186 if( MagSqShift >= sqr(fPreviousSafety) )
00187 {
00188 currentSafety = 0.0 ;
00189 }
00190 else
00191 {
00192 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
00193 }
00194
00195
00196
00197 G4double particleCharge = pParticle->GetCharge() ;
00198 G4double magneticMoment = pParticle->GetMagneticMoment() ;
00199 G4double restMass = pParticleDef->GetPDGMass() ;
00200
00201 fGeometryLimitedStep = false ;
00202
00203
00204
00205
00206
00207
00208
00209
00210 G4FieldManager* fieldMgr=0;
00211 G4bool fieldExertsForce = false ;
00212
00213 G4bool gravityOn = false;
00214 G4bool fieldExists= false;
00215
00216 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
00217 if( fieldMgr != 0 )
00218 {
00219
00220 fieldMgr->ConfigureForTrack( &track );
00221
00222
00223
00224
00225
00226 const G4Field* ptrField= fieldMgr->GetDetectorField();
00227 fieldExists = (ptrField!=0) ;
00228 if( fieldExists )
00229 {
00230 gravityOn= ptrField->IsGravityActive();
00231
00232 if( (particleCharge != 0.0)
00233 || (fUseMagneticMoment && (magneticMoment != 0.0) )
00234 || (gravityOn && (restMass != 0.0) )
00235 )
00236 {
00237 fieldExertsForce = fieldExists;
00238 }
00239 }
00240 }
00241
00242
00243
00244 if( !fieldExertsForce )
00245 {
00246 G4double linearStepLength ;
00247 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
00248 {
00249
00250
00251 geometryStepLength = currentMinimumStep ;
00252 fGeometryLimitedStep = false ;
00253 }
00254 else
00255 {
00256
00257
00258 linearStepLength = fLinearNavigator->ComputeStep( startPosition,
00259 startMomentumDir,
00260 currentMinimumStep,
00261 newSafety) ;
00262
00263
00264 fPreviousSftOrigin = startPosition ;
00265 fPreviousSafety = newSafety ;
00266 fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
00267
00268 currentSafety = newSafety ;
00269
00270 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
00271 if( fGeometryLimitedStep )
00272 {
00273
00274 geometryStepLength = linearStepLength ;
00275 }
00276 else
00277 {
00278
00279 geometryStepLength = currentMinimumStep ;
00280 }
00281 }
00282 fEndPointDistance = geometryStepLength ;
00283
00284
00285
00286 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
00287
00288
00289
00290 fTransportEndMomentumDir = startMomentumDir ;
00291 fTransportEndKineticEnergy = track.GetKineticEnergy() ;
00292 fTransportEndSpin = track.GetPolarization();
00293 fParticleIsLooping = false ;
00294 fMomentumChanged = false ;
00295 fEndGlobalTimeComputed = false ;
00296 }
00297 else
00298 {
00299 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
00300 G4ThreeVector EndUnitMomentum ;
00301 G4double lengthAlongCurve ;
00302
00303 fFieldPropagator->SetChargeMomentumMass( particleCharge,
00304 momentumMagnitude,
00305 restMass ) ;
00306
00307 G4ThreeVector spin = track.GetPolarization() ;
00308 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
00309 track.GetMomentumDirection(),
00310 0.0,
00311 track.GetKineticEnergy(),
00312 restMass,
00313 track.GetVelocity(),
00314 track.GetGlobalTime(),
00315 track.GetProperTime(),
00316 &spin ) ;
00317 if( currentMinimumStep > 0 )
00318 {
00319
00320
00321 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
00322 currentMinimumStep,
00323 currentSafety,
00324 track.GetVolume() ) ;
00325 fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep;
00326 if( fGeometryLimitedStep )
00327 {
00328 geometryStepLength = lengthAlongCurve ;
00329 }
00330 else
00331 {
00332 geometryStepLength = currentMinimumStep ;
00333 }
00334
00335
00336
00337 fPreviousSftOrigin = startPosition ;
00338 fPreviousSafety = currentSafety ;
00339 fpSafetyHelper->SetCurrentSafety( currentSafety, startPosition);
00340 }
00341 else
00342 {
00343 geometryStepLength = lengthAlongCurve= 0.0 ;
00344 fGeometryLimitedStep = false ;
00345 }
00346
00347
00348
00349 fTransportEndPosition = aFieldTrack.GetPosition() ;
00350
00351
00352
00353 fMomentumChanged = true ;
00354 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
00355
00356 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ;
00357
00358 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
00359 {
00360
00361
00362
00363 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
00364 fEndGlobalTimeComputed = true;
00365
00366
00367
00368 }
00369 else
00370 {
00371
00372
00373
00374 fEndGlobalTimeComputed = false;
00375
00376
00377
00378 G4double startEnergy= track.GetKineticEnergy();
00379 G4double endEnergy= fTransportEndKineticEnergy;
00380
00381 static G4int no_inexact_steps=0, no_large_ediff;
00382 G4double absEdiff = std::fabs(startEnergy- endEnergy);
00383 if( absEdiff > perMillion * endEnergy )
00384 {
00385 no_inexact_steps++;
00386
00387 }
00388 if( fVerboseLevel > 1 )
00389 {
00390 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
00391 {
00392 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
00393 no_large_ediff ++;
00394 if( (no_large_ediff% warnModulo) == 0 )
00395 {
00396 no_warnings++;
00397 G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
00398 << " Energy change in Step is above 1^-3 relative value. " << G4endl
00399 << " Relative change in 'tracking' step = "
00400 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
00401 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
00402 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl;
00403 G4cout << " Energy has been corrected -- however, review"
00404 << " field propagation parameters for accuracy." << G4endl;
00405 if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
00406 {
00407 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
00408 << " which determine fractional error per step for integrated quantities. " << G4endl
00409 << " Note also the influence of the permitted number of integration steps."
00410 << G4endl;
00411 }
00412 G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
00413 << " Bad 'endpoint'. Energy change detected"
00414 << " and corrected. "
00415 << " Has occurred already "
00416 << no_large_ediff << " times." << G4endl;
00417 if( no_large_ediff == warnModulo * moduloFactor )
00418 {
00419 warnModulo *= moduloFactor;
00420 }
00421 }
00422 }
00423 }
00424
00425
00426
00427
00428 fTransportEndKineticEnergy= track.GetKineticEnergy();
00429 }
00430
00431 fTransportEndSpin = aFieldTrack.GetSpin();
00432 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
00433 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
00434 }
00435
00436
00437
00438
00439 if( currentMinimumStep == 0.0 )
00440 {
00441 if( currentSafety == 0.0 ) { fGeometryLimitedStep = true; }
00442 }
00443
00444
00445
00446
00447 if( currentSafety < fEndPointDistance )
00448 {
00449 if( particleCharge != 0.0 )
00450 {
00451 G4double endSafety =
00452 fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
00453 currentSafety = endSafety ;
00454 fPreviousSftOrigin = fTransportEndPosition ;
00455 fPreviousSafety = currentSafety ;
00456 fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition);
00457
00458
00459
00460
00461 currentSafety += fEndPointDistance ;
00462
00463 #ifdef G4DEBUG_TRANSPORT
00464 G4cout.precision(12) ;
00465 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
00466 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition
00467 << " and it returned safety= " << endSafety << G4endl ;
00468 G4cout << " Adding endpoint distance " << fEndPointDistance
00469 << " to obtain pseudo-safety= " << currentSafety << G4endl ;
00470 }
00471 else
00472 {
00473 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
00474 G4cout << " Avoiding call to ComputeSafety : " << G4endl;
00475 G4cout << " charge = " << particleCharge << G4endl;
00476 G4cout << " mag moment = " << magneticMoment << G4endl;
00477 #endif
00478 }
00479 }
00480
00481 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
00482
00483 return geometryStepLength ;
00484 }
00485
00487
00488
00489
00490
00491 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track,
00492 const G4Step& stepData )
00493 {
00494 static G4int noCalls=0;
00495 noCalls++;
00496
00497 fParticleChange.Initialize(track) ;
00498
00499
00500
00501 fParticleChange.ProposePosition(fTransportEndPosition) ;
00502 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
00503 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
00504 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
00505
00506 fParticleChange.ProposePolarization(fTransportEndSpin);
00507
00508 G4double deltaTime = 0.0 ;
00509
00510
00511
00512
00513
00514 G4double startTime = track.GetGlobalTime() ;
00515
00516 if (!fEndGlobalTimeComputed)
00517 {
00518
00519
00520 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
00521 G4double stepLength = track.GetStepLength();
00522
00523 deltaTime= 0.0;
00524 if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; }
00525
00526 fCandidateEndGlobalTime = startTime + deltaTime ;
00527 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
00528 }
00529 else
00530 {
00531 deltaTime = fCandidateEndGlobalTime - startTime ;
00532 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
00533 }
00534
00535
00536
00537
00538 G4double restMass = track.GetDynamicParticle()->GetMass() ;
00539 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
00540
00541 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
00542
00543
00544
00545
00546
00547
00548 if ( fParticleIsLooping )
00549 {
00550 G4double endEnergy= fTransportEndKineticEnergy;
00551
00552 if( (endEnergy < fThreshold_Important_Energy)
00553 || (fNoLooperTrials >= fThresholdTrials ) )
00554 {
00555
00556
00557 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
00558
00559
00560 fSumEnergyKilled += endEnergy;
00561 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
00562
00563 #ifdef G4VERBOSE
00564 if( (fVerboseLevel > 1) &&
00565 ( endEnergy > fThreshold_Warning_Energy ) )
00566 {
00567 G4cout << " G4Transportation is killing track that is looping or stuck "
00568 << G4endl
00569 << " This track has " << track.GetKineticEnergy() / MeV
00570 << " MeV energy." << G4endl;
00571 G4cout << " Number of trials = " << fNoLooperTrials
00572 << " No of calls to AlongStepDoIt = " << noCalls
00573 << G4endl;
00574 }
00575 #endif
00576 fNoLooperTrials=0;
00577 }
00578 else
00579 {
00580 fNoLooperTrials ++;
00581 #ifdef G4VERBOSE
00582 if( (fVerboseLevel > 2) )
00583 {
00584 G4cout << " G4Transportation::AlongStepDoIt(): Particle looping - "
00585 << " Number of trials = " << fNoLooperTrials
00586 << " No of calls to = " << noCalls
00587 << G4endl;
00588 }
00589 #endif
00590 }
00591 }
00592 else
00593 {
00594 fNoLooperTrials=0;
00595 }
00596
00597
00598
00599
00600
00601
00602 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
00603 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
00604
00605 return &fParticleChange ;
00606 }
00607
00609
00610
00611
00612
00613
00614 G4double G4Transportation::
00615 PostStepGetPhysicalInteractionLength( const G4Track&,
00616 G4double,
00617 G4ForceCondition* pForceCond )
00618 {
00619 *pForceCond = Forced ;
00620 return DBL_MAX ;
00621 }
00622
00624
00625
00626 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track,
00627 const G4Step& )
00628 {
00629 G4TouchableHandle retCurrentTouchable ;
00630 G4bool isLastStep= false;
00631
00632
00633
00634
00635
00636 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
00637
00638
00639
00640
00641 if(fGeometryLimitedStep)
00642 {
00643
00644
00645
00646
00647 fLinearNavigator->SetGeometricallyLimitedStep() ;
00648 fLinearNavigator->
00649 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
00650 track.GetMomentumDirection(),
00651 fCurrentTouchableHandle,
00652 true ) ;
00653
00654
00655
00656 if( fCurrentTouchableHandle->GetVolume() == 0 )
00657 {
00658 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
00659 }
00660 retCurrentTouchable = fCurrentTouchableHandle ;
00661 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
00662
00663
00664 isLastStep = fLinearNavigator->ExitedMotherVolume()
00665 | fLinearNavigator->EnteredDaughterVolume() ;
00666
00667 #ifdef G4DEBUG_TRANSPORT
00668
00669
00670 G4bool exiting = fLinearNavigator->ExitedMotherVolume();
00671 G4bool entering = fLinearNavigator->EnteredDaughterVolume();
00672 if( ! (exiting || entering) )
00673 {
00674 G4cout << " Transport> : Proposed isLastStep= " << isLastStep
00675 << " Exiting " << fLinearNavigator->ExitedMotherVolume()
00676 << " Entering " << fLinearNavigator->EnteredDaughterVolume()
00677 << G4endl;
00678 }
00679 #endif
00680 }
00681 else
00682 {
00683
00684
00685 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
00686
00687
00688
00689
00690
00691
00692 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
00693 retCurrentTouchable = track.GetTouchableHandle() ;
00694
00695 isLastStep= false;
00696
00697 #ifdef G4DEBUG_TRANSPORT
00698
00699
00700 G4cout << " Transport> Proposed isLastStep= " << isLastStep
00701 << " Geometry did not limit step. " << G4endl;
00702 #endif
00703 }
00704
00705 fParticleChange.ProposeLastStepInVolume(isLastStep);
00706
00707 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
00708 const G4Material* pNewMaterial = 0 ;
00709 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
00710
00711 if( pNewVol != 0 )
00712 {
00713 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
00714 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
00715 }
00716
00717
00718
00719
00720 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
00721 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
00722
00723 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
00724 if( pNewVol != 0 )
00725 {
00726 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
00727 }
00728
00729 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
00730 {
00731
00732
00733 pNewMaterialCutsCouple =
00734 G4ProductionCutsTable::GetProductionCutsTable()
00735 ->GetMaterialCutsCouple(pNewMaterial,
00736 pNewMaterialCutsCouple->GetProductionCuts());
00737 }
00738 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
00739
00740
00741
00742
00743
00744
00745
00746 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
00747
00748 return &fParticleChange ;
00749 }
00750
00751
00752
00753
00754 void
00755 G4Transportation::StartTracking(G4Track* aTrack)
00756 {
00757 G4VProcess::StartTracking(aTrack);
00758
00759
00760
00761
00762
00763
00764 fPreviousSafety = 0.0 ;
00765 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
00766
00767
00768 fNoLooperTrials= 0;
00769
00770
00771
00772
00773
00774
00775 if( DoesGlobalFieldExist() )
00776 {
00777 fFieldPropagator->ClearPropagatorState();
00778
00779
00780
00781
00782
00783 }
00784
00785
00786
00787 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
00788 fieldMgrStore->ClearAllChordFindersState();
00789
00790
00791
00792 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
00793 }