Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // 28 // ------------------------------------------------------------ 29 // GEANT 4 include file implementation 30 // 31 // ------------------------------------------------------------ 32 // 33 // This class is a process responsible for the transportation of 34 // a particle, ie the geometrical propagation that encounters the 35 // geometrical sub-volumes of the detectors. 36 // 37 // It is also tasked with the key role of proposing the "isotropic safety", 38 // which will be used to update the post-step point's safety. 39 // 40 // ======================================================================= 41 // Created: 19 March 1997, J. Apostolakis 42 // ======================================================================= 43 44 #include "G4Transportation.hh" 45 #include "G4TransportationProcessType.hh" 46 #include "G4TransportationLogger.hh" 47 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4ProductionCutsTable.hh" 51 #include "G4ParticleTable.hh" 52 53 #include "G4ChargeState.hh" 54 #include "G4EquationOfMotion.hh" 55 56 #include "G4FieldManagerStore.hh" 57 58 #include "G4Navigator.hh" 59 #include "G4PropagatorInField.hh" 60 #include "G4TransportationManager.hh" 61 62 #include "G4TransportationParameters.hh" 63 64 class G4VSensitiveDetector; 65 66 G4bool G4Transportation::fUseMagneticMoment=false; 67 G4bool G4Transportation::fUseGravity= false; 68 G4bool G4Transportation::fSilenceLooperWarnings= false; 69 70 ////////////////////////////////////////////////////////////////////////// 71 // 72 // Constructor 73 74 G4Transportation::G4Transportation( G4int verbosity, const G4String& aName ) 75 : G4VProcess( aName, fTransportation ), 76 fFieldExertedForce( false ), 77 fPreviousSftOrigin( 0.,0.,0. ), 78 fPreviousSafety( 0.0 ), 79 fEndPointDistance( -1.0 ), 80 fShortStepOptimisation( false ) // Old default: true (=fast short steps) 81 { 82 SetProcessSubType(static_cast<G4int>(TRANSPORTATION)); 83 pParticleChange= &fParticleChange; // Required to conform to G4VProcess 84 SetVerboseLevel(verbosity); 85 86 G4TransportationManager* transportMgr ; 87 88 transportMgr = G4TransportationManager::GetTransportationManager() ; 89 90 fLinearNavigator = transportMgr->GetNavigatorForTracking() ; 91 92 fFieldPropagator = transportMgr->GetPropagatorInField() ; 93 94 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New 95 96 fpLogger = new G4TransportationLogger("G4Transportation", verbosity); 97 98 if( G4TransportationParameters::Exists() ) 99 { 100 auto trParams= G4TransportationParameters::Instance(); 101 102 SetThresholdWarningEnergy( trParams->GetWarningEnergy() ); 103 SetThresholdImportantEnergy( trParams->GetImportantEnergy() ); 104 SetThresholdTrials( trParams->GetNumberOfTrials() ); 105 G4Transportation::fSilenceLooperWarnings= trParams->GetSilenceAllLooperWarnings(); 106 } 107 else { 108 SetHighLooperThresholds(); 109 // Use the old defaults: Warning = 100 MeV, Important = 250 MeV, No Trials = 10; 110 } 111 112 PushThresholdsToLogger(); 113 // Should be done by Set methods in SetHighLooperThresholds -- making sure 114 115 static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0; 116 if ( !pNullTouchableHandle) 117 { 118 pNullTouchableHandle = new G4TouchableHandle; 119 } 120 fCurrentTouchableHandle = *pNullTouchableHandle; 121 // Points to (G4VTouchable*) 0 122 123 124 #ifdef G4VERBOSE 125 if( verboseLevel > 0) 126 { 127 G4cout << " G4Transportation constructor> set fShortStepOptimisation to "; 128 if ( fShortStepOptimisation ) { G4cout << "true" << G4endl; } 129 else { G4cout << "false" << G4endl; } 130 } 131 #endif 132 } 133 134 ////////////////////////////////////////////////////////////////////////// 135 136 G4Transportation::~G4Transportation() 137 { 138 if( fSumEnergyKilled > 0.0 ) 139 { 140 PrintStatistics( G4cout ); 141 } 142 delete fpLogger; 143 } 144 145 ////////////////////////////////////////////////////////////////////////// 146 147 void 148 G4Transportation::PrintStatistics( std::ostream& outStr) const 149 { 150 outStr << " G4Transportation: Statistics for looping particles " << G4endl; 151 if( fSumEnergyKilled > 0.0 || fNumLoopersKilled > 0 ) 152 { 153 outStr << " Sum of energy of looping tracks killed: " 154 << fSumEnergyKilled / CLHEP::MeV << " MeV " 155 << " from " << fNumLoopersKilled << " tracks " << G4endl 156 << " Sum of energy of non-electrons : " 157 << fSumEnergyKilled_NonElectron / CLHEP::MeV << " MeV " 158 << " from " << fNumLoopersKilled_NonElectron << " tracks " 159 << G4endl; 160 outStr << " Max energy of *any type* looper killed: " << fMaxEnergyKilled 161 << " its PDG was " << fMaxEnergyKilledPDG << G4endl; 162 if( fMaxEnergyKilled_NonElectron > 0.0 ) 163 { 164 outStr << " Max energy of non-electron looper killed: " 165 << fMaxEnergyKilled_NonElectron 166 << " its PDG was " << fMaxEnergyKilled_NonElecPDG << G4endl; 167 } 168 if( fMaxEnergySaved > 0.0 ) 169 { 170 outStr << " Max energy of loopers 'saved': " << fMaxEnergySaved << G4endl; 171 outStr << " Sum of energy of loopers 'saved': " 172 << fSumEnergySaved << G4endl; 173 outStr << " Sum of energy of unstable loopers 'saved': " 174 << fSumEnergyUnstableSaved << G4endl; 175 } 176 } 177 else 178 { 179 outStr << " No looping tracks found or killed. " << G4endl; 180 } 181 } 182 183 ////////////////////////////////////////////////////////////////////////// 184 // 185 // Responsibilities: 186 // Find whether the geometry limits the Step, and to what length 187 // Calculate the new value of the safety and return it. 188 // Store the final time, position and momentum. 189 190 G4double G4Transportation::AlongStepGetPhysicalInteractionLength( 191 const G4Track& track, 192 G4double, // previousStepSize 193 G4double currentMinimumStep, G4double& currentSafety, 194 G4GPILSelection* selection) 195 { 196 // Initial actions moved to StartTrack() 197 // -------------------------------------- 198 // Note: in case another process changes touchable handle 199 // it will be necessary to add here (for all steps) 200 // fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 201 202 // GPILSelection is set to defaule value of CandidateForSelection 203 // It is a return value 204 // 205 *selection = CandidateForSelection; 206 207 // Get initial Energy/Momentum of the track 208 // 209 const G4ThreeVector startPosition = track.GetPosition(); 210 const G4ThreeVector startMomentumDir = track.GetMomentumDirection(); 211 212 // The Step Point safety can be limited by other geometries and/or the 213 // assumptions of any process - it's not always the geometrical safety. 214 // We calculate the starting point's isotropic safety here. 215 { 216 const G4double MagSqShift = (startPosition - fPreviousSftOrigin).mag2(); 217 218 if(MagSqShift >= sqr(fPreviousSafety)) 219 currentSafety = 0.0; 220 else 221 currentSafety = fPreviousSafety - std::sqrt(MagSqShift); 222 } 223 224 // Is the particle charged or has it a magnetic moment? 225 // 226 const G4DynamicParticle* pParticle = track.GetDynamicParticle(); 227 228 const G4double particleMass = pParticle->GetMass(); 229 const G4double particleCharge = pParticle->GetCharge(); 230 const G4double kineticEnergy = pParticle->GetKineticEnergy(); 231 232 const G4double magneticMoment = pParticle->GetMagneticMoment(); 233 const G4ThreeVector particleSpin = pParticle->GetPolarization(); 234 235 // There is no need to locate the current volume. It is Done elsewhere: 236 // On track construction 237 // By the tracking, after all AlongStepDoIts, in "Relocation" 238 239 // Check if the particle has a force, EM or gravitational, exerted on it 240 // 241 242 G4bool eligibleEM = 243 (particleCharge != 0.0) || ((magneticMoment != 0.0) && fUseMagneticMoment); 244 G4bool eligibleGrav = (particleMass != 0.0) && fUseGravity; 245 246 fFieldExertedForce = false; 247 248 if(eligibleEM || eligibleGrav) 249 { 250 if(G4FieldManager* fieldMgr = 251 fFieldPropagator->FindAndSetFieldManager(track.GetVolume())) 252 { 253 // User can configure the field Manager for this track 254 fieldMgr->ConfigureForTrack(&track); 255 // Called here to allow a transition from no-field pointer 256 // to finite field (non-zero pointer). 257 258 // If the field manager has no field ptr, the field is zero 259 // by definition ( = there is no field ! ) 260 if(const G4Field* ptrField = fieldMgr->GetDetectorField()) 261 fFieldExertedForce = 262 eligibleEM || (eligibleGrav && ptrField->IsGravityActive()); 263 } 264 } 265 266 G4double geometryStepLength = currentMinimumStep; 267 268 if(currentMinimumStep == 0.0) 269 { 270 fEndPointDistance = 0.0; 271 fGeometryLimitedStep = false; // Old code: = (currentSafety == 0.0); 272 // Changed to avoid problems when setting the step status (now done in AlongStepDoIt) 273 fMomentumChanged = false; 274 fParticleIsLooping = false; 275 fEndGlobalTimeComputed = false; 276 fTransportEndPosition = startPosition; 277 fTransportEndMomentumDir = startMomentumDir; 278 fTransportEndKineticEnergy = kineticEnergy; 279 fTransportEndSpin = particleSpin; 280 } 281 else if(!fFieldExertedForce) 282 { 283 fGeometryLimitedStep = false; 284 if(geometryStepLength > currentSafety || !fShortStepOptimisation) 285 { 286 const G4double linearStepLength = fLinearNavigator->ComputeStep( 287 startPosition, startMomentumDir, currentMinimumStep, currentSafety); 288 289 if(linearStepLength <= currentMinimumStep) 290 { 291 geometryStepLength = linearStepLength; 292 fGeometryLimitedStep = true; 293 } 294 // Remember last safety origin & value. 295 // 296 fPreviousSftOrigin = startPosition; 297 fPreviousSafety = currentSafety; 298 fpSafetyHelper->SetCurrentSafety(currentSafety, startPosition); 299 } 300 301 fEndPointDistance = geometryStepLength; 302 303 fMomentumChanged = false; 304 fParticleIsLooping = false; 305 fEndGlobalTimeComputed = false; 306 fTransportEndPosition = 307 startPosition + geometryStepLength * startMomentumDir; 308 fTransportEndMomentumDir = startMomentumDir; 309 fTransportEndKineticEnergy = kineticEnergy; 310 fTransportEndSpin = particleSpin; 311 } 312 else // A field exerts force 313 { 314 const auto pParticleDef = pParticle->GetDefinition(); 315 const auto particlePDGSpin = pParticleDef->GetPDGSpin(); 316 const auto particlePDGMagM = pParticleDef->GetPDGMagneticMoment(); 317 318 auto equationOfMotion = fFieldPropagator->GetCurrentEquationOfMotion(); 319 320 // The charge can change (dynamic), therefore the use of G4ChargeState 321 // 322 equationOfMotion->SetChargeMomentumMass( 323 G4ChargeState(particleCharge, magneticMoment, particlePDGSpin), 324 pParticle->GetTotalMomentum(), particleMass); 325 326 G4FieldTrack aFieldTrack(startPosition, 327 track.GetGlobalTime(), // Lab. 328 startMomentumDir, kineticEnergy, particleMass, 329 particleCharge, particleSpin, particlePDGMagM, 330 0.0, // Length along track 331 particlePDGSpin); 332 333 // Do the Transport in the field (non recti-linear) 334 // 335 const G4double lengthAlongCurve = fFieldPropagator->ComputeStep( 336 aFieldTrack, currentMinimumStep, currentSafety, track.GetVolume(), 337 kineticEnergy < fThreshold_Important_Energy); 338 339 if(lengthAlongCurve < geometryStepLength) 340 geometryStepLength = lengthAlongCurve; 341 342 // Remember last safety origin & value. 343 // 344 fPreviousSftOrigin = startPosition; 345 fPreviousSafety = currentSafety; 346 fpSafetyHelper->SetCurrentSafety(currentSafety, startPosition); 347 348 fGeometryLimitedStep = fFieldPropagator->IsLastStepInVolume(); 349 // 350 // It is possible that step was reduced in PropagatorInField due to 351 // previous zero steps. To cope with case that reduced step is taken 352 // in full, we must rely on PiF to obtain this value 353 354 G4bool changesEnergy = 355 fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy(); 356 357 fMomentumChanged = true; 358 fParticleIsLooping = fFieldPropagator->IsParticleLooping(); 359 360 fEndGlobalTimeComputed = changesEnergy; 361 fTransportEndPosition = aFieldTrack.GetPosition(); 362 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir(); 363 364 // G4cout << " G4Transport: End of step pMag = " << aFieldTrack.GetMomentum().mag() << G4endl; 365 366 fEndPointDistance = (fTransportEndPosition - startPosition).mag(); 367 368 // Ignore change in energy for fields that conserve energy 369 // This hides the integration error, but gives a better physical answer 370 fTransportEndKineticEnergy = 371 changesEnergy ? aFieldTrack.GetKineticEnergy() : kineticEnergy; 372 fTransportEndSpin = aFieldTrack.GetSpin(); 373 374 if(fEndGlobalTimeComputed) 375 { 376 // If the field can change energy, then the time must be integrated 377 // - so this should have been updated 378 // 379 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight(); 380 381 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); 382 // a cleaner way is to have FieldTrack knowing whether time is updated. 383 } 384 #if defined(G4VERBOSE) || defined(G4DEBUG_TRANSPORT) 385 else 386 { 387 // The energy should be unchanged by field transport, 388 // - so the time changed will be calculated elsewhere 389 // 390 // Check that the integration preserved the energy 391 // - and if not correct this! 392 G4double startEnergy = kineticEnergy; 393 G4double endEnergy = fTransportEndKineticEnergy; 394 395 static G4ThreadLocal G4int no_large_ediff; 396 if(verboseLevel > 1) 397 { 398 if(std::fabs(startEnergy - endEnergy) > perThousand * endEnergy) 399 { 400 static G4ThreadLocal G4int no_warnings = 0, warnModulo = 1, 401 moduloFactor = 10; 402 no_large_ediff++; 403 if((no_large_ediff % warnModulo) == 0) 404 { 405 no_warnings++; 406 std::ostringstream message; 407 message << "Energy change in Step is above 1^-3 relative value. " 408 << G4endl << " Relative change in 'tracking' step = " 409 << std::setw(15) << (endEnergy - startEnergy) / startEnergy 410 << G4endl << " Starting E= " << std::setw(12) 411 << startEnergy / MeV << " MeV " << G4endl 412 << " Ending E= " << std::setw(12) << endEnergy / MeV 413 << " MeV " << G4endl 414 << "Energy has been corrected -- however, review" 415 << " field propagation parameters for accuracy." << G4endl; 416 if((verboseLevel > 2) || (no_warnings < 4) || 417 (no_large_ediff == warnModulo * moduloFactor)) 418 { 419 message << "These include EpsilonStepMax(/Min) in G4FieldManager " 420 << G4endl 421 << "which determine fractional error per step for " 422 "integrated quantities. " 423 << G4endl 424 << "Note also the influence of the permitted number of " 425 "integration steps." 426 << G4endl; 427 } 428 message << "Bad 'endpoint'. Energy change detected and corrected." 429 << G4endl << "Has occurred already " << no_large_ediff 430 << " times."; 431 G4Exception("G4Transportation::AlongStepGetPIL()", "EnergyChange", 432 JustWarning, message); 433 if(no_large_ediff == warnModulo * moduloFactor) 434 { 435 warnModulo *= moduloFactor; 436 } 437 } 438 } 439 } // end of if (verboseLevel) 440 } 441 #endif 442 } 443 444 // Update the safety starting from the end-point, 445 // if it will become negative at the end-point. 446 // 447 if(currentSafety < fEndPointDistance) 448 { 449 if(particleCharge != 0.0) 450 { 451 G4double endSafety = 452 fLinearNavigator->ComputeSafety(fTransportEndPosition); 453 currentSafety = endSafety; 454 fPreviousSftOrigin = fTransportEndPosition; 455 fPreviousSafety = currentSafety; 456 fpSafetyHelper->SetCurrentSafety(currentSafety, fTransportEndPosition); 457 458 // Because the Stepping Manager assumes it is from the start point, 459 // add the StepLength 460 // 461 currentSafety += fEndPointDistance; 462 463 #ifdef G4DEBUG_TRANSPORT 464 G4cout.precision(12); 465 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl; 466 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition 467 << " and it returned safety= " << endSafety << G4endl; 468 G4cout << " Adding endpoint distance " << fEndPointDistance 469 << " to obtain pseudo-safety= " << currentSafety << G4endl; 470 } 471 else 472 { 473 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl; 474 G4cout << " Avoiding call to ComputeSafety : " << G4endl; 475 G4cout << " charge = " << particleCharge << G4endl; 476 G4cout << " mag moment = " << magneticMoment << G4endl; 477 #endif 478 } 479 } 480 481 fFirstStepInVolume = fNewTrack || fLastStepInVolume; 482 fLastStepInVolume = false; 483 fNewTrack = false; 484 485 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume); 486 fParticleChange.ProposeTrueStepLength(geometryStepLength); 487 488 return geometryStepLength; 489 } 490 491 ////////////////////////////////////////////////////////////////////////// 492 // 493 // Initialize ParticleChange (by setting all its members equal 494 // to corresponding members in G4Track) 495 496 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track, 497 const G4Step& stepData ) 498 { 499 #if defined(G4VERBOSE) || defined(G4DEBUG_TRANSPORT) 500 static G4ThreadLocal G4long noCallsASDI=0; 501 noCallsASDI++; 502 #else 503 #define noCallsASDI 0 504 #endif 505 506 if(fGeometryLimitedStep) 507 { 508 stepData.GetPostStepPoint()->SetStepStatus(fGeomBoundary); 509 } 510 511 fParticleChange.Initialize(track) ; 512 513 // Code for specific process 514 // 515 fParticleChange.ProposePosition(fTransportEndPosition) ; 516 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ; 517 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ; 518 fParticleChange.SetMomentumChanged(fMomentumChanged) ; 519 520 fParticleChange.ProposePolarization(fTransportEndSpin); 521 522 G4double deltaTime = 0.0 ; 523 524 // Calculate Lab Time of Flight (ONLY if field Equations used it!) 525 // G4double endTime = fCandidateEndGlobalTime; 526 // G4double delta_time = endTime - startTime; 527 528 G4double startTime = track.GetGlobalTime() ; 529 530 if (!fEndGlobalTimeComputed) 531 { 532 // The time was not integrated .. make the best estimate possible 533 // 534 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity(); 535 G4double stepLength = track.GetStepLength(); 536 537 deltaTime= 0.0; // in case initialVelocity = 0 538 if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; } 539 540 fCandidateEndGlobalTime = startTime + deltaTime ; 541 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ; 542 } 543 else 544 { 545 deltaTime = fCandidateEndGlobalTime - startTime ; 546 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ; 547 } 548 549 550 // Now Correct by Lorentz factor to get delta "proper" Time 551 552 G4double restMass = track.GetDynamicParticle()->GetMass() ; 553 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; 554 555 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; 556 //fParticleChange.ProposeTrueStepLength( track.GetStepLength() ) ; 557 558 // If the particle is caught looping or is stuck (in very difficult 559 // boundaries) in a magnetic field (doing many steps) THEN can kill it ... 560 // 561 if ( fParticleIsLooping ) 562 { 563 G4double endEnergy= fTransportEndKineticEnergy; 564 fNoLooperTrials ++; 565 auto particleType= track.GetDynamicParticle()->GetParticleDefinition(); 566 567 G4bool stable = particleType->GetPDGStable(); 568 G4bool candidateForEnd = (endEnergy < fThreshold_Important_Energy) 569 || (fNoLooperTrials >= fThresholdTrials) ; 570 G4bool unstableAndKillable = !stable && ( fAbandonUnstableTrials != 0); 571 G4bool unstableForEnd = (endEnergy < fThreshold_Important_Energy) 572 && (fNoLooperTrials >= fAbandonUnstableTrials) ; 573 if( (candidateForEnd && stable) || (unstableAndKillable && unstableForEnd) ) 574 { 575 // Kill the looping particle 576 // 577 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 578 G4int particlePDG= particleType->GetPDGEncoding(); 579 const G4int electronPDG= 11; // G4Electron::G4Electron()->GetPDGEncoding(); 580 581 // Simple statistics 582 fSumEnergyKilled += endEnergy; 583 fSumEnerSqKilled += endEnergy * endEnergy; 584 fNumLoopersKilled++; 585 586 if( endEnergy > fMaxEnergyKilled ) { 587 fMaxEnergyKilled = endEnergy; 588 fMaxEnergyKilledPDG = particlePDG; 589 } 590 if( particleType->GetPDGEncoding() != electronPDG ) 591 { 592 fSumEnergyKilled_NonElectron += endEnergy; 593 fSumEnerSqKilled_NonElectron += endEnergy * endEnergy; 594 fNumLoopersKilled_NonElectron++; 595 596 if( endEnergy > fMaxEnergyKilled_NonElectron ) 597 { 598 fMaxEnergyKilled_NonElectron = endEnergy; 599 fMaxEnergyKilled_NonElecPDG = particlePDG; 600 } 601 } 602 603 if( endEnergy > fThreshold_Warning_Energy && ! fSilenceLooperWarnings ) 604 { 605 fpLogger->ReportLoopingTrack( track, stepData, fNoLooperTrials, 606 noCallsASDI, __func__ ); 607 } 608 fNoLooperTrials=0; 609 } 610 else 611 { 612 fMaxEnergySaved = std::max( endEnergy, fMaxEnergySaved); 613 if( fNoLooperTrials == 1 ) { 614 fSumEnergySaved += endEnergy; 615 if ( !stable ) 616 fSumEnergyUnstableSaved += endEnergy; 617 } 618 #ifdef G4VERBOSE 619 if( verboseLevel > 2 && ! fSilenceLooperWarnings ) 620 { 621 G4cout << " " << __func__ 622 << " Particle is looping but is saved ..." << G4endl 623 << " Number of trials = " << fNoLooperTrials << G4endl 624 << " No of calls to = " << noCallsASDI << G4endl; 625 } 626 #endif 627 } 628 } 629 else 630 { 631 fNoLooperTrials=0; 632 } 633 634 // Another (sometimes better way) is to use a user-limit maximum Step size 635 // to alleviate this problem .. 636 637 // Introduce smooth curved trajectories to particle-change 638 // 639 fParticleChange.SetPointerToVectorOfAuxiliaryPoints 640 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() ); 641 642 return &fParticleChange ; 643 } 644 645 ////////////////////////////////////////////////////////////////////////// 646 // 647 // This ensures that the PostStep action is always called, 648 // so that it can do the relocation if it is needed. 649 // 650 651 G4double G4Transportation:: 652 PostStepGetPhysicalInteractionLength( const G4Track&, 653 G4double, // previousStepSize 654 G4ForceCondition* pForceCond ) 655 { 656 fFieldExertedForce = false; // Not known 657 *pForceCond = Forced ; 658 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX 659 } 660 661 ///////////////////////////////////////////////////////////////////////////// 662 663 void G4Transportation::SetTouchableInformation(const G4TouchableHandle& touchable) 664 { 665 const G4VPhysicalVolume* pNewVol = touchable->GetVolume() ; 666 const G4Material* pNewMaterial = 0 ; 667 G4VSensitiveDetector* pNewSensitiveDetector = 0; 668 669 if( pNewVol != 0 ) 670 { 671 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); 672 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector(); 673 } 674 675 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ; 676 fParticleChange.SetSensitiveDetectorInTouchable( pNewSensitiveDetector ) ; 677 // temporarily until Get/Set Material of ParticleChange, 678 // and StepPoint can be made const. 679 680 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; 681 if( pNewVol != 0 ) 682 { 683 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); 684 } 685 686 if ( pNewVol!=0 && pNewMaterialCutsCouple!=0 687 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial ) 688 { 689 // for parametrized volume 690 // 691 pNewMaterialCutsCouple = 692 G4ProductionCutsTable::GetProductionCutsTable() 693 ->GetMaterialCutsCouple(pNewMaterial, 694 pNewMaterialCutsCouple->GetProductionCuts()); 695 } 696 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple ); 697 698 // Set the touchable in ParticleChange 699 // this must always be done because the particle change always 700 // uses this value to overwrite the current touchable pointer. 701 // 702 fParticleChange.SetTouchableHandle(touchable) ; 703 } 704 705 ///////////////////////////////////////////////////////////////////////////// 706 // 707 708 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track, 709 const G4Step& ) 710 { 711 G4TouchableHandle retCurrentTouchable ; // The one to return 712 G4bool isLastStep= false; 713 714 // Initialize ParticleChange (by setting all its members equal 715 // to corresponding members in G4Track) 716 // fParticleChange.Initialize(track) ; // To initialise TouchableChange 717 718 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; 719 720 // If the Step was determined by the volume boundary, 721 // logically relocate the particle 722 723 if(fGeometryLimitedStep) 724 { 725 // fCurrentTouchable will now become the previous touchable, 726 // and what was the previous will be freed. 727 // (Needed because the preStepPoint can point to the previous touchable) 728 729 fLinearNavigator->SetGeometricallyLimitedStep() ; 730 fLinearNavigator-> 731 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), 732 track.GetMomentumDirection(), 733 fCurrentTouchableHandle, 734 true ) ; 735 // Check whether the particle is out of the world volume 736 // If so it has exited and must be killed. 737 // 738 if( fCurrentTouchableHandle->GetVolume() == 0 ) 739 { 740 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 741 } 742 retCurrentTouchable = fCurrentTouchableHandle ; 743 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; 744 745 // Update the Step flag which identifies the Last Step in a volume 746 if( !fFieldExertedForce ) 747 isLastStep = fLinearNavigator->ExitedMotherVolume() 748 || fLinearNavigator->EnteredDaughterVolume() ; 749 else 750 isLastStep = fFieldPropagator->IsLastStepInVolume(); 751 } 752 else // fGeometryLimitedStep is false 753 { 754 // This serves only to move the Navigator's location 755 // 756 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ; 757 758 // The value of the track's current Touchable is retained. 759 // (and it must be correct because we must use it below to 760 // overwrite the (unset) one in particle change) 761 // It must be fCurrentTouchable too ?? 762 // 763 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; 764 retCurrentTouchable = track.GetTouchableHandle() ; 765 766 isLastStep= false; 767 } // endif ( fGeometryLimitedStep ) 768 fLastStepInVolume= isLastStep; 769 770 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume); 771 fParticleChange.ProposeLastStepInVolume(isLastStep); 772 773 SetTouchableInformation(retCurrentTouchable); 774 775 return &fParticleChange ; 776 } 777 778 ///////////////////////////////////////////////////////////////////////////// 779 // New method takes over the responsibility to reset the state of 780 // G4Transportation object at the start of a new track or the resumption 781 // of a suspended track. 782 // 783 784 void 785 G4Transportation::StartTracking(G4Track* aTrack) 786 { 787 G4VProcess::StartTracking(aTrack); 788 fNewTrack= true; 789 fFirstStepInVolume= true; 790 fLastStepInVolume= false; 791 792 // The actions here are those that were taken in AlongStepGPIL 793 // when track.GetCurrentStepNumber()==1 794 795 // Whether field exists should be determined at run level -- TODO 796 G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance(); 797 fAnyFieldExists = fieldMgrStore->size() > 0; 798 799 // reset safety value and center 800 // 801 fPreviousSafety = 0.0 ; 802 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; 803 804 // reset looping counter -- for motion in field 805 fNoLooperTrials= 0; 806 // Must clear this state .. else it depends on last track's value 807 // --> a better solution would set this from state of suspended track TODO ? 808 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. } 809 810 // ChordFinder reset internal state 811 // 812 if( fFieldPropagator && fAnyFieldExists ) 813 { 814 fFieldPropagator->ClearPropagatorState(); 815 // Resets all state of field propagator class (ONLY) including safety 816 // values (in case of overlaps and to wipe for first track). 817 } 818 819 // Make sure to clear the chord finders of all fields (i.e. managers) 820 // 821 fieldMgrStore->ClearAllChordFindersState(); 822 823 // Update the current touchable handle (from the track's) 824 // 825 fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 826 827 // Inform field propagator of new track 828 // 829 fFieldPropagator->PrepareNewTrack(); 830 } 831 832 ///////////////////////////////////////////////////////////////////////////// 833 // 834 835 G4bool G4Transportation::EnableMagneticMoment(G4bool useMoment) 836 { 837 G4bool lastValue= fUseMagneticMoment; 838 fUseMagneticMoment= useMoment; 839 return lastValue; 840 } 841 842 ///////////////////////////////////////////////////////////////////////////// 843 // 844 845 G4bool G4Transportation::EnableGravity(G4bool useGravity) 846 { 847 G4bool lastValue= fUseGravity; 848 fUseGravity= useGravity; 849 return lastValue; 850 } 851 852 ///////////////////////////////////////////////////////////////////////////// 853 // 854 // Supress (or not) warnings about 'looping' particles 855 856 void G4Transportation::SetSilenceLooperWarnings( G4bool val) 857 { 858 fSilenceLooperWarnings= val; // Flag to *Supress* all 'looper' warnings 859 } 860 861 ///////////////////////////////////////////////////////////////////////////// 862 // 863 G4bool G4Transportation::GetSilenceLooperWarnings() 864 { 865 return fSilenceLooperWarnings; 866 } 867 868 869 ///////////////////////////////////////////////////////////////////////////// 870 // 871 void G4Transportation::SetHighLooperThresholds() 872 { 873 // Restores the old high values -- potentially appropriate for energy-frontier 874 // HEP experiments. 875 // Caution: All tracks with E < 100 MeV that are found to loop are 876 SetThresholdWarningEnergy( 100.0 * CLHEP::MeV ); // Warn above this energy 877 SetThresholdImportantEnergy( 250.0 * CLHEP::MeV ); // Extra trial above this En 878 879 G4int maxTrials = 10; 880 SetThresholdTrials( maxTrials ); 881 882 PushThresholdsToLogger(); // Again, to be sure 883 if( verboseLevel ) ReportLooperThresholds(); 884 } 885 886 ///////////////////////////////////////////////////////////////////////////// 887 void G4Transportation::SetLowLooperThresholds() // Values for low-E applications 888 { 889 // These values were the default in Geant4 10.5 - beta 890 SetThresholdWarningEnergy( 1.0 * CLHEP::keV ); // Warn above this En 891 SetThresholdImportantEnergy( 1.0 * CLHEP::MeV ); // Extra trials above it 892 893 G4int maxTrials = 30; // A new value - was 10 894 SetThresholdTrials( maxTrials ); 895 896 PushThresholdsToLogger(); // Again, to be sure 897 if( verboseLevel ) ReportLooperThresholds(); 898 } 899 900 ///////////////////////////////////////////////////////////////////////////// 901 // 902 void 903 G4Transportation::ReportMissingLogger( const char* methodName ) 904 { 905 const char* message= "Logger object missing from G4Transportation object"; 906 G4String classAndMethod= G4String("G4Transportation") + G4String( methodName ); 907 G4Exception(classAndMethod, "Missing Logger", JustWarning, message); 908 } 909 910 911 ///////////////////////////////////////////////////////////////////////////// 912 // 913 void 914 G4Transportation::ReportLooperThresholds() 915 { 916 PushThresholdsToLogger(); // To be absolutely certain they are in sync 917 fpLogger->ReportLooperThresholds("G4Transportation"); 918 } 919 920 ///////////////////////////////////////////////////////////////////////////// 921 // 922 void G4Transportation::ProcessDescription(std::ostream& outStr) const 923 924 // StreamInfo(std::ostream& out, const G4ParticleDefinition& part, G4bool rst) const 925 926 { 927 G4String indent = " "; // : ""); 928 G4long oldPrec= outStr.precision(6); 929 // outStr << std::setprecision(6); 930 outStr << G4endl << indent << GetProcessName() << ": "; 931 932 outStr << " Parameters for looping particles: " << G4endl 933 << " warning-E = " << fThreshold_Warning_Energy / CLHEP::MeV << " MeV " << G4endl 934 << " important E = " << fThreshold_Important_Energy / CLHEP::MeV << " MeV " << G4endl 935 << " thresholdTrials " << fThresholdTrials << G4endl; 936 outStr.precision(oldPrec); 937 } 938