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