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