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