| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
28 // ------------------------------------------------------------
29 // GEANT 4 class implementation
30 //
31 // =======================================================================
32 // Created: 19 March 1997, J. Apostolakis
33 // =======================================================================
34
35 #include "G4CoupledTransportation.hh"
36 #include "G4TransportationProcessType.hh"
37 #include "G4TransportationLogger.hh"
38
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ProductionCutsTable.hh"
42 #include "G4ParticleTable.hh"
43 #include "G4ChordFinder.hh"
44 #include "G4Field.hh"
45 #include "G4FieldTrack.hh"
46 #include "G4FieldManagerStore.hh"
47 #include "G4PathFinder.hh"
48
49 #include "G4PropagatorInField.hh"
50 #include "G4TransportationManager.hh"
51
52 class G4VSensitiveDetector;
53
54 G4bool G4CoupledTransportation::fSignifyStepInAnyVolume= false;
55 // This mode must apply to all threads
56
57 //////////////////////////////////////////////////////////////////////////
58 //
59 // Constructor
60
61 G4CoupledTransportation::G4CoupledTransportation( G4int verbosity )
62 : G4Transportation( verbosity, "CoupledTransportation" ),
63 fPreviousMassSafety( 0.0 ),
64 fPreviousFullSafety( 0.0 ),
65 fMassGeometryLimitedStep( false ),
66 fFirstStepInMassVolume( true )
67 {
68 SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION));
69 // SetVerboseLevel is called in the constructor of G4Transportation
70
71 if( verboseLevel > 0 )
72 {
73 G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
74 G4cout << " Verbose level is " << verboseLevel << G4endl;
75 G4cout << " Reports First/Last in "
76 << (fSignifyStepInAnyVolume ? " any " : " mass " )
77 << " geometry " << G4endl;
78 }
79 fPathFinder= G4PathFinder::GetInstance();
80 }
81
82 //////////////////////////////////////////////////////////////////////////
83
84 G4CoupledTransportation::~G4CoupledTransportation()
85 {
86 }
87
88 //////////////////////////////////////////////////////////////////////////
89 //
90 // Responsibilities:
91 // Find whether the geometry limits the Step, and to what length
92 // Calculate the new value of the safety and return it.
93 // Store the final time, position and momentum.
94
95 G4double G4CoupledTransportation::
96 AlongStepGetPhysicalInteractionLength( const G4Track& track,
97 G4double, // previousStepSize
98 G4double currentMinimumStep,
99 G4double& proposedSafetyForStart,
100 G4GPILSelection* selection )
101 {
102 G4double geometryStepLength;
103 G4double startFullSafety= 0.0; // estimated safety for start point (all geometries)
104 G4double safetyProposal= -1.0; // local copy of proposal
105
106 G4ThreeVector EndUnitMomentum ;
107 G4double lengthAlongCurve = 0.0 ;
108
109 fParticleIsLooping = false ;
110
111 // Initial actions moved to StartTrack()
112 // --------------------------------------
113 // Note: in case another process changes touchable handle
114 // it will be necessary to add here (for all steps)
115 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
116
117 // GPILSelection is set to defaule value of CandidateForSelection
118 // It is a return value
119 //
120 *selection = CandidateForSelection ;
121
122 fFirstStepInMassVolume = fNewTrack || fMassGeometryLimitedStep ;
123 fFirstStepInVolume = fNewTrack || fGeometryLimitedStep ;
124
125 #ifdef G4DEBUG_TRANSPORT
126 G4cout << " CoupledTransport::AlongStep GPIL: "
127 << " 1st-step: any= " <<fFirstStepInVolume << " ( geom= "
128 << fGeometryLimitedStep << " ) "
129 << " mass= " << fFirstStepInMassVolume << " ( geom= "
130 << fMassGeometryLimitedStep << " ) "
131 << " newTrack= " << fNewTrack << G4endl;
132 #endif
133
134 // fLastStepInVolume= false;
135 fNewTrack = false;
136
137 // Get initial Energy/Momentum of the track
138 //
139 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
140 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
141 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
142 G4ThreeVector startPosition = track.GetPosition() ;
143 G4VPhysicalVolume* currentVolume= track.GetVolume();
144
145 #ifdef G4DEBUG_TRANSPORT
146 if( verboseLevel > 1 )
147 {
148 G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume "
149 << currentVolume->GetName() << G4endl;
150 }
151 #endif
152 // G4double theTime = track.GetGlobalTime() ;
153
154 // The Step Point safety can be limited by other geometries and/or the
155 // assumptions of any process - it's not always the geometrical safety.
156 // We calculate the starting point's isotropic safety here.
157 //
158 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
159 G4double MagSqShift = OriginShift.mag2() ;
160 startFullSafety= 0.0;
161
162 // Recall that FullSafety <= MassSafety
163 // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
164 if( MagSqShift < sqr(fPreviousFullSafety) )
165 {
166 G4double mag_shift= std::sqrt(MagSqShift);
167 startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
168 // Need to be consistent between full safety with Mass safety
169 // in order reproduce results in simple case
170 // --> use same calculation method
171
172 // Only compute full safety if massSafety > 0. Else it remains 0
173 // startFullSafety = fPathFinder->ComputeSafety( startPosition );
174 }
175
176 // Is the particle charged or has it a magnetic moment?
177 //
178 G4double particleCharge = pParticle->GetCharge() ;
179 G4double magneticMoment = pParticle->GetMagneticMoment() ;
180 G4double restMass = pParticle->GetMass() ;
181
182 fMassGeometryLimitedStep = false ; // Set default - alex
183 fGeometryLimitedStep = false;
184
185 // There is no need to locate the current volume. It is Done elsewhere:
186 // On track construction
187 // By the tracking, after all AlongStepDoIts, in "Relocation"
188
189 // Check if the particle has a force, EM or gravitational, exerted on it
190 //
191 G4FieldManager* fieldMgr= nullptr;
192 G4bool fieldExertsForce = false ;
193
194 const G4Field* ptrField= nullptr;
195
196 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
197 G4bool eligibleEM = (particleCharge != 0.0)
198 || ( fUseMagneticMoment && (magneticMoment != 0.0) );
199 G4bool eligibleGrav = fUseGravity && (restMass != 0.0) ;
200
201 if( (fieldMgr!=nullptr) && (eligibleEM||eligibleGrav) )
202 {
203 // Message the field Manager, to configure it for this track
204 //
205 fieldMgr->ConfigureForTrack( &track );
206
207 // The above call can transition from a null field-ptr oto a finite field.
208 // If the field manager has no field ptr, the field is zero
209 // by definition ( = there is no field ! )
210 //
211 ptrField= fieldMgr->GetDetectorField();
212
213 if( ptrField != nullptr)
214 {
215 fieldExertsForce = eligibleEM
216 || ( eligibleGrav && ptrField->IsGravityActive() );
217 }
218 }
219 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
220
221 if( fieldExertsForce )
222 {
223 auto equationOfMotion= fFieldPropagator->GetCurrentEquationOfMotion();
224
225 G4ChargeState chargeState(particleCharge, // Charge can change (dynamic)
226 magneticMoment,
227 pParticleDef->GetPDGSpin() );
228 if( equationOfMotion )
229 {
230 equationOfMotion->SetChargeMomentumMass( chargeState,
231 momentumMagnitude,
232 restMass );
233 }
234 #ifdef G4DEBUG_TRANSPORT
235 else
236 {
237 G4cerr << " ERROR in G4CoupledTransportation> "
238 << "Cannot find valid Equation of motion: " << G4endl;
239 << " Unable to pass Charge, Momentum and Mass " << G4endl;
240 }
241 #endif
242 }
243
244 G4ThreeVector polarizationVec = track.GetPolarization() ;
245 G4FieldTrack aFieldTrack = G4FieldTrack(startPosition,
246 track.GetGlobalTime(), // Lab.
247 track.GetMomentumDirection(),
248 track.GetKineticEnergy(),
249 restMass,
250 particleCharge,
251 polarizationVec,
252 pParticleDef->GetPDGMagneticMoment(),
253 0.0, // Length along track
254 pParticleDef->GetPDGSpin() ) ;
255 G4int stepNo= track.GetCurrentStepNumber();
256
257 ELimited limitedStep;
258 G4FieldTrack endTrackState('a'); // Default values
259
260 fMassGeometryLimitedStep = false ; // default
261 fGeometryLimitedStep = false;
262 if( currentMinimumStep > 0 )
263 {
264 G4double newMassSafety= 0.0; // temp. for recalculation
265
266 // Do the Transport in the field (non recti-linear)
267 //
268 lengthAlongCurve = fPathFinder->ComputeStep( aFieldTrack,
269 currentMinimumStep,
270 G4TransportationManager::kMassNavigatorId,
271 stepNo,
272 newMassSafety,
273 limitedStep,
274 endTrackState,
275 currentVolume ) ;
276
277 G4double newFullSafety= fPathFinder->GetCurrentSafety();
278 // this was estimated already in step above
279
280 if( limitedStep == kUnique || limitedStep == kSharedTransport )
281 {
282 fMassGeometryLimitedStep = true ;
283 }
284
285 fGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0);
286
287 #ifdef G4DEBUG_TRANSPORT
288 if( fMassGeometryLimitedStep && !fGeometryLimitedStep )
289 {
290 std::ostringstream message;
291 message << " ERROR in determining geometries limiting the step" << G4endl;
292 message << " Limiting: mass=" << fMassGeometryLimitedStep
293 << " any= " << fGeometryLimitedStep << G4endl;
294 message << "Incompatible conditions - by which geometries was it limited ?";
295 G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()",
296 "PathFinderConfused", FatalException, message);
297 }
298 #endif
299
300 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep);
301
302 // Momentum: Magnitude and direction can be changed too now ...
303 //
304 fMomentumChanged = true ;
305 fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
306
307 // Remember last safety origin & value.
308 fPreviousSftOrigin = startPosition ;
309 fPreviousMassSafety = newMassSafety ;
310 fPreviousFullSafety = newFullSafety ;
311 // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
312
313 #ifdef G4DEBUG_TRANSPORT
314 if( verboseLevel > 1 )
315 {
316 G4cout << "G4Transport:CompStep> "
317 << " called the pathfinder for a new step at " << startPosition
318 << " and obtained step = " << lengthAlongCurve << G4endl;
319 G4cout << " New safety (preStep) = " << newMassSafety << G4endl;
320 }
321 #endif
322
323 // Store as best estimate value
324 startFullSafety = newFullSafety ;
325
326 // Get the End-Position and End-Momentum (Dir-ection)
327 fTransportEndPosition = endTrackState.GetPosition() ;
328 fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ;
329 }
330 else
331 {
332 geometryStepLength = lengthAlongCurve= 0.0 ;
333 fMomentumChanged = false ;
334 // fMassGeometryLimitedStep = false ; // --- ???
335 // fGeometryLimitedStep = true;
336 fTransportEndMomentumDir = track.GetMomentumDirection();
337 fTransportEndKineticEnergy = track.GetKineticEnergy();
338
339 fTransportEndPosition = startPosition;
340
341 endTrackState= aFieldTrack; // Ensures that time is updated
342 }
343 // G4FieldTrack aTrackState(endTrackState);
344
345 if( !fieldExertsForce )
346 {
347 fParticleIsLooping = false ;
348 fMomentumChanged = false ;
349 fEndGlobalTimeComputed = false ;
350 }
351 else
352 {
353 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
354
355 #ifdef G4DEBUG_TRANSPORT
356 if( verboseLevel > 1 )
357 {
358 G4cout << " G4CT::CS End Position = "
359 << fTransportEndPosition << G4endl;
360 G4cout << " G4CT::CS End Direction = "
361 << fTransportEndMomentumDir << G4endl;
362 }
363 #endif
364 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
365 {
366 // If the field can change energy, then the time must be integrated
367 // - so this should have been updated
368 //
369 fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight();
370 fEndGlobalTimeComputed = true;
371
372 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
373 // a cleaner way is to have FieldTrack knowing whether time
374 // is updated
375 }
376 else
377 {
378 // The energy should be unchanged by field transport,
379 // - so the time changed will be calculated elsewhere
380 //
381 fEndGlobalTimeComputed = false;
382
383 #ifdef G4VERBOSE
384 // Check that the integration preserved the energy
385 // - and if not correct this!
386 G4double startEnergy= track.GetKineticEnergy();
387 G4double endEnergy= fTransportEndKineticEnergy;
388
389 G4double absEdiff = std::fabs(startEnergy- endEnergy);
390 if( (verboseLevel > 1) && ( absEdiff > perThousand * endEnergy) )
391 {
392 ReportInexactEnergy(startEnergy, endEnergy);
393 } // end of if (verboseLevel)
394 #endif
395 // Correct the energy for fields that conserve it
396 // This - hides the integration error
397 // - but gives a better physical answer
398 fTransportEndKineticEnergy= track.GetKineticEnergy();
399 }
400 }
401
402 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
403 fTransportEndSpin = endTrackState.GetSpin();
404
405 // Calculate the safety
406
407 safetyProposal= startFullSafety; // used to be startMassSafety
408 // Changed to accomodate processes that cannot update the safety
409
410 // Update safety for the end-point, if becomes negative at the end-point.
411
412 if( (startFullSafety < fEndPointDistance )
413 && ( particleCharge != 0.0 ) ) // Only needed to prepare for MSC
414 // && !fGeometryLimitedStep ) // To-Try: No safety update if at boundary
415 {
416 G4double endFullSafety =
417 fPathFinder->ComputeSafety( fTransportEndPosition);
418 // Expected mission -- only mass geometry's safety
419 // fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
420 // Yet discrete processes only have poststep -- and this cannot
421 // currently revise the safety
422 // ==> so we use the all-geometry safety as a precaution
423
424 fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition);
425 // Pushing safety to Helper avoids recalculation at this point
426
427 G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value
428 G4double endMassSafety= fPathFinder->ObtainSafety( G4TransportationManager::kMassNavigatorId, centerPt);
429 // Retrieves the mass value from PathFinder (it calculated it)
430
431 fPreviousMassSafety = endMassSafety ;
432 fPreviousFullSafety = endFullSafety;
433 fPreviousSftOrigin = fTransportEndPosition ;
434
435 // The convention (Stepping Manager's) is safety from the start point
436 //
437 safetyProposal = endFullSafety + fEndPointDistance;
438 // --> was endMassSafety
439 // Changed to accomodate processes that cannot update the safety
440
441 #ifdef G4DEBUG_TRANSPORT
442 G4int prec= G4cout.precision(12) ;
443 G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
444 G4cout << " Revised Safety at endpoint " << fTransportEndPosition
445 << " give safety values: Mass= " << endMassSafety
446 << " All= " << endFullSafety << G4endl ;
447 G4cout << " Adding endpoint distance " << fEndPointDistance
448 << " to obtain pseudo-safety= " << safetyProposal << G4endl ;
449 G4cout.precision(prec);
450 }
451 else
452 {
453 G4int prec= G4cout.precision(12) ;
454 G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
455 G4cout << " Quick Safety estimate at endpoint "
456 << fTransportEndPosition
457 << " gives safety endpoint value = "
458 << startFullSafety - fEndPointDistance
459 << " using start-point value " << startFullSafety
460 << " and endpointDistance " << fEndPointDistance << G4endl;
461 G4cout.precision(prec);
462 #endif
463 }
464
465 proposedSafetyForStart= safetyProposal;
466 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
467
468 return geometryStepLength ;
469 }
470
471 /////////////////////////////////////////////////////////////////////////////
472
473 void G4CoupledTransportation::
474 ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector,
475 const G4String& Quantity )
476 {
477 G4ThreeVector moveVec = ( NewVector - OldVector );
478
479 G4cerr << G4endl
480 << "**************************************************************"
481 << G4endl;
482 G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
483 << " and value from Track in PostStepDoIt. " << G4endl
484 << "Change of " << Quantity << " is " << moveVec.mag() / mm
485 << " mm long, "
486 << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
487 << "Endpoint of ComputeStep was " << OldVector
488 << " and current position to locate is " << NewVector << G4endl;
489 }
490
491 /////////////////////////////////////////////////////////////////////////////
492
493 G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track,
494 const G4Step& )
495 {
496 G4TouchableHandle retCurrentTouchable ; // The one to return
497
498 // Initialize ParticleChange (by setting all its members equal
499 // to corresponding members in G4Track)
500 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
501
502 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
503
504 if( fSignifyStepInAnyVolume )
505 {
506 fParticleChange.ProposeFirstStepInVolume( fFirstStepInVolume );
507 }
508 else
509 {
510 fParticleChange.ProposeFirstStepInVolume( fFirstStepInMassVolume );
511 }
512
513 // Check that the end position and direction are preserved
514 // since call to AlongStepDoIt
515
516 #ifdef G4DEBUG_TRANSPORT
517 if( ( verboseLevel > 0 )
518 && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) )
519 {
520 ReportMove( track.GetPosition(), fTransportEndPosition,
521 "End of Step Position" );
522 G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl;
523 }
524
525 // If the Step was determined by the volume boundary, relocate the particle
526 // The pathFinder will know that the geometry limited the step (!?)
527
528 if( verboseLevel > 0 )
529 {
530 G4cout << " Calling PathFinder::Locate() from "
531 << " G4CoupledTransportation::PostStepDoIt " << G4endl;
532 G4cout << " fGeometryLimitedStep is " << fGeometryLimitedStep << G4endl;
533 }
534 #endif
535
536 if(fGeometryLimitedStep)
537 {
538 fPathFinder->Locate( track.GetPosition(),
539 track.GetMomentumDirection(),
540 true);
541
542 // fCurrentTouchable will now become the previous touchable,
543 // and what was the previous will be freed.
544 // (Needed because the preStepPoint can point to the previous touchable)
545
546 fCurrentTouchableHandle=
547 fPathFinder->CreateTouchableHandle( G4TransportationManager::kMassNavigatorId );
548
549 #ifdef G4DEBUG_TRANSPORT
550 if( verboseLevel > 1 )
551 {
552 G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume();
553 G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = "
554 << vol;
555 if( vol ) { G4cout << "Name=" << vol->GetName(); }
556 G4cout << G4endl;
557 }
558 #endif
559
560 // Check whether the particle is out of the world volume
561 // If so it has exited and must be killed.
562 //
563 if( fCurrentTouchableHandle->GetVolume() == 0 )
564 {
565 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
566 }
567 retCurrentTouchable = fCurrentTouchableHandle ;
568 // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
569 }
570 else // fGeometryLimitedStep is false
571 {
572 #ifdef G4DEBUG_TRANSPORT
573 if( verboseLevel > 1 )
574 {
575 G4cout << "G4CoupledTransportation::PostStepDoIt -- "
576 << " fGeometryLimitedStep = " << fGeometryLimitedStep
577 << " must be false " << G4endl;
578 }
579 #endif
580 // This serves only to move each of the Navigator's location
581 //
582 // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
583
584 fPathFinder->ReLocate( track.GetPosition() );
585 // track.GetMomentumDirection() );
586
587 // Keep the value of the track's current Touchable is retained,
588 // and use it to overwrite the (unset) one in particle change.
589 // Expect this must be fCurrentTouchable too
590 // - could it be different, eg at the start of a step ?
591 //
592 retCurrentTouchable = track.GetTouchableHandle() ;
593 // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
594 } // endif ( fGeometryLimitedStep )
595
596 #ifdef G4DEBUG_NAVIGATION
597 G4cout << " CoupledTransport::AlongStep GPIL: "
598 << " last-step: any= " << fGeometryLimitedStep << " . ..... x . "
599 << " mass= " << fMassGeometryLimitedStep << G4endl;
600 #endif
601
602 if( fSignifyStepInAnyVolume )
603 fParticleChange.ProposeLastStepInVolume(fGeometryLimitedStep);
604 else
605 fParticleChange.ProposeLastStepInVolume(fMassGeometryLimitedStep);
606
607 SetTouchableInformation(retCurrentTouchable);
608
609 return &fParticleChange ;
610 }
611
612 /////////////////////////////////////////////////////////////////////////////
613 // New method takes over the responsibility to reset the state of
614 // G4CoupledTransportation object:
615 // - at the start of a new track, and
616 // - on the resumption of a suspended track.
617 //
618 void
619 G4CoupledTransportation::StartTracking(G4Track* aTrack)
620 {
621 G4Transportation::StartTracking(aTrack);
622
623 G4ThreeVector position = aTrack->GetPosition();
624 G4ThreeVector direction = aTrack->GetMomentumDirection();
625
626 fPathFinder->PrepareNewTrack( position, direction);
627 // This implies a call to fPathFinder->Locate( position, direction );
628
629 // reset safety value and center
630 //
631 fPreviousMassSafety = 0.0 ;
632 fPreviousFullSafety = 0.0 ;
633 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
634 }
635
636 /////////////////////////////////////////////////////////////////////////////
637
638 void
639 G4CoupledTransportation::EndTracking()
640 {
641 G4TransportationManager::GetTransportationManager()->InactivateAll();
642 fPathFinder->EndTrack();
643 // Resets TransportationManager to use ordinary Navigator
644 }
645
646 /////////////////////////////////////////////////////////////////////////////
647
648 void
649 G4CoupledTransportation::
650 ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
651 {
652 static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
653 moduloFactor= 10, no_large_ediff= 0;
654
655 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
656 {
657 no_large_ediff ++;
658 if( (no_large_ediff% warnModulo) == 0 )
659 {
660 no_warnings++;
661 std::ostringstream message;
662 message << "Energy change in Step is above 1^-3 relative value. "
663 << G4endl
664 << " Relative change in 'tracking' step = "
665 << std::setw(15) << (endEnergy-startEnergy)/startEnergy
666 << G4endl
667 << " Starting E= " << std::setw(12) << startEnergy / MeV
668 << " MeV " << G4endl
669 << " Ending E= " << std::setw(12) << endEnergy / MeV
670 << " MeV " << G4endl
671 << "Energy has been corrected -- however, review"
672 << " field propagation parameters for accuracy." << G4endl;
673 if ( (verboseLevel > 2 ) || (no_warnings<4)
674 || (no_large_ediff == warnModulo * moduloFactor) )
675 {
676 message << "These include EpsilonStepMax(/Min) in G4FieldManager,"
677 << G4endl
678 << "which determine fractional error per step for integrated quantities."
679 << G4endl
680 << "Note also the influence of the permitted number of integration steps."
681 << G4endl;
682 }
683 message << "Bad 'endpoint'. Energy change detected and corrected."
684 << G4endl
685 << "Has occurred already " << no_large_ediff << " times.";
686 G4Exception("G4CoupledTransportation::AlongStepGetPIL()",
687 "EnergyChange", JustWarning, message);
688 if( no_large_ediff == warnModulo * moduloFactor )
689 {
690 warnModulo *= moduloFactor;
691 }
692 }
693 }
694 }
695