Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4CoupledTransportation.cc
이 파일의 문서화 페이지로 가기
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 // $Id: G4CoupledTransportation.cc 110805 2018-06-15 06:52:15Z gcosmo $
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class implementation
31 // =======================================================================
32 // Modified:
33 // 10 Jan 2015, M.Kelsey: Use G4DynamicParticle mass, NOT PDGMass
34 // 13 May 2006, J. Apostolakis: Revised for parallel navigation (PathFinder)
35 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
36 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
37 // 21 June 2003, J.Apostolakis: Calling field manager with
38 // track, to enable it to configure its accuracy
39 // 13 May 2003, J.Apostolakis: Zero field areas now taken into
40 // account correclty in all cases (thanks to W Pokorski).
41 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
42 // correction for spin tracking
43 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack
44 // 22 Sept 2000, V.Grichine: update of Kinetic Energy
45 // Created: 19 March 1997, J. Apostolakis
46 // =======================================================================
47 
51 
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4ProductionCutsTable.hh"
55 #include "G4ParticleTable.hh"
56 #include "G4ChordFinder.hh"
57 #include "G4Field.hh"
58 #include "G4FieldTrack.hh"
59 #include "G4FieldManagerStore.hh"
60 
61 #include "G4Transportation.hh" // In order to use fUseMagneticMoment
62 
64 
66 // This mode must apply to all threads
67 
70 //
71 // Constructor
72 
74  : G4VProcess( G4String("CoupledTransportation"), fTransportation ),
75  fTransportEndPosition(0.0, 0.0, 0.0),
76  fTransportEndMomentumDir(0.0, 0.0, 0.0),
77  fTransportEndKineticEnergy(0.0),
78  fTransportEndSpin(0.0, 0.0, 0.0), // fTransportEndPolarization(0.0, 0.0, 0.0),
79  fMomentumChanged(false),
80  fEndGlobalTimeComputed(false),
81  fCandidateEndGlobalTime(0.0),
82  fParticleIsLooping( false ),
83  fNewTrack( true ),
84  fPreviousSftOrigin( 0.,0.,0. ),
85  fPreviousMassSafety( 0.0 ),
86  fPreviousFullSafety( 0.0 ),
87  fMassGeometryLimitedStep( false ),
88  fAnyGeometryLimitedStep( false ),
89  fEndpointDistance( -1.0 ),
90  fThreshold_Warning_Energy( 1.0 * CLHEP::kiloelectronvolt ), // Old value: 100 * MeV ),
91  fThreshold_Important_Energy( 1.0 * MeV ), // Old value: 250 * MeV ),
92  fThresholdTrials( 10 ),
93  fNoLooperTrials( 0 ),
94  fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
95  fFirstStepInMassVolume( true ),
96  fFirstStepInAnyVolume( true )
97 {
98  SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION));
99  SetVerboseLevel(verbosity);
100 
101  G4TransportationManager* transportMgr ;
102 
104 
105  fMassNavigator = transportMgr->GetNavigatorForTracking() ;
106  fFieldPropagator = transportMgr->GetPropagatorInField() ;
107  // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
108  fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator );
109  if( verboseLevel > 0 )
110  {
111  G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
112  G4cout << " Verbose level is " << verboseLevel << G4endl;
113  G4cout << " Navigator Id obtained in G4CoupledTransportation constructor "
114  << fNavigatorId << G4endl;
115  G4cout << " Reports First/Last in "
116  << (fSignifyStepInAnyVolume ? " any " : " mass " ) << " geometry " << G4endl;
117  }
119  fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
120 
121  fpLogger = new G4TransportationLogger("G4Transportation", verbosity);
123  GetThresholdTrials() );
124 
125  // Following assignment is to fix small memory leak from simple use of 'new'
126  static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
127  if ( !pNullTouchableHandle) { pNullTouchableHandle = new G4TouchableHandle; }
128  fCurrentTouchableHandle = *pNullTouchableHandle;
129  // Points to (G4VTouchable*) 0
130 
131  G4FieldManager *globalFieldMgr= transportMgr->GetFieldManager();
132  fGlobalFieldExists= globalFieldMgr ? globalFieldMgr->GetDetectorField() : 0 ;
133 }
134 
136 
138 {
139  delete fpLogger;
140 
141  if( (verboseLevel > 0) || (fSumEnergyKilled > 0.0 ) )
142  {
143  G4cout << " G4CoupledTransportation: Statistics for looping particles "
144  << G4endl;
145  G4cout << " Sum of energy of loopers killed: "
146  << fSumEnergyKilled / MeV << " MeV " << G4endl;
147  G4cout << " Max energy of loopers killed: "
148  << fMaxEnergyKilled / MeV << " MeV " << G4endl;
149  }
150 }
151 
153 //
154 // Responsibilities:
155 // Find whether the geometry limits the Step, and to what length
156 // Calculate the new value of the safety and return it.
157 // Store the final time, position and momentum.
158 
161  G4double, // previousStepSize
162  G4double currentMinimumStep,
163  G4double& proposedSafetyForStart,
164  G4GPILSelection* selection )
165 {
166  G4double geometryStepLength;
167  G4double startMassSafety= 0.0; // estimated safety for start point (mass geometry)
168  G4double startFullSafety= 0.0; // estimated safety for start point (all geometries)
169  G4double safetyProposal= -1.0; // local copy of proposal
170 
171  G4ThreeVector EndUnitMomentum ;
172  G4double lengthAlongCurve=0.0 ;
173 
174  fParticleIsLooping = false ;
175 
176  // Initial actions moved to StartTrack()
177  // --------------------------------------
178  // Note: in case another process changes touchable handle
179  // it will be necessary to add here (for all steps)
180  // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
181 
182  // GPILSelection is set to defaule value of CandidateForSelection
183  // It is a return value
184  //
185  *selection = CandidateForSelection ;
186 
189 
190 #ifdef G4DEBUG_TRANSPORT
191  G4cout << " CoupledTransport::AlongStep GPIL: "
192  << " 1st-step: any= " <<fFirstStepInAnyVolume << " ( geom= " << fAnyGeometryLimitedStep << " ) "
193  << " mass= " << fFirstStepInMassVolume << " ( geom= " << fMassGeometryLimitedStep << " ) "
194  << " newTrack= " << fNewTrack << G4endl;
195  // G4cout << " fLastStep-in-Vol= " << fLastStepInVolume << G4endl;
196 #endif
197 
198  // fLastStepInVolume= false;
199  fNewTrack = false;
200 
201  // Get initial Energy/Momentum of the track
202  //
203  const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
204  const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
205  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
206  G4ThreeVector startPosition = track.GetPosition() ;
207  G4VPhysicalVolume* currentVolume= track.GetVolume();
208 
209 #ifdef G4DEBUG_TRANSPORT
210  if( verboseLevel > 1 )
211  {
212  G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume "
213  << currentVolume->GetName() << G4endl;
214  }
215 #endif
216  // G4double theTime = track.GetGlobalTime() ;
217 
218  // The Step Point safety can be limited by other geometries and/or the
219  // assumptions of any process - it's not always the geometrical safety.
220  // We calculate the starting point's isotropic safety here.
221  //
222  G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
223  G4double MagSqShift = OriginShift.mag2() ;
224  startMassSafety = 0.0;
225  startFullSafety= 0.0;
226 
227  // Recall that FullSafety <= MassSafety
228  // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
229  if( MagSqShift < sqr(fPreviousFullSafety) ) // Revision proposed by Alex H, 2 Oct 07
230  {
231  G4double mag_shift= std::sqrt(MagSqShift);
232  startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0);
233  startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
234  // Need to be consistent between full safety with Mass safety
235  // in order reproduce results in simple case --> use same calculation method
236 
237  // Only compute full safety if massSafety > 0. Else it remains 0
238  // startFullSafety = fPathFinder->ComputeSafety( startPosition );
239  }
240 
241  // Is the particle charged or has it a magnetic moment?
242  //
243  G4double particleCharge = pParticle->GetCharge() ;
244  G4double magneticMoment = pParticle->GetMagneticMoment() ;
245  G4double restMass = pParticle->GetMass() ;
246 
247  fMassGeometryLimitedStep = false ; // Set default - alex
248  fAnyGeometryLimitedStep = false;
249 
250  // fEndGlobalTimeComputed = false ;
251 
252  // There is no need to locate the current volume. It is Done elsewhere:
253  // On track construction
254  // By the tracking, after all AlongStepDoIts, in "Relocation"
255 
256  // Check if the particle has a force, EM or gravitational, exerted on it
257  //
258  G4FieldManager* fieldMgr=0;
259  G4bool fieldExertsForce = false ;
260 
261  G4bool gravityOn = false;
262  const G4Field* ptrField= 0;
263 
264  fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
265  if( fieldMgr != 0 )
266  {
267  // Message the field Manager, to configure it for this track
268  fieldMgr->ConfigureForTrack( &track );
269  // Here it can transition from a null field-ptr to a finite field
270 
271  // If the field manager has no field ptr, the field is zero
272  // by definition ( = there is no field ! )
273  ptrField= fieldMgr->GetDetectorField();
274 
275  if( ptrField != 0)
276  {
277  gravityOn= ptrField->IsGravityActive();
278  if( (particleCharge != 0.0)
279  || (fUseMagneticMoment && (magneticMoment != 0.0) )
280  || (gravityOn && (restMass != 0.0)) )
281  {
282  fieldExertsForce = true;
283  }
284  }
285  }
286  G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
287 
288  if( fieldExertsForce )
289  {
290  auto equationOfMotion= fFieldPropagator->GetCurrentEquationOfMotion();
291 
292  G4ChargeState chargeState(particleCharge, // Charge can change (dynamic)
293  magneticMoment,
294  pParticleDef->GetPDGSpin() );
295  // For insurance, could set it again
296  // chargeState.SetPDGSpin( pParticleDef->GetPDGSpin() );
297 
298  if( equationOfMotion )
299  {
300  equationOfMotion->SetChargeMomentumMass( chargeState,
301  momentumMagnitude,
302  restMass );
303  }
304 #ifdef G4DEBUG_TRANSPORT
305  else
306  {
307  G4cerr << " ERROR in G4CoupledTransportation> "
308  << "Cannot find valid Equation of motion: " << G4endl;
309  << " Unable to pass Charge, Momentum and Mass " << G4endl;
310  }
311 #endif
312  }
313 
314  G4ThreeVector polarizationVec = track.GetPolarization() ;
315  G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
316  track.GetGlobalTime(), // Lab.
317  track.GetMomentumDirection(),
318  track.GetKineticEnergy(),
319  restMass,
320  particleCharge,
321  polarizationVec,
322  pParticleDef->GetPDGMagneticMoment(),
323  0.0, // Length along track
324  pParticleDef->GetPDGSpin() ) ;
325  G4int stepNo= track.GetCurrentStepNumber();
326 
327  ELimited limitedStep;
328  G4FieldTrack endTrackState('a'); // Default values
329 
330  fMassGeometryLimitedStep = false ; // default
331  fAnyGeometryLimitedStep = false ;
332  if( currentMinimumStep > 0 )
333  {
334  G4double newMassSafety= 0.0; // temp. for recalculation
335 
336  // Do the Transport in the field (non recti-linear)
337  //
338  lengthAlongCurve = fPathFinder->ComputeStep( aFieldTrack,
339  currentMinimumStep,
340  fNavigatorId,
341  stepNo,
342  newMassSafety,
343  limitedStep,
344  endTrackState,
345  currentVolume ) ;
346 
347  G4double newFullSafety= fPathFinder->GetCurrentSafety();
348  // this was estimated already in step above
349 
350  if( limitedStep == kUnique || limitedStep == kSharedTransport )
351  {
352  fMassGeometryLimitedStep = true ;
353  }
354 
356 
357 #ifdef G4DEBUG_TRANSPORT
359  {
360  std::ostringstream message;
361  message << " ERROR in determining geometries limiting the step" << G4endl;
362  message << " Limiting: mass=" << fMassGeometryLimitedStep
363  << " any= " << fAnyGeometryLimitedStep << G4endl;
364  message << "Incompatible conditions - by which geometries was it limited ?"<<G4endl;
365  G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()",
366  "PathFinderConfused", FatalException, message);
367  }
368 #endif
369 
370  geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep);
371 
372  // Momentum: Magnitude and direction can be changed too now ...
373  //
374  fMomentumChanged = true ;
375  fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
376 
377  // Remember last safety origin & value.
378  fPreviousSftOrigin = startPosition ;
379  fPreviousMassSafety = newMassSafety ;
380  fPreviousFullSafety = newFullSafety ;
381  // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
382 
383 #ifdef G4DEBUG_TRANSPORT
384  if( verboseLevel > 1 )
385  {
386  G4cout << "G4Transport:CompStep> "
387  << " called the pathfinder for a new step at " << startPosition
388  << " and obtained step = " << lengthAlongCurve << G4endl;
389  G4cout << " New safety (preStep) = " << newMassSafety
390  << " versus precalculated = " << startMassSafety << G4endl;
391  }
392 #endif
393 
394  // Store as best estimate value
395  startMassSafety = newMassSafety ;
396  startFullSafety = newFullSafety ;
397 
398  // Get the End-Position and End-Momentum (Dir-ection)
399  fTransportEndPosition = endTrackState.GetPosition() ;
400  fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ;
401  }
402  else
403  {
404  geometryStepLength = lengthAlongCurve= 0.0 ;
405  fMomentumChanged = false ;
406  // fMassGeometryLimitedStep = false ; // --- ???
407  // fAnyGeometryLimitedStep = true;
410 
411  fTransportEndPosition = startPosition;
412 
413  endTrackState= aFieldTrack; // Ensures that time is updated
414 
415  // If the step length requested is 0, and we are on a boundary
416  // then a boundary will also limit the step.
417  if( startMassSafety == 0.0 )
418  {
419  fMassGeometryLimitedStep = true ;
421  }
422  // TODO: Add explicit logical status for being at a boundary
423  }
424  // G4FieldTrack aTrackState(endTrackState);
425 
426  if( !fieldExertsForce )
427  {
428  fParticleIsLooping = false ;
429  fMomentumChanged = false ;
430  fEndGlobalTimeComputed = false ;
431  }
432  else
433  {
434 
435 #ifdef G4DEBUG_TRANSPORT
436  if( verboseLevel > 1 )
437  {
438  G4cout << " G4CT::CS End Position = "
440  G4cout << " G4CT::CS End Direction = "
442  }
443 #endif
445  {
446  // If the field can change energy, then the time must be integrated
447  // - so this should have been updated
448  //
450  fEndGlobalTimeComputed = true;
451 
452  // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
453  // a cleaner way is to have FieldTrack knowing whether time is updated
454  }
455  else
456  {
457  // The energy should be unchanged by field transport,
458  // - so the time changed will be calculated elsewhere
459  //
460  fEndGlobalTimeComputed = false;
461 
462  // Check that the integration preserved the energy
463  // - and if not correct this!
464  G4double startEnergy= track.GetKineticEnergy();
466 
467  static G4ThreadLocal G4int no_inexact_steps=0; // , no_large_ediff;
468  G4double absEdiff = std::fabs(startEnergy- endEnergy);
469  if( absEdiff > perMillion * endEnergy )
470  {
471  no_inexact_steps++;
472  // Possible statistics keeping here ...
473  }
474 #ifdef G4VERBOSE
475  if( (verboseLevel > 1) && ( absEdiff > perThousand * endEnergy) )
476  {
477  ReportInexactEnergy(startEnergy, endEnergy);
478  } // end of if (verboseLevel)
479 #endif
480  // Correct the energy for fields that conserve it
481  // This - hides the integration error
482  // - but gives a better physical answer
484  }
485  }
486 
487  fEndpointDistance = (fTransportEndPosition - startPosition).mag() ;
489 
490  fTransportEndSpin = endTrackState.GetSpin();
491 
492  // Calculate the safety
493 
494  safetyProposal= startFullSafety; // used to be startMassSafety
495  // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06
496 
497  // Update safety for the end-point, if becomes negative at the end-point.
498 
499  if( (startFullSafety < fEndpointDistance )
500  && ( particleCharge != 0.0 ) ) // Only needed to prepare for Mult Scat.
501  // && !fAnyGeometryLimitedStep ) // To-Try: No safety update if at boundary
502  {
503  G4double endFullSafety =
505  // Expected mission -- only mass geometry's safety
506  // fMassNavigator->ComputeSafety( fTransportEndPosition) ;
507  // Yet discrete processes only have poststep -- and this cannot
508  // currently revise the safety
509  // ==> so we use the all-geometry safety as a precaution
510 
512  // Pushing safety to Helper avoids recalculation at this point
513 
514  G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value
515  G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt);
516  // Retrieves the mass value from PathFinder (it calculated it)
517 
518  fPreviousMassSafety = endMassSafety ;
519  fPreviousFullSafety = endFullSafety;
520  fPreviousSftOrigin = fTransportEndPosition ;
521 
522  // The convention (Stepping Manager's) is safety from the start point
523  //
524  safetyProposal = endFullSafety + fEndpointDistance;
525  // --> was endMassSafety
526  // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06
527 
528 #ifdef G4DEBUG_TRANSPORT
529  G4int prec= G4cout.precision(12) ;
530  G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
531  G4cout << " Revised Safety at endpoint " << fTransportEndPosition
532  << " give safety values: Mass= " << endMassSafety
533  << " All= " << endFullSafety << G4endl ;
534  G4cout << " Adding endpoint distance " << fEndpointDistance
535  << " to obtain pseudo-safety= " << safetyProposal << G4endl ;
536  G4cout.precision(prec);
537  }
538  else
539  {
540  G4int prec= G4cout.precision(12) ;
541  G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
542  G4cout << " Quick Safety estimate at endpoint " << fTransportEndPosition
543  << " gives safety endpoint value = " << startFullSafety - fEndpointDistance
544  << " using start-point value " << startFullSafety
545  << " and endpointDistance " << fEndpointDistance << G4endl;
546  G4cout.precision(prec);
547 #endif
548  }
549 
550  proposedSafetyForStart= safetyProposal;
551  fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
552 
553  return geometryStepLength ;
554 }
555 
557 
560  const G4Step& stepData )
561 {
562  static G4ThreadLocal G4long noCallsCT_ASDI=0;
563  const char *methodName= "AlongStepDoIt";
564 
565  noCallsCT_ASDI++;
566 
567  fParticleChange.Initialize(track) ;
568  // sets all its members to the value of corresponding members in G4Track
569 
570  // Code specific for Transport
571  //
573  // G4cout << " G4CoupledTransportation::AlongStepDoIt"
574  // << " proposes position = " << fTransportEndPosition
575  // << " and end momentum direction = " << fTransportEndMomentumDir << G4endl;
579 
581 
582  G4double deltaTime = 0.0 ;
583 
584  // Calculate Lab Time of Flight (ONLY if field Equations used it!)
585  // G4double endTime = fCandidateEndGlobalTime;
586  // G4double delta_time = endTime - startTime;
587 
588  G4double startTime = track.GetGlobalTime() ;
589 
591  {
592  G4double finalInverseVel= DBL_MAX, initialInverseVel=DBL_MAX;
593 
594  // The time was not integrated .. make the best estimate possible
595  //
596  G4double finalVelocity = track.GetVelocity() ;
597  if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; }
598  G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
599  if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; }
600  G4double stepLength = track.GetStepLength() ;
601 
602  if (finalVelocity > 0.0)
603  {
604  // deltaTime = stepLength/finalVelocity ;
605  G4double meanInverseVelocity = 0.5 * ( initialInverseVel + finalInverseVel );
606  deltaTime = stepLength * meanInverseVelocity ;
607  // G4cout << " dt = s * mean(1/v) , with " << " s = " << stepLength
608  // << " mean(1/v)= " << meanInverseVelocity << G4endl;
609  }
610  else
611  {
612  deltaTime = stepLength * initialInverseVel ;
613  // G4cout << " dt = s / initV " << " s = " << stepLength
614  // << " 1 / initV= " << initialInverseVel << G4endl;
615  } // Could do with better estimate for final step (finalVelocity = 0) ?
616 
617  fCandidateEndGlobalTime = startTime + deltaTime ;
618  fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
619 
620  // G4cout << " Calculated global time from start = " << startTime << " and "
621  // << " delta time = " << deltaTime << G4endl;
622  }
623  else
624  {
625  deltaTime = fCandidateEndGlobalTime - startTime ;
627  // G4cout << " Calculated global time from candidate end time = "
628  // << fCandidateEndGlobalTime << " and start time = " << startTime << G4endl;
629  }
630 
631  // G4cout << " G4CoupledTransportation::AlongStepDoIt "
632  // << " flag whether computed time = " << fEndGlobalTimeComputed << " and "
633  // << " is proposes end time " << fCandidateEndGlobalTime << G4endl;
634 
635  // Now Correct by Lorentz factor to get "proper" deltaTime
636 
637  G4double restMass = track.GetDynamicParticle()->GetMass() ;
638  G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
639 
640  fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
641  // fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
642 
643  // If the particle is caught looping or is stuck (in very difficult
644  // boundaries) in a magnetic field (doing many steps) THEN this kills it ...
645  //
646  if ( fParticleIsLooping )
647  {
649 
650  if( (endEnergy < fThreshold_Important_Energy)
652  {
653  // Kill the looping particle
654  //
656 
657  // 'Bare' statistics
658  fSumEnergyKilled += endEnergy;
659  if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
660 
661  if( endEnergy > fThreshold_Warning_Energy )
662  {
663  fpLogger->ReportLoopingTrack( track, stepData, fNoLooperTrials, noCallsCT_ASDI, methodName );
664  }
665  fNoLooperTrials=0;
666  }
667  else
668  {
669  fNoLooperTrials ++;
670 #ifdef G4VERBOSE
671  if( verboseLevel > 2 )
672  {
673  G4cout << " ** G4CoupledTransportation::AlongStepDoIt(): Particle looping - " << G4endl
674  << " Number of trials (this track) = " << fNoLooperTrials
675  << G4endl
676  << " Steps by this track: " << track.GetCurrentStepNumber()
677  << G4endl
678  << " Total no of calls to this method (all tracks) = "
679  << noCallsCT_ASDI << G4endl;
680  }
681 #endif
682  }
683  }
684  else
685  {
686  fNoLooperTrials=0;
687  }
688 
689  // Another (sometimes better way) is to use a user-limit maximum Step size
690  // to alleviate this problem ..
691 
692  // Add smooth curved trajectories to particle-change
693  //
694  // fParticleChange.SetPointerToVectorOfAuxiliaryPoints
695  // (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
696 
697  return &fParticleChange ;
698 }
699 
701 //
702 // This ensures that the PostStep action is always called,
703 // so that it can do the relocation if it is needed.
704 //
705 
708  G4double, // previousStepSize
709  G4ForceCondition* pForceCond )
710 {
711  // Must act as PostStep action -- to relocate particle
712  *pForceCond = Forced ;
713  return DBL_MAX ;
714 }
715 
717 
719 ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector,
720  const G4String& Quantity )
721 {
722  G4ThreeVector moveVec = ( NewVector - OldVector );
723 
724  G4cerr << G4endl
725  << "**************************************************************"
726  << G4endl;
727  G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
728  << " and value from Track in PostStepDoIt. " << G4endl
729  << "Change of " << Quantity << " is " << moveVec.mag() / mm << " mm long, "
730  << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
731  << "Endpoint of ComputeStep was " << OldVector
732  << " and current position to locate is " << NewVector << G4endl;
733 }
734 
736 
738  const G4Step& )
739 {
740  G4TouchableHandle retCurrentTouchable ; // The one to return
741 
742  // Initialize ParticleChange (by setting all its members equal
743  // to corresponding members in G4Track)
744  // fParticleChange.Initialize(track) ; // To initialise TouchableChange
745 
747 
749  {
751  }
752  else
753  {
755  }
756 
757  // Check that the end position and direction are preserved
758  // since call to AlongStepDoIt
759 
760 #ifdef G4DEBUG_TRANSPORT
761  if( ( verboseLevel > 0 )
762  && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) )
763  {
764  ReportMove( track.GetPosition(), fTransportEndPosition, "End of Step Position" );
765  G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl;
766  }
767 
768  // If the Step was determined by the volume boundary, relocate the particle
769  // The pathFinder will know that the geometry limited the step (!?)
770 
771  if( verboseLevel > 0 )
772  {
773  G4cout << " Calling PathFinder::Locate() from "
774  << " G4CoupledTransportation::PostStepDoIt " << G4endl;
775  G4cout << " fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep << G4endl;
776 
777  }
778 #endif
779 
781  {
782  fPathFinder->Locate( track.GetPosition(),
783  track.GetMomentumDirection(),
784  true);
785 
786  // fCurrentTouchable will now become the previous touchable,
787  // and what was the previous will be freed.
788  // (Needed because the preStepPoint can point to the previous touchable)
789 
792 
793 #ifdef G4DEBUG_TRANSPORT
794  if( verboseLevel > 0 )
795  {
796  G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = "
797  << fNavigatorId << G4endl;
798  }
799  if( verboseLevel > 1 )
800  {
802  G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = "
803  << vol;
804  if( vol ) { G4cout << "Name=" << vol->GetName(); }
805  G4cout << G4endl;
806  }
807 #endif
808 
809  // Check whether the particle is out of the world volume
810  // If so it has exited and must be killed.
811  //
812  if( fCurrentTouchableHandle->GetVolume() == 0 )
813  {
815  }
816  retCurrentTouchable = fCurrentTouchableHandle ;
817  // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
818  }
819  else // fAnyGeometryLimitedStep is false
820  {
821 #ifdef G4DEBUG_TRANSPORT
822  if( verboseLevel > 1 )
823  {
824  G4cout << "G4CoupledTransportation::PostStepDoIt -- "
825  << " fAnyGeometryLimitedStep = " << fAnyGeometryLimitedStep
826  << " must be false " << G4endl;
827  }
828 #endif
829  // This serves only to move each of the Navigator's location
830  //
831  // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
832 
833  fPathFinder->ReLocate( track.GetPosition() );
834  // track.GetMomentumDirection() );
835 
836  // Keep the value of the track's current Touchable is retained,
837  // and use it to overwrite the (unset) one in particle change.
838  // Expect this must be fCurrentTouchable too
839  // - could it be different, eg at the start of a step ?
840  //
841  retCurrentTouchable = track.GetTouchableHandle() ;
842  // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
843  } // endif ( fAnyGeometryLimitedStep )
844 
845 #ifdef G4DEBUG_NAVIGATION
846  G4cout << " CoupledTransport::AlongStep GPIL: "
847  << " last-step: any= " << fAnyGeometryLimitedStep << " . ..... x . "
848  << " mass= " << fMassGeometryLimitedStep
849  << G4endl;
850 #endif
851 
854  else
856 
857  const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
858  const G4Material* pNewMaterial = 0 ;
859  const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
860 
861  if( pNewVol != 0 )
862  {
863  pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
864  pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
865  }
866 
867  // ( const_cast<G4Material *> pNewMaterial ) ;
868  // ( const_cast<G4VSensitiveDetetor *> pNewSensitiveDetector) ;
869 
872  // "temporarily" until Get/Set Material of ParticleChange,
873  // and StepPoint can be made const.
874 
875  const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
876  if( pNewVol != 0 )
877  {
878  pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
879  if( pNewMaterialCutsCouple!=0
880  && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
881  {
882  // for parametrized volume
883  //
884  pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()
885  ->GetMaterialCutsCouple(pNewMaterial,
886  pNewMaterialCutsCouple->GetProductionCuts());
887  }
888  }
889  fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
890 
891  // Must always set the touchable in ParticleChange, whether relocated or not
892  //
893  fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
894 
895  return &fParticleChange ;
896 }
897 
899 // New method takes over the responsibility to reset the state of
900 // G4CoupledTransportation object:
901 // - at the start of a new track, and
902 // - on the resumption of a suspended track.
903 //
904 void
906 {
907 
908  G4TransportationManager* transportMgr =
910 
911  // G4VProcess::StartTracking(aTrack);
912  fNewTrack= true;
913 
914  // The 'initialising' actions
915  // once taken in AlongStepGPIL -- if ( track.GetCurrentStepNumber()==1 )
916 
917  // fStartedNewTrack= true;
918 
919  fMassNavigator = transportMgr->GetNavigatorForTracking() ;
921 
922  G4ThreeVector position = aTrack->GetPosition();
923  G4ThreeVector direction = aTrack->GetMomentumDirection();
924 
925  fPathFinder->PrepareNewTrack( position, direction);
926  // This implies a call to fPathFinder->Locate( position, direction );
927 
928  // Global field, if any, must exist before tracking is started
930 
931  // reset safety value and center
932  //
933  fPreviousMassSafety = 0.0 ;
934  fPreviousFullSafety = 0.0 ;
935  fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
936 
937  // reset looping counter -- for motion in field
938  fNoLooperTrials= 0;
939 
940  // Must clear this state .. else it depends on last track's value
941  // --> a better solution would set this from state of suspended track TODO ?
942  // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
943 
944  // ChordFinder reset internal state
945  //
946  if( fGlobalFieldExists )
947  {
949  // Resets safety values, in case of overlaps.
950 
952  if( chordF ) { chordF->ResetStepEstimate(); }
953  }
954 
955  // Clear the chord finders of all fields (ie managers) derived objects
956  //
958  fieldMgrStore->ClearAllChordFindersState();
959 
960 #ifdef G4DEBUG_TRANSPORT
961  if( verboseLevel > 1 )
962  {
963  G4cout << " Returning touchable handle " << fCurrentTouchableHandle << G4endl;
964  }
965 #endif
966 
967  // Update the current touchable handle (from the track's)
968  //
970 }
971 
973 
974 void
976 {
978  fPathFinder->EndTrack(); // Resets TransportationManager to use ordinary Navigator
979 }
980 
982 
983 void
985 ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
986 {
987  static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
988  moduloFactor= 10, no_large_ediff= 0;
989 
990  if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
991  {
992  no_large_ediff ++;
993  if( (no_large_ediff% warnModulo) == 0 )
994  {
995  no_warnings++;
996  G4cout << "WARNING - G4CoupledTransportation::AlongStepGetPIL() "
997  << " Energy change in Step is above 1^-3 relative value. " << G4endl
998  << " Relative change in 'tracking' step = "
999  << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
1000  << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
1001  << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl;
1002  G4cout << " Energy has been corrected -- however, review"
1003  << " field propagation parameters for accuracy." << G4endl;
1004  if( (verboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
1005  {
1006  G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
1007  << " which determine fractional error per step for integrated quantities. " << G4endl
1008  << " Note also the influence of the permitted number of integration steps."
1009  << G4endl;
1010  }
1011  G4cerr << "ERROR - G4CoupledTransportation::AlongStepGetPIL()" << G4endl
1012  << " Bad 'endpoint'. Energy change detected"
1013  << " and corrected. "
1014  << " Has occurred already "
1015  << no_large_ediff << " times." << G4endl;
1016  if( no_large_ediff == warnModulo * moduloFactor )
1017  {
1018  warnModulo *= moduloFactor;
1019  }
1020  }
1021  }
1022 }
1023 
1025 
1027 {
1028  G4bool lastValue= fUseMagneticMoment;
1029  fUseMagneticMoment= useMoment;
1031  return lastValue;
1032 }
G4double GetLocalTime() const
G4double GetKineticEnergy() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const double prec
Definition: RanecuEngine.cc:58
virtual void ConfigureForTrack(const G4Track *)
CLHEP::Hep3Vector G4ThreeVector
G4SafetyHelper * GetSafetyHelper() const
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
G4double ObtainSafety(G4int navId, G4ThreeVector &globalCenterPoint)
void ProposeGlobalTime(G4double t)
void ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
G4LogicalVolume * GetLogicalVolume() const
G4bool DoesFieldChangeEnergy() const
static constexpr double MeV
Definition: G4SIunits.hh:214
G4VSensitiveDetector * GetSensitiveDetector() const
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
static G4FieldManagerStore * GetInstance()
void SetThresholds(G4double newEnWarn, G4double importantEnergy, G4int newMaxTrials)
G4StepPoint * GetPreStepPoint() const
static constexpr double mm
Definition: G4SIunits.hh:115
G4double GetThresholdImportantEnergy() const
G4double GetStepLength() const
#define G4endl
Definition: G4ios.hh:61
G4double GetCurrentSafety() const
const G4ThreeVector & GetMomentumDirection() const
G4ReferenceCountedHandle< G4VTouchable > G4TouchableHandle
void ReportMove(G4ThreeVector OldVector, G4ThreeVector NewVector, const G4String &Quantity)
G4EquationOfMotion * GetCurrentEquationOfMotion()
void SetSensitiveDetectorInTouchable(G4VSensitiveDetector *fSensitiveDetector)
void message(RunManager *runmanager)
Definition: ts_scorers.cc:72
G4TransportationLogger * fpLogger
G4FieldManager * GetCurrentFieldManager()
G4ParticleChangeForTransport fParticleChange
const G4TouchableHandle & GetTouchableHandle() const
const G4ThreeVector & GetPolarization() const
static constexpr double perMillion
Definition: G4SIunits.hh:334
void ReLocate(const G4ThreeVector &position)
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
Definition: G4VTouchable.cc:44
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &stepData)
G4Navigator * GetNavigatorForTracking() const
G4PropagatorInField * fFieldPropagator
G4int GetCurrentStepNumber() const
G4double GetThresholdWarningEnergy() const
#define G4ThreadLocal
Definition: tls.hh:69
G4Material * GetMaterial() const
G4double GetVelocity() const
G4double GetProperTime() const
virtual void Initialize(const G4Track &)
G4double ComputeStep(const G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4int navigatorId, G4int stepNo, G4double &pNewSafety, ELimited &limitedStep, G4FieldTrack &EndState, G4VPhysicalVolume *currentVolume)
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:440
G4ThreeVector GetSpin() const
G4ThreeVector GetPosition() const
void ProposeLocalTime(G4double t)
G4double GetPDGSpin() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
G4ParticleDefinition * GetDefinition() const
G4FieldManager * GetFieldManager() const
void ProposeLastStepInVolume(G4bool flag)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
#define fPreviousSftOrigin
void ResetStepEstimate()
G4double GetKineticEnergy() const
long G4long
Definition: G4Types.hh:80
void ProposePosition(G4double x, G4double y, G4double z)
G4double GetMagneticMoment() const
G4GPILSelection
static G4bool fUseMagneticMoment
G4TouchableHandle CreateTouchableHandle(G4int navId) const
G4double GetGlobalTime() const
static constexpr double kiloelectronvolt
Definition: G4SIunits.hh:207
G4double GetCharge() const
G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
#define fNewTrack
void StartTracking(G4Track *aTrack)
const G4ThreeVector & GetPosition() const
const G4Field * GetDetectorField() const
double mag2() const
Definition: G4Step.hh:76
void ProposeEnergy(G4double finalEnergy)
void ProposeTrueStepLength(G4double truePathLength)
G4double GetVelocity() const
unsigned int GetNumberGeometriesLimitingStep() const
G4GLOB_DLL std::ostream G4cerr
G4double GetPDGMagneticMoment() const
static G4TransportationManager * GetTransportationManager()
void Locate(const G4ThreeVector &position, const G4ThreeVector &direction, G4bool relativeSearch=true)
static G4ProductionCutsTable * GetProductionCutsTable()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
double mag() const
ELimited
int G4int
Definition: G4Types.hh:78
void ProposeProperTime(G4double finalProperTime)
G4double PostStepGetPhysicalInteractionLength(const G4Track &, G4double previousStepSize, G4ForceCondition *pForceCond)
G4int verboseLevel
Definition: G4VProcess.hh:371
static G4bool EnableUseMagneticMoment(G4bool useMoment=true)
void ProposeFirstStepInVolume(G4bool flag)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void PrepareNewTrack(const G4ThreeVector &position, const G4ThreeVector &direction, G4VPhysicalVolume *massStartVol=0)
G4int ActivateNavigator(G4Navigator *aNavigator)
const G4ThreeVector & GetMomentumDir() const
void SetMaterialInTouchable(G4Material *fMaterial)
G4bool IsParticleLooping() const
G4TouchableHandle fCurrentTouchableHandle
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4ForceCondition
G4int GetThresholdTrials() const
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4ProductionCuts * GetProductionCuts() const
G4PropagatorInField * GetPropagatorInField() const
G4ChordFinder * GetChordFinder()
void ReportLoopingTrack(const G4Track &track, const G4Step &stepInfo, G4int numTrials, G4long noCalls, const char *methodName) const
const G4Material * GetMaterial() const
G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
G4double GetMass() const
T sqr(const T &x)
Definition: templates.hh:145
G4double GetTotalMomentum() const
void SetTouchableHandle(const G4TouchableHandle &fTouchable)
void SetMaterialCutsCoupleInTouchable(const G4MaterialCutsCouple *fMaterialCutsCouple)
G4TrackStatus GetTrackStatus() const
G4CoupledTransportation(G4int verbosityLevel=0)
const G4ThreeVector & GetMomentumDirection() const
static G4PathFinder * GetInstance()
Definition: G4PathFinder.cc:57
#define DBL_MAX
Definition: templates.hh:83
void ProposeTrackStatus(G4TrackStatus status)
G4double GetLabTimeOfFlight() const
const G4DynamicParticle * GetDynamicParticle() const
static constexpr double perThousand
Definition: G4SIunits.hh:333
void SetMomentumChanged(G4bool b)
G4double ComputeSafety(const G4ThreeVector &globalPoint)
const G4String & GetName() const
G4bool IsGravityActive() const
Definition: G4Field.hh:98
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VPhysicalVolume * GetVolume() const
G4double GetTotalEnergy() const