Geant4 Cross Reference |
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 // class G4ITNavigator1 Implementation 28 // 29 // Original author: Paul Kent, July 95/96 30 // 31 // G4ITNavigator1 is a duplicate version of G4Navigator starting from Geant4.9.5 32 // initially written by Paul Kent and colleagues. 33 // The only difference resides in the way the information is saved and managed 34 // 35 // History: 36 // - Created. Paul Kent, Jul 95/96 37 // - Zero step protections J.A. / G.C., Nov 2004 38 // - Added check mode G. Cosmo, Mar 2004 39 // - Made Navigator Abstract G. Cosmo, Nov 2003 40 // - G4ITNavigator1 created M.K., Nov 2012 41 // -------------------------------------------------------------------- 42 43 #include <iomanip> 44 45 #include "G4ITNavigator1.hh" 46 //#include "G4ITNavigator.hh" 47 #include "G4ios.hh" 48 #include "G4SystemOfUnits.hh" 49 #include "G4GeometryTolerance.hh" 50 #include "G4VPhysicalVolume.hh" 51 52 #define G4DEBUG_NAVIGATION 1 53 #include "G4VoxelSafety.hh" 54 55 // ******************************************************************** 56 // Constructor 57 // ******************************************************************** 58 // 59 G4ITNavigator1::G4ITNavigator1() 60 { 61 fActive= false; 62 fLastTriedStepComputation= false; 63 64 ResetStackAndState(); 65 // Initialises also all 66 // - exit / entry flags 67 // - flags & variables for exit normals 68 // - zero step counters 69 // - blocked volume 70 71 fActionThreshold_NoZeroSteps = 10; 72 fAbandonThreshold_NoZeroSteps = 25; 73 74 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 75 fregularNav.SetNormalNavigation( &fnormalNav ); 76 77 fStepEndPoint = G4ThreeVector( kInfinity, kInfinity, kInfinity ); 78 fLastStepEndPointLocal = G4ThreeVector( kInfinity, kInfinity, kInfinity ); 79 80 fpVoxelSafety= new G4VoxelSafety(); 81 fpSaveState = nullptr; 82 83 // this->SetVerboseLevel(3); 84 // this->CheckMode(true); 85 } 86 87 // !> 88 89 G4ITNavigator1::G4SaveNavigatorState::G4SaveNavigatorState() 90 { 91 sWasLimitedByGeometry = false; 92 sEntering = false; 93 sExiting = false; 94 sLocatedOnEdge = false; 95 sLastStepWasZero = 0; 96 sEnteredDaughter = false; 97 sExitedMother = false; 98 sPushed = false; 99 100 sValidExitNormal = false; 101 sExitNormal = G4ThreeVector(0,0,0); 102 103 sPreviousSftOrigin = G4ThreeVector(0,0,0); 104 sPreviousSafety = 0.0; 105 106 sNumberZeroSteps = 0; 107 108 spBlockedPhysicalVolume = nullptr; 109 sBlockedReplicaNo = -1; 110 111 sLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 ); 112 sLocatedOutsideWorld = false; 113 } 114 115 // <! 116 117 // ******************************************************************** 118 // Destructor 119 // ******************************************************************** 120 // 121 G4ITNavigator1::~G4ITNavigator1() 122 { delete fpVoxelSafety; } 123 124 // ******************************************************************** 125 // ResetHierarchyAndLocate 126 // ******************************************************************** 127 // 128 G4VPhysicalVolume* 129 G4ITNavigator1::ResetHierarchyAndLocate(const G4ThreeVector &p, 130 const G4ThreeVector &direction, 131 const G4TouchableHistory &h) 132 { 133 ResetState(); 134 fHistory = *h.GetHistory(); 135 SetupHierarchy(); 136 fLastTriedStepComputation= false; // Redundant, but best 137 return LocateGlobalPointAndSetup(p, &direction, true, false); 138 } 139 140 // ******************************************************************** 141 // LocateGlobalPointAndSetup 142 // 143 // Locate the point in the hierarchy return 0 if outside 144 // The direction is required 145 // - if on an edge shared by more than two surfaces 146 // (to resolve likely looping in tracking) 147 // - at initial location of a particle 148 // (to resolve potential ambiguity at boundary) 149 // 150 // Flags on exit: (comments to be completed) 151 // fEntering - True if entering `daughter' volume (or replica) 152 // whether daughter of last mother directly 153 // or daughter of that volume's ancestor. 154 // ******************************************************************** 155 // 156 G4VPhysicalVolume* 157 G4ITNavigator1::LocateGlobalPointAndSetup( const G4ThreeVector& globalPoint, 158 const G4ThreeVector* pGlobalDirection, 159 const G4bool relativeSearch, 160 const G4bool ignoreDirection ) 161 { 162 G4bool notKnownContained=true, noResult; 163 G4VPhysicalVolume *targetPhysical; 164 G4LogicalVolume *targetLogical; 165 G4VSolid *targetSolid=nullptr; 166 G4ThreeVector localPoint, globalDirection; 167 EInside insideCode; 168 169 G4bool considerDirection = (!ignoreDirection) || fLocatedOnEdge; 170 171 fLastTriedStepComputation= false; 172 fChangedGrandMotherRefFrame= false; // For local exit normal 173 174 if( considerDirection && pGlobalDirection != nullptr ) 175 { 176 globalDirection=*pGlobalDirection; 177 } 178 179 180 #ifdef G4VERBOSE 181 if( fVerbose > 2 ) 182 { 183 G4long oldcoutPrec = G4cout.precision(8); 184 G4cout << "*** G4ITNavigator1::LocateGlobalPointAndSetup: ***" << G4endl; 185 G4cout << " Called with arguments: " << G4endl 186 << " Globalpoint = " << globalPoint << G4endl 187 << " RelativeSearch = " << relativeSearch << G4endl; 188 if( fVerbose == 4 ) 189 { 190 G4cout << " ----- Upon entering:" << G4endl; 191 PrintState(); 192 } 193 G4cout.precision(oldcoutPrec); 194 } 195 #endif 196 197 if ( !relativeSearch ) 198 { 199 ResetStackAndState(); 200 } 201 else 202 { 203 if ( fWasLimitedByGeometry ) 204 { 205 fWasLimitedByGeometry = false; 206 fEnteredDaughter = fEntering; // Remember 207 fExitedMother = fExiting; // Remember 208 if ( fExiting ) 209 { 210 if ( fHistory.GetDepth() != 0u ) 211 { 212 fBlockedPhysicalVolume = fHistory.GetTopVolume(); 213 fBlockedReplicaNo = fHistory.GetTopReplicaNo(); 214 fHistory.BackLevel(); 215 } 216 else 217 { 218 fLastLocatedPointLocal = localPoint; 219 fLocatedOutsideWorld = true; 220 return nullptr; // Have exited world volume 221 } 222 // A fix for the case where a volume is "entered" at an edge 223 // and a coincident surface exists outside it. 224 // - This stops it from exiting further volumes and cycling 225 // - However ReplicaNavigator treats this case itself 226 // 227 if ( fLocatedOnEdge && (VolumeType(fBlockedPhysicalVolume)!=kReplica )) 228 { 229 fExiting= false; 230 } 231 } 232 else 233 if ( fEntering ) 234 { 235 switch (VolumeType(fBlockedPhysicalVolume)) 236 { 237 case kNormal: 238 fHistory.NewLevel(fBlockedPhysicalVolume, kNormal, 239 fBlockedPhysicalVolume->GetCopyNo()); 240 break; 241 case kReplica: 242 freplicaNav.ComputeTransformation(fBlockedReplicaNo, 243 fBlockedPhysicalVolume); 244 fHistory.NewLevel(fBlockedPhysicalVolume, kReplica, 245 fBlockedReplicaNo); 246 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo); 247 break; 248 case kParameterised: 249 if( fBlockedPhysicalVolume->GetRegularStructureId() == 0 ) 250 { 251 G4VSolid *pSolid; 252 G4VPVParameterisation *pParam; 253 G4TouchableHistory parentTouchable( fHistory ); 254 pParam = fBlockedPhysicalVolume->GetParameterisation(); 255 pSolid = pParam->ComputeSolid(fBlockedReplicaNo, 256 fBlockedPhysicalVolume); 257 pSolid->ComputeDimensions(pParam, fBlockedReplicaNo, 258 fBlockedPhysicalVolume); 259 pParam->ComputeTransformation(fBlockedReplicaNo, 260 fBlockedPhysicalVolume); 261 fHistory.NewLevel(fBlockedPhysicalVolume, kParameterised, 262 fBlockedReplicaNo); 263 fBlockedPhysicalVolume->SetCopyNo(fBlockedReplicaNo); 264 // 265 // Set the correct solid and material in Logical Volume 266 // 267 G4LogicalVolume *pLogical; 268 pLogical = fBlockedPhysicalVolume->GetLogicalVolume(); 269 pLogical->SetSolid( pSolid ); 270 pLogical->UpdateMaterial(pParam -> 271 ComputeMaterial(fBlockedReplicaNo, 272 fBlockedPhysicalVolume, 273 &parentTouchable)); 274 } 275 break; 276 case kExternal: 277 G4Exception("G4ITNavigator1::LocateGlobalPointAndSetup()", 278 "GeomNav0001", FatalException, 279 "Not applicable for external volumes."); 280 break; 281 } 282 fEntering = false; 283 fBlockedPhysicalVolume = nullptr; 284 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint); 285 notKnownContained = false; 286 } 287 } 288 else 289 { 290 fBlockedPhysicalVolume = nullptr; 291 fEntering = false; 292 fEnteredDaughter = false; // Full Step was not taken, did not enter 293 fExiting = false; 294 fExitedMother = false; // Full Step was not taken, did not exit 295 } 296 } 297 // 298 // Search from top of history up through geometry until 299 // containing volume found: 300 // If on 301 // o OUTSIDE - Back up level, not/no longer exiting volumes 302 // o SURFACE and EXITING - Back up level, setting new blocking no.s 303 // else 304 // o containing volume found 305 // 306 G4int noLevelsExited=0 ; 307 308 while (notKnownContained) 309 { 310 if ( fHistory.GetTopVolumeType()!=kReplica ) 311 { 312 targetSolid = fHistory.GetTopVolume()->GetLogicalVolume()->GetSolid(); 313 localPoint = fHistory.GetTopTransform().TransformPoint(globalPoint); 314 insideCode = targetSolid->Inside(localPoint); 315 #ifdef G4VERBOSE 316 if(( fVerbose == 1 ) && ( fCheck )) 317 { 318 G4String solidResponse = "-kInside-"; 319 if (insideCode == kOutside) 320 solidResponse = "-kOutside-"; 321 else if (insideCode == kSurface) 322 solidResponse = "-kSurface-"; 323 G4cout << "*** G4ITNavigator1::LocateGlobalPointAndSetup(): ***" << G4endl 324 << " Invoked Inside() for solid: " << targetSolid->GetName() 325 << ". Solid replied: " << solidResponse << G4endl 326 << " For local point p: " << localPoint << G4endl; 327 } 328 #endif 329 } 330 else 331 { 332 insideCode = freplicaNav.BackLocate(fHistory, globalPoint, localPoint, 333 fExiting, notKnownContained); 334 // !CARE! if notKnownContained returns false then the point is within 335 // the containing placement volume of the replica(s). If insidecode 336 // will result in the history being backed up one level, then the 337 // local point returned is the point in the system of this new level 338 } 339 340 341 if ( insideCode==kOutside ) 342 { 343 noLevelsExited++; 344 if ( fHistory.GetDepth() != 0u ) 345 { 346 fBlockedPhysicalVolume = fHistory.GetTopVolume(); 347 fBlockedReplicaNo = fHistory.GetTopReplicaNo(); 348 fHistory.BackLevel(); 349 fExiting = false; 350 351 if( noLevelsExited > 1 ) 352 { 353 // The first transformation was done by the sub-navigator 354 // 355 const G4RotationMatrix* mRot = fBlockedPhysicalVolume->GetRotation(); 356 if( mRot != nullptr ) 357 { 358 fGrandMotherExitNormal *= (*mRot).inverse(); 359 fChangedGrandMotherRefFrame= true; 360 } 361 } 362 } 363 else 364 { 365 fLastLocatedPointLocal = localPoint; 366 fLocatedOutsideWorld = true; 367 // No extra transformation for ExitNormal - is in frame of Top Volume 368 return nullptr; // Have exited world volume 369 } 370 } 371 else 372 if ( insideCode==kSurface ) 373 { 374 G4bool isExiting = fExiting; 375 if( (!fExiting)&&considerDirection ) 376 { 377 // Figure out whether we are exiting this level's volume 378 // by using the direction 379 // 380 G4bool directionExiting = false; 381 G4ThreeVector localDirection = 382 fHistory.GetTopTransform().TransformAxis(globalDirection); 383 384 // Make sure localPoint in correct reference frame 385 // ( Was it already correct ? How ? ) 386 // 387 localPoint= fHistory.GetTopTransform().TransformPoint(globalPoint); 388 if ( fHistory.GetTopVolumeType()!=kReplica ) 389 { 390 G4ThreeVector normal = targetSolid->SurfaceNormal(localPoint); 391 directionExiting = normal.dot(localDirection) > 0.0; 392 isExiting = isExiting || directionExiting; 393 } 394 } 395 if( isExiting ) 396 { 397 noLevelsExited++; 398 if ( fHistory.GetDepth() != 0u ) 399 { 400 fBlockedPhysicalVolume = fHistory.GetTopVolume(); 401 fBlockedReplicaNo = fHistory.GetTopReplicaNo(); 402 fHistory.BackLevel(); 403 // 404 // Still on surface but exited volume not necessarily convex 405 // 406 fValidExitNormal = false; 407 408 if( noLevelsExited > 1 ) 409 { 410 // The first transformation was done by the sub-navigator 411 // 412 const G4RotationMatrix* mRot=fBlockedPhysicalVolume->GetRotation(); 413 if( mRot != nullptr ) 414 { 415 fGrandMotherExitNormal *= (*mRot).inverse(); 416 fChangedGrandMotherRefFrame= true; 417 } 418 } 419 } 420 else 421 { 422 fLastLocatedPointLocal = localPoint; 423 fLocatedOutsideWorld = true; 424 // No extra transformation for ExitNormal, is in frame of Top Vol 425 return nullptr; // Have exited world volume 426 } 427 } 428 else 429 { 430 notKnownContained=false; 431 } 432 } 433 else 434 { 435 notKnownContained=false; 436 } 437 } // END while (notKnownContained) 438 // 439 // Search downwards until deepest containing volume found, 440 // blocking fBlockedPhysicalVolume/BlockedReplicaNum 441 // 442 // 3 Cases: 443 // 444 // o Parameterised daughters 445 // =>Must be one G4PVParameterised daughter & voxels 446 // o Positioned daughters & voxels 447 // o Positioned daughters & no voxels 448 449 noResult = true; // noResult should be renamed to 450 // something like enteredLevel, as that is its meaning. 451 do 452 { 453 // Determine `type' of current mother volume 454 // 455 targetPhysical = fHistory.GetTopVolume(); 456 if (targetPhysical == nullptr) { break; } 457 targetLogical = targetPhysical->GetLogicalVolume(); 458 switch( CharacteriseDaughters(targetLogical) ) 459 { 460 case kNormal: 461 if ( targetLogical->GetVoxelHeader() != nullptr ) // use optimised navigation 462 { 463 noResult = fvoxelNav.LevelLocate(fHistory, 464 fBlockedPhysicalVolume, 465 fBlockedReplicaNo, 466 globalPoint, 467 pGlobalDirection, 468 considerDirection, 469 localPoint); 470 } 471 else // do not use optimised navigation 472 { 473 noResult = fnormalNav.LevelLocate(fHistory, 474 fBlockedPhysicalVolume, 475 fBlockedReplicaNo, 476 globalPoint, 477 pGlobalDirection, 478 considerDirection, 479 localPoint); 480 } 481 break; 482 case kReplica: 483 noResult = freplicaNav.LevelLocate(fHistory, 484 fBlockedPhysicalVolume, 485 fBlockedReplicaNo, 486 globalPoint, 487 pGlobalDirection, 488 considerDirection, 489 localPoint); 490 break; 491 case kParameterised: 492 if( GetDaughtersRegularStructureId(targetLogical) != 1 ) 493 { 494 noResult = fparamNav.LevelLocate(fHistory, 495 fBlockedPhysicalVolume, 496 fBlockedReplicaNo, 497 globalPoint, 498 pGlobalDirection, 499 considerDirection, 500 localPoint); 501 } 502 else // Regular structure 503 { 504 noResult = fregularNav.LevelLocate(fHistory, 505 fBlockedPhysicalVolume, 506 fBlockedReplicaNo, 507 globalPoint, 508 pGlobalDirection, 509 considerDirection, 510 localPoint); 511 } 512 break; 513 case kExternal: 514 G4Exception("G4ITNavigator1::LocateGlobalPointAndSetup()", 515 "GeomNav0001", FatalException, 516 "Not applicable for external volumes."); 517 break; 518 } 519 520 // LevelLocate returns true if it finds a daughter volume 521 // in which globalPoint is inside (or on the surface). 522 523 if ( noResult ) 524 { 525 // Entering a daughter after ascending 526 // 527 // The blocked volume is no longer valid - it was for another level 528 // 529 fBlockedPhysicalVolume = nullptr; 530 fBlockedReplicaNo = -1; 531 532 // fEntering should be false -- else blockedVolume is assumed good. 533 // fEnteredDaughter is used for ExitNormal 534 // 535 fEntering = false; 536 fEnteredDaughter = true; 537 538 if( fExitedMother ) 539 { 540 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume(); 541 const G4RotationMatrix* mRot = enteredPhysical->GetRotation(); 542 if( mRot != nullptr ) 543 { 544 fGrandMotherExitNormal *= (*mRot).inverse(); 545 } 546 } 547 548 #ifdef G4DEBUG_NAVIGATION 549 if( fVerbose > 2 ) 550 { 551 G4VPhysicalVolume* enteredPhysical = fHistory.GetTopVolume(); 552 G4cout << "*** G4ITNavigator1::LocateGlobalPointAndSetup() ***" << G4endl; 553 G4cout << " Entering volume: " << enteredPhysical->GetName() 554 << G4endl; 555 } 556 #endif 557 } 558 } while (noResult); 559 560 fLastLocatedPointLocal = localPoint; 561 562 #ifdef G4VERBOSE 563 if( fVerbose >= 4 ) 564 { 565 G4long oldcoutPrec = G4cout.precision(8); 566 G4String curPhysVol_Name("None"); 567 if (targetPhysical != nullptr) { curPhysVol_Name = targetPhysical->GetName(); } 568 G4cout << " Return value = new volume = " << curPhysVol_Name << G4endl; 569 G4cout << " ----- Upon exiting:" << G4endl; 570 PrintState(); 571 if( fVerbose == 5 ) 572 { 573 G4cout << "Upon exiting LocateGlobalPointAndSetup():" << G4endl; 574 G4cout << " History = " << G4endl << fHistory << G4endl << G4endl; 575 } 576 G4cout.precision(oldcoutPrec); 577 } 578 #endif 579 580 fLocatedOutsideWorld= false; 581 582 return targetPhysical; 583 } 584 585 // ******************************************************************** 586 // LocateGlobalPointWithinVolume 587 // 588 // -> the state information of this Navigator and its subNavigators 589 // is updated in order to start the next step at pGlobalpoint 590 // -> no check is performed whether pGlobalpoint is inside the 591 // original volume (this must be the case). 592 // 593 // Note: a direction could be added to the arguments, to aid in future 594 // optional checking (via the old code below, flagged by OLD_LOCATE). 595 // [ This would be done only in verbose mode ] 596 // ******************************************************************** 597 // 598 void 599 G4ITNavigator1::LocateGlobalPointWithinVolume(const G4ThreeVector& pGlobalpoint) 600 { 601 fLastLocatedPointLocal = ComputeLocalPoint(pGlobalpoint); 602 fLastTriedStepComputation= false; 603 fChangedGrandMotherRefFrame= false; // Frame for Exit Normal 604 605 #ifdef G4DEBUG_NAVIGATION 606 if( fVerbose > 2 ) 607 { 608 G4cout << "Entering LocateGlobalWithinVolume(): History = " << G4endl; 609 G4cout << fHistory << G4endl; 610 } 611 #endif 612 613 // For the case of Voxel (or Parameterised) volume the respective 614 // Navigator must be messaged to update its voxel information etc 615 616 // Update the state of the Sub Navigators 617 // - in particular any voxel information they store/cache 618 // 619 G4VPhysicalVolume* motherPhysical = fHistory.GetTopVolume(); 620 G4LogicalVolume* motherLogical = motherPhysical->GetLogicalVolume(); 621 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader(); 622 623 if ( fHistory.GetTopVolumeType()!=kReplica ) 624 { 625 switch( CharacteriseDaughters(motherLogical) ) 626 { 627 case kNormal: 628 if ( pVoxelHeader != nullptr ) 629 { 630 fvoxelNav.VoxelLocate( pVoxelHeader, fLastLocatedPointLocal ); 631 } 632 break; 633 case kParameterised: 634 if( GetDaughtersRegularStructureId(motherLogical) != 1 ) 635 { 636 // Resets state & returns voxel node 637 // 638 fparamNav.ParamVoxelLocate( pVoxelHeader, fLastLocatedPointLocal ); 639 } 640 break; 641 case kReplica: 642 G4Exception("G4ITNavigator1::LocateGlobalPointWithinVolume()", 643 "GeomNav0001", FatalException, 644 "Not applicable for replicated volumes."); 645 break; 646 case kExternal: 647 G4Exception("G4ITNavigator1::LocateGlobalPointWithinVolume()", 648 "GeomNav0001", FatalException, 649 "Not applicable for external volumes."); 650 break; 651 } 652 } 653 654 // Reset the state variables 655 // - which would have been affected 656 // by the 'equivalent' call to LocateGlobalPointAndSetup 657 // - who's values have been invalidated by the 'move'. 658 // 659 fBlockedPhysicalVolume = nullptr; 660 fBlockedReplicaNo = -1; 661 fEntering = false; 662 fEnteredDaughter = false; // Boundary not encountered, did not enter 663 fExiting = false; 664 fExitedMother = false; // Boundary not encountered, did not exit 665 } 666 667 // !> 668 G4ITNavigatorState_Lock1* G4ITNavigator1::GetNavigatorState() 669 { 670 SetSavedState(); 671 return fpSaveState; 672 } 673 674 void G4ITNavigator1::SetNavigatorState(G4ITNavigatorState_Lock1* navState) 675 { 676 fpSaveState = (G4SaveNavigatorState*) navState; 677 if(navState != nullptr) RestoreSavedState(); 678 } 679 680 void G4ITNavigator1::NewNavigatorState() 681 { 682 fpSaveState = new G4SaveNavigatorState(); 683 ResetState(); 684 } 685 686 687 // ******************************************************************** 688 // SetSavedState 689 // 690 // Save the state, in case this is a parasitic call 691 // Save fValidExitNormal, fExitNormal, fExiting, fEntering, 692 // fBlockedPhysicalVolume, fBlockedReplicaNo, fLastStepWasZero; 693 // ******************************************************************** 694 // 695 void G4ITNavigator1::SetSavedState() 696 { 697 // !> 698 // This check can be avoid if instead, at every first step of a track, 699 // the IT tracking uses NewNavigatorSate 700 // The normal tracking would just call once NewNavigatorState() before tracking 701 702 // if(fpSaveState == 0) 703 // fpSaveState = new G4SaveNavigatorState; 704 // <! 705 706 // fSaveExitNormal = fExitNormal; 707 fpSaveState->sExitNormal = fExitNormal; 708 fpSaveState->sValidExitNormal = fValidExitNormal; 709 fpSaveState->sExiting = fExiting; 710 fpSaveState->sEntering = fEntering; 711 712 fpSaveState->spBlockedPhysicalVolume = fBlockedPhysicalVolume; 713 fpSaveState->sBlockedReplicaNo = fBlockedReplicaNo, 714 715 fpSaveState->sLastStepWasZero = static_cast<G4int>(fLastStepWasZero); 716 717 // !> 718 fpSaveState->sPreviousSftOrigin = fPreviousSftOrigin; 719 fpSaveState->sPreviousSafety = fPreviousSafety; 720 fpSaveState->sNumberZeroSteps = fNumberZeroSteps; 721 fpSaveState->sLocatedOnEdge = fLocatedOnEdge; 722 fpSaveState->sWasLimitedByGeometry= fWasLimitedByGeometry; 723 fpSaveState->sPushed=fPushed; 724 fpSaveState->sNumberZeroSteps=fNumberZeroSteps; 725 fpSaveState->sEnteredDaughter = fEnteredDaughter; 726 fpSaveState->sExitedMother = fExitedMother; 727 728 fpSaveState->sLastLocatedPointLocal = fLastLocatedPointLocal; 729 fpSaveState->sLocatedOutsideWorld = fLocatedOutsideWorld; 730 // <! 731 } 732 733 // ******************************************************************** 734 // RestoreSavedState 735 // 736 // Restore the state (in Compute Step), in case this is a parasitic call 737 // ******************************************************************** 738 // 739 void G4ITNavigator1::RestoreSavedState() 740 { 741 fExitNormal = fpSaveState->sExitNormal; 742 fValidExitNormal = fpSaveState->sValidExitNormal; 743 fExiting = fpSaveState->sExiting; 744 fEntering = fpSaveState->sEntering; 745 746 fBlockedPhysicalVolume = fpSaveState->spBlockedPhysicalVolume; 747 fBlockedReplicaNo = fpSaveState->sBlockedReplicaNo, 748 749 fLastStepWasZero = (fpSaveState->sLastStepWasZero != 0); 750 751 // !> 752 fPreviousSftOrigin = fpSaveState->sPreviousSftOrigin ; 753 fPreviousSafety = fpSaveState->sPreviousSafety ; 754 fNumberZeroSteps = fpSaveState->sNumberZeroSteps ; 755 fLocatedOnEdge = fpSaveState->sLocatedOnEdge ; 756 fWasLimitedByGeometry = fpSaveState->sWasLimitedByGeometry; 757 fPushed = fpSaveState->sPushed; 758 fNumberZeroSteps = fpSaveState->sNumberZeroSteps; 759 fEnteredDaughter= fpSaveState->sEnteredDaughter ; 760 fExitedMother = fpSaveState->sExitedMother ; 761 762 fLastLocatedPointLocal = fpSaveState->sLastLocatedPointLocal ; 763 fLocatedOutsideWorld = fpSaveState->sLocatedOutsideWorld; 764 // <! 765 } 766 // <! 767 768 // ******************************************************************** 769 // ComputeStep 770 // 771 // Computes the next geometric Step: intersections with current 772 // mother and `daughter' volumes. 773 // 774 // NOTE: 775 // 776 // Flags on entry: 777 // -------------- 778 // fValidExitNormal - Normal of exited volume is valid (convex, not a 779 // coincident boundary) 780 // fExitNormal - Surface normal of exited volume 781 // fExiting - True if have exited solid 782 // 783 // fBlockedPhysicalVolume - Ptr to exited volume (or 0) 784 // fBlockedReplicaNo - Replication no of exited volume 785 // fLastStepWasZero - True if last Step size was zero. 786 // 787 // Flags on exit: 788 // ------------- 789 // fValidExitNormal - True if surface normal of exited volume is valid 790 // fExitNormal - Surface normal of exited volume rotated to mothers 791 // reference system 792 // fExiting - True if exiting mother 793 // fEntering - True if entering `daughter' volume (or replica) 794 // fBlockedPhysicalVolume - Ptr to candidate (entered) volume 795 // fBlockedReplicaNo - Replication no of candidate (entered) volume 796 // fLastStepWasZero - True if this Step size was zero. 797 // ******************************************************************** 798 // 799 G4double G4ITNavigator1::ComputeStep( const G4ThreeVector &pGlobalpoint, 800 const G4ThreeVector &pDirection, 801 const G4double pCurrentProposedStepLength, 802 G4double &pNewSafety) 803 { 804 G4ThreeVector localDirection = ComputeLocalAxis(pDirection); 805 G4double Step = kInfinity; 806 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume(); 807 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume(); 808 809 // All state relating to exiting normals must be reset 810 // 811 fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.); 812 // Reset value - to erase its memory 813 fChangedGrandMotherRefFrame= false; 814 // Reset - used for local exit normal 815 fGrandMotherExitNormal= G4ThreeVector( 0., 0., 0.); 816 fCalculatedExitNormal = false; 817 // Reset for new step 818 819 static G4ThreadLocal G4int sNavCScalls=0; 820 sNavCScalls++; 821 822 fLastTriedStepComputation= true; 823 824 #ifdef G4VERBOSE 825 if( fVerbose > 0 ) 826 { 827 G4cout << "*** G4ITNavigator1::ComputeStep: ***" << G4endl; 828 G4cout << " Volume = " << motherPhysical->GetName() 829 << " - Proposed step length = " << pCurrentProposedStepLength 830 << G4endl; 831 #ifdef G4DEBUG_NAVIGATION 832 if( fVerbose >= 2 ) 833 { 834 G4cout << " Called with the arguments: " << G4endl 835 << " Globalpoint = " << std::setw(25) << pGlobalpoint << G4endl 836 << " Direction = " << std::setw(25) << pDirection << G4endl; 837 if( fVerbose >= 4 ) 838 { 839 G4cout << " ---- Upon entering : State" << G4endl; 840 PrintState(); 841 } 842 } 843 #endif 844 } 845 #endif 846 847 G4ThreeVector newLocalPoint = ComputeLocalPoint(pGlobalpoint); 848 if( newLocalPoint != fLastLocatedPointLocal ) 849 { 850 // Check whether the relocation is within safety 851 // 852 G4ThreeVector oldLocalPoint = fLastLocatedPointLocal; 853 G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2(); 854 855 if ( moveLenSq >= kCarTolerance*kCarTolerance ) 856 { 857 #ifdef G4VERBOSE 858 ComputeStepLog(pGlobalpoint, moveLenSq); 859 #endif 860 // Relocate the point within the same volume 861 // 862 LocateGlobalPointWithinVolume( pGlobalpoint ); 863 fLastTriedStepComputation= true; // Ensure that this is set again !! 864 } 865 } 866 if ( fHistory.GetTopVolumeType()!=kReplica ) 867 { 868 switch( CharacteriseDaughters(motherLogical) ) 869 { 870 case kNormal: 871 if ( motherLogical->GetVoxelHeader() != nullptr ) 872 { 873 Step = fvoxelNav.ComputeStep(fLastLocatedPointLocal, 874 localDirection, 875 pCurrentProposedStepLength, 876 pNewSafety, 877 fHistory, 878 fValidExitNormal, 879 fExitNormal, 880 fExiting, 881 fEntering, 882 &fBlockedPhysicalVolume, 883 fBlockedReplicaNo); 884 885 } 886 else 887 { 888 if( motherPhysical->GetRegularStructureId() == 0 ) 889 { 890 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal, 891 localDirection, 892 pCurrentProposedStepLength, 893 pNewSafety, 894 fHistory, 895 fValidExitNormal, 896 fExitNormal, 897 fExiting, 898 fEntering, 899 &fBlockedPhysicalVolume, 900 fBlockedReplicaNo); 901 } 902 else // Regular (non-voxelised) structure 903 { 904 LocateGlobalPointAndSetup( pGlobalpoint, &pDirection, true, true ); 905 fLastTriedStepComputation= true; // Ensure that this is set again !! 906 // 907 // if physical process limits the step, the voxel will not be the 908 // one given by ComputeStepSkippingEqualMaterials() and the local 909 // point will be wrongly calculated. 910 911 // There is a problem: when msc limits the step and the point is 912 // assigned wrongly to phantom in previous step (while it is out 913 // of the container volume). Then LocateGlobalPointAndSetup() has 914 // reset the history topvolume to world. 915 // 916 if(fHistory.GetTopVolume()->GetRegularStructureId() == 0 ) 917 { 918 G4Exception("G4ITNavigator1::ComputeStep()", 919 "GeomNav1001", JustWarning, 920 "Point is relocated in voxels, while it should be outside!"); 921 Step = fnormalNav.ComputeStep(fLastLocatedPointLocal, 922 localDirection, 923 pCurrentProposedStepLength, 924 pNewSafety, 925 fHistory, 926 fValidExitNormal, 927 fExitNormal, 928 fExiting, 929 fEntering, 930 &fBlockedPhysicalVolume, 931 fBlockedReplicaNo); 932 } 933 else 934 { 935 Step = fregularNav. 936 ComputeStepSkippingEqualMaterials(fLastLocatedPointLocal, 937 localDirection, 938 pCurrentProposedStepLength, 939 pNewSafety, 940 fHistory, 941 fValidExitNormal, 942 fExitNormal, 943 fExiting, 944 fEntering, 945 &fBlockedPhysicalVolume, 946 fBlockedReplicaNo, 947 motherPhysical); 948 } 949 } 950 } 951 break; 952 case kParameterised: 953 if( GetDaughtersRegularStructureId(motherLogical) != 1 ) 954 { 955 Step = fparamNav.ComputeStep(fLastLocatedPointLocal, 956 localDirection, 957 pCurrentProposedStepLength, 958 pNewSafety, 959 fHistory, 960 fValidExitNormal, 961 fExitNormal, 962 fExiting, 963 fEntering, 964 &fBlockedPhysicalVolume, 965 fBlockedReplicaNo); 966 } 967 else // Regular structure 968 { 969 Step = fregularNav.ComputeStep(fLastLocatedPointLocal, 970 localDirection, 971 pCurrentProposedStepLength, 972 pNewSafety, 973 fHistory, 974 fValidExitNormal, 975 fExitNormal, 976 fExiting, 977 fEntering, 978 &fBlockedPhysicalVolume, 979 fBlockedReplicaNo); 980 } 981 break; 982 case kReplica: 983 G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0001", 984 FatalException, "Not applicable for replicated volumes."); 985 break; 986 case kExternal: 987 G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0001", 988 FatalException, "Not applicable for external volumes."); 989 break; 990 } 991 } 992 else 993 { 994 // In the case of a replica, it must handle the exiting 995 // edge/corner problem by itself 996 // 997 G4bool exitingReplica = fExitedMother; 998 G4bool calculatedExitNormal= false; 999 1000 Step = freplicaNav.ComputeStep(pGlobalpoint, 1001 pDirection, 1002 fLastLocatedPointLocal, 1003 localDirection, 1004 pCurrentProposedStepLength, 1005 pNewSafety, 1006 fHistory, 1007 fValidExitNormal, 1008 calculatedExitNormal, 1009 fExitNormal, 1010 exitingReplica, 1011 fEntering, 1012 &fBlockedPhysicalVolume, 1013 fBlockedReplicaNo); 1014 fExiting= exitingReplica; 1015 fCalculatedExitNormal= calculatedExitNormal; 1016 } 1017 1018 // Remember last safety origin & value. 1019 // 1020 fPreviousSftOrigin = pGlobalpoint; 1021 fPreviousSafety = pNewSafety; 1022 1023 // Count zero steps - one can occur due to changing momentum at a boundary 1024 // - one, two (or a few) can occur at common edges between 1025 // volumes 1026 // - more than two is likely a problem in the geometry 1027 // description or the Navigation 1028 1029 // Rule of thumb: likely at an Edge if two consecutive steps are zero, 1030 // because at least two candidate volumes must have been 1031 // checked 1032 // 1033 fLocatedOnEdge = fLastStepWasZero && (Step==0.0); 1034 fLastStepWasZero = (Step==0.0); 1035 if (fPushed) { fPushed = fLastStepWasZero; } 1036 1037 // Handle large number of consecutive zero steps 1038 // 1039 if ( fLastStepWasZero ) 1040 { 1041 fNumberZeroSteps++; 1042 #ifdef G4DEBUG_NAVIGATION 1043 if( fNumberZeroSteps > 1 ) 1044 { 1045 G4cout << "G4ITNavigator1::ComputeStep(): another zero step, # " 1046 << fNumberZeroSteps 1047 << " at " << pGlobalpoint 1048 << " in volume " << motherPhysical->GetName() 1049 << " nav-comp-step calls # " << sNavCScalls 1050 << G4endl; 1051 } 1052 #endif 1053 if( fNumberZeroSteps > fActionThreshold_NoZeroSteps-1 ) 1054 { 1055 // Act to recover this stuck track. Pushing it along direction 1056 // 1057 Step += 100*kCarTolerance; 1058 #ifdef G4VERBOSE 1059 if ((!fPushed) && (fWarnPush)) 1060 { 1061 std::ostringstream message; 1062 message << "Track stuck or not moving." << G4endl 1063 << " Track stuck, not moving for " 1064 << fNumberZeroSteps << " steps" << G4endl 1065 << " in volume -" << motherPhysical->GetName() 1066 << "- at point " << pGlobalpoint << G4endl 1067 << " direction: " << pDirection << "." << G4endl 1068 << " Potential geometry or navigation problem !" 1069 << G4endl 1070 << " Trying pushing it of " << Step << " mm ..."; 1071 G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav1002", 1072 JustWarning, message, "Potential overlap in geometry!"); 1073 } 1074 #endif 1075 fPushed = true; 1076 } 1077 if( fNumberZeroSteps > fAbandonThreshold_NoZeroSteps-1 ) 1078 { 1079 // Must kill this stuck track 1080 // 1081 std::ostringstream message; 1082 message << "Stuck Track: potential geometry or navigation problem." 1083 << G4endl 1084 << " Track stuck, not moving for " 1085 << fNumberZeroSteps << " steps" << G4endl 1086 << " in volume -" << motherPhysical->GetName() 1087 << "- at point " << pGlobalpoint << G4endl 1088 << " direction: " << pDirection << "."; 1089 motherPhysical->CheckOverlaps(5000, 0.0); 1090 G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0003", 1091 EventMustBeAborted, message); 1092 } 1093 } 1094 else 1095 { 1096 if (!fPushed) fNumberZeroSteps = 0; 1097 } 1098 1099 fEnteredDaughter = fEntering; // I expect to enter a volume in this Step 1100 fExitedMother = fExiting; 1101 1102 fStepEndPoint = pGlobalpoint + Step * pDirection; 1103 fLastStepEndPointLocal = fLastLocatedPointLocal + Step * localDirection; 1104 1105 if( fExiting ) 1106 { 1107 #ifdef G4DEBUG_NAVIGATION 1108 if( fVerbose > 2 ) 1109 { 1110 G4cout << " At G4Nav CompStep End - if(exiting) - fExiting= " << fExiting 1111 << " fValidExitNormal = " << fValidExitNormal << G4endl; 1112 G4cout << " fExitNormal= " << fExitNormal << G4endl; 1113 } 1114 #endif 1115 1116 if(fValidExitNormal || fCalculatedExitNormal) 1117 { 1118 if ( fHistory.GetTopVolumeType()!=kReplica ) 1119 { 1120 // Convention: fExitNormal is in the 'grand-mother' coordinate system 1121 // 1122 fGrandMotherExitNormal= fExitNormal; 1123 fCalculatedExitNormal= true; 1124 } 1125 else 1126 { 1127 fGrandMotherExitNormal = fExitNormal; 1128 } 1129 } 1130 else 1131 { 1132 // We must calculate the normal anyway (in order to have it if requested) 1133 // 1134 G4ThreeVector finalLocalPoint = 1135 fLastLocatedPointLocal + localDirection*Step; 1136 1137 if ( fHistory.GetTopVolumeType()!=kReplica ) 1138 { 1139 // Find normal in the 'mother' coordinate system 1140 // 1141 G4ThreeVector exitNormalMotherFrame= 1142 motherLogical->GetSolid()->SurfaceNormal(finalLocalPoint); 1143 1144 // Transform it to the 'grand-mother' coordinate system 1145 // 1146 const G4RotationMatrix* mRot = motherPhysical->GetRotation(); 1147 if( mRot != nullptr ) 1148 { 1149 fChangedGrandMotherRefFrame= true; 1150 fGrandMotherExitNormal = (*mRot).inverse() * exitNormalMotherFrame; 1151 } 1152 else 1153 { 1154 fGrandMotherExitNormal = exitNormalMotherFrame; 1155 } 1156 1157 // Do not set fValidExitNormal -- this signifies 1158 // that the solid is convex! 1159 // 1160 fCalculatedExitNormal= true; 1161 } 1162 else 1163 { 1164 fCalculatedExitNormal = false; 1165 // 1166 // Nothing can be done at this stage currently - to solve this 1167 // Replica Navigation must have calculated the normal for this case 1168 // already. 1169 // Cases: mother is not convex, and exit is at previous replica level 1170 1171 #ifdef G4DEBUG_NAVIGATION 1172 G4ExceptionDescription desc; 1173 1174 desc << "Problem in ComputeStep: Replica Navigation did not provide" 1175 << " valid exit Normal. " << G4endl; 1176 desc << " Do not know how calculate it in this case." << G4endl; 1177 desc << " Location = " << finalLocalPoint << G4endl; 1178 desc << " Volume name = " << motherPhysical->GetName() 1179 << " copy/replica No = " << motherPhysical->GetCopyNo() << G4endl; 1180 G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav0003", 1181 JustWarning, desc, "Normal not available for exiting."); 1182 #endif 1183 } 1184 } 1185 1186 // Now transform it to the global reference frame !! 1187 // 1188 if( fValidExitNormal || fCalculatedExitNormal ) 1189 { 1190 auto depth = (G4int)fHistory.GetDepth(); 1191 if( depth > 0 ) 1192 { 1193 G4AffineTransform GrandMotherToGlobalTransf = 1194 fHistory.GetTransform(depth-1).Inverse(); 1195 fExitNormalGlobalFrame = 1196 GrandMotherToGlobalTransf.TransformAxis( fGrandMotherExitNormal ); 1197 } 1198 else 1199 { 1200 fExitNormalGlobalFrame= fGrandMotherExitNormal; 1201 } 1202 } 1203 else 1204 { 1205 fExitNormalGlobalFrame= G4ThreeVector( 0., 0., 0.); 1206 } 1207 } 1208 fStepEndPoint= pGlobalpoint+Step*pDirection; 1209 1210 if( (Step == pCurrentProposedStepLength) && (!fExiting) && (!fEntering) ) 1211 { 1212 // This if Step is not really limited by the geometry. 1213 // The Navigator is obliged to return "infinity" 1214 // 1215 Step = kInfinity; 1216 } 1217 1218 #ifdef G4VERBOSE 1219 if( fVerbose > 1 ) 1220 { 1221 if( fVerbose >= 4 ) 1222 { 1223 G4cout << " ----- Upon exiting :" << G4endl; 1224 PrintState(); 1225 } 1226 G4cout << " Returned step= " << Step; 1227 if( fVerbose > 5 ) G4cout << G4endl; 1228 if( Step == kInfinity ) 1229 { 1230 G4cout << " Requested step= " << pCurrentProposedStepLength ; 1231 if( fVerbose > 5) G4cout << G4endl; 1232 } 1233 G4cout << " Safety = " << pNewSafety << G4endl; 1234 } 1235 #endif 1236 1237 return Step; 1238 } 1239 1240 // ******************************************************************** 1241 // CheckNextStep 1242 // 1243 // Compute the step without altering the navigator state 1244 // ******************************************************************** 1245 // 1246 G4double G4ITNavigator1::CheckNextStep( const G4ThreeVector& pGlobalpoint, 1247 const G4ThreeVector& pDirection, 1248 const G4double pCurrentProposedStepLength, 1249 G4double& pNewSafety) 1250 { 1251 G4double step; 1252 1253 // Save the state, for this parasitic call 1254 // 1255 SetSavedState(); 1256 1257 step = ComputeStep ( pGlobalpoint, 1258 pDirection, 1259 pCurrentProposedStepLength, 1260 pNewSafety ); 1261 1262 // If a parasitic call, then attempt to restore the key parts of the state 1263 // 1264 RestoreSavedState(); 1265 1266 return step; 1267 } 1268 1269 // ******************************************************************** 1270 // ResetState 1271 // 1272 // Resets stack and minimum of navigator state `machine' 1273 // ******************************************************************** 1274 // 1275 void G4ITNavigator1::ResetState() 1276 { 1277 fWasLimitedByGeometry = false; 1278 fEntering = false; 1279 fExiting = false; 1280 fLocatedOnEdge = false; 1281 fLastStepWasZero = false; 1282 fEnteredDaughter = false; 1283 fExitedMother = false; 1284 fPushed = false; 1285 1286 fValidExitNormal = false; 1287 fChangedGrandMotherRefFrame= false; 1288 fCalculatedExitNormal = false; 1289 1290 fExitNormal = G4ThreeVector(0,0,0); 1291 fGrandMotherExitNormal = G4ThreeVector(0,0,0); 1292 fExitNormalGlobalFrame = G4ThreeVector(0,0,0); 1293 1294 fPreviousSftOrigin = G4ThreeVector(0,0,0); 1295 fPreviousSafety = 0.0; 1296 1297 fNumberZeroSteps = 0; 1298 1299 fBlockedPhysicalVolume = nullptr; 1300 fBlockedReplicaNo = -1; 1301 1302 fLastLocatedPointLocal = G4ThreeVector( kInfinity, -kInfinity, 0.0 ); 1303 fLocatedOutsideWorld = false; 1304 } 1305 1306 // ******************************************************************** 1307 // SetupHierarchy 1308 // 1309 // Renavigates & resets hierarchy described by current history 1310 // o Reset volumes 1311 // o Recompute transforms and/or solids of replicated/parameterised volumes 1312 // ******************************************************************** 1313 // 1314 void G4ITNavigator1::SetupHierarchy() 1315 { 1316 G4int i; 1317 const auto cdepth = (G4int)fHistory.GetDepth(); 1318 G4VPhysicalVolume *current; 1319 G4VSolid *pSolid; 1320 G4VPVParameterisation *pParam; 1321 1322 for ( i=1; i<=cdepth; i++ ) 1323 { 1324 current = fHistory.GetVolume(i); 1325 switch ( fHistory.GetVolumeType(i) ) 1326 { 1327 case kNormal: 1328 break; 1329 case kReplica: 1330 freplicaNav.ComputeTransformation(fHistory.GetReplicaNo(i), current); 1331 break; 1332 case kParameterised: { 1333 G4int replicaNo; 1334 pParam = current->GetParameterisation(); 1335 replicaNo = fHistory.GetReplicaNo(i); 1336 pSolid = pParam->ComputeSolid(replicaNo, current); 1337 1338 // Set up dimensions & transform in solid/physical volume 1339 // 1340 pSolid->ComputeDimensions(pParam, replicaNo, current); 1341 pParam->ComputeTransformation(replicaNo, current); 1342 1343 G4TouchableHistory touchable( fHistory ); 1344 touchable.MoveUpHistory(); // move up to the parent level 1345 1346 // Set up the correct solid and material in Logical Volume 1347 // 1348 G4LogicalVolume *pLogical = current->GetLogicalVolume(); 1349 pLogical->SetSolid( pSolid ); 1350 pLogical->UpdateMaterial( pParam -> 1351 ComputeMaterial(replicaNo, current, &touchable) ); 1352 break; 1353 } 1354 case kExternal: 1355 G4Exception("G4ITNavigator1::SetupHierarchy()", 1356 "GeomNav0001", FatalException, 1357 "Not applicable for external volumes."); 1358 break; 1359 } 1360 } 1361 } 1362 1363 // ******************************************************************** 1364 // GetLocalExitNormal 1365 // 1366 // Obtains the Normal vector to a surface (in local coordinates) 1367 // pointing out of previous volume and into current volume 1368 // ******************************************************************** 1369 // 1370 G4ThreeVector G4ITNavigator1::GetLocalExitNormal( G4bool* valid ) 1371 { 1372 G4ThreeVector ExitNormal(0.,0.,0.); 1373 G4VSolid *currentSolid=nullptr; 1374 G4LogicalVolume *candidateLogical; 1375 if ( fLastTriedStepComputation ) 1376 { 1377 // use fLastLocatedPointLocal and next candidate volume 1378 // 1379 G4ThreeVector nextSolidExitNormal(0.,0.,0.); 1380 1381 if( fEntering && (fBlockedPhysicalVolume!=nullptr) ) 1382 { 1383 candidateLogical= fBlockedPhysicalVolume->GetLogicalVolume(); 1384 if( candidateLogical != nullptr ) 1385 { 1386 // fLastStepEndPointLocal is in the coordinates of the mother 1387 // we need it in the daughter's coordinate system. 1388 1389 // The following code should also work in case of Replica 1390 { 1391 // First transform fLastLocatedPointLocal to the new daughter 1392 // coordinates 1393 // 1394 G4AffineTransform MotherToDaughterTransform= 1395 GetMotherToDaughterTransform( fBlockedPhysicalVolume, 1396 fBlockedReplicaNo, 1397 VolumeType(fBlockedPhysicalVolume) ); 1398 G4ThreeVector daughterPointOwnLocal= 1399 MotherToDaughterTransform.TransformPoint( fLastStepEndPointLocal ); 1400 1401 // OK if it is a parameterised volume 1402 // 1403 EInside inSideIt; 1404 G4bool onSurface; 1405 G4double safety= -1.0; 1406 currentSolid= candidateLogical->GetSolid(); 1407 inSideIt = currentSolid->Inside(daughterPointOwnLocal); 1408 onSurface = (inSideIt == kSurface); 1409 if( ! onSurface ) 1410 { 1411 if( inSideIt == kOutside ) 1412 { 1413 safety = (currentSolid->DistanceToIn(daughterPointOwnLocal)); 1414 onSurface = safety < 100.0 * kCarTolerance; 1415 } 1416 else if (inSideIt == kInside ) 1417 { 1418 safety = (currentSolid->DistanceToOut(daughterPointOwnLocal)); 1419 onSurface = safety < 100.0 * kCarTolerance; 1420 } 1421 } 1422 1423 if( onSurface ) 1424 { 1425 nextSolidExitNormal = 1426 currentSolid->SurfaceNormal(daughterPointOwnLocal); 1427 1428 // Entering the solid ==> opposite 1429 // 1430 ExitNormal = -nextSolidExitNormal; 1431 fCalculatedExitNormal= true; 1432 } 1433 else 1434 { 1435 #ifdef G4VERBOSE 1436 if(( fVerbose == 1 ) && ( fCheck )) 1437 { 1438 std::ostringstream message; 1439 message << "Point not on surface ! " << G4endl 1440 << " Point = " 1441 << daughterPointOwnLocal << G4endl 1442 << " Physical volume = " 1443 << fBlockedPhysicalVolume->GetName() << G4endl 1444 << " Logical volume = " 1445 << candidateLogical->GetName() << G4endl 1446 << " Solid = " << currentSolid->GetName() 1447 << " Type = " 1448 << currentSolid->GetEntityType() << G4endl 1449 << *currentSolid << G4endl; 1450 if( inSideIt == kOutside ) 1451 { 1452 message << "Point is Outside. " << G4endl 1453 << " Safety (from outside) = " << safety << G4endl; 1454 } 1455 else // if( inSideIt == kInside ) 1456 { 1457 message << "Point is Inside. " << G4endl 1458 << " Safety (from inside) = " << safety << G4endl; 1459 } 1460 G4Exception("G4ITNavigator1::GetLocalExitNormal()", "GeomNav1001", 1461 JustWarning, message); 1462 } 1463 #endif 1464 } 1465 *valid = onSurface; // was =true; 1466 } 1467 } 1468 } 1469 else if ( fExiting ) 1470 { 1471 ExitNormal = fGrandMotherExitNormal; 1472 *valid = true; 1473 fCalculatedExitNormal= true; // Should be true already 1474 } 1475 else // i.e. ( fBlockedPhysicalVolume == 0 ) 1476 { 1477 *valid = false; 1478 G4Exception("G4ITNavigator1::GetLocalExitNormal()", 1479 "GeomNav0003", JustWarning, 1480 "Incorrect call to GetLocalSurfaceNormal." ); 1481 } 1482 } 1483 else // ( ! fLastTriedStepComputation ) ie. last call was to Locate 1484 { 1485 if ( EnteredDaughterVolume() ) 1486 { 1487 G4VSolid* daughterSolid =fHistory.GetTopVolume()->GetLogicalVolume() 1488 ->GetSolid(); 1489 ExitNormal= -(daughterSolid->SurfaceNormal(fLastLocatedPointLocal)); 1490 if( std::fabs(ExitNormal.mag2()-1.0 ) > CLHEP::perMillion ) 1491 { 1492 G4ExceptionDescription desc; 1493 desc << " Parameters of solid: " << *daughterSolid 1494 << " Point for surface = " << fLastLocatedPointLocal << std::endl; 1495 G4Exception("G4ITNavigator1::GetLocalExitNormal()", 1496 "GeomNav0003", FatalException, desc, 1497 "Surface Normal returned by Solid is not a Unit Vector." ); 1498 } 1499 fCalculatedExitNormal= true; 1500 *valid = true; 1501 } 1502 else 1503 { 1504 if( fExitedMother ) 1505 { 1506 ExitNormal = fGrandMotherExitNormal; 1507 *valid = true; 1508 fCalculatedExitNormal= true; 1509 } 1510 else // We are not at a boundary. ExitNormal remains (0,0,0) 1511 { 1512 *valid = false; 1513 fCalculatedExitNormal= false; 1514 G4ExceptionDescription message; 1515 message << "Function called when *NOT* at a Boundary." << G4endl; 1516 G4Exception("G4ITNavigator1::GetLocalExitNormal()", 1517 "GeomNav0003", JustWarning, message); 1518 } 1519 } 1520 } 1521 return ExitNormal; 1522 } 1523 1524 // ******************************************************************** 1525 // GetMotherToDaughterTransform 1526 // 1527 // Obtains the mother to daughter affine transformation 1528 // ******************************************************************** 1529 // 1530 G4AffineTransform 1531 G4ITNavigator1::GetMotherToDaughterTransform( G4VPhysicalVolume *pEnteringPhysVol, // not Const 1532 G4int enteringReplicaNo, 1533 EVolume enteringVolumeType ) 1534 { 1535 switch (enteringVolumeType) 1536 { 1537 case kNormal: // Nothing is needed to prepare the transformation 1538 break; // It is stored already in the physical volume (placement) 1539 case kReplica: // Sets the transform in the Replica - tbc 1540 G4Exception("G4ITNavigator1::GetMotherToDaughterTransform()", 1541 "GeomNav0001", FatalException, 1542 "Method NOT Implemented yet for replica volumes."); 1543 break; 1544 case kParameterised: 1545 if( pEnteringPhysVol->GetRegularStructureId() == 0 ) 1546 { 1547 G4VPVParameterisation *pParam = 1548 pEnteringPhysVol->GetParameterisation(); 1549 G4VSolid* pSolid = 1550 pParam->ComputeSolid(enteringReplicaNo, pEnteringPhysVol); 1551 pSolid->ComputeDimensions(pParam, enteringReplicaNo, pEnteringPhysVol); 1552 1553 // Sets the transform in the Parameterisation 1554 // 1555 pParam->ComputeTransformation(enteringReplicaNo, pEnteringPhysVol); 1556 1557 // Set the correct solid and material in Logical Volume 1558 // 1559 G4LogicalVolume* pLogical = pEnteringPhysVol->GetLogicalVolume(); 1560 pLogical->SetSolid( pSolid ); 1561 } 1562 break; 1563 case kExternal: 1564 G4Exception("G4ITNavigator1::GetMotherToDaughterTransform()", 1565 "GeomNav0001", FatalException, 1566 "Not applicable for external volumes."); 1567 break; 1568 } 1569 return G4AffineTransform(pEnteringPhysVol->GetRotation(), 1570 pEnteringPhysVol->GetTranslation()).Invert(); 1571 } 1572 1573 // ******************************************************************** 1574 // GetLocalExitNormalAndCheck 1575 // 1576 // Obtains the Normal vector to a surface (in local coordinates) 1577 // pointing out of previous volume and into current volume, and 1578 // checks the current point against expected 'local' value. 1579 // ******************************************************************** 1580 // 1581 G4ThreeVector G4ITNavigator1:: 1582 GetLocalExitNormalAndCheck( 1583 #ifdef G4DEBUG_NAVIGATION 1584 const G4ThreeVector& ExpectedBoundaryPointGlobal, 1585 #else 1586 const G4ThreeVector&, 1587 #endif 1588 G4bool* pValid) 1589 { 1590 #ifdef G4DEBUG_NAVIGATION 1591 // Check Current point against expected 'local' value 1592 // 1593 if ( fLastTriedStepComputation ) 1594 { 1595 G4ThreeVector ExpectedBoundaryPointLocal; 1596 1597 const G4AffineTransform& GlobalToLocal= GetGlobalToLocalTransform(); 1598 ExpectedBoundaryPointLocal = 1599 GlobalToLocal.TransformPoint( ExpectedBoundaryPointGlobal ); 1600 1601 // Add here: Comparison against expected position, 1602 // i.e. the endpoint of ComputeStep 1603 } 1604 #endif 1605 1606 return GetLocalExitNormal( pValid); 1607 } 1608 1609 // ******************************************************************** 1610 // GetGlobalExitNormal 1611 // 1612 // Obtains the Normal vector to a surface (in global coordinates) 1613 // pointing out of previous volume and into current volume 1614 // ******************************************************************** 1615 // 1616 G4ThreeVector 1617 G4ITNavigator1::GetGlobalExitNormal(const G4ThreeVector& IntersectPointGlobal, 1618 G4bool* pNormalCalculated) 1619 { 1620 G4bool validNormal; 1621 G4ThreeVector localNormal, globalNormal; 1622 1623 if( fLastTriedStepComputation && fExiting ) 1624 { 1625 // This was computed in ComputeStep -- and only on arrival at boundary 1626 // 1627 globalNormal = fExitNormalGlobalFrame; 1628 *pNormalCalculated = true; // ComputeStep always computes it if Exiting 1629 // (fExiting==true) 1630 } 1631 else 1632 { 1633 localNormal = GetLocalExitNormalAndCheck(IntersectPointGlobal,&validNormal); 1634 *pNormalCalculated = fCalculatedExitNormal; 1635 1636 #ifdef G4DEBUG_NAVIGATION 1637 if( (!validNormal) && !fCalculatedExitNormal) 1638 { 1639 G4ExceptionDescription edN; 1640 edN << " Calculated = " << fCalculatedExitNormal << G4endl; 1641 edN << " Entering= " << fEntering << G4endl; 1642 G4int oldVerbose= this->GetVerboseLevel(); 1643 this->SetVerboseLevel(4); 1644 edN << " State of Navigator: " << G4endl; 1645 edN << *this << G4endl; 1646 this->SetVerboseLevel( oldVerbose ); 1647 1648 G4Exception("G4ITNavigator1::GetGlobalExitNormal()", 1649 "GeomNav0003", JustWarning, edN, 1650 "LocalExitNormalAndCheck() did not calculate Normal."); 1651 } 1652 #endif 1653 1654 G4double localMag2= localNormal.mag2(); 1655 if( validNormal && (std::fabs(localMag2-1.0)) > CLHEP::perMillion ) 1656 { 1657 G4ExceptionDescription edN; 1658 1659 edN << "G4ITNavigator1::GetGlobalExitNormal: " 1660 << " Using Local Normal - from call to GetLocalExitNormalAndCheck. " 1661 << G4endl 1662 << " Local Exit Normal = " << localNormal << " || = " 1663 << std::sqrt(localMag2) << G4endl 1664 << " Global Exit Normal = " << globalNormal << " || = " 1665 << globalNormal.mag() << G4endl; 1666 edN << " Calculated It = " << fCalculatedExitNormal << G4endl; 1667 1668 G4Exception("G4ITNavigator1::GetGlobalExitNormal()", 1669 "GeomNav0003",JustWarning, edN, 1670 "Value obtained from new local *solid* is incorrect."); 1671 localNormal = localNormal.unit(); // Should we correct it ?? 1672 } 1673 G4AffineTransform localToGlobal = GetLocalToGlobalTransform(); 1674 globalNormal = localToGlobal.TransformAxis( localNormal ); 1675 } 1676 1677 #ifdef G4DEBUG_NAVIGATION 1678 // Temporary extra checks 1679 if( fLastTriedStepComputation && fExiting) 1680 { 1681 localNormal = GetLocalExitNormalAndCheck( IntersectPointGlobal, &validNormal); 1682 *pNormalCalculated = fCalculatedExitNormal; 1683 1684 G4AffineTransform localToGlobal = GetLocalToGlobalTransform(); 1685 globalNormal = localToGlobal.TransformAxis( localNormal ); 1686 1687 // Check the value computed against fExitNormalGlobalFrame 1688 G4ThreeVector diffNorm = globalNormal - fExitNormalGlobalFrame; 1689 if( diffNorm.mag2() > perMillion*CLHEP::perMillion) 1690 { 1691 G4ExceptionDescription edDfn; 1692 edDfn << "Found difference in normals in case of exiting mother " 1693 << "- when Get is called after ComputingStep " << G4endl; 1694 edDfn << " Magnitude of diff = " << diffNorm.mag() << G4endl; 1695 edDfn << " Normal stored (Global) = " << fExitNormalGlobalFrame 1696 << G4endl; 1697 edDfn << " Global Computed from Local = " << globalNormal << G4endl; 1698 G4Exception("G4ITNavigator1::GetGlobalExitNormal()", "GeomNav0003", 1699 JustWarning, edDfn); 1700 } 1701 } 1702 #endif 1703 1704 return globalNormal; 1705 } 1706 1707 // To make the new Voxel Safety the default, uncomment the next line 1708 #define G4NEW_SAFETY 1 1709 1710 // ******************************************************************** 1711 // ComputeSafety 1712 // 1713 // It assumes that it will be 1714 // i) called at the Point in the same volume as the EndPoint of the 1715 // ComputeStep. 1716 // ii) after (or at the end of) ComputeStep OR after the relocation. 1717 // ******************************************************************** 1718 // 1719 G4double G4ITNavigator1::ComputeSafety( const G4ThreeVector &pGlobalpoint, 1720 const G4double pMaxLength, 1721 const G4bool keepState) 1722 { 1723 G4double newSafety = 0.0; 1724 1725 #ifdef G4DEBUG_NAVIGATION 1726 G4long oldcoutPrec = G4cout.precision(8); 1727 if( fVerbose > 0 ) 1728 { 1729 G4cout << "*** G4ITNavigator1::ComputeSafety: ***" << G4endl 1730 << " Called at point: " << pGlobalpoint << G4endl; 1731 1732 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume(); 1733 G4cout << " Volume = " << motherPhysical->GetName() 1734 << " - Maximum length = " << pMaxLength << G4endl; 1735 if( fVerbose >= 4 ) 1736 { 1737 G4cout << " ----- Upon entering Compute Safety:" << G4endl; 1738 PrintState(); 1739 } 1740 } 1741 #endif 1742 1743 if (keepState) { SetSavedState(); } 1744 1745 G4double distEndpointSq = (pGlobalpoint-fStepEndPoint).mag2(); 1746 G4bool stayedOnEndpoint = distEndpointSq < kCarTolerance*kCarTolerance; 1747 G4bool endpointOnSurface = fEnteredDaughter || fExitedMother; 1748 1749 if( !(endpointOnSurface && stayedOnEndpoint) ) 1750 { 1751 // Pseudo-relocate to this point (updates voxel information only) 1752 // 1753 LocateGlobalPointWithinVolume( pGlobalpoint ); 1754 // --->> DANGER: Side effects on sub-navigator voxel information <<--- 1755 // Could be replaced again by 'granular' calls to sub-navigator 1756 // locates (similar side-effects, but faster. 1757 // Solutions: 1758 // 1) Re-locate (to where?) 1759 // 2) Insure that the methods using (G4ComputeStep?) 1760 // does a relocation (if information is disturbed only ?) 1761 1762 #ifdef G4DEBUG_NAVIGATION 1763 if( fVerbose >= 2 ) 1764 { 1765 G4cout << " G4ITNavigator1::ComputeSafety() relocates-in-volume to point: " 1766 << pGlobalpoint << G4endl; 1767 } 1768 #endif 1769 G4VPhysicalVolume *motherPhysical = fHistory.GetTopVolume(); 1770 G4LogicalVolume *motherLogical = motherPhysical->GetLogicalVolume(); 1771 G4SmartVoxelHeader* pVoxelHeader = motherLogical->GetVoxelHeader(); 1772 G4ThreeVector localPoint = ComputeLocalPoint(pGlobalpoint); 1773 1774 if ( fHistory.GetTopVolumeType()!=kReplica ) 1775 { 1776 switch(CharacteriseDaughters(motherLogical)) 1777 { 1778 case kNormal: 1779 if ( pVoxelHeader != nullptr ) 1780 { 1781 #ifdef G4NEW_SAFETY 1782 G4double safetyTwo = fpVoxelSafety->ComputeSafety(localPoint, 1783 *motherPhysical, pMaxLength); 1784 newSafety= safetyTwo; // Faster and best available 1785 #else 1786 G4double safetyOldVoxel; 1787 safetyOldVoxel = 1788 fvoxelNav.ComputeSafety(localPoint,fHistory,pMaxLength); 1789 newSafety= safetyOldVoxel; 1790 #endif 1791 } 1792 else 1793 { 1794 newSafety=fnormalNav.ComputeSafety(localPoint,fHistory,pMaxLength); 1795 } 1796 break; 1797 case kParameterised: 1798 if( GetDaughtersRegularStructureId(motherLogical) != 1 ) 1799 { 1800 newSafety=fparamNav.ComputeSafety(localPoint,fHistory,pMaxLength); 1801 } 1802 else // Regular structure 1803 { 1804 newSafety=fregularNav.ComputeSafety(localPoint,fHistory,pMaxLength); 1805 } 1806 break; 1807 case kReplica: 1808 G4Exception("G4ITNavigator1::ComputeSafety()", "GeomNav0001", 1809 FatalException, "Not applicable for replicated volumes."); 1810 break; 1811 case kExternal: 1812 G4Exception("G4ITNavigator1::ComputeSafety()", 1813 "GeomNav0001", FatalException, 1814 "Not applicable for external volumes."); 1815 break; 1816 } 1817 } 1818 else 1819 { 1820 newSafety = freplicaNav.ComputeSafety(pGlobalpoint, localPoint, 1821 fHistory, pMaxLength); 1822 } 1823 } 1824 else // if( endpointOnSurface && stayedOnEndpoint ) 1825 { 1826 #ifdef G4DEBUG_NAVIGATION 1827 if( fVerbose >= 2 ) 1828 { 1829 G4cout << " G4ITNavigator1::ComputeSafety() finds that point - " 1830 << pGlobalpoint << " - is on surface " << G4endl; 1831 if( fEnteredDaughter ) { G4cout << " entered new daughter volume"; } 1832 if( fExitedMother ) { G4cout << " and exited previous volume."; } 1833 G4cout << G4endl; 1834 G4cout << " EndPoint was = " << fStepEndPoint << G4endl; 1835 } 1836 #endif 1837 newSafety = 0.0; 1838 } 1839 1840 // Remember last safety origin & value 1841 // 1842 fPreviousSftOrigin = pGlobalpoint; 1843 fPreviousSafety = newSafety; 1844 1845 if (keepState) { RestoreSavedState(); } 1846 1847 #ifdef G4DEBUG_NAVIGATION 1848 if( fVerbose > 1 ) 1849 { 1850 G4cout << " ---- Exiting ComputeSafety " << G4endl; 1851 if( fVerbose > 2 ) { PrintState(); } 1852 G4cout << " Returned value of Safety = " << newSafety << G4endl; 1853 } 1854 G4cout.precision(oldcoutPrec); 1855 #endif 1856 1857 return newSafety; 1858 } 1859 1860 // ******************************************************************** 1861 // CreateTouchableHistoryHandle 1862 // ******************************************************************** 1863 // 1864 G4TouchableHandle G4ITNavigator1::CreateTouchableHistoryHandle() const 1865 { 1866 return G4TouchableHandle( CreateTouchableHistory() ); 1867 } 1868 1869 // ******************************************************************** 1870 // PrintState 1871 // ******************************************************************** 1872 // 1873 void G4ITNavigator1::PrintState() const 1874 { 1875 G4long oldcoutPrec = G4cout.precision(4); 1876 if( fVerbose == 4 ) 1877 { 1878 G4cout << "The current state of G4ITNavigator1 is: " << G4endl; 1879 G4cout << " ValidExitNormal= " << fValidExitNormal << G4endl 1880 << " ExitNormal = " << fExitNormal << G4endl 1881 << " Exiting = " << fExiting << G4endl 1882 << " Entering = " << fEntering << G4endl 1883 << " BlockedPhysicalVolume= " ; 1884 if (fBlockedPhysicalVolume==nullptr) 1885 G4cout << "None"; 1886 else 1887 G4cout << fBlockedPhysicalVolume->GetName(); 1888 G4cout << G4endl 1889 << " BlockedReplicaNo = " << fBlockedReplicaNo << G4endl 1890 << " LastStepWasZero = " << fLastStepWasZero << G4endl 1891 << G4endl; 1892 } 1893 if( ( 1 < fVerbose) && (fVerbose < 4) ) 1894 { 1895 G4cout << G4endl; // Make sure to line up 1896 G4cout << std::setw(30) << " ExitNormal " << " " 1897 << std::setw( 5) << " Valid " << " " 1898 << std::setw( 9) << " Exiting " << " " 1899 << std::setw( 9) << " Entering" << " " 1900 << std::setw(15) << " Blocked:Volume " << " " 1901 << std::setw( 9) << " ReplicaNo" << " " 1902 << std::setw( 8) << " LastStepZero " << " " 1903 << G4endl; 1904 G4cout << "( " << std::setw(7) << fExitNormal.x() 1905 << ", " << std::setw(7) << fExitNormal.y() 1906 << ", " << std::setw(7) << fExitNormal.z() << " ) " 1907 << std::setw( 5) << fValidExitNormal << " " 1908 << std::setw( 9) << fExiting << " " 1909 << std::setw( 9) << fEntering << " "; 1910 if ( fBlockedPhysicalVolume==nullptr ) 1911 { 1912 G4cout << std::setw(15) << "None"; 1913 } 1914 else 1915 { 1916 G4cout << std::setw(15)<< fBlockedPhysicalVolume->GetName(); 1917 } 1918 G4cout << std::setw( 9) << fBlockedReplicaNo << " " 1919 << std::setw( 8) << fLastStepWasZero << " " 1920 << G4endl; 1921 } 1922 if( fVerbose > 2 ) 1923 { 1924 G4cout.precision(8); 1925 G4cout << " Current Localpoint = " << fLastLocatedPointLocal << G4endl; 1926 G4cout << " PreviousSftOrigin = " << fPreviousSftOrigin << G4endl; 1927 G4cout << " PreviousSafety = " << fPreviousSafety << G4endl; 1928 } 1929 G4cout.precision(oldcoutPrec); 1930 } 1931 1932 // ******************************************************************** 1933 // ComputeStepLog 1934 // ******************************************************************** 1935 // 1936 void G4ITNavigator1::ComputeStepLog(const G4ThreeVector& pGlobalpoint, 1937 G4double moveLenSq) const 1938 { 1939 // The following checks only make sense if the move is larger 1940 // than the tolerance. 1941 1942 static const G4double fAccuracyForWarning = kCarTolerance, 1943 fAccuracyForException = 1000*kCarTolerance; 1944 1945 G4ThreeVector OriginalGlobalpoint = fHistory.GetTopTransform().Inverse(). 1946 TransformPoint(fLastLocatedPointLocal); 1947 1948 G4double shiftOriginSafSq = (fPreviousSftOrigin-pGlobalpoint).mag2(); 1949 1950 // Check that the starting point of this step is 1951 // within the isotropic safety sphere of the last point 1952 // to a accuracy/precision given by fAccuracyForWarning. 1953 // If so give warning. 1954 // If it fails by more than fAccuracyForException exit with error. 1955 // 1956 if( shiftOriginSafSq >= sqr(fPreviousSafety) ) 1957 { 1958 G4double shiftOrigin = std::sqrt(shiftOriginSafSq); 1959 G4double diffShiftSaf = shiftOrigin - fPreviousSafety; 1960 1961 if( diffShiftSaf > fAccuracyForWarning ) 1962 { 1963 G4long oldcoutPrec= G4cout.precision(8); 1964 G4long oldcerrPrec= G4cerr.precision(10); 1965 std::ostringstream message, suggestion; 1966 message << "Accuracy error or slightly inaccurate position shift." 1967 << G4endl 1968 << " The Step's starting point has moved " 1969 << std::sqrt(moveLenSq)/mm << " mm " << G4endl 1970 << " since the last call to a Locate method." << G4endl 1971 << " This has resulted in moving " 1972 << shiftOrigin/mm << " mm " 1973 << " from the last point at which the safety " 1974 << " was calculated " << G4endl 1975 << " which is more than the computed safety= " 1976 << fPreviousSafety/mm << " mm at that point." << G4endl 1977 << " This difference is " 1978 << diffShiftSaf/mm << " mm." << G4endl 1979 << " The tolerated accuracy is " 1980 << fAccuracyForException/mm << " mm."; 1981 1982 suggestion << " "; 1983 static G4ThreadLocal G4int warnNow = 0; 1984 if( ((++warnNow % 100) == 1) ) 1985 { 1986 message << G4endl 1987 << " This problem can be due to either " << G4endl 1988 << " - a process that has proposed a displacement" 1989 << " larger than the current safety , or" << G4endl 1990 << " - inaccuracy in the computation of the safety"; 1991 suggestion << "We suggest that you " << G4endl 1992 << " - find i) what particle is being tracked, and " 1993 << " ii) through what part of your geometry " << G4endl 1994 << " for example by re-running this event with " 1995 << G4endl 1996 << " /tracking/verbose 1 " << G4endl 1997 << " - check which processes you declare for" 1998 << " this particle (and look at non-standard ones)" 1999 << G4endl 2000 << " - in case, create a detailed logfile" 2001 << " of this event using:" << G4endl 2002 << " /tracking/verbose 6 "; 2003 } 2004 G4Exception("G4ITNavigator1::ComputeStep()", 2005 "GeomNav1002", JustWarning, 2006 message, G4String(suggestion.str())); 2007 G4cout.precision(oldcoutPrec); 2008 G4cerr.precision(oldcerrPrec); 2009 } 2010 #ifdef G4DEBUG_NAVIGATION 2011 else 2012 { 2013 G4cerr << "WARNING - G4ITNavigator1::ComputeStep()" << G4endl 2014 << " The Step's starting point has moved " 2015 << std::sqrt(moveLenSq) << "," << G4endl 2016 << " which has taken it to the limit of" 2017 << " the current safety. " << G4endl; 2018 } 2019 #endif 2020 } 2021 G4double safetyPlus = fPreviousSafety + fAccuracyForException; 2022 if ( shiftOriginSafSq > sqr(safetyPlus) ) 2023 { 2024 std::ostringstream message; 2025 message << "May lead to a crash or unreliable results." << G4endl 2026 << " Position has shifted considerably without" 2027 << " notifying the navigator !" << G4endl 2028 << " Tolerated safety: " << safetyPlus << G4endl 2029 << " Computed shift : " << shiftOriginSafSq; 2030 G4Exception("G4ITNavigator1::ComputeStep()", "GeomNav1002", 2031 JustWarning, message); 2032 } 2033 } 2034 2035 // ******************************************************************** 2036 // Operator << 2037 // ******************************************************************** 2038 // 2039 std::ostream& operator << (std::ostream &os,const G4ITNavigator1 &n) 2040 { 2041 // Old version did only the following: 2042 // os << "Current History: " << G4endl << n.fHistory; 2043 // Old behaviour is recovered for fVerbose = 0 2044 2045 // Adapted from G4ITNavigator1::PrintState() const 2046 2047 G4long oldcoutPrec = os.precision(4); 2048 if( n.fVerbose >= 4 ) 2049 { 2050 os << "The current state of G4ITNavigator1 is: " << G4endl; 2051 os << " ValidExitNormal= " << n.fValidExitNormal << G4endl 2052 << " ExitNormal = " << n.fExitNormal << G4endl 2053 << " Exiting = " << n.fExiting << G4endl 2054 << " Entering = " << n.fEntering << G4endl 2055 << " BlockedPhysicalVolume= " ; 2056 if (n.fBlockedPhysicalVolume==nullptr) 2057 os << "None"; 2058 else 2059 os << n.fBlockedPhysicalVolume->GetName(); 2060 os << G4endl 2061 << " BlockedReplicaNo = " << n.fBlockedReplicaNo << G4endl 2062 << " LastStepWasZero = " << n.fLastStepWasZero << G4endl 2063 << G4endl; 2064 } 2065 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) ) 2066 { 2067 os << G4endl; // Make sure to line up 2068 os << std::setw(30) << " ExitNormal " << " " 2069 << std::setw( 5) << " Valid " << " " 2070 << std::setw( 9) << " Exiting " << " " 2071 << std::setw( 9) << " Entering" << " " 2072 << std::setw(15) << " Blocked:Volume " << " " 2073 << std::setw( 9) << " ReplicaNo" << " " 2074 << std::setw( 8) << " LastStepZero " << " " 2075 << G4endl; 2076 os << "( " << std::setw(7) << n.fExitNormal.x() 2077 << ", " << std::setw(7) << n.fExitNormal.y() 2078 << ", " << std::setw(7) << n.fExitNormal.z() << " ) " 2079 << std::setw( 5) << n.fValidExitNormal << " " 2080 << std::setw( 9) << n.fExiting << " " 2081 << std::setw( 9) << n.fEntering << " "; 2082 if ( n.fBlockedPhysicalVolume==nullptr ) 2083 { os << std::setw(15) << "None"; } 2084 else 2085 { os << std::setw(15)<< n.fBlockedPhysicalVolume->GetName(); } 2086 os << std::setw( 9) << n.fBlockedReplicaNo << " " 2087 << std::setw( 8) << n.fLastStepWasZero << " " 2088 << G4endl; 2089 } 2090 if( n.fVerbose > 2 ) 2091 { 2092 os.precision(8); 2093 os << " Current Localpoint = " << n.fLastLocatedPointLocal << G4endl; 2094 os << " PreviousSftOrigin = " << n.fPreviousSftOrigin << G4endl; 2095 os << " PreviousSafety = " << n.fPreviousSafety << G4endl; 2096 } 2097 if( n.fVerbose > 3 || n.fVerbose == 0 ) 2098 { 2099 os << "Current History: " << G4endl << n.fHistory; 2100 } 2101 2102 os.precision(oldcoutPrec); 2103 return os; 2104 } 2105