source: trunk/source/geometry/navigation/src/G4PathFinder.cc@ 1354

Last change on this file since 1354 was 1347, checked in by garnier, 15 years ago

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