Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PathFinder.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: G4PathFinder.cc 90009 2015-05-08 07:42:39Z gcosmo $
28 // GEANT4 tag $ Name: $
29 //
30 // class G4PathFinder Implementation
31 //
32 // Original author: John Apostolakis, April 2006
33 //
34 // --------------------------------------------------------------------
35 
36 #include <iomanip>
37 
38 #include "G4PathFinder.hh"
39 
40 #include "G4SystemOfUnits.hh"
41 #include "G4GeometryTolerance.hh"
42 #include "G4Navigator.hh"
43 #include "G4PropagatorInField.hh"
45 #include "G4MultiNavigator.hh"
46 #include "G4SafetyHelper.hh"
47 
48 // Initialise the static instance of the singleton
49 //
50 G4ThreadLocal G4PathFinder* G4PathFinder::fpPathFinder=0;
51 
52 // ----------------------------------------------------------------------------
53 // GetInstance()
54 //
55 // Retrieve the static instance of the singleton
56 //
58 {
59  if( ! fpPathFinder )
60  {
61  fpPathFinder = new G4PathFinder;
62  }
63  return fpPathFinder;
64 }
65 
66 // ----------------------------------------------------------------------------
67 // Constructor
68 //
70  : fEndState( G4ThreeVector(), G4ThreeVector(), 0., 0., 0., 0., 0.),
71  fFieldExertedForce(false),
72  fRelocatedPoint(true),
73  fLastStepNo(-1), fCurrentStepNo(-1),
74  fVerboseLevel(0)
75 {
76  fpMultiNavigator= new G4MultiNavigator();
77 
79  fpFieldPropagator = fpTransportManager->GetPropagatorInField();
80 
82 
83  fNoActiveNavigators= 0;
85  fLastLocatedPosition= Big3Vector;
86  fSafetyLocation= Big3Vector;
87  fPreSafetyLocation= Big3Vector;
88  fPreStepLocation= Big3Vector;
89 
90  fPreSafetyMinValue= -1.0;
91  fMinSafety_PreStepPt= -1.0;
92  fMinSafety_atSafLocation= -1.0;
93  fMinStep= -1.0;
94  fTrueMinStep= -1.0;
95  fPreStepCenterRenewed= false;
96  fNewTrack= false;
97  fNoGeometriesLimiting= 0;
98 
99  for( G4int num=0; num< fMaxNav; ++num )
100  {
101  fpNavigator[num] = 0;
102  fLimitTruth[num] = false;
103  fLimitedStep[num] = kUndefLimited;
104  fCurrentStepSize[num] = -1.0;
105  fLocatedVolume[num] = 0;
106  fPreSafetyValues[num]= -1.0;
107  fCurrentPreStepSafety[num] = -1.0;
108  fNewSafetyComputed[num]= -1.0;
109  }
110 }
111 
112 // ----------------------------------------------------------------------------
113 // Destructor
114 //
116 {
117  delete fpMultiNavigator;
118  if (fpPathFinder) { delete fpPathFinder; fpPathFinder=0; }
119 }
120 
121 // ----------------------------------------------------------------------------
122 //
123 void
125 {
126  G4Navigator *navigatorForPropagation=0, *massNavigator=0;
127 
128  massNavigator= fpTransportManager->GetNavigatorForTracking();
129  if( enableChoice )
130  {
131  navigatorForPropagation= fpMultiNavigator;
132 
133  // Enable SafetyHelper to use PF
134  //
135  fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(true);
136  }
137  else
138  {
139  navigatorForPropagation= massNavigator;
140 
141  // Disable SafetyHelper to use PF
142  //
143  fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(false);
144  }
145  fpFieldPropagator->SetNavigatorForPropagating(navigatorForPropagation);
146 }
147 
148 // ----------------------------------------------------------------------------
149 //
150 G4double
151 G4PathFinder::ComputeStep( const G4FieldTrack &InitialFieldTrack,
152  G4double proposedStepLength,
153  G4int navigatorNo,
154  G4int stepNo, // find next step
155  G4double &pNewSafety, // for this geom
156  ELimited &limitedStep,
157  G4FieldTrack &EndState,
158  G4VPhysicalVolume* currentVolume)
159 {
160  G4double possibleStep= -1.0;
161 
162 #ifdef G4DEBUG_PATHFINDER
163  if( fVerboseLevel > 2 )
164  {
165  G4cout << " -------------------------" << G4endl;
166  G4cout << " G4PathFinder::ComputeStep - entered " << G4endl;
167  G4cout << " - stepNo = " << std::setw(4) << stepNo << " "
168  << " navigatorId = " << std::setw(2) << navigatorNo << " "
169  << " proposed step len = " << proposedStepLength << " " << G4endl;
170  G4cout << " PF::ComputeStep requested step "
171  << " from " << InitialFieldTrack.GetPosition()
172  << " dir " << InitialFieldTrack.GetMomentumDirection() << G4endl;
173  }
174 #endif
175 #ifdef G4VERBOSE
176  if( navigatorNo >= fNoActiveNavigators )
177  {
178  std::ostringstream message;
179  message << "Bad Navigator ID !" << G4endl
180  << " Requested Navigator ID = " << navigatorNo << G4endl
181  << " Number of active navigators = " << fNoActiveNavigators;
182  G4Exception("G4PathFinder::ComputeStep()", "GeomNav0002",
183  FatalException, message);
184  }
185 #endif
186 
187  if( fNewTrack || (stepNo != fLastStepNo) )
188  {
189  // This is a new track or a new step, so we must make the step
190  // ( else we can simply retrieve its results for this Navigator Id )
191 
192  G4FieldTrack currentState= InitialFieldTrack;
193 
194  fCurrentStepNo = stepNo;
195 
196  // Check whether a process shifted the position
197  // since the last step -- by physics processes
198  //
199  G4ThreeVector newPosition = InitialFieldTrack.GetPosition();
200  G4ThreeVector moveVector= newPosition - fLastLocatedPosition;
201  G4double moveLenSq= moveVector.mag2();
202  if( moveLenSq > kCarTolerance * kCarTolerance )
203  {
204  G4ThreeVector newDirection = InitialFieldTrack.GetMomentumDirection();
205 #ifdef G4DEBUG_PATHFINDER
206  if( fVerboseLevel > 2 )
207  {
208  G4double moveLen= std::sqrt( moveLenSq );
209  G4cout << " G4PathFinder::ComputeStep : Point moved since last step "
210  << " -- at step # = " << stepNo << G4endl
211  << " by " << moveLen << " to " << newPosition << G4endl;
212  }
213 #endif
214  MovePoint(); // Unintentional changed -- ????
215 
216  // Relocate to cope with this move -- else could abort !?
217  //
218  Locate( newPosition, newDirection );
219  }
220 
221  // Check whether the particle have an (EM) field force exerting upon it
222  //
223  G4double particleCharge= currentState.GetCharge();
224 
225  G4FieldManager* fieldMgr=0;
226  G4bool fieldExertsForce = false ;
227  if( (particleCharge != 0.0) )
228  {
229  fieldMgr= fpFieldPropagator->FindAndSetFieldManager( currentVolume );
230 
231  // Protect for case where field manager has no field (= field is zero)
232  //
233  fieldExertsForce = (fieldMgr != 0)
234  && (fieldMgr->GetDetectorField() != 0);
235  }
236  fFieldExertedForce = fieldExertsForce; // Store for use in later calls
237  // referring to this 'step'.
238 
239  fNoGeometriesLimiting= -1; // At start of track, no process limited step
240  if( fieldExertsForce )
241  {
242  DoNextCurvedStep( currentState, proposedStepLength, currentVolume );
243  //--------------
244  }else{
245  DoNextLinearStep( currentState, proposedStepLength );
246  //--------------
247  }
248  fLastStepNo= stepNo;
249 
250 #ifdef G4DEBUG_PATHFINDER
251  if ( (fNoGeometriesLimiting < 0)
252  || (fNoGeometriesLimiting > fNoActiveNavigators) )
253  {
254  std::ostringstream message;
255  message << "Number of geometries limiting the step not set." << G4endl
256  << " Number of geometries limiting step = "
257  << fNoGeometriesLimiting;
258  G4Exception("G4PathFinder::ComputeStep()",
259  "GeomNav0002", FatalException, message);
260  }
261 #endif
262  }
263 #ifdef G4DEBUG_PATHFINDER
264  else
265  {
266  const G4double checkTolerance = 1.0e-9;
267  if( proposedStepLength < fTrueMinStep * ( 1.0 + checkTolerance) ) // For 2nd+ geometry
268  {
269  std::ostringstream message;
270  message.precision( 12 );
271  message << "Problem in step size request." << G4endl
272  << " Being requested to make a step which is shorter"
273  << " than the minimum Step " << G4endl
274  << " already computed for any Navigator/geometry during"
275  << " this tracking-step: " << G4endl
276  << " This could happen due to an error in process ordering."
277  << G4endl
278  << " Check that all physics processes are registered"
279  << " before all processes with a navigator/geometry."
280  << G4endl
281  << " If using pre-packaged physics list and/or"
282  << " functionality, please report this error."
283  << G4endl << G4endl
284  << " Additional information for problem: " << G4endl
285  << " Steps request/proposed = " << proposedStepLength
286  << G4endl
287  << " MinimumStep (true) = " << fTrueMinStep
288  << G4endl
289  << " MinimumStep (navraw) = " << fMinStep
290  << G4endl
291  << " Navigator raw return value" << G4endl
292  << " Requested step now = " << proposedStepLength
293  << G4endl
294  << " Difference min-req (absolute) = "
295  << fTrueMinStep-proposedStepLength << G4endl
296  << " Relative (to max of two) = "
297  << (fTrueMinStep-proposedStepLength)
298  / std::max(proposedStepLength, fTrueMinStep) << G4endl
299  << " -- Step info> stepNo= " << stepNo
300  << " last= " << fLastStepNo
301  << " newTr= " << fNewTrack << G4endl;
302  G4Exception("G4PathFinder::ComputeStep()",
303  "GeomNav0003", FatalException, message);
304  }
305  else
306  {
307  // This is neither a new track nor a new step -- just another
308  // client accessing information for the current track, step
309  // We will simply retrieve the results of the synchronous
310  // stepping for this Navigator Id below.
311  //
312  if( fVerboseLevel > 1 )
313  {
314  G4cout << " G4P::CS -> Not calling DoNextLinearStep: "
315  << " stepNo= " << stepNo << " last= " << fLastStepNo
316  << " new= " << fNewTrack << " Step already done" << G4endl;
317  }
318  }
319  }
320 #endif
321 
322  fNewTrack= false;
323 
324  // Prepare the information to return
325 
326  pNewSafety = fCurrentPreStepSafety[ navigatorNo ];
327  limitedStep = fLimitedStep[ navigatorNo ];
328  fRelocatedPoint= false;
329 
330  possibleStep= std::min(proposedStepLength, fCurrentStepSize[ navigatorNo ]);
331  EndState = fEndState; // now corrected for smaller step, if needed
332 
333 #ifdef G4DEBUG_PATHFINDER
334  if( fVerboseLevel > 0 )
335  {
336  G4cout << " G4PathFinder::ComputeStep returns "
337  << fCurrentStepSize[ navigatorNo ]
338  << " for Navigator " << navigatorNo
339  << " Limited step = " << limitedStep
340  << " Safety(mm) = " << pNewSafety / mm
341  << G4endl;
342  }
343 #endif
344 
345  return possibleStep;
346 }
347 
348 // ----------------------------------------------------------------------
349 
350 void
352  const G4ThreeVector& direction,
353  G4VPhysicalVolume* massStartVol)
354 {
355  // Key purposes:
356  // - Check and cache set of active navigators
357  // - Reset state for new track
358 
359  G4int num=0;
360 
362  // Switch PropagatorInField to use MultiNavigator
363 
364  fpTransportManager->GetSafetyHelper()->InitialiseHelper();
365  // Reinitialise state of safety helper -- avoid problems with overlaps
366 
367  fNewTrack= true;
368  this->MovePoint(); // Signal further that the last status is wiped
369 
370  fpFieldPropagator->PrepareNewTrack(); // Inform field propagator of new track
371 
372  // Message the G4NavigatorPanel / Dispatcher to find active navigators
373  //
374  std::vector<G4Navigator*>::iterator pNavigatorIter;
375 
376  fNoActiveNavigators= fpTransportManager-> GetNoActiveNavigators();
377  if( fNoActiveNavigators > fMaxNav )
378  {
379  std::ostringstream message;
380  message << "Too many active Navigators / worlds." << G4endl
381  << " Transportation Manager has "
382  << fNoActiveNavigators << " active navigators." << G4endl
383  << " This is more than the number allowed = "
384  << fMaxNav << " !";
385  G4Exception("G4PathFinder::PrepareNewTrack()", "GeomNav0002",
386  FatalException, message);
387  }
388 
389  fpMultiNavigator->PrepareNavigators();
390  //------------------------------------
391 
392  pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
393  for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
394  {
395  // Keep information in C-array ... for creating touchables - at least
396 
397  fpNavigator[num] = *pNavigatorIter;
398  fLimitTruth[num] = false;
399  fLimitedStep[num] = kDoNot;
400  fCurrentStepSize[num] = 0.0;
401  fLocatedVolume[num] = 0;
402  }
403  fNoGeometriesLimiting= 0; // At start of track, no process limited step
404 
405  // In case of one geometry, the tracking will have done the locating!!
406 
407  if( fNoActiveNavigators > 1 )
408  {
409  Locate( position, direction, false );
410  }
411  else
412  {
413  // Update state -- depending on the tracking's call to Mass Navigator
414 
415  fLastLocatedPosition= position;
416  fLocatedVolume[0]= massStartVol; // This information must be given
417  // by transportation
418  fLimitedStep[0] = kDoNot;
419  fCurrentStepSize[0] = 0.0;
420  }
421 
422  // Reset Safety Information -- as in case of overlaps this can cause
423  // inconsistencies ...
424  //
425  fMinSafety_PreStepPt= fPreSafetyMinValue= fMinSafety_atSafLocation= 0.0;
426 
427  for( num=0; num< fNoActiveNavigators; ++num )
428  {
429  fPreSafetyValues[num]= 0.0;
430  fNewSafetyComputed[num]= 0.0;
431  fCurrentPreStepSafety[num] = 0.0;
432  }
433 
434  // The first location for each Navigator must be non-relative
435  // or else call ResetStackAndState() for each Navigator
436 
437  fRelocatedPoint= false;
438 }
439 
440 
442  // Signal end of tracking of current track.
443  // Reset TransportationManager to use 'ordinary' Navigator
444  // Reset internal state, if needed
445 {
446  EnableParallelNavigation(false); // Else it will be continue to be used
447 }
448 
449 void G4PathFinder::ReportMove( const G4ThreeVector& OldVector,
450  const G4ThreeVector& NewVector,
451  const G4String& Quantity ) const
452 {
453  G4ThreeVector moveVec = ( NewVector - OldVector );
454 
455  G4int prc= G4cerr.precision(12);
456  std::ostringstream message;
457  message << "Endpoint moved between value returned by ComputeStep()"
458  << " and call to Locate(). " << G4endl
459  << " Change of " << Quantity << " is "
460  << moveVec.mag() / mm << " mm long" << G4endl
461  << " and its vector is "
462  << (1.0/mm) * moveVec << " mm " << G4endl
463  << " Endpoint of ComputeStep() was " << OldVector << G4endl
464  << " and current position to locate is " << NewVector;
465  G4Exception("G4PathFinder::ReportMove()", "GeomNav1002",
466  JustWarning, message);
467  G4cerr.precision(prc);
468 }
469 
470 void
472  const G4ThreeVector& direction,
473  G4bool relative)
474 {
475  // Locate the point in each geometry
476 
477  std::vector<G4Navigator*>::iterator pNavIter=
478  fpTransportManager->GetActiveNavigatorsIterator();
479 
480  G4ThreeVector lastEndPosition= fEndState.GetPosition();
481  G4ThreeVector moveVec = (position - lastEndPosition );
482  G4double moveLenSq= moveVec.mag2();
483  if( (!fNewTrack) && (!fRelocatedPoint)
484  && ( moveLenSq> 10*kCarTolerance*kCarTolerance ) )
485  {
486  ReportMove( position, lastEndPosition, "Position" );
487  }
488  fLastLocatedPosition= position;
489 
490 #ifdef G4DEBUG_PATHFINDER
491  if( fVerboseLevel > 2 )
492  {
493  G4cout << G4endl;
494  G4cout << " G4PathFinder::Locate : entered " << G4endl;
495  G4cout << " -------------------- -------" << G4endl;
496  G4cout << " Locating at position " << position
497  << " with direction " << direction
498  << " relative= " << relative << G4endl;
499  if ( (fVerboseLevel > 1) || ( moveLenSq > 0.0) )
500  {
501  G4cout << " lastEndPosition = " << lastEndPosition
502  << " moveVec = " << moveVec
503  << " newTr = " << fNewTrack
504  << " relocated = " << fRelocatedPoint << G4endl;
505  }
506 
507  G4cout << " Located at " << position ;
508  if( fNoActiveNavigators > 1 ) { G4cout << G4endl; }
509  }
510 #endif
511 
512  for ( G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
513  {
514  // ... who limited the step ....
515 
516  if( fLimitTruth[num] ) { (*pNavIter)->SetGeometricallyLimitedStep(); }
517 
518  G4VPhysicalVolume *pLocated=
519  (*pNavIter)->LocateGlobalPointAndSetup( position, &direction,
520  relative,
521  false);
522  // Set the state related to the location
523  //
524  fLocatedVolume[num] = pLocated;
525 
526  // Clear state related to the step
527  //
528  fLimitedStep[num] = kDoNot;
529  fCurrentStepSize[num] = 0.0;
530 
531 #ifdef G4DEBUG_PATHFINDER
532  if( fVerboseLevel > 2 )
533  {
534  G4cout << " - In world " << num << " geomLimStep= " << fLimitTruth[num]
535  << " gives volume= " << pLocated ;
536  if( pLocated )
537  {
538  G4cout << " name = '" << pLocated->GetName() << "'";
539  G4cout << " - CopyNo= " << pLocated->GetCopyNo();
540  }
541  G4cout << G4endl;
542  }
543 #endif
544  }
545 
546  fRelocatedPoint= false;
547 }
548 
550 {
551  // Locate the point in each geometry
552 
553  std::vector<G4Navigator*>::iterator pNavIter=
554  fpTransportManager->GetActiveNavigatorsIterator();
555 
556 #ifdef G4DEBUG_PATHFINDER
557 
558  // Check that this relocation does not violate safety
559  // - at endpoint (computed from start point) AND
560  // - at last safety location (likely just called)
561 
562  G4ThreeVector lastEndPosition= fEndState.GetPosition();
563 
564  // Calculate end-point safety ...
565  //
566  G4double DistanceStartEnd= (lastEndPosition - fPreStepLocation).mag();
567  G4double endPointSafety_raw = fMinSafety_PreStepPt - DistanceStartEnd;
568  G4double endPointSafety_Est1 = std::max( 0.0, endPointSafety_raw );
569 
570  // ... and check move from endpoint against this endpoint safety
571  //
572  G4ThreeVector moveVecEndPos = position - lastEndPosition;
573  G4double moveLenEndPosSq = moveVecEndPos.mag2();
574 
575  // Check that move from endpoint of last step is within safety
576  // -- or check against last location or relocation ??
577  //
578  G4ThreeVector moveVecSafety= position - fSafetyLocation;
579  G4double moveLenSafSq= moveVecSafety.mag2();
580 
581  G4double distCheckEnd_sq= ( moveLenEndPosSq - endPointSafety_Est1
582  *endPointSafety_Est1 );
583  G4double distCheckSaf_sq= ( moveLenSafSq - fMinSafety_atSafLocation
584  *fMinSafety_atSafLocation );
585 
586  G4bool longMoveEnd = distCheckEnd_sq > 0.0;
587  G4bool longMoveSaf = distCheckSaf_sq > 0.0;
588 
589  G4double revisedSafety= 0.0;
590 
591  if( (!fNewTrack) && ( longMoveEnd && longMoveSaf ) )
592  {
593  // Recompute ComputeSafety for end position
594  //
595  revisedSafety= ComputeSafety(lastEndPosition);
596 
597  const G4double kRadTolerance =
599  const G4double cErrorTolerance=1e-12;
600  // Maximum relative error from roundoff of arithmetic
601 
602  G4double distCheckRevisedEnd= moveLenEndPosSq-revisedSafety*revisedSafety;
603 
604  G4bool longMoveRevisedEnd= ( distCheckRevisedEnd > 0. ) ;
605 
606  G4double moveMinusSafety= 0.0;
607  G4double moveLenEndPosition= std::sqrt( moveLenEndPosSq );
608  moveMinusSafety = moveLenEndPosition - revisedSafety;
609 
610  if ( longMoveRevisedEnd && (moveMinusSafety > 0.0 )
611  && ( revisedSafety > 0.0 ) )
612  {
613  // Take into account possibility of roundoff error causing
614  // this apparent move further than safety
615 
616  if( fVerboseLevel > 0 )
617  {
618  G4cout << " G4PF:Relocate> Ratio to revised safety is "
619  << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
620  }
621 
622  G4double absMoveMinusSafety= std::fabs(moveMinusSafety);
623  G4bool smallRatio= absMoveMinusSafety < kRadTolerance * revisedSafety ;
624  G4double maxCoordPos = std::max(
625  std::max( std::fabs(position.x()),
626  std::fabs(position.y())),
627  std::fabs(position.z()) );
628  G4bool smallValue= absMoveMinusSafety < cErrorTolerance * maxCoordPos;
629  if( ! (smallRatio || smallValue) )
630  {
631  G4cout << " G4PF:Relocate> Ratio to revised safety is "
632  << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
633  G4cout << " Difference of move and safety is not very small."
634  << G4endl;
635  }
636  else
637  {
638  moveMinusSafety = 0.0;
639  longMoveRevisedEnd = false; // Numerical issue -- not too long!
640 
641  G4cout << " Difference of move & safety is very small in magnitude, "
642  << absMoveMinusSafety << G4endl;
643  if( smallRatio )
644  {
645  G4cout << " ratio to safety " << revisedSafety
646  << " is " << absMoveMinusSafety / revisedSafety
647  << "smaller than " << kRadTolerance << " of safety ";
648  }
649  else
650  {
651  G4cout << " as fraction " << absMoveMinusSafety / maxCoordPos
652  << " of position vector max-coord " << maxCoordPos
653  << " smaller than " << cErrorTolerance ;
654  }
655  G4cout << " -- reset moveMinusSafety to "
656  << moveMinusSafety << G4endl;
657  }
658  }
659 
660  if ( longMoveEnd && longMoveSaf
661  && longMoveRevisedEnd && (moveMinusSafety>0.0) )
662  {
663  G4int oldPrec= G4cout.precision(9);
664  std::ostringstream message;
665  message << "ReLocation is further than end-safety value." << G4endl
666  << " Moved from last endpoint by " << moveLenEndPosition
667  << " compared to end safety (from preStep point) = "
668  << endPointSafety_Est1 << G4endl
669  << " --> last PreSafety Location was " << fPreSafetyLocation
670  << G4endl
671  << " safety value = " << fPreSafetyMinValue << G4endl
672  << " --> last PreStep Location was " << fPreStepLocation
673  << G4endl
674  << " safety value = " << fMinSafety_PreStepPt << G4endl
675  << " --> last EndStep Location was " << lastEndPosition
676  << G4endl
677  << " safety value = " << endPointSafety_Est1
678  << " raw-value = " << endPointSafety_raw << G4endl
679  << " --> Calling again at this endpoint, we get "
680  << revisedSafety << " as safety value." << G4endl
681  << " --> last position for safety " << fSafetyLocation
682  << G4endl
683  << " its safety value = " << fMinSafety_atSafLocation
684  << G4endl
685  << " move from safety location = "
686  << std::sqrt(moveLenSafSq) << G4endl
687  << " again= " << moveVecSafety.mag() << G4endl
688  << " safety - Move-from-end= "
689  << revisedSafety - moveLenEndPosition
690  << " (negative is Bad.)" << G4endl
691  << " Debug: distCheckRevisedEnd = "
692  << distCheckRevisedEnd;
693  ReportMove( lastEndPosition, position, "Position" );
694  G4Exception("G4PathFinder::ReLocate", "GeomNav0003",
695  FatalException, message);
696  G4cout.precision(oldPrec);
697  }
698  }
699 
700  if( fVerboseLevel > 2 )
701  {
702  G4cout << G4endl;
703  G4cout << " G4PathFinder::ReLocate : entered " << G4endl;
704  G4cout << " ---------------------- -------" << G4endl;
705  G4cout << " *Re*Locating at position " << position << G4endl;
706  // << " with direction " << direction
707  // << " relative= " << relative << G4endl;
708  if ( (fVerboseLevel > -1) || ( moveLenEndPosSq > 0.0) )
709  {
710  G4cout << " lastEndPosition = " << lastEndPosition
711  << " moveVec from step-end = " << moveVecEndPos
712  << " is new Track = " << fNewTrack
713  << " relocated = " << fRelocatedPoint << G4endl;
714  }
715  }
716 #endif // G4DEBUG_PATHFINDER
717 
718  for ( G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
719  {
720  // ... none limited the step
721 
722  (*pNavIter)->LocateGlobalPointWithinVolume( position );
723 
724  // Clear state related to the step
725  //
726  fLimitedStep[num] = kDoNot;
727  fCurrentStepSize[num] = 0.0;
728  fLimitTruth[num] = false;
729  }
730 
731  fLastLocatedPosition= position;
732  fRelocatedPoint= false;
733 
734 #ifdef G4DEBUG_PATHFINDER
735  if( fVerboseLevel > 2 )
736  {
737  G4cout << " G4PathFinder::ReLocate : exiting "
738  << " at position " << fLastLocatedPosition << G4endl << G4endl;
739  }
740 #endif
741 }
742 
743 // -----------------------------------------------------------------------------
744 
746 {
747  // Recompute safety for the relevant point
748 
749  G4double minSafety= kInfinity;
750 
751  std::vector<G4Navigator*>::iterator pNavigatorIter;
752  pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
753 
754  for( G4int num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num )
755  {
756  G4double safety = (*pNavigatorIter)->ComputeSafety( position, DBL_MAX, true );
757  if( safety < minSafety ) { minSafety = safety; }
758  fNewSafetyComputed[num]= safety;
759  }
760 
761  fSafetyLocation= position;
762  fMinSafety_atSafLocation = minSafety;
763 
764 #ifdef G4DEBUG_PATHFINDER
765  if( fVerboseLevel > 1 )
766  {
767  G4cout << " G4PathFinder::ComputeSafety - returns "
768  << minSafety << " at location " << position << G4endl;
769  }
770 #endif
771  return minSafety;
772 }
773 
774 
775 // -----------------------------------------------------------------------------
776 
779 {
780 #ifdef G4DEBUG_PATHFINDER
781  if( fVerboseLevel > 2 )
782  {
783  G4cout << "G4PathFinder::CreateTouchableHandle : navId = "
784  << navId << " -- " << GetNavigator(navId) << G4endl;
785  }
786 #endif
787 
788  G4TouchableHistory* touchHist;
789  touchHist= GetNavigator(navId) -> CreateTouchableHistory();
790 
791  G4VPhysicalVolume* locatedVolume= fLocatedVolume[navId];
792  if( locatedVolume == 0 )
793  {
794  // Workaround to ensure that the touchable is fixed !! // TODO: fix
795 
796  touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() );
797  }
798 
799 #ifdef G4DEBUG_PATHFINDER
800  if( fVerboseLevel > 2 )
801  {
802  G4String VolumeName("None");
803  if( locatedVolume ) { VolumeName= locatedVolume->GetName(); }
804  G4cout << " Touchable History created at address " << touchHist
805  << " volume = " << locatedVolume << " name= " << VolumeName
806  << G4endl;
807  }
808 #endif
809 
810  return G4TouchableHandle(touchHist);
811 }
812 
813 G4double
815  G4double proposedStepLength )
816 {
817  std::vector<G4Navigator*>::iterator pNavigatorIter;
818  G4double safety= 0.0, step=0.0;
819  G4double minSafety= kInfinity, minStep;
820 
821  const G4int IdTransport= 0; // Id of Mass Navigator !!
822  G4int num=0;
823 
824 #ifdef G4DEBUG_PATHFINDER
825  if( fVerboseLevel > 2 )
826  {
827  G4cout << " G4PathFinder::DoNextLinearStep : entered " << G4endl;
828  G4cout << " Input field track= " << initialState << G4endl;
829  G4cout << " Requested step= " << proposedStepLength << G4endl;
830  }
831 #endif
832 
833  G4ThreeVector initialPosition= initialState.GetPosition();
834  G4ThreeVector initialDirection= initialState.GetMomentumDirection();
835 
836  G4ThreeVector OriginShift = initialPosition - fPreSafetyLocation;
837  G4double MagSqShift = OriginShift.mag2() ;
838  G4double MagShift; // Only given value if it larger than minimum safety
839 
840  // Potential optimisation using Maximum Value of safety!
841  // if( MagSqShift >= sqr(fPreSafetyMaxValue ) ){
842  // MagShift= kInfinity; // Not a useful value -- all will not use/ignore
843  // else
844  // MagShift= std::sqrt(MagSqShift) ;
845 
846  MagShift= std::sqrt(MagSqShift) ;
847 
848 #ifdef G4PATHFINDER_OPTIMISATION
849 
850  G4double fullSafety; // For all geometries, for prestep point
851 
852  if( MagSqShift >= sqr(fPreSafetyMinValue ) )
853  {
854  fullSafety = 0.0 ;
855  }
856  else
857  {
858  fullSafety = fPreSafetyMinValue - MagShift;
859  }
860  if( proposedStepLength < fullSafety )
861  {
862  // Move is smaller than all safeties
863  // -> so we do not have to move the safety center
864 
865  fPreStepCenterRenewed= false;
866 
867  for( num=0; num< fNoActiveNavigators; ++num )
868  {
869  fCurrentStepSize[num]= kInfinity;
870  safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
871  minSafety= std::min( safety, minSafety );
872  fCurrentPreStepSafety[num]= safety;
873  }
874  minStep= kInfinity;
875 
876 #ifdef G4DEBUG_PATHFINDER
877  if( fVerboseLevel > 2 )
878  {
879  G4cout << "G4PathFinder::DoNextLinearStep : Quick Stepping. " << G4endl
880  << " proposedStepLength " << proposedStepLength
881  << " < (full) safety = " << fullSafety
882  << " at " << initialPosition
883  << G4endl;
884  }
885 #endif
886  }
887  else
888 #endif // End of G4PATHFINDER_OPTIMISATION 1
889  {
890  // Move is larger than at least one of the safeties
891  // -> so we must move the safety center!
892 
893  fPreStepCenterRenewed= true;
894  pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator();
895 
896  minStep= kInfinity; // Not proposedStepLength;
897 
898  for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
899  {
900  safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
901 
902 #ifdef G4PATHFINDER_OPTIMISATION
903  if( proposedStepLength <= safety ) // Should be just < safety ?
904  {
905  // The Step is guaranteed to be taken
906 
907  step= kInfinity; // ComputeStep Would return this
908 
909 #ifdef G4DEBUG_PATHFINDER
910  G4cout.precision(8);
911  G4cout << "PathFinder::ComputeStep> small proposed step = "
912  << proposedStepLength
913  << " <= safety = " << safety << " for nav " << num
914  << " Step fully taken. " << G4endl;
915 #endif
916  }
917  else
918 #endif // End of G4PATHFINDER_OPTIMISATION 2
919  {
920 #ifdef G4DEBUG_PATHFINDER
921  G4double previousSafety= safety;
922 #endif
923  step= (*pNavigatorIter)->ComputeStep( initialPosition,
924  initialDirection,
925  proposedStepLength,
926  safety );
927  minStep = std::min( step, minStep);
928 
929  // TODO: consider whether/how to reduce the proposed step
930  // to the latest minStep value - to reduce calculations
931 
932 #ifdef G4DEBUG_PATHFINDER
933  if( fVerboseLevel > 0)
934  {
935  G4cout.precision(8);
936  G4cout << "PathFinder::ComputeStep> long proposed step = "
937  << proposedStepLength
938  << " > safety = " << previousSafety
939  << " for nav " << num
940  << " . New safety = " << safety << " step= " << step
941  << G4endl;
942  }
943 #endif
944  }
945  fCurrentStepSize[num] = step;
946 
947  // Save safety value, must be done for all geometries "together"
948  // (even if not recomputed using call to ComputeStep)
949  // since they share the fPreSafetyLocation
950 
951  fPreSafetyValues[num]= safety;
952  fCurrentPreStepSafety[num]= safety;
953 
954  minSafety= std::min( safety, minSafety );
955 
956 #ifdef G4DEBUG_PATHFINDER
957  if( fVerboseLevel > 2 )
958  {
959  G4cout << "G4PathFinder::DoNextLinearStep : Navigator ["
960  << num << "] -- step size " << step << G4endl;
961  }
962 #endif
963  }
964 
965  // Only change these when safety is recalculated
966  // it is good/relevant only for safety calculations
967 
968  fPreSafetyLocation= initialPosition;
969  fPreSafetyMinValue= minSafety;
970  } // end of else for if( proposedStepLength <= fullSafety)
971 
972  // For use in Relocation, need PreStep point location, min-safety
973  //
974  fPreStepLocation= initialPosition;
975  fMinSafety_PreStepPt= minSafety;
976 
977  fMinStep= minStep;
978 
979  if( fMinStep == kInfinity )
980  {
981  minStep = proposedStepLength; // Use this below for endpoint !!
982  }
983  fTrueMinStep = minStep;
984 
985  // Set the EndState
986 
987  G4ThreeVector endPosition;
988 
989  fEndState= initialState;
990  endPosition= initialPosition + minStep * initialDirection ;
991 
992 #ifdef G4DEBUG_PATHFINDER
993  if( fVerboseLevel > 1 )
994  {
995  G4cout << "G4PathFinder::DoNextLinearStep : "
996  << " initialPosition = " << initialPosition
997  << " and endPosition = " << endPosition<< G4endl;
998  }
999 #endif
1000 
1001  fEndState.SetPosition( endPosition );
1002  fEndState.SetProperTimeOfFlight( -1.000 ); // Not defined YET
1003 
1004  if( fNoActiveNavigators == 1 )
1005  {
1006  G4bool transportLimited = (fMinStep!= kInfinity);
1007  fLimitTruth[IdTransport] = transportLimited;
1008  fLimitedStep[IdTransport] = transportLimited ? kUnique : kDoNot;
1009 
1010  // Set fNoGeometriesLimiting - as WhichLimited does
1011  fNoGeometriesLimiting = transportLimited ? 1 : 0;
1012  }
1013  else
1014  {
1015  WhichLimited();
1016  }
1017 
1018 #ifdef G4DEBUG_PATHFINDER
1019  if( fVerboseLevel > 2 )
1020  {
1021  G4cout << " G4PathFinder::DoNextLinearStep : exits returning "
1022  << minStep << G4endl;
1023  G4cout << " Endpoint values = " << fEndState << G4endl;
1024  G4cout << G4endl;
1025  }
1026 #endif
1027 
1028  return minStep;
1029 }
1030 
1032 {
1033  // Flag which processes limited the step
1034 
1035  G4int num=-1, last=-1;
1036  G4int noLimited=0;
1037  ELimited shared= kSharedOther;
1038 
1039  const G4int IdTransport= 0; // Id of Mass Navigator !!
1040 
1041  // Assume that [IdTransport] is Mass / Transport
1042  //
1043  G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep)
1044  && ( fMinStep!= kInfinity) ;
1045 
1046  if( transportLimited ) {
1047  shared= kSharedTransport;
1048  }
1049 
1050  for ( num= 0; num < fNoActiveNavigators; num++ )
1051  {
1052  G4bool limitedStep;
1053 
1054  G4double step= fCurrentStepSize[num];
1055 
1056  limitedStep = ( std::fabs(step - fMinStep) < kCarTolerance )
1057  && ( step != kInfinity);
1058 
1059  fLimitTruth[ num ] = limitedStep;
1060  if( limitedStep )
1061  {
1062  noLimited++;
1063  fLimitedStep[num] = shared;
1064  last= num;
1065  }
1066  else
1067  {
1068  fLimitedStep[num] = kDoNot;
1069  }
1070  }
1071  fNoGeometriesLimiting= noLimited; // Save # processes limiting step
1072 
1073  if( (last > -1) && (noLimited == 1 ) )
1074  {
1075  fLimitedStep[ last ] = kUnique;
1076  }
1077 
1078 #ifdef G4DEBUG_PATHFINDER
1079  if( fVerboseLevel > 1 )
1080  {
1081  PrintLimited(); // --> for tracing
1082  if( fVerboseLevel > 4 ) {
1083  G4cout << " G4PathFinder::WhichLimited - exiting. " << G4endl;
1084  }
1085  }
1086 #endif
1087 }
1088 
1090 {
1091  // Report results -- for checking
1092 
1093  G4cout << "G4PathFinder::PrintLimited reports: " ;
1094  G4cout << " Minimum step (true)= " << fTrueMinStep
1095  << " reported min = " << fMinStep
1096  << G4endl;
1097  if( (fCurrentStepNo <= 2) || (fVerboseLevel>=2) )
1098  {
1099  G4cout << std::setw(5) << " Step#" << " "
1100  << std::setw(5) << " NavId" << " "
1101  << std::setw(12) << " step-size " << " "
1102  << std::setw(12) << " raw-size " << " "
1103  << std::setw(12) << " pre-safety " << " "
1104  << std::setw(15) << " Limited / flag" << " "
1105  << std::setw(15) << " World " << " "
1106  << G4endl;
1107  }
1108  G4int num;
1109  for ( num= 0; num < fNoActiveNavigators; num++ )
1110  {
1111  G4double rawStep = fCurrentStepSize[num];
1112  G4double stepLen = fCurrentStepSize[num];
1113  if( stepLen > fTrueMinStep )
1114  {
1115  stepLen = fTrueMinStep; // did not limit (went as far as asked)
1116  }
1117  G4int oldPrec= G4cout.precision(9);
1118 
1119  G4cout << std::setw(5) << fCurrentStepNo << " "
1120  << std::setw(5) << num << " "
1121  << std::setw(12) << stepLen << " "
1122  << std::setw(12) << rawStep << " "
1123  << std::setw(12) << fCurrentPreStepSafety[num] << " "
1124  << std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " ";
1125  G4String limitedStr= LimitedString(fLimitedStep[num]);
1126  G4cout << " " << std::setw(15) << limitedStr << " ";
1127  G4cout.precision(oldPrec);
1128 
1129  G4Navigator *pNav= GetNavigator( num );
1130  G4String WorldName( "Not-Set" );
1131  if (pNav)
1132  {
1133  G4VPhysicalVolume *pWorld= pNav->GetWorldVolume();
1134  if( pWorld )
1135  {
1136  WorldName = pWorld->GetName();
1137  }
1138  }
1139  G4cout << " " << WorldName ;
1140  G4cout << G4endl;
1141  }
1142 
1143  if( fVerboseLevel > 4 )
1144  {
1145  G4cout << " G4PathFinder::PrintLimited - exiting. " << G4endl;
1146  }
1147 }
1148 
1149 G4double
1151  G4double proposedStepLength,
1152  G4VPhysicalVolume* pCurrentPhysicalVolume )
1153 {
1154  const G4double toleratedRelativeError= 1.0e-10;
1155  G4double minStep= kInfinity, newSafety=0.0;
1156  G4int numNav;
1157  G4FieldTrack fieldTrack= initialState;
1158  G4ThreeVector startPoint= initialState.GetPosition();
1159 
1160 
1161  G4EquationOfMotion* equationOfMotion =
1162  (fpFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper())
1163  ->GetEquationOfMotion();
1164 
1165  equationOfMotion->SetChargeMomentumMass( *(initialState.GetChargeState()),
1166  initialState.GetMomentum().mag2(),
1167  initialState.GetRestMass() );
1168 
1169 #ifdef G4DEBUG_PATHFINDER
1170  G4int prc= G4cout.precision(9);
1171  if( fVerboseLevel > 2 )
1172  {
1173  G4cout << " G4PathFinder::DoNextCurvedStep ****** " << G4endl;
1174  G4cout << " Initial value of field track is " << fieldTrack
1175  << " and proposed step= " << proposedStepLength << G4endl;
1176  }
1177 #endif
1178 
1179  fPreStepCenterRenewed= true; // Always update PreSafety with PreStep point
1180 
1181  if( fNoActiveNavigators > 1 )
1182  {
1183  // Calculate the safety values before making the step
1184 
1185  G4double minSafety= kInfinity, safety;
1186  for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1187  {
1188  safety= fpNavigator[numNav]->ComputeSafety( startPoint, DBL_MAX, false );
1189  fPreSafetyValues[numNav]= safety;
1190  fCurrentPreStepSafety[numNav]= safety;
1191  minSafety = std::min( safety, minSafety );
1192  }
1193 
1194  // Save safety value, related position
1195 
1196  fPreSafetyLocation= startPoint;
1197  fPreSafetyMinValue= minSafety;
1198  fPreStepLocation= startPoint;
1199  fMinSafety_PreStepPt= minSafety;
1200  }
1201 
1202  // Allow Propagator In Field to do the hard work, calling G4MultiNavigator
1203  //
1204  minStep= fpFieldPropagator->ComputeStep( fieldTrack,
1205  proposedStepLength,
1206  newSafety,
1207  pCurrentPhysicalVolume );
1208 
1209  // fieldTrack now contains the endpoint information
1210  //
1211  fEndState= fieldTrack;
1212  fMinStep= minStep;
1213  fTrueMinStep = std::min( minStep, proposedStepLength );
1214 
1215  if( fNoActiveNavigators== 1 )
1216  {
1217  // Update the 'PreSafety' sphere - as any ComputeStep was called
1218  // (must be done anyway in field)
1219 
1220  fPreSafetyValues[0]= newSafety;
1221  fPreSafetyLocation= startPoint;
1222  fPreSafetyMinValue= newSafety;
1223 
1224  // Update the current 'PreStep' point's values - mandatory
1225  //
1226  fCurrentPreStepSafety[0]= newSafety;
1227  fPreStepLocation= startPoint;
1228  fMinSafety_PreStepPt= newSafety;
1229  }
1230 
1231 #ifdef G4DEBUG_PATHFINDER
1232  if( fVerboseLevel > 2 )
1233  {
1234  G4cout << "G4PathFinder::DoNextCurvedStep : " << G4endl
1235  << " initialState = " << initialState << G4endl
1236  << " and endState = " << fEndState << G4endl;
1237  G4cout << "G4PathFinder::DoNextCurvedStep : "
1238  << " minStep = " << minStep
1239  << " proposedStepLength " << proposedStepLength
1240  << " safety = " << newSafety << G4endl;
1241  }
1242 #endif
1243  G4double currentStepSize; // = 0.0;
1244  if( minStep < proposedStepLength ) // if == , then a boundary found at end ??
1245  {
1246  // Recover the remaining information from MultiNavigator
1247  // especially regarding which Navigator limited the step
1248 
1249  G4int noLimited= 0; // No geometries limiting step
1250  for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1251  {
1252  G4double finalStep, lastPreSafety=0.0, minStepLast;
1253  ELimited didLimit;
1254  G4bool limited;
1255 
1256  finalStep= fpMultiNavigator->ObtainFinalStep( numNav, lastPreSafety,
1257  minStepLast, didLimit );
1258 
1259  // Calculate the step for this geometry, using the
1260  // final step (the only one which can differ.)
1261 
1262  currentStepSize = fTrueMinStep;
1263  G4double diffStep= 0.0;
1264  if( (minStepLast != kInfinity) )
1265  {
1266  diffStep = (finalStep-minStepLast);
1267  if ( std::abs(diffStep) <= toleratedRelativeError * finalStep )
1268  {
1269  diffStep = 0.0;
1270  }
1271  currentStepSize += diffStep;
1272  }
1273  fCurrentStepSize[numNav] = currentStepSize;
1274 
1275  // TODO: could refine the way to obtain safeties for > 1 geometries
1276  // - for pre step safety
1277  // notify MultiNavigator about new set of sub-steps
1278  // allow it to return this value in ObtainFinalStep
1279  // instead of lastPreSafety (or as well?)
1280  // - for final step start (available)
1281  // get final Step start from MultiNavigator
1282  // and corresponding safety values
1283  // and/or ALSO calculate ComputeSafety at endpoint
1284  // endSafety= fpNavigator[numNav]->ComputeSafety( endPoint );
1285 
1286  fLimitedStep[numNav] = didLimit;
1287  fLimitTruth[numNav] = limited = (didLimit != kDoNot );
1288  if( limited ) { noLimited++; }
1289 
1290 #ifdef G4DEBUG_PATHFINDER
1291  G4bool StepError= (currentStepSize < 0)
1292  || ( (minStepLast != kInfinity) && (diffStep < 0) ) ;
1293  if( StepError || (fVerboseLevel > 2) )
1294  {
1295  G4String limitedString= LimitedString( fLimitedStep[numNav] );
1296 
1297  G4cout << " G4PathFinder::ComputeStep. Geometry " << numNav
1298  << " step= " << fCurrentStepSize[numNav]
1299  << " from final-step= " << finalStep
1300  << " fTrueMinStep= " << fTrueMinStep
1301  << " minStepLast= " << minStepLast
1302  << " limited = " << (fLimitTruth[numNav] ? "YES" : " NO")
1303  << " ";
1304  G4cout << " status = " << limitedString << " #= " << didLimit
1305  << G4endl;
1306 
1307  if( StepError )
1308  {
1309  std::ostringstream message;
1310  message << "Incorrect calculation of step size for one navigator"
1311  << G4endl
1312  << " currentStepSize = " << currentStepSize
1313  << ", diffStep= " << diffStep << G4endl
1314  << "ERROR in computing step size for this navigator.";
1315  G4Exception("G4PathFinder::DoNextCurvedStep",
1316  "GeomNav0003", FatalException, message);
1317  }
1318  }
1319 #endif
1320  } // for num Navigators
1321 
1322  fNoGeometriesLimiting= noLimited; // Save # processes limiting step
1323  }
1324  else if ( (minStep == proposedStepLength)
1325  || (minStep == kInfinity)
1326  || ( std::abs(minStep-proposedStepLength)
1327  < toleratedRelativeError * proposedStepLength ) )
1328  {
1329  // In case the step was not limited, use default responses
1330  // --> all Navigators
1331  // Also avoid problems in case of PathFinder using safety to optimise
1332  // - it is possible that the Navigators were not called
1333  // if the safety was already satisfactory.
1334  // (In that case calling ObtainFinalStep gives invalid results.)
1335 
1336  currentStepSize= minStep;
1337  for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1338  {
1339  fCurrentStepSize[numNav] = minStep;
1340  // Safety for endpoint ?? // Can eventuall improve it -- see TODO above
1341  fLimitedStep[numNav] = kDoNot;
1342  fLimitTruth[numNav] = false;
1343  }
1344  fNoGeometriesLimiting= 0; // Save # processes limiting step
1345  }
1346  else // (minStep > proposedStepLength) and not (minStep == kInfinity)
1347  {
1348  std::ostringstream message;
1349  message << "Incorrect calculation of step size for one navigator." << G4endl
1350  << " currentStepSize = " << minStep << " is larger than "
1351  << " proposed StepSize = " << proposedStepLength << ".";
1352  G4Exception("G4PathFinder::DoNextCurvedStep()",
1353  "GeomNav0003", FatalException, message);
1354  }
1355 
1356 #ifdef G4DEBUG_PATHFINDER
1357  if( fVerboseLevel > 2 )
1358  {
1359  G4cout << " Exiting G4PathFinder::DoNextCurvedStep " << G4endl;
1360  PrintLimited();
1361  }
1362  G4cout.precision(prc);
1363 #endif
1364 
1365  return minStep;
1366 }
1367 
1368 
1370  const G4ThreeVector &pGlobalPoint,
1371  const G4ThreeVector &pDirection,
1372  const G4double aProposedMove,
1373  G4double *prDistance,
1374  G4double *prNewSafety)const
1375 {
1376  G4bool retval= true;
1377 
1378  if( fNoActiveNavigators > 0 )
1379  {
1380  // Calculate the safety values before making the step
1381 
1382  G4double minSafety= kInfinity;
1383  G4double minMove= kInfinity;
1384  int numNav;
1385  for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
1386  {
1387  G4double distance, safety;
1388  G4bool moveIsOK;
1389  moveIsOK= fpNavigator[numNav]->RecheckDistanceToCurrentBoundary(
1390  pGlobalPoint,
1391  pDirection,
1392  aProposedMove,
1393  &distance,
1394  &safety);
1395  minSafety = std::min( safety, minSafety );
1396  minMove = std::min( distance, minMove );
1397  // The first surface encountered will determine it
1398  // - even if it is at a negative distance.
1399  retval &= moveIsOK;
1400  }
1401 
1402  *prDistance= minMove;
1403  if( prNewSafety ) *prNewSafety= minSafety;
1404 
1405  }else{
1406  retval= false;
1407  }
1408 
1409  return retval;
1410 }
1411 
1412 
1413 
1415 {
1416  static G4String StrDoNot("DoNot"),
1417  StrUnique("Unique"),
1418  StrUndefined("Undefined"),
1419  StrSharedTransport("SharedTransport"),
1420  StrSharedOther("SharedOther");
1421 
1422  G4String* limitedStr;
1423  switch ( lim )
1424  {
1425  case kDoNot: limitedStr= &StrDoNot; break;
1426  case kUnique: limitedStr = &StrUnique; break;
1427  case kSharedTransport: limitedStr= &StrSharedTransport; break;
1428  case kSharedOther: limitedStr = &StrSharedOther; break;
1429  default: limitedStr = &StrUndefined; break;
1430  }
1431  return *limitedStr;
1432 }
1433 
1435 {
1436  fPreSafetyLocation= fSafetyLocation;
1437  fPreSafetyMinValue= fMinSafety_atSafLocation;
1438  for( G4int nav=0; nav < fNoActiveNavigators; ++nav )
1439  {
1440  fPreSafetyValues[nav]= fNewSafetyComputed[nav];
1441  }
1442 }
void PrepareNewTrack(const G4ThreeVector &position, const G4ThreeVector &direction, G4VPhysicalVolume *massStartVol=0)
G4String & LimitedString(ELimited lim)
static G4PathFinder * GetInstance()
Definition: G4PathFinder.cc:57
void SetPosition(G4ThreeVector nPos)
void Locate(const G4ThreeVector &position, const G4ThreeVector &direction, G4bool relativeSearch=true)
void SetNavigatorForPropagating(G4Navigator *SimpleOrMultiNavigator)
G4SafetyHelper * GetSafetyHelper() const
void InitialiseHelper()
void EnableParallelNavigation(G4bool enableChoice=true)
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
virtual void SetChargeMomentumMass(G4ChargeState particleCharge, G4double MomentumXc, G4double MassXc2)=0
double x() const
const G4MagIntegratorStepper * GetStepper() const
void ReLocate(const G4ThreeVector &position)
ELimited
G4double GetSurfaceTolerance() const
G4Navigator * GetNavigatorForTracking() const
G4TouchableHandle CreateTouchableHandle(G4int navId) const
void WhichLimited()
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
#define fFieldExertedForce
void MovePoint()
double z() const
G4double DoNextCurvedStep(const G4FieldTrack &FieldTrack, G4double proposedStepLength, G4VPhysicalVolume *pCurrentPhysVolume)
G4ReferenceCountedHandle< G4VTouchable > G4TouchableHandle
#define position
Definition: xmlparse.cc:622
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
G4double GetRadialTolerance() const
bool G4bool
Definition: G4Types.hh:79
void EnableParallelNavigation(G4bool parallel)
G4double GetRestMass() const
G4double ComputeStep(const G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4int navigatorId, G4int stepNo, G4double &pNewSafety, ELimited &limitedStep, G4FieldTrack &EndState, G4VPhysicalVolume *currentVolume)
virtual G4bool RecheckDistanceToCurrentBoundary(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double CurrentProposedStepLength, G4double *prDistance, G4double *prNewSafety=0) const
G4double DoNextLinearStep(const G4FieldTrack &FieldTrack, G4double proposedStepLength)
G4double ComputeSafety(const G4ThreeVector &globalPoint)
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
void PushPostSafetyToPreSafety()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
G4bool RecheckDistanceToCurrentBoundary(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double *prDistance, G4double *prNewSafety=0) const
G4double GetCharge() const
#define fCurrentStepNo
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4Navigator * >::iterator GetActiveNavigatorsIterator()
G4ChordFinder * GetChordFinder()
double y() const
G4double ObtainFinalStep(G4int navigatorId, G4double &pNewSafety, G4double &minStepLast, ELimited &limitedStep)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
#define G4endl
Definition: G4ios.hh:61
void UpdateYourself(G4VPhysicalVolume *pPhysVol, const G4NavigationHistory *history=0)
#define fLastStepNo
const G4NavigationHistory * GetHistory() const
T sqr(const T &x)
Definition: templates.hh:145
#define fEndState
G4ThreeVector GetMomentum() const
double G4double
Definition: G4Types.hh:76
void PrintLimited()
const G4Field * GetDetectorField() const
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
#define fRelocatedPoint
G4MagInt_Driver * GetIntegrationDriver()
G4PropagatorInField * GetPropagatorInField() const
double mag() const
#define DBL_MAX
Definition: templates.hh:83
const G4ChargeState * GetChargeState() const
void ReportMove(const G4ThreeVector &OldV, const G4ThreeVector &NewV, const G4String &Quantity) const
G4VPhysicalVolume * GetWorldVolume() const
G4Navigator * GetNavigator(G4int n) const
static G4GeometryTolerance * GetInstance()
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GetMomentumDirection() const
void SetProperTimeOfFlight(G4double tofProper)