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