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