00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <iomanip>
00037
00038 #include "G4PathFinder.hh"
00039
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4GeometryTolerance.hh"
00042 #include "G4Navigator.hh"
00043 #include "G4PropagatorInField.hh"
00044 #include "G4TransportationManager.hh"
00045 #include "G4MultiNavigator.hh"
00046 #include "G4SafetyHelper.hh"
00047
00048
00049
00050 G4PathFinder* G4PathFinder::fpPathFinder=0;
00051
00052
00053
00054
00055
00056
00057 G4PathFinder* G4PathFinder::GetInstance()
00058 {
00059 static G4PathFinder theInstance;
00060 if( ! fpPathFinder )
00061 {
00062 fpPathFinder = &theInstance;
00063 }
00064 return fpPathFinder;
00065 }
00066
00067
00068
00069
00070 G4PathFinder::G4PathFinder()
00071 : fEndState( G4ThreeVector(), G4ThreeVector(), 0., 0., 0., 0., 0.),
00072 fFieldExertedForce(false),
00073 fRelocatedPoint(true),
00074 fLastStepNo(-1), fCurrentStepNo(-1),
00075 fVerboseLevel(0)
00076 {
00077 fpMultiNavigator= new G4MultiNavigator();
00078
00079 fpTransportManager= G4TransportationManager::GetTransportationManager();
00080 fpFieldPropagator = fpTransportManager->GetPropagatorInField();
00081
00082 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00083
00084 fNoActiveNavigators= 0;
00085 G4ThreeVector Big3Vector( kInfinity, kInfinity, kInfinity );
00086 fLastLocatedPosition= Big3Vector;
00087 fSafetyLocation= Big3Vector;
00088 fPreSafetyLocation= Big3Vector;
00089 fPreStepLocation= Big3Vector;
00090
00091 fPreSafetyMinValue= -1.0;
00092 fMinSafety_PreStepPt= -1.0;
00093 fMinSafety_atSafLocation= -1.0;
00094 fMinStep= -1.0;
00095 fTrueMinStep= -1.0;
00096 fPreStepCenterRenewed= false;
00097 fNewTrack= false;
00098 fNoGeometriesLimiting= 0;
00099
00100 for( register int num=0; num< fMaxNav; ++num )
00101 {
00102 fpNavigator[num] = 0;
00103 fLimitTruth[num] = false;
00104 fLimitedStep[num] = kUndefLimited;
00105 fCurrentStepSize[num] = -1.0;
00106 fLocatedVolume[num] = 0;
00107 fPreSafetyValues[num]= -1.0;
00108 fCurrentPreStepSafety[num] = -1.0;
00109 fNewSafetyComputed[num]= -1.0;
00110 }
00111 }
00112
00113
00114
00115
00116 G4PathFinder::~G4PathFinder()
00117 {
00118 delete fpMultiNavigator;
00119 }
00120
00121
00122
00123 void
00124 G4PathFinder::EnableParallelNavigation(G4bool enableChoice)
00125 {
00126 G4Navigator *navigatorForPropagation=0, *massNavigator=0;
00127
00128 massNavigator= fpTransportManager->GetNavigatorForTracking();
00129 if( enableChoice )
00130 {
00131 navigatorForPropagation= fpMultiNavigator;
00132
00133
00134
00135 fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(true);
00136 }
00137 else
00138 {
00139 navigatorForPropagation= massNavigator;
00140
00141
00142
00143 fpTransportManager->GetSafetyHelper()->EnableParallelNavigation(false);
00144 }
00145 fpFieldPropagator->SetNavigatorForPropagating(navigatorForPropagation);
00146 }
00147
00148
00149
00150 G4double
00151 G4PathFinder::ComputeStep( const G4FieldTrack &InitialFieldTrack,
00152 G4double proposedStepLength,
00153 G4int navigatorNo,
00154 G4int stepNo,
00155 G4double &pNewSafety,
00156 ELimited &limitedStep,
00157 G4FieldTrack &EndState,
00158 G4VPhysicalVolume* currentVolume)
00159 {
00160 G4double possibleStep= -1.0;
00161
00162 #ifdef G4DEBUG_PATHFINDER
00163 if( fVerboseLevel > 2 )
00164 {
00165 G4cout << " -------------------------" << G4endl;
00166 G4cout << " G4PathFinder::ComputeStep - entered " << G4endl;
00167 G4cout << " - stepNo = " << std::setw(4) << stepNo << " "
00168 << " navigatorId = " << std::setw(2) << navigatorNo << " "
00169 << " proposed step len = " << proposedStepLength << " " << G4endl;
00170 G4cout << " PF::ComputeStep requested step "
00171 << " from " << InitialFieldTrack.GetPosition()
00172 << " dir " << InitialFieldTrack.GetMomentumDirection() << G4endl;
00173 }
00174 #endif
00175 #ifdef G4VERBOSE
00176 if( navigatorNo >= fNoActiveNavigators )
00177 {
00178 std::ostringstream message;
00179 message << "Bad Navigator ID !" << G4endl
00180 << " Requested Navigator ID = " << navigatorNo << G4endl
00181 << " Number of active navigators = " << fNoActiveNavigators;
00182 G4Exception("G4PathFinder::ComputeStep()", "GeomNav0002",
00183 FatalException, message);
00184 }
00185 #endif
00186
00187 if( fNewTrack || (stepNo != fLastStepNo) )
00188 {
00189
00190
00191
00192 G4FieldTrack currentState= InitialFieldTrack;
00193
00194 fCurrentStepNo = stepNo;
00195
00196
00197
00198
00199 G4ThreeVector newPosition = InitialFieldTrack.GetPosition();
00200 G4ThreeVector moveVector= newPosition - fLastLocatedPosition;
00201 G4double moveLenSq= moveVector.mag2();
00202 if( moveLenSq > kCarTolerance * kCarTolerance )
00203 {
00204 G4ThreeVector newDirection = InitialFieldTrack.GetMomentumDirection();
00205 #ifdef G4DEBUG_PATHFINDER
00206 if( fVerboseLevel > 2 )
00207 {
00208 G4double moveLen= std::sqrt( moveLenSq );
00209 G4cout << " G4PathFinder::ComputeStep : Point moved since last step "
00210 << " -- at step # = " << stepNo << G4endl
00211 << " by " << moveLen << " to " << newPosition << G4endl;
00212 }
00213 #endif
00214 MovePoint();
00215
00216
00217
00218 Locate( newPosition, newDirection );
00219 }
00220
00221
00222
00223 G4double particleCharge= currentState.GetCharge();
00224
00225 G4FieldManager* fieldMgr=0;
00226 G4bool fieldExertsForce = false ;
00227 if( (particleCharge != 0.0) )
00228 {
00229 fieldMgr= fpFieldPropagator->FindAndSetFieldManager( currentVolume );
00230
00231
00232
00233 fieldExertsForce = (fieldMgr != 0)
00234 && (fieldMgr->GetDetectorField() != 0);
00235 }
00236 fFieldExertedForce = fieldExertsForce;
00237
00238
00239 fNoGeometriesLimiting= -1;
00240 if( fieldExertsForce )
00241 {
00242 DoNextCurvedStep( currentState, proposedStepLength, currentVolume );
00243
00244 }else{
00245 DoNextLinearStep( currentState, proposedStepLength );
00246
00247 }
00248 fLastStepNo= stepNo;
00249
00250 #ifdef G4DEBUG_PATHFINDER
00251 if ( (fNoGeometriesLimiting < 0)
00252 || (fNoGeometriesLimiting > fNoActiveNavigators) )
00253 {
00254 std::ostringstream message;
00255 message << "Number of geometries limiting the step not set." << G4endl
00256 << " Number of geometries limiting step = "
00257 << fNoGeometriesLimiting;
00258 G4Exception("G4PathFinder::ComputeStep()",
00259 "GeomNav0002", FatalException, message);
00260 }
00261 #endif
00262 }
00263 #ifdef G4DEBUG_PATHFINDER
00264 else
00265 {
00266 if( proposedStepLength < fTrueMinStep )
00267 {
00268 std::ostringstream message;
00269 message << "Problem in step size request." << G4endl
00270 << " Error can be caused by incorrect process ordering."
00271 << " Being requested to make a step which is shorter"
00272 << " than the minimum Step " << G4endl
00273 << " already computed for any Navigator/geometry during"
00274 << " this tracking-step: " << G4endl
00275 << " This can happen due to an error in process ordering."
00276 << G4endl
00277 << " Check that all physics processes are registered"
00278 << G4endl
00279 << " before all processes with a navigator/geometry."
00280 << G4endl
00281 << " If using pre-packaged physics list and/or"
00282 << G4endl
00283 << " functionality, please report this error."
00284 << G4endl << G4endl
00285 << " Additional information for problem: " << G4endl
00286 << " Steps request/proposed = " << proposedStepLength
00287 << G4endl
00288 << " MinimumStep (true) = " << fTrueMinStep
00289 << G4endl
00290 << " MinimumStep (navraw) = " << fMinStep
00291 << G4endl
00292 << " Navigator raw return value" << G4endl
00293 << " Requested step now = " << proposedStepLength
00294 << G4endl
00295 << " Difference min-req = "
00296 << fTrueMinStep-proposedStepLength << G4endl
00297 << " -- Step info> stepNo= " << stepNo
00298 << " last= " << fLastStepNo
00299 << " newTr= " << fNewTrack << G4endl;
00300 G4Exception("G4PathFinder::ComputeStep()",
00301 "GeomNav0003", FatalException, message);
00302 }
00303 else
00304 {
00305
00306
00307
00308
00309
00310 if( fVerboseLevel > 1 )
00311 {
00312 G4cout << " G4P::CS -> Not calling DoNextLinearStep: "
00313 << " stepNo= " << stepNo << " last= " << fLastStepNo
00314 << " new= " << fNewTrack << " Step already done" << G4endl;
00315 }
00316 }
00317 }
00318 #endif
00319
00320 fNewTrack= false;
00321
00322
00323
00324 pNewSafety = fCurrentPreStepSafety[ navigatorNo ];
00325 limitedStep = fLimitedStep[ navigatorNo ];
00326 fRelocatedPoint= false;
00327
00328 possibleStep= std::min(proposedStepLength, fCurrentStepSize[ navigatorNo ]);
00329 EndState = fEndState;
00330
00331 #ifdef G4DEBUG_PATHFINDER
00332 if( fVerboseLevel > 0 )
00333 {
00334 G4cout << " G4PathFinder::ComputeStep returns "
00335 << fCurrentStepSize[ navigatorNo ]
00336 << " for Navigator " << navigatorNo
00337 << " Limited step = " << limitedStep
00338 << " Safety(mm) = " << pNewSafety / mm
00339 << G4endl;
00340 }
00341 #endif
00342
00343 return possibleStep;
00344 }
00345
00346
00347
00348 void
00349 G4PathFinder::PrepareNewTrack( const G4ThreeVector& position,
00350 const G4ThreeVector& direction,
00351 G4VPhysicalVolume* massStartVol)
00352 {
00353
00354
00355
00356
00357 G4int num=0;
00358
00359 EnableParallelNavigation(true);
00360
00361
00362 fpTransportManager->GetSafetyHelper()->InitialiseHelper();
00363
00364
00365 fNewTrack= true;
00366 this->MovePoint();
00367
00368
00369
00370 std::vector<G4Navigator*>::iterator pNavigatorIter;
00371
00372 fNoActiveNavigators= fpTransportManager-> GetNoActiveNavigators();
00373 if( fNoActiveNavigators > fMaxNav )
00374 {
00375 std::ostringstream message;
00376 message << "Too many active Navigators / worlds." << G4endl
00377 << " Transportation Manager has "
00378 << fNoActiveNavigators << " active navigators." << G4endl
00379 << " This is more than the number allowed = "
00380 << fMaxNav << " !";
00381 G4Exception("G4PathFinder::PrepareNewTrack()", "GeomNav0002",
00382 FatalException, message);
00383 }
00384
00385 fpMultiNavigator->PrepareNavigators();
00386
00387
00388 pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
00389 for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
00390 {
00391
00392
00393 fpNavigator[num] = *pNavigatorIter;
00394 fLimitTruth[num] = false;
00395 fLimitedStep[num] = kDoNot;
00396 fCurrentStepSize[num] = 0.0;
00397 fLocatedVolume[num] = 0;
00398 }
00399 fNoGeometriesLimiting= 0;
00400
00401
00402
00403 if( fNoActiveNavigators > 1 )
00404 {
00405 Locate( position, direction, false );
00406 }
00407 else
00408 {
00409
00410
00411 fLastLocatedPosition= position;
00412 fLocatedVolume[0]= massStartVol;
00413
00414 fLimitedStep[0] = kDoNot;
00415 fCurrentStepSize[0] = 0.0;
00416 }
00417
00418
00419
00420
00421 fMinSafety_PreStepPt= fPreSafetyMinValue= fMinSafety_atSafLocation= 0.0;
00422
00423 for( num=0; num< fNoActiveNavigators; ++num )
00424 {
00425 fPreSafetyValues[num]= 0.0;
00426 fNewSafetyComputed[num]= 0.0;
00427 fCurrentPreStepSafety[num] = 0.0;
00428 }
00429
00430
00431
00432
00433 fRelocatedPoint= false;
00434 }
00435
00436 void G4PathFinder::ReportMove( const G4ThreeVector& OldVector,
00437 const G4ThreeVector& NewVector,
00438 const G4String& Quantity ) const
00439 {
00440 G4ThreeVector moveVec = ( NewVector - OldVector );
00441
00442 G4int prc= G4cerr.precision(12);
00443 std::ostringstream message;
00444 message << "Endpoint moved between value returned by ComputeStep()"
00445 << " and call to Locate(). " << G4endl
00446 << " Change of " << Quantity << " is "
00447 << moveVec.mag() / mm << " mm long" << G4endl
00448 << " and its vector is "
00449 << (1.0/mm) * moveVec << " mm " << G4endl
00450 << " Endpoint of ComputeStep() was " << OldVector << G4endl
00451 << " and current position to locate is " << NewVector;
00452 G4Exception("G4PathFinder::ReportMove()", "GeomNav1002",
00453 JustWarning, message);
00454 G4cerr.precision(prc);
00455 }
00456
00457 void
00458 G4PathFinder::Locate( const G4ThreeVector& position,
00459 const G4ThreeVector& direction,
00460 G4bool relative)
00461 {
00462
00463
00464 std::vector<G4Navigator*>::iterator pNavIter=
00465 fpTransportManager->GetActiveNavigatorsIterator();
00466
00467 G4ThreeVector lastEndPosition= fEndState.GetPosition();
00468 G4ThreeVector moveVec = (position - lastEndPosition );
00469 G4double moveLenSq= moveVec.mag2();
00470 if( (!fNewTrack) && (!fRelocatedPoint)
00471 && ( moveLenSq> kCarTolerance*kCarTolerance ) )
00472 {
00473 ReportMove( position, lastEndPosition, "Position" );
00474 }
00475 fLastLocatedPosition= position;
00476
00477 #ifdef G4DEBUG_PATHFINDER
00478 if( fVerboseLevel > 2 )
00479 {
00480 G4cout << G4endl;
00481 G4cout << " G4PathFinder::Locate : entered " << G4endl;
00482 G4cout << " -------------------- -------" << G4endl;
00483 G4cout << " Locating at position " << position
00484 << " with direction " << direction
00485 << " relative= " << relative << G4endl;
00486 if ( (fVerboseLevel > 1) || ( moveLenSq > 0.0) )
00487 {
00488 G4cout << " lastEndPosition = " << lastEndPosition
00489 << " moveVec = " << moveVec
00490 << " newTr = " << fNewTrack
00491 << " relocated = " << fRelocatedPoint << G4endl;
00492 }
00493
00494 G4cout << " Located at " << position ;
00495 if( fNoActiveNavigators > 1 ) { G4cout << G4endl; }
00496 }
00497 #endif
00498
00499 for ( register G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
00500 {
00501
00502
00503 if( fLimitTruth[num] ) { (*pNavIter)->SetGeometricallyLimitedStep(); }
00504
00505 G4VPhysicalVolume *pLocated=
00506 (*pNavIter)->LocateGlobalPointAndSetup( position, &direction,
00507 relative,
00508 false);
00509
00510
00511 fLocatedVolume[num] = pLocated;
00512
00513
00514
00515 fLimitedStep[num] = kDoNot;
00516 fCurrentStepSize[num] = 0.0;
00517
00518 #ifdef G4DEBUG_PATHFINDER
00519 if( fVerboseLevel > 2 )
00520 {
00521 G4cout << " - In world " << num << " geomLimStep= " << fLimitTruth[num]
00522 << " gives volume= " << pLocated ;
00523 if( pLocated )
00524 {
00525 G4cout << " name = '" << pLocated->GetName() << "'";
00526 G4cout << " - CopyNo= " << pLocated->GetCopyNo();
00527 }
00528 G4cout << G4endl;
00529 }
00530 #endif
00531 }
00532
00533 fRelocatedPoint= false;
00534 }
00535
00536 void G4PathFinder::ReLocate( const G4ThreeVector& position )
00537 {
00538
00539
00540 std::vector<G4Navigator*>::iterator pNavIter=
00541 fpTransportManager->GetActiveNavigatorsIterator();
00542
00543 #ifdef G4DEBUG_PATHFINDER
00544
00545
00546
00547
00548
00549 G4ThreeVector lastEndPosition= fEndState.GetPosition();
00550
00551
00552
00553 G4double DistanceStartEnd= (lastEndPosition - fPreStepLocation).mag();
00554 G4double endPointSafety_raw = fMinSafety_PreStepPt - DistanceStartEnd;
00555 G4double endPointSafety_Est1 = std::max( 0.0, endPointSafety_raw );
00556
00557
00558
00559 G4ThreeVector moveVecEndPos = position - lastEndPosition;
00560 G4double moveLenEndPosSq = moveVecEndPos.mag2();
00561
00562
00563
00564
00565 G4ThreeVector moveVecSafety= position - fSafetyLocation;
00566 G4double moveLenSafSq= moveVecSafety.mag2();
00567
00568 G4double distCheckEnd_sq= ( moveLenEndPosSq - endPointSafety_Est1
00569 *endPointSafety_Est1 );
00570 G4double distCheckSaf_sq= ( moveLenSafSq - fMinSafety_atSafLocation
00571 *fMinSafety_atSafLocation );
00572
00573 G4bool longMoveEnd = distCheckEnd_sq > 0.0;
00574 G4bool longMoveSaf = distCheckSaf_sq > 0.0;
00575
00576 G4double revisedSafety= 0.0;
00577
00578 if( (!fNewTrack) && ( longMoveEnd && longMoveSaf ) )
00579 {
00580
00581
00582 revisedSafety= ComputeSafety(lastEndPosition);
00583
00584 const G4double kRadTolerance =
00585 G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00586 const G4double cErrorTolerance=1e-12;
00587
00588
00589 G4double distCheckRevisedEnd= moveLenEndPosSq-revisedSafety*revisedSafety;
00590
00591 G4bool longMoveRevisedEnd= ( distCheckRevisedEnd > 0. ) ;
00592
00593 G4double moveMinusSafety= 0.0;
00594 G4double moveLenEndPosition= std::sqrt( moveLenEndPosSq );
00595 moveMinusSafety = moveLenEndPosition - revisedSafety;
00596
00597 if ( longMoveRevisedEnd && (moveMinusSafety > 0.0 )
00598 && ( revisedSafety > 0.0 ) )
00599 {
00600
00601
00602
00603 if( fVerboseLevel > 0 )
00604 {
00605 G4cout << " G4PF:Relocate> Ratio to revised safety is "
00606 << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
00607 }
00608
00609 G4double absMoveMinusSafety= std::fabs(moveMinusSafety);
00610 G4bool smallRatio= absMoveMinusSafety < kRadTolerance * revisedSafety ;
00611 G4double maxCoordPos = std::max(
00612 std::max( std::fabs(position.x()),
00613 std::fabs(position.y())),
00614 std::fabs(position.z()) );
00615 G4bool smallValue= absMoveMinusSafety < cErrorTolerance * maxCoordPos;
00616 if( ! (smallRatio || smallValue) )
00617 {
00618 G4cout << " G4PF:Relocate> Ratio to revised safety is "
00619 << std::fabs(moveMinusSafety)/revisedSafety << G4endl;
00620 G4cout << " Difference of move and safety is not very small."
00621 << G4endl;
00622 }
00623 else
00624 {
00625 moveMinusSafety = 0.0;
00626 longMoveRevisedEnd = false;
00627
00628 G4cout << " Difference of move & safety is very small in magnitude, "
00629 << absMoveMinusSafety << G4endl;
00630 if( smallRatio )
00631 {
00632 G4cout << " ratio to safety " << revisedSafety
00633 << " is " << absMoveMinusSafety / revisedSafety
00634 << "smaller than " << kRadTolerance << " of safety ";
00635 }
00636 else
00637 {
00638 G4cout << " as fraction " << absMoveMinusSafety / maxCoordPos
00639 << " of position vector max-coord " << maxCoordPos
00640 << " smaller than " << cErrorTolerance ;
00641 }
00642 G4cout << " -- reset moveMinusSafety to "
00643 << moveMinusSafety << G4endl;
00644 }
00645 }
00646
00647 if ( longMoveEnd && longMoveSaf
00648 && longMoveRevisedEnd && (moveMinusSafety>0.0) )
00649 {
00650 G4int oldPrec= G4cout.precision(9);
00651 std::ostringstream message;
00652 message << "ReLocation is further than end-safety value." << G4endl
00653 << " Moved from last endpoint by " << moveLenEndPosition
00654 << " compared to end safety (from preStep point) = "
00655 << endPointSafety_Est1 << G4endl
00656 << " --> last PreSafety Location was " << fPreSafetyLocation
00657 << G4endl
00658 << " safety value = " << fPreSafetyMinValue << G4endl
00659 << " --> last PreStep Location was " << fPreStepLocation
00660 << G4endl
00661 << " safety value = " << fMinSafety_PreStepPt << G4endl
00662 << " --> last EndStep Location was " << lastEndPosition
00663 << G4endl
00664 << " safety value = " << endPointSafety_Est1
00665 << " raw-value = " << endPointSafety_raw << G4endl
00666 << " --> Calling again at this endpoint, we get "
00667 << revisedSafety << " as safety value." << G4endl
00668 << " --> last position for safety " << fSafetyLocation
00669 << G4endl
00670 << " its safety value = " << fMinSafety_atSafLocation
00671 << G4endl
00672 << " move from safety location = "
00673 << std::sqrt(moveLenSafSq) << G4endl
00674 << " again= " << moveVecSafety.mag() << G4endl
00675 << " safety - Move-from-end= "
00676 << revisedSafety - moveLenEndPosition
00677 << " (negative is Bad.)" << G4endl
00678 << " Debug: distCheckRevisedEnd = "
00679 << distCheckRevisedEnd;
00680 ReportMove( lastEndPosition, position, "Position" );
00681 G4Exception("G4PathFinder::ReLocate", "GeomNav0003",
00682 FatalException, message);
00683 G4cout.precision(oldPrec);
00684 }
00685 }
00686
00687 if( fVerboseLevel > 2 )
00688 {
00689 G4cout << G4endl;
00690 G4cout << " G4PathFinder::ReLocate : entered " << G4endl;
00691 G4cout << " ---------------------- -------" << G4endl;
00692 G4cout << " *Re*Locating at position " << position << G4endl;
00693
00694
00695 if ( (fVerboseLevel > -1) || ( moveLenEndPosSq > 0.0) )
00696 {
00697 G4cout << " lastEndPosition = " << lastEndPosition
00698 << " moveVec from step-end = " << moveVecEndPos
00699 << " is new Track = " << fNewTrack
00700 << " relocated = " << fRelocatedPoint << G4endl;
00701 }
00702 }
00703 #endif // G4DEBUG_PATHFINDER
00704
00705 for ( register G4int num=0; num< fNoActiveNavigators ; ++pNavIter,++num )
00706 {
00707
00708
00709 (*pNavIter)->LocateGlobalPointWithinVolume( position );
00710
00711
00712
00713 fLimitedStep[num] = kDoNot;
00714 fCurrentStepSize[num] = 0.0;
00715 fLimitTruth[num] = false;
00716 }
00717
00718 fLastLocatedPosition= position;
00719 fRelocatedPoint= false;
00720
00721 #ifdef G4DEBUG_PATHFINDER
00722 if( fVerboseLevel > 2 )
00723 {
00724 G4cout << " G4PathFinder::ReLocate : exiting "
00725 << " at position " << fLastLocatedPosition << G4endl << G4endl;
00726 }
00727 #endif
00728 }
00729
00730
00731
00732 G4double G4PathFinder::ComputeSafety( const G4ThreeVector& position )
00733 {
00734
00735
00736 G4double minSafety= kInfinity;
00737
00738 std::vector<G4Navigator*>::iterator pNavigatorIter;
00739 pNavigatorIter= fpTransportManager->GetActiveNavigatorsIterator();
00740
00741 for( register G4int num=0; num<fNoActiveNavigators; ++pNavigatorIter,++num )
00742 {
00743 G4double safety = (*pNavigatorIter)->ComputeSafety( position,true );
00744 if( safety < minSafety ) { minSafety = safety; }
00745 fNewSafetyComputed[num]= safety;
00746 }
00747
00748 fSafetyLocation= position;
00749 fMinSafety_atSafLocation = minSafety;
00750
00751 #ifdef G4DEBUG_PATHFINDER
00752 if( fVerboseLevel > 1 )
00753 {
00754 G4cout << " G4PathFinder::ComputeSafety - returns "
00755 << minSafety << " at location " << position << G4endl;
00756 }
00757 #endif
00758 return minSafety;
00759 }
00760
00761
00762
00763
00764 G4TouchableHandle
00765 G4PathFinder::CreateTouchableHandle( G4int navId ) const
00766 {
00767 #ifdef G4DEBUG_PATHFINDER
00768 if( fVerboseLevel > 2 )
00769 {
00770 G4cout << "G4PathFinder::CreateTouchableHandle : navId = "
00771 << navId << " -- " << GetNavigator(navId) << G4endl;
00772 }
00773 #endif
00774
00775 G4TouchableHistory* touchHist;
00776 touchHist= GetNavigator(navId) -> CreateTouchableHistory();
00777
00778 G4VPhysicalVolume* locatedVolume= fLocatedVolume[navId];
00779 if( locatedVolume == 0 )
00780 {
00781
00782
00783 touchHist->UpdateYourself( locatedVolume, touchHist->GetHistory() );
00784 }
00785
00786 #ifdef G4DEBUG_PATHFINDER
00787 if( fVerboseLevel > 2 )
00788 {
00789 G4String VolumeName("None");
00790 if( locatedVolume ) { VolumeName= locatedVolume->GetName(); }
00791 G4cout << " Touchable History created at address " << touchHist
00792 << " volume = " << locatedVolume << " name= " << VolumeName
00793 << G4endl;
00794 }
00795 #endif
00796
00797 return G4TouchableHandle(touchHist);
00798 }
00799
00800 G4double
00801 G4PathFinder::DoNextLinearStep( const G4FieldTrack &initialState,
00802 G4double proposedStepLength )
00803 {
00804 std::vector<G4Navigator*>::iterator pNavigatorIter;
00805 G4double safety= 0.0, step=0.0;
00806 G4double minSafety= kInfinity, minStep;
00807
00808 const G4int IdTransport= 0;
00809 register G4int num=0;
00810
00811 #ifdef G4DEBUG_PATHFINDER
00812 if( fVerboseLevel > 2 )
00813 {
00814 G4cout << " G4PathFinder::DoNextLinearStep : entered " << G4endl;
00815 G4cout << " Input field track= " << initialState << G4endl;
00816 G4cout << " Requested step= " << proposedStepLength << G4endl;
00817 }
00818 #endif
00819
00820 G4ThreeVector initialPosition= initialState.GetPosition();
00821 G4ThreeVector initialDirection= initialState.GetMomentumDirection();
00822
00823 G4ThreeVector OriginShift = initialPosition - fPreSafetyLocation;
00824 G4double MagSqShift = OriginShift.mag2() ;
00825 G4double MagShift;
00826
00827
00828
00829
00830
00831
00832
00833 MagShift= std::sqrt(MagSqShift) ;
00834
00835 #ifdef G4PATHFINDER_OPTIMISATION
00836
00837 G4double fullSafety;
00838
00839 if( MagSqShift >= sqr(fPreSafetyMinValue ) )
00840 {
00841 fullSafety = 0.0 ;
00842 }
00843 else
00844 {
00845 fullSafety = fPreSafetyMinValue - MagShift;
00846 }
00847 if( proposedStepLength < fullSafety )
00848 {
00849
00850
00851
00852 fPreStepCenterRenewed= false;
00853
00854 for( num=0; num< fNoActiveNavigators; ++num )
00855 {
00856 fCurrentStepSize[num]= kInfinity;
00857 safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
00858 minSafety= std::min( safety, minSafety );
00859 fCurrentPreStepSafety[num]= safety;
00860 }
00861 minStep= kInfinity;
00862
00863 #ifdef G4DEBUG_PATHFINDER
00864 if( fVerboseLevel > 2 )
00865 {
00866 G4cout << "G4PathFinder::DoNextLinearStep : Quick Stepping. " << G4endl
00867 << " proposedStepLength " << proposedStepLength
00868 << " < (full) safety = " << fullSafety
00869 << " at " << initialPosition
00870 << G4endl;
00871 }
00872 #endif
00873 }
00874 else
00875 #endif // End of G4PATHFINDER_OPTIMISATION 1
00876 {
00877
00878
00879
00880 fPreStepCenterRenewed= true;
00881 pNavigatorIter= fpTransportManager-> GetActiveNavigatorsIterator();
00882
00883 minStep= kInfinity;
00884
00885 for( num=0; num< fNoActiveNavigators; ++pNavigatorIter,++num )
00886 {
00887 safety = std::max( 0.0, fPreSafetyValues[num] - MagShift);
00888
00889 #ifdef G4PATHFINDER_OPTIMISATION
00890 if( proposedStepLength <= safety )
00891 {
00892
00893
00894 step= kInfinity;
00895
00896 #ifdef G4DEBUG_PATHFINDER
00897 G4cout.precision(8);
00898 G4cout << "PathFinder::ComputeStep> small proposed step = "
00899 << proposedStepLength
00900 << " <= safety = " << safety << " for nav " << num
00901 << " Step fully taken. " << G4endl;
00902 #endif
00903 }
00904 else
00905 #endif // End of G4PATHFINDER_OPTIMISATION 2
00906 {
00907 #ifdef G4DEBUG_PATHFINDER
00908 G4double previousSafety= safety;
00909 #endif
00910 step= (*pNavigatorIter)->ComputeStep( initialPosition,
00911 initialDirection,
00912 proposedStepLength,
00913 safety );
00914 minStep = std::min( step, minStep);
00915
00916
00917
00918
00919 #ifdef G4DEBUG_PATHFINDER
00920 if( fVerboseLevel > 0)
00921 {
00922 G4cout.precision(8);
00923 G4cout << "PathFinder::ComputeStep> long proposed step = "
00924 << proposedStepLength
00925 << " > safety = " << previousSafety
00926 << " for nav " << num
00927 << " . New safety = " << safety << " step= " << step
00928 << G4endl;
00929 }
00930 #endif
00931 }
00932 fCurrentStepSize[num] = step;
00933
00934
00935
00936
00937
00938 fPreSafetyValues[num]= safety;
00939 fCurrentPreStepSafety[num]= safety;
00940
00941 minSafety= std::min( safety, minSafety );
00942
00943 #ifdef G4DEBUG_PATHFINDER
00944 if( fVerboseLevel > 2 )
00945 {
00946 G4cout << "G4PathFinder::DoNextLinearStep : Navigator ["
00947 << num << "] -- step size " << step << G4endl;
00948 }
00949 #endif
00950 }
00951
00952
00953
00954
00955 fPreSafetyLocation= initialPosition;
00956 fPreSafetyMinValue= minSafety;
00957 }
00958
00959
00960
00961 fPreStepLocation= initialPosition;
00962 fMinSafety_PreStepPt= minSafety;
00963
00964 fMinStep= minStep;
00965
00966 if( fMinStep == kInfinity )
00967 {
00968 minStep = proposedStepLength;
00969 }
00970 fTrueMinStep = minStep;
00971
00972
00973
00974 G4ThreeVector endPosition;
00975
00976 fEndState= initialState;
00977 endPosition= initialPosition + minStep * initialDirection ;
00978
00979 #ifdef G4DEBUG_PATHFINDER
00980 if( fVerboseLevel > 1 )
00981 {
00982 G4cout << "G4PathFinder::DoNextLinearStep : "
00983 << " initialPosition = " << initialPosition
00984 << " and endPosition = " << endPosition<< G4endl;
00985 }
00986 #endif
00987
00988 fEndState.SetPosition( endPosition );
00989 fEndState.SetProperTimeOfFlight( -1.000 );
00990
00991 if( fNoActiveNavigators == 1 )
00992 {
00993 G4bool transportLimited = (fMinStep!= kInfinity);
00994 fLimitTruth[IdTransport] = transportLimited;
00995 fLimitedStep[IdTransport] = transportLimited ? kUnique : kDoNot;
00996
00997
00998 fNoGeometriesLimiting = transportLimited ? 1 : 0;
00999 }
01000 else
01001 {
01002 WhichLimited();
01003 }
01004
01005 #ifdef G4DEBUG_PATHFINDER
01006 if( fVerboseLevel > 2 )
01007 {
01008 G4cout << " G4PathFinder::DoNextLinearStep : exits returning "
01009 << minStep << G4endl;
01010 G4cout << " Endpoint values = " << fEndState << G4endl;
01011 G4cout << G4endl;
01012 }
01013 #endif
01014
01015 return minStep;
01016 }
01017
01018 void G4PathFinder::WhichLimited()
01019 {
01020
01021
01022 G4int num=-1, last=-1;
01023 G4int noLimited=0;
01024 ELimited shared= kSharedOther;
01025
01026 const G4int IdTransport= 0;
01027
01028
01029
01030 G4bool transportLimited = (fCurrentStepSize[IdTransport] == fMinStep)
01031 && ( fMinStep!= kInfinity) ;
01032
01033 if( transportLimited ) {
01034 shared= kSharedTransport;
01035 }
01036
01037 for ( num= 0; num < fNoActiveNavigators; num++ )
01038 {
01039 G4bool limitedStep;
01040
01041 G4double step= fCurrentStepSize[num];
01042
01043 limitedStep = ( std::fabs(step - fMinStep) < kCarTolerance )
01044 && ( step != kInfinity);
01045
01046 fLimitTruth[ num ] = limitedStep;
01047 if( limitedStep )
01048 {
01049 noLimited++;
01050 fLimitedStep[num] = shared;
01051 last= num;
01052 }
01053 else
01054 {
01055 fLimitedStep[num] = kDoNot;
01056 }
01057 }
01058 fNoGeometriesLimiting= noLimited;
01059
01060 if( (last > -1) && (noLimited == 1 ) )
01061 {
01062 fLimitedStep[ last ] = kUnique;
01063 }
01064
01065 #ifdef G4DEBUG_PATHFINDER
01066 if( fVerboseLevel > 1 )
01067 {
01068 PrintLimited();
01069 if( fVerboseLevel > 4 ) {
01070 G4cout << " G4PathFinder::WhichLimited - exiting. " << G4endl;
01071 }
01072 }
01073 #endif
01074 }
01075
01076 void G4PathFinder::PrintLimited()
01077 {
01078
01079
01080 G4cout << "G4PathFinder::PrintLimited reports: " ;
01081 G4cout << " Minimum step (true)= " << fTrueMinStep
01082 << " reported min = " << fMinStep
01083 << G4endl;
01084 if( (fCurrentStepNo <= 2) || (fVerboseLevel>=2) )
01085 {
01086 G4cout << std::setw(5) << " Step#" << " "
01087 << std::setw(5) << " NavId" << " "
01088 << std::setw(12) << " step-size " << " "
01089 << std::setw(12) << " raw-size " << " "
01090 << std::setw(12) << " pre-safety " << " "
01091 << std::setw(15) << " Limited / flag" << " "
01092 << std::setw(15) << " World " << " "
01093 << G4endl;
01094 }
01095 G4int num;
01096 for ( num= 0; num < fNoActiveNavigators; num++ )
01097 {
01098 G4double rawStep = fCurrentStepSize[num];
01099 G4double stepLen = fCurrentStepSize[num];
01100 if( stepLen > fTrueMinStep )
01101 {
01102 stepLen = fTrueMinStep;
01103 }
01104 G4int oldPrec= G4cout.precision(9);
01105
01106 G4cout << std::setw(5) << fCurrentStepNo << " "
01107 << std::setw(5) << num << " "
01108 << std::setw(12) << stepLen << " "
01109 << std::setw(12) << rawStep << " "
01110 << std::setw(12) << fCurrentPreStepSafety[num] << " "
01111 << std::setw(5) << (fLimitTruth[num] ? "YES" : " NO") << " ";
01112 G4String limitedStr= LimitedString(fLimitedStep[num]);
01113 G4cout << " " << std::setw(15) << limitedStr << " ";
01114 G4cout.precision(oldPrec);
01115
01116 G4Navigator *pNav= GetNavigator( num );
01117 G4String WorldName( "Not-Set" );
01118 if (pNav)
01119 {
01120 G4VPhysicalVolume *pWorld= pNav->GetWorldVolume();
01121 if( pWorld )
01122 {
01123 WorldName = pWorld->GetName();
01124 }
01125 }
01126 G4cout << " " << WorldName ;
01127 G4cout << G4endl;
01128 }
01129
01130 if( fVerboseLevel > 4 )
01131 {
01132 G4cout << " G4PathFinder::PrintLimited - exiting. " << G4endl;
01133 }
01134 }
01135
01136 G4double
01137 G4PathFinder::DoNextCurvedStep( const G4FieldTrack &initialState,
01138 G4double proposedStepLength,
01139 G4VPhysicalVolume* pCurrentPhysicalVolume )
01140 {
01141 const G4double toleratedRelativeError= 1.0e-10;
01142 G4double minStep= kInfinity, newSafety=0.0;
01143 G4int numNav;
01144 G4FieldTrack fieldTrack= initialState;
01145 G4ThreeVector startPoint= initialState.GetPosition();
01146
01147 #ifdef G4DEBUG_PATHFINDER
01148 G4int prc= G4cout.precision(9);
01149 if( fVerboseLevel > 2 )
01150 {
01151 G4cout << " G4PathFinder::DoNextCurvedStep ****** " << G4endl;
01152 G4cout << " Initial value of field track is " << fieldTrack
01153 << " and proposed step= " << proposedStepLength << G4endl;
01154 }
01155 #endif
01156
01157 fPreStepCenterRenewed= true;
01158
01159 if( fNoActiveNavigators > 1 )
01160 {
01161
01162
01163 G4double minSafety= kInfinity, safety;
01164 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
01165 {
01166 safety= fpNavigator[numNav]->ComputeSafety( startPoint, false );
01167 fPreSafetyValues[numNav]= safety;
01168 fCurrentPreStepSafety[numNav]= safety;
01169 minSafety = std::min( safety, minSafety );
01170 }
01171
01172
01173
01174 fPreSafetyLocation= startPoint;
01175 fPreSafetyMinValue= minSafety;
01176 fPreStepLocation= startPoint;
01177 fMinSafety_PreStepPt= minSafety;
01178 }
01179
01180
01181
01182 minStep= fpFieldPropagator->ComputeStep( fieldTrack,
01183 proposedStepLength,
01184 newSafety,
01185 pCurrentPhysicalVolume );
01186
01187
01188
01189 fEndState= fieldTrack;
01190 fMinStep= minStep;
01191 fTrueMinStep = std::min( minStep, proposedStepLength );
01192
01193 if( fNoActiveNavigators== 1 )
01194 {
01195
01196
01197
01198 fPreSafetyValues[0]= newSafety;
01199 fPreSafetyLocation= startPoint;
01200 fPreSafetyMinValue= newSafety;
01201
01202
01203
01204 fCurrentPreStepSafety[0]= newSafety;
01205 fPreStepLocation= startPoint;
01206 fMinSafety_PreStepPt= newSafety;
01207 }
01208
01209 #ifdef G4DEBUG_PATHFINDER
01210 if( fVerboseLevel > 2 )
01211 {
01212 G4cout << "G4PathFinder::DoNextCurvedStep : " << G4endl
01213 << " initialState = " << initialState << G4endl
01214 << " and endState = " << fEndState << G4endl;
01215 G4cout << "G4PathFinder::DoNextCurvedStep : "
01216 << " minStep = " << minStep
01217 << " proposedStepLength " << proposedStepLength
01218 << " safety = " << newSafety << G4endl;
01219 }
01220 #endif
01221 G4double currentStepSize;
01222 if( minStep < proposedStepLength )
01223 {
01224
01225
01226
01227 G4int noLimited= 0;
01228 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
01229 {
01230 G4double finalStep, lastPreSafety=0.0, minStepLast;
01231 ELimited didLimit;
01232 G4bool limited;
01233
01234 finalStep= fpMultiNavigator->ObtainFinalStep( numNav, lastPreSafety,
01235 minStepLast, didLimit );
01236
01237
01238
01239
01240 currentStepSize = fTrueMinStep;
01241 G4double diffStep= 0.0;
01242 if( (minStepLast != kInfinity) )
01243 {
01244 diffStep = (finalStep-minStepLast);
01245 if ( std::abs(diffStep) <= toleratedRelativeError * finalStep )
01246 {
01247 diffStep = 0.0;
01248 }
01249 currentStepSize += diffStep;
01250 }
01251 fCurrentStepSize[numNav] = currentStepSize;
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264 fLimitedStep[numNav] = didLimit;
01265 fLimitTruth[numNav] = limited = (didLimit != kDoNot );
01266 if( limited ) { noLimited++; }
01267
01268 #ifdef G4DEBUG_PATHFINDER
01269 G4bool StepError= (currentStepSize < 0)
01270 || ( (minStepLast != kInfinity) && (diffStep < 0) ) ;
01271 if( StepError || (fVerboseLevel > 2) )
01272 {
01273 G4String limitedString= LimitedString( fLimitedStep[numNav] );
01274
01275 G4cout << " G4PathFinder::ComputeStep. Geometry " << numNav
01276 << " step= " << fCurrentStepSize[numNav]
01277 << " from final-step= " << finalStep
01278 << " fTrueMinStep= " << fTrueMinStep
01279 << " minStepLast= " << minStepLast
01280 << " limited = " << (fLimitTruth[numNav] ? "YES" : " NO")
01281 << " ";
01282 G4cout << " status = " << limitedString << " #= " << didLimit
01283 << G4endl;
01284
01285 if( StepError )
01286 {
01287 std::ostringstream message;
01288 message << "Incorrect calculation of step size for one navigator"
01289 << G4endl
01290 << " currentStepSize = " << currentStepSize
01291 << ", diffStep= " << diffStep << G4endl
01292 << "ERROR in computing step size for this navigator.";
01293 G4Exception("G4PathFinder::DoNextCurvedStep",
01294 "GeomNav0003", FatalException, message);
01295 }
01296 }
01297 #endif
01298 }
01299
01300 fNoGeometriesLimiting= noLimited;
01301 }
01302 else if ( (minStep == proposedStepLength)
01303 || (minStep == kInfinity)
01304 || ( std::abs(minStep-proposedStepLength)
01305 < toleratedRelativeError * proposedStepLength ) )
01306 {
01307
01308
01309
01310
01311
01312
01313
01314 currentStepSize= minStep;
01315 for( numNav=0; numNav < fNoActiveNavigators; ++numNav )
01316 {
01317 fCurrentStepSize[numNav] = minStep;
01318
01319 fLimitedStep[numNav] = kDoNot;
01320 fLimitTruth[numNav] = false;
01321 }
01322 fNoGeometriesLimiting= 0;
01323 }
01324 else
01325 {
01326 std::ostringstream message;
01327 message << "Incorrect calculation of step size for one navigator." << G4endl
01328 << " currentStepSize = " << minStep << " is larger than "
01329 << " proposed StepSize = " << proposedStepLength << ".";
01330 G4Exception("G4PathFinder::DoNextCurvedStep()",
01331 "GeomNav0003", FatalException, message);
01332 }
01333
01334 #ifdef G4DEBUG_PATHFINDER
01335 if( fVerboseLevel > 2 )
01336 {
01337 G4cout << " Exiting G4PathFinder::DoNextCurvedStep " << G4endl;
01338 PrintLimited();
01339 }
01340 G4cout.precision(prc);
01341 #endif
01342
01343 return minStep;
01344 }
01345
01346 G4String& G4PathFinder::LimitedString( ELimited lim )
01347 {
01348 static G4String StrDoNot("DoNot"),
01349 StrUnique("Unique"),
01350 StrUndefined("Undefined"),
01351 StrSharedTransport("SharedTransport"),
01352 StrSharedOther("SharedOther");
01353
01354 G4String* limitedStr;
01355 switch ( lim )
01356 {
01357 case kDoNot: limitedStr= &StrDoNot; break;
01358 case kUnique: limitedStr = &StrUnique; break;
01359 case kSharedTransport: limitedStr= &StrSharedTransport; break;
01360 case kSharedOther: limitedStr = &StrSharedOther; break;
01361 default: limitedStr = &StrUndefined; break;
01362 }
01363 return *limitedStr;
01364 }
01365
01366 void G4PathFinder::PushPostSafetyToPreSafety()
01367 {
01368 fPreSafetyLocation= fSafetyLocation;
01369 fPreSafetyMinValue= fMinSafety_atSafLocation;
01370 for( register G4int nav=0; nav < fNoActiveNavigators; ++nav )
01371 {
01372 fPreSafetyValues[nav]= fNewSafetyComputed[nav];
01373 }
01374 }