Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // class G4ITNavigator2 Implementation 28 // 29 // Original author: Paul Kent, July 95/96 30 // 31 // G4ITNavigator2 is a duplicate version of G4 32 // initially written by Paul Kent and colleagu 33 // The only difference resides in the way the 34 // 35 // History: 36 // - Created. 37 // - Zero step protections 38 // - Added check mode 39 // - Made Navigator Abstract 40 // - G4ITNavigator2 created 41 // ------------------------------------------- 42 43 #include <iomanip> 44 45 //#include "G4ITNavigator2.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 #include "Randomize.hh" 55 #include "G4VoxelLimits.hh" 56 57 #define fHistory fpNavigatorState->fHistory 58 #define fLastTriedStepComputation fpNavigatorS 59 #define fEnteredDaughter fpNavigatorState->fEn 60 #define fExitedMother fpNavigatorState->fExite 61 #define fWasLimitedByGeometry fpNavigatorState 62 #define fStepEndPoint fpNavigatorState->fStepE 63 #define fLastStepEndPointLocal fpNavigatorStat 64 #define fPushed fpNavigatorState->fPushed 65 #define fLastTriedStepComputation fpNavigatorS 66 #define fEntering fpNavigatorState->fEntering 67 #define fExiting fpNavigatorState->fExiting 68 #define fBlockedPhysicalVolume fpNavigatorStat 69 #define fBlockedReplicaNo fpNavigatorState->fB 70 #define fLastLocatedPointLocal fpNavigatorStat 71 #define fLocatedOutsideWorld fpNavigatorState- 72 #define fValidExitNormal fpNavigatorState->fVa 73 #define fExitNormal fpNavigatorState->fExitNor 74 #define fGrandMotherExitNormal fpNavigatorStat 75 #define fChangedGrandMotherRefFrame fpNavigato 76 #define fExitNormalGlobalFrame fpNavigatorStat 77 #define fCalculatedExitNormal fpNavigatorState 78 #define fLastStepWasZero fpNavigatorState->fLa 79 #define fLocatedOnEdge fpNavigatorState->fLoca 80 #define fNumberZeroSteps fpNavigatorState->fNu 81 #define fPreviousSftOrigin fpNavigatorState->f 82 #define fPreviousSafety fpNavigatorState->fPre 83 84 //-------------------------------------------- 85 86 namespace BoundingBox 87 { 88 enum Boundary 89 { 90 kMin, 91 kMax 92 }; 93 } 94 95 // ******************************************* 96 // Constructor 97 // ******************************************* 98 // 99 G4ITNavigator2::G4ITNavigator2() 100 { 101 fActive= false; 102 fActionThreshold_NoZeroSteps = 1000; 103 fAbandonThreshold_NoZeroSteps = 2500; 104 kCarTolerance = G4GeometryTolerance::GetInst 105 fregularNav.SetNormalNavigation( &fnormalNav 106 107 fpNavigatorState = nullptr; 108 // fSaveState = 0; 109 fpVoxelSafety= new G4VoxelSafety(); 110 111 // this->SetVerboseLevel(3); 112 // this->CheckMode(true); 113 } 114 115 // ******************************************* 116 // Destructor 117 // ******************************************* 118 // 119 G4ITNavigator2::~G4ITNavigator2() 120 { delete fpVoxelSafety; } 121 122 // ******************************************* 123 // ResetHierarchyAndLocate 124 // ******************************************* 125 // 126 G4VPhysicalVolume* 127 G4ITNavigator2::ResetHierarchyAndLocate(const 128 const G4T 129 const G4T 130 { 131 // ResetState(); 132 fHistory = *h.GetHistory(); 133 SetupHierarchy(); 134 fLastTriedStepComputation= false; // Redund 135 return LocateGlobalPointAndSetup(p, &directi 136 } 137 138 // ******************************************* 139 // LocateGlobalPointAndSetup 140 // 141 // Locate the point in the hierarchy return 0 142 // The direction is required 143 // - if on an edge shared by more than two 144 // (to resolve likely looping in tracking 145 // - at initial location of a particle 146 // (to resolve potential ambiguity at bou 147 // 148 // Flags on exit: (comments to be completed) 149 // fEntering - True if entering `daugh 150 // whether daughter of las 151 // or daughter of that vol 152 // ******************************************* 153 // 154 G4VPhysicalVolume* 155 G4ITNavigator2::LocateGlobalPointAndSetup( con 156 const 157 const 158 const 159 { 160 CheckNavigatorStateIsValid(); 161 G4bool notKnownContained=true, noResult; 162 G4VPhysicalVolume *targetPhysical; 163 G4LogicalVolume *targetLogical; 164 G4VSolid *targetSolid=nullptr; 165 G4ThreeVector localPoint, globalDirection; 166 EInside insideCode; 167 168 G4bool considerDirection = (!ignoreDirection 169 170 fLastTriedStepComputation= false; 171 fChangedGrandMotherRefFrame= false; // For 172 173 if( considerDirection && pGlobalDirection != 174 { 175 globalDirection=*pGlobalDirection; 176 } 177 178 179 #ifdef G4VERBOSE 180 if( fVerbose > 2 ) 181 { 182 G4long oldcoutPrec = G4cout.precision(8); 183 G4cout << "*** G4ITNavigator2::LocateGloba 184 G4cout << " Called with arguments: " << 185 << " Globalpoint = " << globalPo 186 << " RelativeSearch = " << relat 187 if( fVerbose == 4 ) 188 { 189 G4cout << " ----- Upon entering:" << 190 PrintState(); 191 } 192 G4cout.precision(oldcoutPrec); 193 } 194 #endif 195 196 G4int noLevelsExited=0 ; 197 198 if ( !relativeSearch ) 199 { 200 ResetStackAndState(); 201 } 202 else 203 { 204 if ( fWasLimitedByGeometry ) 205 { 206 fWasLimitedByGeometry = false; 207 fEnteredDaughter = fEntering; // Remem 208 fExitedMother = fExiting; // Remem 209 if ( fExiting ) 210 { 211 noLevelsExited++; // count this first 212 213 if ( fHistory.GetDepth() ) 214 { 215 fBlockedPhysicalVolume = fHistory.Ge 216 fBlockedReplicaNo = fHistory.GetTopR 217 fHistory.BackLevel(); 218 } 219 else 220 { 221 fLastLocatedPointLocal = localPoint; 222 fLocatedOutsideWorld = true; 223 fBlockedPhysicalVolume = nullptr; 224 fBlockedReplicaNo = -1; 225 fEntering = false; // No 226 fEnteredDaughter = false; 227 fExitedMother = true; // ?? 228 229 return nullptr; // Have ex 230 } 231 // A fix for the case where a volume i 232 // and a coincident surface exists out 233 // - This stops it from exiting furth 234 // - However ReplicaNavigator treats 235 // 236 // assert( fBlockedPhysicalVolume!=0 ) 237 238 // Expect to be on edge => on surface 239 // 240 if ( fLocatedOnEdge && (VolumeType(fBl 241 { 242 fExiting= false; 243 // Consider effect on Exit Normal !? 244 } 245 } 246 else 247 if ( fEntering ) 248 { 249 // assert( fBlockedPhysicalVolume!=0 250 251 switch (VolumeType(fBlockedPhysicalV 252 { 253 case kNormal: 254 fHistory.NewLevel(fBlockedPhysic 255 fBlockedPhysic 256 break; 257 case kReplica: 258 freplicaNav.ComputeTransformatio 259 260 fHistory.NewLevel(fBlockedPhysic 261 fBlockedReplic 262 fBlockedPhysicalVolume->SetCopyN 263 break; 264 case kParameterised: 265 if( fBlockedPhysicalVolume->GetR 266 { 267 G4VSolid *pSolid; 268 G4VPVParameterisation *pParam; 269 G4TouchableHistory parentTouch 270 pParam = fBlockedPhysicalVolum 271 pSolid = pParam->ComputeSolid( 272 273 pSolid->ComputeDimensions(pPar 274 fBlo 275 pParam->ComputeTransformation( 276 277 fHistory.NewLevel(fBlockedPhys 278 fBlockedRepl 279 fBlockedPhysicalVolume->SetCop 280 // 281 // Set the correct solid and m 282 // 283 G4LogicalVolume *pLogical; 284 pLogical = fBlockedPhysicalVol 285 pLogical->SetSolid( pSolid ); 286 pLogical->UpdateMaterial(pPara 287 ComputeMaterial(fBlockedRepl 288 fBlockedPhys 289 &parentTouch 290 } 291 break; 292 case kExternal: 293 G4Exception("G4ITNavigator2::Loca 294 "GeomNav0001", FatalE 295 "Not applicable for e 296 break; 297 } 298 fEntering = false; 299 fBlockedPhysicalVolume = nullptr; 300 localPoint = fHistory.GetTopTransfor 301 notKnownContained = false; 302 } 303 } 304 else 305 { 306 fBlockedPhysicalVolume = nullptr; 307 fEntering = false; 308 fEnteredDaughter = false; // Full Step 309 fExiting = false; 310 fExitedMother = false; // Full Step 311 } 312 } 313 // 314 // Search from top of history up through geo 315 // containing volume found: 316 // If on 317 // o OUTSIDE - Back up level, not/no longer 318 // o SURFACE and EXITING - Back up level, se 319 // else 320 // o containing volume found 321 // 322 323 while (notKnownContained) 324 { 325 if ( fHistory.GetTopVolumeType()!=kReplica 326 { 327 targetSolid = fHistory.GetTopVolume()->G 328 localPoint = fHistory.GetTopTransform(). 329 insideCode = targetSolid->Inside(localPo 330 #ifdef G4VERBOSE 331 if(( fVerbose == 1 ) && ( fCheck )) 332 { 333 G4String solidResponse = "-kInside-"; 334 if (insideCode == kOutside) 335 solidResponse = "-kOutside-"; 336 else if (insideCode == kSurface) 337 solidResponse = "-kSurface-"; 338 G4cout << "*** G4ITNavigator2::Locate 339 << " Invoked Inside() for s 340 << ". Solid replied: " << soli 341 << " For local point p: " < 342 } 343 #endif 344 } 345 else 346 { 347 insideCode = freplicaNav.BackLocate(fHis 348 fExi 349 // !CARE! if notKnownContained returns f 350 // the containing placement volume of th 351 // will result in the history being back 352 // local point returned is the point in 353 } 354 355 356 if ( insideCode==kOutside ) 357 { 358 noLevelsExited++; 359 if ( fHistory.GetDepth() ) 360 { 361 fBlockedPhysicalVolume = fHistory.GetT 362 fBlockedReplicaNo = fHistory.GetTopRep 363 fHistory.BackLevel(); 364 fExiting = false; 365 366 if( noLevelsExited > 1 ) 367 { 368 // The first transformation was done 369 // 370 const G4RotationMatrix* mRot = fBloc 371 if( mRot != nullptr ) 372 { 373 fGrandMotherExitNormal *= (*mRot). 374 fChangedGrandMotherRefFrame= true; 375 } 376 } 377 } 378 else 379 { 380 fLastLocatedPointLocal = localPoint; 381 fLocatedOutsideWorld = true; 382 // No extra transformation for ExitN 383 return nullptr; // Have exited 384 } 385 } 386 else 387 if ( insideCode==kSurface ) 388 { 389 G4bool isExiting = fExiting; 390 if( (!fExiting)&&considerDirection ) 391 { 392 // Figure out whether we are exiting 393 // by using the direction 394 // 395 G4bool directionExiting = false; 396 G4ThreeVector localDirection = 397 fHistory.GetTopTransform().Trans 398 399 // Make sure localPoint in correct r 400 // ( Was it already correct ? Ho 401 // 402 localPoint= fHistory.GetTopTransform 403 if ( fHistory.GetTopVolumeType()!=kR 404 { 405 G4ThreeVector normal = targetSolid 406 directionExiting = normal.dot(loca 407 isExiting = isExiting || direction 408 } 409 } 410 if( isExiting ) 411 { 412 noLevelsExited++; 413 if ( fHistory.GetDepth() ) 414 { 415 fBlockedPhysicalVolume = fHistory. 416 fBlockedReplicaNo = fHistory.GetTo 417 fHistory.BackLevel(); 418 // 419 // Still on surface but exited vol 420 // 421 fValidExitNormal = false; 422 423 if( noLevelsExited > 1 ) 424 { 425 // The first transformation was 426 // 427 const G4RotationMatrix* mRot = 428 fBlockedPhysicalVolume->Ge 429 if( mRot != nullptr ) 430 { 431 fGrandMotherExitNormal *= (*mR 432 fChangedGrandMotherRefFrame= t 433 } 434 } 435 } 436 else 437 { 438 fLastLocatedPointLocal = localPoin 439 fLocatedOutsideWorld = true; 440 // No extra transformation for E 441 return nullptr; // Have e 442 } 443 } 444 else 445 { 446 notKnownContained=false; 447 } 448 } 449 else 450 { 451 notKnownContained=false; 452 } 453 } // END while (notKnownContained) 454 // 455 // Search downwards until deepest containing 456 // blocking fBlockedPhysicalVolume/BlockedRe 457 // 458 // 3 Cases: 459 // 460 // o Parameterised daughters 461 // =>Must be one G4PVParameterised daughte 462 // o Positioned daughters & voxels 463 // o Positioned daughters & no voxels 464 465 noResult = true; // noResult should be rena 466 // something like enteredL 467 do 468 { 469 // Determine `type' of current mother volu 470 // 471 targetPhysical = fHistory.GetTopVolume(); 472 if (targetPhysical == nullptr) { break; } 473 targetLogical = targetPhysical->GetLogical 474 switch( CharacteriseDaughters(targetLogica 475 { 476 case kNormal: 477 if ( targetLogical->GetVoxelHeader() ! 478 { 479 noResult = fvoxelNav.LevelLocate(fHi 480 fBl 481 fBl 482 glo 483 pGl 484 con 485 loc 486 } 487 else // do not u 488 { 489 noResult = fnormalNav.LevelLocate(fH 490 fB 491 fB 492 gl 493 pG 494 co 495 lo 496 } 497 break; 498 case kReplica: 499 noResult = freplicaNav.LevelLocate(fHi 500 fBl 501 fBl 502 glo 503 pGl 504 con 505 loc 506 break; 507 case kParameterised: 508 if( GetDaughtersRegularStructureId(tar 509 { 510 noResult = fparamNav.LevelLocate(fHi 511 fBl 512 fBl 513 glo 514 pGl 515 con 516 loc 517 } 518 else // Regular structure 519 { 520 noResult = fregularNav.LevelLocate(f 521 f 522 f 523 g 524 p 525 c 526 l 527 } 528 break; 529 case kExternal: 530 G4Exception("G4ITNavigator2::LocateGlo 531 "GeomNav0001", FatalExcept 532 "Not applicable for extern 533 break; 534 } 535 536 // LevelLocate returns true if it finds a 537 // in which globalPoint is inside (or on t 538 539 if ( noResult ) 540 { 541 // Entering a daughter after ascending 542 // 543 // The blocked volume is no longer valid 544 // 545 fBlockedPhysicalVolume = nullptr; 546 fBlockedReplicaNo = -1; 547 548 // fEntering should be false -- else blo 549 // fEnteredDaughter is used for ExitNorm 550 // 551 fEntering = false; 552 fEnteredDaughter = true; 553 554 if( fExitedMother ) 555 { 556 G4VPhysicalVolume* enteredPhysical = f 557 const G4RotationMatrix* mRot = entered 558 if( mRot != nullptr ) 559 { 560 // Go deeper, i.e. move 'down' in th 561 // Apply direct rotation, not invers 562 // 563 fGrandMotherExitNormal *= (*mRot); 564 fChangedGrandMotherRefFrame= true; 565 } 566 } 567 568 #ifdef G4DEBUG_NAVIGATION 569 if( fVerbose > 2 ) 570 { 571 G4VPhysicalVolume* enteredPhysical = 572 G4cout << "*** G4ITNavigator2::Locate 573 G4cout << " Entering volume: " << 574 << G4endl; 575 } 576 #endif 577 } 578 } while (noResult); 579 580 fLastLocatedPointLocal = localPoint; 581 582 #ifdef G4VERBOSE 583 if( fVerbose >= 4 ) 584 { 585 G4long oldcoutPrec = G4cout.precision(8); 586 G4String curPhysVol_Name("None"); 587 if (targetPhysical != nullptr) { curPhysV 588 G4cout << " Return value = new volume = 589 G4cout << " ----- Upon exiting:" << G4e 590 PrintState(); 591 if( fVerbose >= 5 ) 592 { 593 G4cout << "Upon exiting LocateGlobalPoin 594 G4cout << " History = " << G4endl << 595 } 596 G4cout.precision(oldcoutPrec); 597 } 598 #endif 599 600 fLocatedOutsideWorld= false; 601 602 return targetPhysical; 603 } 604 605 // ******************************************* 606 // LocateGlobalPointWithinVolume 607 // 608 // -> the state information of this Navigator 609 // is updated in order to start the next st 610 // -> no check is performed whether pGlobalpoi 611 // original volume (this must be the case). 612 // 613 // Note: a direction could be added to the arg 614 // optional checking (via the old code b 615 // [ This would be done only in verbose 616 // ******************************************* 617 // 618 void 619 G4ITNavigator2::LocateGlobalPointWithinVolume( 620 { 621 CheckNavigatorStateIsValid(); 622 623 #ifdef G4DEBUG_NAVIGATION 624 // Check: Either step was not limited by a 625 // or else the full step is no long 626 assert( !fWasLimitedByGeometry ); 627 #endif 628 629 fLastLocatedPointLocal = ComputeLocalPoint( 630 fLastTriedStepComputation= false; 631 fChangedGrandMotherRefFrame= false; // Fr 632 633 #ifdef G4DEBUG_NAVIGATION 634 if( fVerbose > 2 ) 635 { 636 G4cout << "Entering LocateGlobalWithinVol 637 G4cout << fHistory << G4endl; 638 } 639 #endif 640 641 // For the case of Voxel (or Parameterised) 642 // Navigator must be messaged to update its 643 644 // Update the state of the Sub Navigators 645 // - in particular any voxel information th 646 // 647 G4VPhysicalVolume* motherPhysical = fHisto 648 G4LogicalVolume* motherLogical = mother 649 G4SmartVoxelHeader* pVoxelHeader = mother 650 651 if ( fHistory.GetTopVolumeType()!=kReplica 652 { 653 switch( CharacteriseDaughters(motherLogic 654 { 655 case kNormal: 656 if ( pVoxelHeader != nullptr ) 657 { 658 fvoxelNav.VoxelLocate( pVoxelHeader 659 } 660 break; 661 case kParameterised: 662 if( GetDaughtersRegularStructureId(mo 663 { 664 // Resets state & returns voxel nod 665 // 666 fparamNav.ParamVoxelLocate( pVoxelH 667 } 668 break; 669 case kReplica: 670 G4Exception("G4ITNavigator2::LocateGl 671 "GeomNav0001", FatalExcep 672 "Not applicable for repli 673 break; 674 case kExternal: 675 G4Exception("G4ITNavigator2::LocateGl 676 "GeomNav0001", FatalExcep 677 "Not applicable for exter 678 break; 679 } 680 } 681 682 // Reset the state variables 683 // - which would have been affected 684 // by the 'equivalent' call to LocateGl 685 // - who's values have been invalidated b 686 // 687 fBlockedPhysicalVolume = nullptr; 688 fBlockedReplicaNo = -1; 689 fEntering = false; 690 fEnteredDaughter = false; // Boundary not 691 fExiting = false; 692 fExitedMother = false; // Boundary not 693 } 694 695 // !> 696 697 void G4ITNavigator2::CheckNavigatorState() con 698 { 699 if(fpNavigatorState == nullptr) 700 { 701 G4ExceptionDescription exceptionDescri 702 exceptionDescription << "The navigator 703 exceptionDescription << "Either NewNav 704 exceptionDescription << "or the provid 705 706 G4Exception("G4ITNavigator::CheckNavig 707 "NavigatorStateNotValid",F 708 return; 709 } 710 } 711 712 G4ITNavigatorState_Lock2* G4ITNavigator2::GetN 713 { 714 return fpNavigatorState; 715 } 716 717 void G4ITNavigator2::SetNavigatorState(G4ITNav 718 { 719 fpNavigatorState = (G4NavigatorState*) nav 720 if(fpNavigatorState != nullptr) SetupHiera 721 // fpSaveState = (G4SaveNavigatorState*) na 722 // if(navState) RestoreSavedState(); 723 } 724 725 void G4ITNavigator2::ResetNavigatorState() 726 { 727 fpNavigatorState = nullptr; 728 } 729 730 void G4ITNavigator2::NewNavigatorState() 731 { 732 fpNavigatorState = new G4NavigatorState(); 733 if(fTopPhysical == nullptr) 734 { 735 G4ExceptionDescription exceptionDescri 736 exceptionDescription << "No World Volu 737 738 G4Exception("G4ITNavigator::NewNavigat 739 "NoWorldVolume",FatalExcep 740 return; 741 } 742 743 fHistory.SetFirstEntry(fTopPhysical ); 744 SetupHierarchy(); 745 } 746 747 void G4ITNavigator2::NewNavigatorState(const G 748 { 749 fpNavigatorState = new G4NavigatorState(); 750 if(fTopPhysical == nullptr) 751 { 752 G4ExceptionDescription exceptionDescri 753 exceptionDescription << "No World Volu 754 755 G4Exception("G4ITNavigator::NewNavigat 756 "NoWorldVolume",FatalExcep 757 return; 758 } 759 760 fHistory = *h.GetHistory(); 761 fLastTriedStepComputation= false; // Redu 762 SetupHierarchy(); 763 } 764 765 G4VPhysicalVolume* G4ITNavigator2::NewNavigato 766 767 { 768 fpNavigatorState = new G4NavigatorState(); 769 770 if(fTopPhysical == nullptr) 771 { 772 G4ExceptionDescription exceptionDescri 773 exceptionDescription << "No World Volu 774 775 G4Exception("G4ITNavigator::NewNavigat 776 "NoWorldVolume",FatalExcep 777 return nullptr; 778 } 779 780 fHistory.SetFirstEntry(fTopPhysical ); 781 SetupHierarchy(); 782 return LocateGlobalPointAndSetup(p, &direc 783 } 784 785 // ******************************************* 786 // SetSavedState 787 // 788 // Save the state, in case this is a parasitic 789 // Save fValidExitNormal, fExitNormal, fExitin 790 // fBlockedPhysicalVolume, fBlockedReplic 791 // ******************************************* 792 // 793 /* 794 void G4ITNavigator2::SetSavedState() 795 { 796 fSaveState = *fpNavigatorState; 797 } 798 */ 799 800 /* 801 void G4ITNavigator2::SetSavedState() 802 { 803 // !> 804 // This check can be avoid if instead, at 805 // the IT tracking uses NewNavigatorSate 806 // The normal tracking would just call onc 807 808 // if(fpSaveState == 0) 809 // fpSaveState = new G4SaveNavigatorSta 810 // <! 811 812 // fSaveExitNormal = fExitNormal; 813 fpSaveState->sExitNormal = fExitNormal; 814 fpSaveState->sValidExitNormal = fValidExitNo 815 fpSaveState->sExiting = fExiting; 816 fpSaveState->sEntering = fEntering; 817 818 fpSaveState->spBlockedPhysicalVolume = fBloc 819 fpSaveState->sBlockedReplicaNo = fBlockedRep 820 821 fpSaveState->sLastStepWasZero = fLastStepWas 822 823 // !> 824 fpSaveState->sPreviousSftOrigin = fPreviousS 825 fpSaveState->sPreviousSafety = fPreviousSafe 826 fpSaveState->sNumberZeroSteps = fNumberZeroS 827 fpSaveState->sLocatedOnEdge = fLocatedOnEdge 828 fpSaveState->sWasLimitedByGeometry= fWasLimi 829 fpSaveState->sPushed=fPushed; 830 fpSaveState->sNumberZeroSteps=fNumberZeroSte 831 fpSaveState->sEnteredDaughter = fEnteredDaug 832 fpSaveState->sExitedMother = fExitedMother; 833 834 fpSaveState->sLastLocatedPointLocal = fLastL 835 fpSaveState->sLocatedOutsideWorld = fLocated 836 // <! 837 } 838 */ 839 // ******************************************* 840 // RestoreSavedState 841 // 842 // Restore the state (in Compute Step), in cas 843 // ******************************************* 844 // 845 /* 846 void G4ITNavigator2::RestoreSavedState() 847 { 848 *fpNavigatorState = fSaveState; 849 } 850 */ 851 /* 852 void G4ITNavigator2::RestoreSavedState() 853 { 854 fExitNormal = fpSaveState->sExitNormal; 855 fValidExitNormal = fpSaveState->sValidExitNo 856 fExiting = fpSaveState->sExiting; 857 fEntering = fpSaveState->sEntering; 858 859 fBlockedPhysicalVolume = fpSaveState->spBloc 860 fBlockedReplicaNo = fpSaveState->sBlockedRep 861 862 fLastStepWasZero = fpSaveState->sLastStepWas 863 864 // !> 865 fPreviousSftOrigin = fpSaveState->sPreviousS 866 fPreviousSafety = fpSaveState->sPreviousSafe 867 fNumberZeroSteps = fpSaveState->sNumberZeroS 868 fLocatedOnEdge = fpSaveState->sLocatedOnEdge 869 fWasLimitedByGeometry = fpSaveState->sWasLim 870 fPushed = fpSaveState->sPushed; 871 fNumberZeroSteps = fpSaveState->sNumberZeroS 872 fEnteredDaughter= fpSaveState->sEnteredDaugh 873 fExitedMother = fpSaveState->sExitedMother ; 874 875 fLastLocatedPointLocal = fpSaveState->sLastL 876 fLocatedOutsideWorld = fpSaveState->sLocated 877 // <! 878 } 879 */ 880 // <! 881 882 // ******************************************* 883 // ComputeStep 884 // 885 // Computes the next geometric Step: intersect 886 // mother and `daughter' volumes. 887 // 888 // NOTE: 889 // 890 // Flags on entry: 891 // -------------- 892 // fValidExitNormal - Normal of exited volume 893 // coincident boundary) 894 // fExitNormal - Surface normal of exite 895 // fExiting - True if have exited sol 896 // 897 // fBlockedPhysicalVolume - Ptr to exited volu 898 // fBlockedReplicaNo - Replication no of exite 899 // fLastStepWasZero - True if last Step size 900 // 901 // Flags on exit: 902 // ------------- 903 // fValidExitNormal - True if surface normal 904 // fExitNormal - Surface normal of exite 905 // reference system 906 // fExiting - True if exiting mother 907 // fEntering - True if entering `daugh 908 // fBlockedPhysicalVolume - Ptr to candidate ( 909 // fBlockedReplicaNo - Replication no of candi 910 // fLastStepWasZero - True if this Step size 911 // ******************************************* 912 // 913 G4double G4ITNavigator2::ComputeStep( const G4 914 const G4Thr 915 const G4dou 916 G4dou 917 { 918 CheckNavigatorStateIsValid(); 919 G4ThreeVector localDirection = ComputeLocalA 920 G4double Step = kInfinity; 921 G4VPhysicalVolume *motherPhysical = fHistor 922 G4LogicalVolume *motherLogical = motherPhysi 923 924 // All state relating to exiting normals mus 925 // 926 fExitNormalGlobalFrame= G4ThreeVector( 0., 0 927 // Reset value - to erase its memory 928 fChangedGrandMotherRefFrame= false; 929 // Reset - used for local exit normal 930 fGrandMotherExitNormal= G4ThreeVector( 0., 0 931 fCalculatedExitNormal = false; 932 // Reset for new step 933 934 fLastTriedStepComputation= true; 935 936 #ifdef G4VERBOSE 937 if( fVerbose > 0 ) 938 { 939 G4cout << "*** G4ITNavigator2::ComputeStep 940 G4cout << " Volume = " << motherPhysica 941 << " - Proposed step length = " << 942 << G4endl; 943 #ifdef G4DEBUG_NAVIGATION 944 if( fVerbose >= 2 ) 945 { 946 G4cout << " Called with the arguments: 947 << " Globalpoint = " << std::set 948 << " Direction = " << std::set 949 if( fVerbose >= 4 ) 950 { 951 G4cout << " ---- Upon entering : Stat 952 PrintState(); 953 } 954 } 955 #endif 956 } 957 #endif 958 959 G4ThreeVector newLocalPoint = ComputeLocalPo 960 if( newLocalPoint != fLastLocatedPointLocal 961 { 962 // Check whether the relocation is within 963 // 964 G4ThreeVector oldLocalPoint = fLastLocated 965 G4double moveLenSq = (newLocalPoint-oldLoc 966 967 if ( moveLenSq >= kCarTolerance*kCarTolera 968 { 969 #ifdef G4VERBOSE 970 ComputeStepLog(pGlobalpoint, moveLenSq); 971 #endif 972 // Relocate the point within the same vo 973 // 974 LocateGlobalPointWithinVolume( pGlobalpo 975 fLastTriedStepComputation= true; // 976 } 977 } 978 if ( fHistory.GetTopVolumeType()!=kReplica ) 979 { 980 switch( CharacteriseDaughters(motherLogica 981 { 982 case kNormal: 983 if ( motherLogical->GetVoxelHeader() ! 984 { 985 LocateGlobalPointWithinVolume(pGloba 986 Step = fvoxelNav.ComputeStep(fLastLo 987 localDi 988 pCurren 989 pNewSaf 990 fHistor 991 fValidE 992 fExitNo 993 fExitin 994 fEnteri 995 &fBlock 996 fBlocke 997 998 } 999 else 1000 { 1001 if( motherPhysical->GetRegularStruc 1002 { 1003 Step = fnormalNav.ComputeStep(fLa 1004 loc 1005 pCu 1006 pNe 1007 fHi 1008 fVa 1009 fEx 1010 fEx 1011 fEn 1012 &fB 1013 fBl 1014 } 1015 else // Regular (non-voxelised) st 1016 { 1017 LocateGlobalPointAndSetup( pGloba 1018 fLastTriedStepComputation= true; 1019 // 1020 // if physical process limits the 1021 // one given by ComputeStepSkippi 1022 // point will be wrongly calculat 1023 1024 // There is a problem: when msc l 1025 // assigned wrongly to phantom in 1026 // of the container volume). Then 1027 // reset the history topvolume to 1028 // 1029 if(fHistory.GetTopVolume()->GetRe 1030 { 1031 G4Exception("G4ITNavigator2::Co 1032 "GeomNav1001", Just 1033 "Point is relocated in voxels 1034 Step = fnormalNav.ComputeStep(f 1035 l 1036 p 1037 p 1038 f 1039 f 1040 f 1041 f 1042 f 1043 & 1044 f 1045 } 1046 else 1047 { 1048 Step = fregularNav. 1049 ComputeStepSkippingEqualMa 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 } 1062 } 1063 } 1064 break; 1065 case kParameterised: 1066 if( GetDaughtersRegularStructureId(mo 1067 { 1068 Step = fparamNav.ComputeStep(fLastL 1069 localD 1070 pCurre 1071 pNewSa 1072 fHisto 1073 fValid 1074 fExitN 1075 fExiti 1076 fEnter 1077 &fBloc 1078 fBlock 1079 } 1080 else // Regular structure 1081 { 1082 Step = fregularNav.ComputeStep(fLas 1083 loca 1084 pCur 1085 pNew 1086 fHis 1087 fVal 1088 fExi 1089 fExi 1090 fEnt 1091 &fBl 1092 fBlo 1093 } 1094 break; 1095 case kReplica: 1096 G4Exception("G4ITNavigator2::ComputeS 1097 FatalException, "Not appl 1098 break; 1099 case kExternal: 1100 G4Exception("G4ITNavigator2::ComputeS 1101 FatalException, "Not appl 1102 break; 1103 } 1104 } 1105 else 1106 { 1107 // In the case of a replica, it must hand 1108 // edge/corner problem by itself 1109 // 1110 G4bool exitingReplica = fExitedMother; 1111 G4bool calculatedExitNormal; 1112 Step = freplicaNav.ComputeStep(pGlobalpoi 1113 pDirection 1114 fLastLocat 1115 localDirec 1116 pCurrentPr 1117 pNewSafety 1118 fHistory, 1119 fValidExit 1120 calculated 1121 fExitNorma 1122 exitingRep 1123 fEntering, 1124 &fBlockedP 1125 fBlockedRe 1126 fExiting= exitingReplica; 1127 fCalculatedExitNormal= calculatedExitNorm 1128 } 1129 1130 1131 // G4cout << " !!!! Step = " << Step << G4en 1132 1133 // Remember last safety origin & value. 1134 // 1135 fPreviousSftOrigin = pGlobalpoint; 1136 fPreviousSafety = pNewSafety; 1137 1138 // Count zero steps - one can occur due to 1139 // - one, two (or a few) c 1140 // volumes 1141 // - more than two is like 1142 // description or the Na 1143 1144 // Rule of thumb: likely at an Edge if two 1145 // because at least two cand 1146 // checked 1147 // 1148 fLocatedOnEdge = fLastStepWasZero && (Ste 1149 fLastStepWasZero = (Step==0.0); 1150 if (fPushed) { fPushed = fLastStepWasZero; 1151 1152 // Handle large number of consecutive zero 1153 // 1154 if ( fLastStepWasZero ) 1155 { 1156 fNumberZeroSteps++; 1157 #ifdef G4DEBUG_NAVIGATION 1158 if( fNumberZeroSteps > 1 ) 1159 { 1160 G4cout << "G4ITNavigator2::ComputeStep 1161 << fNumberZeroSteps 1162 << " at " << pGlobalpoint 1163 << " in volume " << motherPhysi 1164 << G4endl; 1165 } 1166 #endif 1167 if( fNumberZeroSteps > fActionThreshold_N 1168 { 1169 // Act to recover this stuck track. Pu 1170 // 1171 Step += 100*kCarTolerance; 1172 #ifdef G4VERBOSE 1173 if ((!fPushed) && (fWarnPush)) 1174 { 1175 std::ostringstream message; 1176 message << "Track stuck or not movin 1177 << " Track stuck, n 1178 << fNumberZeroSteps << " ste 1179 << " in volume -" < 1180 << "- at point " << pGlobalp 1181 << " direction: " < 1182 << " Potential geom 1183 << G4endl 1184 << " Trying pushing 1185 G4Exception("G4ITNavigator2::Compute 1186 JustWarning, message, "P 1187 } 1188 #endif 1189 fPushed = true; 1190 } 1191 if( fNumberZeroSteps > fAbandonThreshold_ 1192 { 1193 // Must kill this stuck track 1194 // 1195 std::ostringstream message; 1196 message << "Stuck Track: potential geom 1197 << G4endl 1198 << " Track stuck, not mo 1199 << fNumberZeroSteps << " steps" 1200 << " in volume -" << mot 1201 << "- at point " << pGlobalpoin 1202 << " direction: " << pDi 1203 motherPhysical->CheckOverlaps(5000, 0.0 1204 G4Exception("G4ITNavigator2::ComputeSte 1205 EventMustBeAborted, message 1206 } 1207 } 1208 else 1209 { 1210 if (!fPushed) fNumberZeroSteps = 0; 1211 } 1212 1213 fEnteredDaughter = fEntering; // I expect 1214 fExitedMother = fExiting; 1215 1216 fStepEndPoint = pGlobalpoint 1217 + std::min(Step,pCurrentPropo 1218 fLastStepEndPointLocal = fLastLocatedPointL 1219 1220 if( fExiting ) 1221 { 1222 #ifdef G4DEBUG_NAVIGATION 1223 if( fVerbose > 2 ) 1224 { 1225 G4cout << " At G4Nav CompStep End - if( 1226 << " fValidExitNormal = " << fVa 1227 G4cout << " fExitNormal= " << fExitNorm 1228 } 1229 #endif 1230 1231 if(fValidExitNormal || fCalculatedExitNor 1232 { 1233 if ( fHistory.GetTopVolumeType()!=kRep 1234 { 1235 // Convention: fExitNormal is in the 1236 // 1237 fGrandMotherExitNormal= fExitNormal; 1238 fCalculatedExitNormal= true; 1239 } 1240 else 1241 { 1242 fGrandMotherExitNormal = fExitNormal; 1243 } 1244 } 1245 else 1246 { 1247 // We must calculate the normal anyway 1248 // 1249 G4ThreeVector finalLocalPoint = 1250 fLastLocatedPointLocal + localDir 1251 1252 if ( fHistory.GetTopVolumeType()!=kRep 1253 { 1254 // Find normal in the 'mother' coordi 1255 // 1256 G4ThreeVector exitNormalMotherFrame= 1257 motherLogical->GetSolid()->Surface 1258 1259 // Transform it to the 'grand-mother' 1260 // 1261 const G4RotationMatrix* mRot = mother 1262 if( mRot != nullptr ) 1263 { 1264 fChangedGrandMotherRefFrame= true; 1265 fGrandMotherExitNormal = (*mRot).in 1266 } 1267 else 1268 { 1269 fGrandMotherExitNormal = exitNormal 1270 } 1271 1272 // Do not set fValidExitNormal -- thi 1273 // that the solid is convex! 1274 // 1275 fCalculatedExitNormal= true; 1276 } 1277 else 1278 { 1279 fCalculatedExitNormal = false; 1280 // 1281 // Nothing can be done at this stage 1282 // Replica Navigation must have calcu 1283 // already. 1284 // Cases: mother is not convex, and e 1285 1286 #ifdef G4DEBUG_NAVIGATION 1287 G4ExceptionDescription desc; 1288 1289 desc << "Problem in ComputeStep: Rep 1290 << " valid exit Normal. " << G4e 1291 desc << " Do not know how calculate i 1292 desc << " Location = " << finalLo 1293 desc << " Volume name = " << motherP 1294 << " copy/replica No = " << mot 1295 G4Exception("G4ITNavigator2::ComputeS 1296 JustWarning, desc, "Norma 1297 #endif 1298 } 1299 } 1300 1301 // Now transform it to the global referen 1302 // 1303 if( fValidExitNormal || fCalculatedExitNo 1304 { 1305 auto depth= (G4int)fHistory.GetDepth() 1306 if( depth > 0 ) 1307 { 1308 G4AffineTransform GrandMotherToGlobal 1309 fHistory.GetTransform(depth-1).Inve 1310 fExitNormalGlobalFrame = 1311 GrandMotherToGlobalTransf.Transform 1312 } 1313 else 1314 { 1315 fExitNormalGlobalFrame= fGrandMotherE 1316 } 1317 } 1318 else 1319 { 1320 fExitNormalGlobalFrame= G4ThreeVector( 1321 } 1322 } 1323 1324 if( (Step == pCurrentProposedStepLength) && 1325 { 1326 // This if Step is not really limited by 1327 // The Navigator is obliged to return "in 1328 // 1329 Step = kInfinity; 1330 } 1331 1332 #ifdef G4VERBOSE 1333 if( fVerbose > 1 ) 1334 { 1335 if( fVerbose >= 4 ) 1336 { 1337 G4cout << " ----- Upon exiting :" << 1338 PrintState(); 1339 } 1340 G4cout << " Returned step= " << Step; 1341 if( fVerbose > 5 ) G4cout << G4endl; 1342 if( Step == kInfinity ) 1343 { 1344 G4cout << " Requested step= " << pCurr 1345 if( fVerbose > 5) G4cout << G4endl; 1346 } 1347 G4cout << " Safety = " << pNewSafety << 1348 } 1349 #endif 1350 1351 return Step; 1352 } 1353 1354 // ****************************************** 1355 // CheckNextStep 1356 // 1357 // Compute the step without altering the navi 1358 // ****************************************** 1359 // 1360 G4double G4ITNavigator2::CheckNextStep( const 1361 const 1362 const 1363 1364 { 1365 CheckNavigatorStateIsValid(); 1366 G4double step; 1367 1368 // Save the state, for this parasitic call 1369 // 1370 //SetSavedState(); 1371 // G4SaveNavigatorState savedState(fpNavigat 1372 G4NavigatorState savedState(*fpNavigatorSta 1373 1374 step = ComputeStep ( pGlobalpoint, 1375 pDirection, 1376 pCurrentProposedStepLe 1377 pNewSafety ); 1378 1379 // If a parasitic call, then attempt to res 1380 // 1381 *fpNavigatorState = savedState; 1382 //RestoreSavedState(); 1383 // NOTE: the state of the current subnaviga 1384 // ***> TODO: restore subnavigator state 1385 // if( last_located) Need 1386 // if( last_computed step) Need 1387 1388 1389 return step; 1390 } 1391 1392 // ****************************************** 1393 // ResetState 1394 // 1395 // Resets stack and minimum of navigator stat 1396 // ****************************************** 1397 // 1398 void G4ITNavigator2::ResetState() 1399 { 1400 fWasLimitedByGeometry = false; 1401 fEntering = false; 1402 fExiting = false; 1403 fLocatedOnEdge = false; 1404 fLastStepWasZero = false; 1405 fEnteredDaughter = false; 1406 fExitedMother = false; 1407 fPushed = false; 1408 1409 fValidExitNormal = false; 1410 fChangedGrandMotherRefFrame= false; 1411 fCalculatedExitNormal = false; 1412 1413 fExitNormal = G4ThreeVector(0,0, 1414 fGrandMotherExitNormal = G4ThreeVector(0,0, 1415 fExitNormalGlobalFrame = G4ThreeVector(0,0, 1416 1417 fPreviousSftOrigin = G4ThreeVector(0,0, 1418 fPreviousSafety = 0.0; 1419 1420 fNumberZeroSteps = 0; 1421 1422 fBlockedPhysicalVolume = nullptr; 1423 fBlockedReplicaNo = -1; 1424 1425 fLastLocatedPointLocal = G4ThreeVector( kIn 1426 fLocatedOutsideWorld = false; 1427 } 1428 1429 // ****************************************** 1430 // SetupHierarchy 1431 // 1432 // Renavigates & resets hierarchy described b 1433 // o Reset volumes 1434 // o Recompute transforms and/or solids of re 1435 // ****************************************** 1436 // 1437 void G4ITNavigator2::SetupHierarchy() 1438 { 1439 G4int i; 1440 const auto cdepth = (G4int)fHistory.GetDep 1441 G4VPhysicalVolume *current; 1442 G4VSolid *pSolid; 1443 G4VPVParameterisation *pParam; 1444 1445 for ( i=1; i<=cdepth; i++ ) 1446 { 1447 current = fHistory.GetVolume(i); 1448 switch ( fHistory.GetVolumeType(i) ) 1449 { 1450 case kNormal: 1451 break; 1452 case kReplica: 1453 freplicaNav.ComputeTransformation(fHi 1454 break; 1455 case kParameterised: 1456 { 1457 G4int replicaNo; 1458 pParam = current->GetParameterisation 1459 replicaNo = fHistory.GetReplicaNo(i); 1460 pSolid = pParam->ComputeSolid(replica 1461 1462 // Set up dimensions & transform in s 1463 // 1464 pSolid->ComputeDimensions(pParam, rep 1465 pParam->ComputeTransformation(replica 1466 1467 G4TouchableHistory *pTouchable= nullp 1468 if( pParam->IsNested() ) 1469 { 1470 pTouchable= new G4TouchableHistory 1471 pTouchable->MoveUpHistory(); // Mo 1472 // Adequate only if Nested at th 1473 // To extend to other cases: 1474 // pTouchable->MoveUpHistory(cdept 1475 // Move to the parent level of * 1476 // Could replace this line and c 1477 // c-tor for History(levels to d 1478 } 1479 // Set up the correct solid and mater 1480 // 1481 G4LogicalVolume *pLogical = current-> 1482 pLogical->SetSolid( pSolid ); 1483 pLogical->UpdateMaterial( pParam -> 1484 ComputeMaterial(replicaNo, current, 1485 delete pTouchable; 1486 } 1487 break; 1488 1489 case kExternal: 1490 G4Exception("G4ITNavigator2::SetupHie 1491 "GeomNav0001", FatalExcep 1492 "Not applicable for exter 1493 break; 1494 1495 } 1496 } 1497 } 1498 1499 // ****************************************** 1500 // GetLocalExitNormal 1501 // 1502 // Obtains the Normal vector to a surface (in 1503 // pointing out of previous volume and into c 1504 // ****************************************** 1505 // 1506 G4ThreeVector G4ITNavigator2::GetLocalExitNor 1507 { 1508 CheckNavigatorStateIsValid(); 1509 G4ThreeVector ExitNormal(0.,0.,0.); 1510 G4VSolid *currentSolid=nullptr; 1511 G4LogicalVolume *candidateLogical; 1512 if ( fLastTriedStepComputation ) 1513 { 1514 // use fLastLocatedPointLocal and next ca 1515 // 1516 G4ThreeVector nextSolidExitNormal(0.,0.,0 1517 1518 if( fEntering && (fBlockedPhysicalVolume! 1519 { 1520 candidateLogical= fBlockedPhysicalVolum 1521 if( candidateLogical != nullptr ) 1522 { 1523 // fLastStepEndPointLocal is in the c 1524 // we need it in the daughter's coord 1525 1526 // The following code should also wor 1527 { 1528 // First transform fLastLocatedPoin 1529 // coordinates 1530 // 1531 G4AffineTransform MotherToDaughterT 1532 GetMotherToDaughterTransform( fBl 1533 fBl 1534 Vol 1535 G4ThreeVector daughterPointOwnLocal 1536 MotherToDaughterTransform.Transfo 1537 1538 // OK if it is a parameterised volu 1539 // 1540 EInside inSideIt; 1541 G4bool onSurface; 1542 G4double safety= -1.0; 1543 currentSolid= candidateLogical->Get 1544 inSideIt = currentSolid->Inside(d 1545 onSurface = (inSideIt == kSurface) 1546 if( ! onSurface ) 1547 { 1548 if( inSideIt == kOutside ) 1549 { 1550 safety = (currentSolid->Distanc 1551 onSurface = safety < 100.0 * kC 1552 } 1553 else if (inSideIt == kInside ) 1554 { 1555 safety = (currentSolid->Distanc 1556 onSurface = safety < 100.0 * kC 1557 } 1558 } 1559 1560 if( onSurface ) 1561 { 1562 nextSolidExitNormal = 1563 currentSolid->SurfaceNormal(dau 1564 1565 // Entering the solid ==> opposit 1566 // 1567 ExitNormal = -nextSolidExitNormal 1568 fCalculatedExitNormal= true; 1569 } 1570 else 1571 { 1572 #ifdef G4VERBOSE 1573 if(( fVerbose == 1 ) && ( fCheck 1574 { 1575 std::ostringstream message; 1576 message << "Point not on surfac 1577 << " Point = 1578 << daughterPointOwnLoca 1579 << " Physical volume = 1580 << fBlockedPhysicalVolu 1581 << " Logical volume = 1582 << candidateLogical->Ge 1583 << " Solid = 1584 << " Type = 1585 << currentSolid->GetEnt 1586 << *currentSolid << G4e 1587 if( inSideIt == kOutside ) 1588 { 1589 message << "Point is Outside. 1590 << " Safety (from ou 1591 } 1592 else // if( inSideIt == kInside 1593 { 1594 message << "Point is Inside. 1595 << " Safety (from in 1596 } 1597 G4Exception("G4ITNavigator2::Ge 1598 JustWarning, messag 1599 } 1600 #endif 1601 } 1602 *valid = onSurface; // was =tru 1603 } 1604 } 1605 } 1606 else if ( fExiting ) 1607 { 1608 ExitNormal = fGrandMotherExitNormal; 1609 *valid = true; 1610 fCalculatedExitNormal= true; // Should 1611 } 1612 else // i.e. ( fBlockedPhysicalVolume = 1613 { 1614 *valid = false; 1615 G4Exception("G4ITNavigator2::GetLocalEx 1616 "GeomNav0003", JustWarning, 1617 "Incorrect call to GetLocal 1618 } 1619 } 1620 else // ( ! fLastTriedStepComputation ) ie 1621 { 1622 if ( EnteredDaughterVolume() ) 1623 { 1624 G4VSolid* daughterSolid =fHistory.GetTo 1625 1626 ExitNormal= -(daughterSolid->SurfaceNor 1627 if( std::fabs(ExitNormal.mag2()-1.0 ) > 1628 { 1629 G4ExceptionDescription desc; 1630 desc << " Parameters of solid: " << * 1631 << " Point for surface = " << fL 1632 G4Exception("G4ITNavigator2::GetLocal 1633 "GeomNav0003", FatalExcep 1634 "Surface Normal returned 1635 } 1636 fCalculatedExitNormal= true; 1637 *valid = true; 1638 } 1639 else 1640 { 1641 if( fExitedMother ) 1642 { 1643 ExitNormal = fGrandMotherExitNormal; 1644 *valid = true; 1645 fCalculatedExitNormal= true; 1646 } 1647 else // We are not at a boundary. Exit 1648 { 1649 *valid = false; 1650 fCalculatedExitNormal= false; 1651 G4ExceptionDescription message; 1652 message << "Function called when *NOT 1653 G4Exception("G4ITNavigator2::GetLocal 1654 "GeomNav0003", JustWarnin 1655 } 1656 } 1657 } 1658 return ExitNormal; 1659 } 1660 1661 // ****************************************** 1662 // GetMotherToDaughterTransform 1663 // 1664 // Obtains the mother to daughter affine tran 1665 // ****************************************** 1666 // 1667 G4AffineTransform 1668 G4ITNavigator2::GetMotherToDaughterTransform( 1669 G4 1670 EV 1671 { 1672 CheckNavigatorStateIsValid(); 1673 switch (enteringVolumeType) 1674 { 1675 case kNormal: // Nothing is needed to pr 1676 break; // It is stored already in 1677 case kReplica: // Sets the transform in t 1678 G4Exception("G4ITNavigator2::GetMotherT 1679 "GeomNav0001", FatalExcepti 1680 "Method NOT Implemented yet 1681 break; 1682 case kParameterised: 1683 if( pEnteringPhysVol->GetRegularStructu 1684 { 1685 G4VPVParameterisation *pParam = 1686 pEnteringPhysVol->GetParameterisati 1687 G4VSolid* pSolid = 1688 pParam->ComputeSolid(enteringReplic 1689 pSolid->ComputeDimensions(pParam, ent 1690 1691 // Sets the transform in the Paramete 1692 // 1693 pParam->ComputeTransformation(enterin 1694 1695 // Set the correct solid and material 1696 // 1697 G4LogicalVolume* pLogical = pEntering 1698 pLogical->SetSolid( pSolid ); 1699 } 1700 break; 1701 case kExternal: 1702 G4Exception("G4ITNavigator2::GetMothe 1703 "GeomNav0001", FatalExcep 1704 "Not applicable for exter 1705 break; 1706 } 1707 return G4AffineTransform(pEnteringPhysVol-> 1708 pEnteringPhysVol-> 1709 } 1710 1711 // ****************************************** 1712 // GetLocalExitNormalAndCheck 1713 // 1714 // Obtains the Normal vector to a surface (in 1715 // pointing out of previous volume and into c 1716 // checks the current point against expected 1717 // ****************************************** 1718 // 1719 G4ThreeVector G4ITNavigator2:: 1720 GetLocalExitNormalAndCheck( 1721 #ifdef G4DEBUG_NAVIGATION 1722 const G4ThreeVecto 1723 #else 1724 const G4ThreeVecto 1725 #endif 1726 G4bool* 1727 { 1728 CheckNavigatorStateIsValid(); 1729 #ifdef G4DEBUG_NAVIGATION 1730 // Check Current point against expected 'lo 1731 // 1732 if ( fLastTriedStepComputation ) 1733 { 1734 G4ThreeVector ExpectedBoundaryPointLocal; 1735 1736 const G4AffineTransform& GlobalToLocal= G 1737 ExpectedBoundaryPointLocal = 1738 GlobalToLocal.TransformPoint( ExpectedB 1739 1740 // Add here: Comparison against expected 1741 // i.e. the endpoint of Comput 1742 } 1743 #endif 1744 1745 return GetLocalExitNormal( pValid); 1746 } 1747 1748 // ****************************************** 1749 // GetGlobalExitNormal 1750 // 1751 // Obtains the Normal vector to a surface (in 1752 // pointing out of previous volume and into c 1753 // ****************************************** 1754 // 1755 G4ThreeVector 1756 G4ITNavigator2::GetGlobalExitNormal(const G4T 1757 G4bool 1758 { 1759 CheckNavigatorStateIsValid(); 1760 G4bool validNormal; 1761 G4ThreeVector localNormal, globalNormal; 1762 G4bool usingStored; 1763 1764 usingStored= 1765 fCalculatedExitNormal && 1766 ( ( fLastTriedStepComputation && fExiti 1767 || 1768 (!fLastTriedStepComputation 1769 && (IntersectPointGlobal-fStepEndPo 1770 < (10.0*kCarTolerance*kCarTolera 1771 ) // Calculated it 'just' before & t 1772 // but it did not move position 1773 ); 1774 1775 if( usingStored ) 1776 { 1777 // This was computed in ComputeStep -- an 1778 // 1779 globalNormal = fExitNormalGlobalFrame; 1780 G4double normMag2 = globalNormal.mag2(); 1781 if( std::fabs ( normMag2 - 1.0 ) < perMil 1782 { 1783 *pNormalCalculated = true; // ComputeS 1784 // (fExitin 1785 } 1786 else 1787 { 1788 G4ExceptionDescription message; 1789 1790 message << " ERROR> Expected normal-gl 1791 << " - but |normal| = " << s 1792 << " - and |normal|^ = " << 1793 << " which differs from 1.0 by 1794 << " n = " << fExitNormalGlo 1795 message << "========================== 1796 << G4endl; 1797 G4int oldVerbose = fVerbose; 1798 fVerbose=4; 1799 message << " State of Navigator: " < 1800 message << *this << G4endl; 1801 fVerbose = oldVerbose; 1802 message << "========================== 1803 << G4endl; 1804 1805 G4Exception("G4ITNavigator2::GetGlobal 1806 "GeomNav0003",JustWarning, 1807 "Value obtained from stored glo 1808 1809 // (Re)Compute it now -- as either it 1810 // 1811 localNormal = GetLocalExitNormalAndChe 1812 1813 *pNormalCalculated = fCalculatedExitNo 1814 1815 G4AffineTransform localToGlobal = GetL 1816 globalNormal = localToGlobal.Transform 1817 } 1818 } 1819 else 1820 { 1821 localNormal = GetLocalExitNormalAndCheck( 1822 *pNormalCalculated = fCalculatedExitNorma 1823 1824 #ifdef G4DEBUG_NAVIGATION 1825 usingStored= false; 1826 1827 if( (!validNormal) && !fCalculatedExitNor 1828 { 1829 G4ExceptionDescription edN; 1830 edN << " Calculated = " << fCalculated 1831 edN << " Entering= " << fEntering << 1832 G4int oldVerbose= this->GetVerboseLevel 1833 this->SetVerboseLevel(4); 1834 edN << " State of Navigator: " << G4e 1835 edN << *this << G4endl; 1836 this->SetVerboseLevel( oldVerbose ); 1837 1838 G4Exception("G4ITNavigator2::GetGlobalE 1839 "GeomNav0003", JustWarning, 1840 "LocalExitNormalAndCheck() 1841 } 1842 #endif 1843 1844 G4double localMag2= localNormal.mag2(); 1845 if( validNormal && (std::fabs(localMag2- 1846 { 1847 G4ExceptionDescription edN; 1848 1849 edN << "G4ITNavigator2::GetGlobalExitN 1850 << " Using Local Normal - from ca 1851 << G4endl 1852 << " Local Exit Normal : " << " 1853 << " vec = " << localNormal << G4e 1854 << " Global Exit Normal : " << " 1855 << " vec = " << globalNormal << G4 1856 edN << " Calculated It = " << fC 1857 1858 G4Exception("G4ITNavigator2::GetGlobal 1859 "GeomNav0003",JustWarning, 1860 "Value obtained from new l 1861 localNormal = localNormal.unit(); // S 1862 } 1863 G4AffineTransform localToGlobal = GetLoc 1864 globalNormal = localToGlobal.TransformAx 1865 } 1866 1867 #ifdef G4DEBUG_NAVIGATION 1868 if( usingStored ) 1869 { 1870 G4ThreeVector globalNormAgn; 1871 1872 localNormal= GetLocalExitNormalAndCheck(I 1873 1874 G4AffineTransform localToGlobal = GetLoca 1875 globalNormAgn = localToGlobal.TransformAx 1876 1877 // Check the value computed against fExit 1878 G4ThreeVector diffNorm = globalNormAgn - 1879 if( diffNorm.mag2() > perMillion*CLHEP::p 1880 { 1881 G4ExceptionDescription edDfn; 1882 edDfn << "Found difference in normals i 1883 << "- when Get is called after Co 1884 edDfn << " Magnitude of diff = " 1885 edDfn << " Normal stored (Global) 1886 << G4endl; 1887 edDfn << " Global Computed from Local 1888 G4Exception("G4ITNavigator::GetGlobalEx 1889 JustWarning, edDfn); 1890 } 1891 } 1892 #endif 1893 1894 return globalNormal; 1895 } 1896 1897 // To make the new Voxel Safety the default, 1898 #define G4NEW_SAFETY 1 1899 1900 // ****************************************** 1901 // ComputeSafety 1902 // 1903 // It assumes that it will be 1904 // i) called at the Point in the same volume 1905 // ComputeStep. 1906 // ii) after (or at the end of) ComputeStep O 1907 // ****************************************** 1908 // 1909 G4double G4ITNavigator2::ComputeSafety( const 1910 const G4 1911 const G4 1912 { 1913 CheckNavigatorStateIsValid(); 1914 G4double newSafety = 0.0; 1915 1916 #ifdef G4DEBUG_NAVIGATION 1917 G4int oldcoutPrec = G4cout.precision(8); 1918 if( fVerbose > 0 ) 1919 { 1920 G4cout << "*** G4ITNavigator2::ComputeSaf 1921 << " Called at point: " << pGlo 1922 1923 G4VPhysicalVolume *motherPhysical = fHis 1924 G4cout << " Volume = " << motherPhysic 1925 << " - Maximum length = " << pMaxL 1926 if( fVerbose >= 4 ) 1927 { 1928 G4cout << " ----- Upon entering Com 1929 PrintState(); 1930 } 1931 } 1932 #endif 1933 1934 G4SaveNavigatorState* savedState(nullptr); 1935 1936 G4double distEndpointSq = (pGlobalpoint-fSt 1937 G4bool stayedOnEndpoint = distEndpointSq 1938 G4bool endpointOnSurface = fEnteredDaught 1939 1940 if( endpointOnSurface && stayedOnEndpoint ) 1941 { 1942 #ifdef G4DEBUG_NAVIGATION 1943 if( fVerbose >= 2 ) 1944 { 1945 G4cout << " G4Navigator::ComputeSa 1946 << pGlobalpoint << " - is on surface 1947 if( fEnteredDaughter ) { G4cout << " 1948 if( fExitedMother ) { G4cout << " 1949 G4cout << G4endl; 1950 G4cout << " EndPoint was = " << fStep 1951 } 1952 #endif 1953 newSafety = 0.0; 1954 // return newSafety; 1955 } 1956 else // if( !(endpointOnSurface && stayedOn 1957 { 1958 1959 if (keepState) 1960 { 1961 // SetSavedState(); 1962 savedState = new G4SaveNavigatorState(fpN 1963 } 1964 1965 // Pseudo-relocate to this point (updates 1966 // 1967 LocateGlobalPointWithinVolume( pGlobalpoi 1968 // --->> DANGER: Side effects on sub-na 1969 // Could be replaced again by 'gr 1970 // locates (similar side-effects, 1971 // Solutions: 1972 // 1) Re-locate (to where?) 1973 // 2) Insure that the methods us 1974 // does a relocation (if info 1975 1976 #ifdef G4DEBUG_NAVIGATION 1977 if( fVerbose >= 2 ) 1978 { 1979 G4cout << " G4ITNavigator2::ComputeSaf 1980 << pGlobalpoint << G4endl; 1981 } 1982 #endif 1983 G4VPhysicalVolume *motherPhysical = fHist 1984 G4LogicalVolume *motherLogical = motherPh 1985 G4SmartVoxelHeader* pVoxelHeader = mother 1986 G4ThreeVector localPoint = ComputeLocalPo 1987 1988 if ( fHistory.GetTopVolumeType()!=kReplic 1989 { 1990 switch(CharacteriseDaughters(motherLogi 1991 { 1992 case kNormal: 1993 if ( pVoxelHeader != nullptr ) 1994 { 1995 #ifdef G4NEW_SAFETY 1996 G4double safetyTwo = fpVoxelSafet 1997 *m 1998 newSafety= safetyTwo; // Faster 1999 #else 2000 G4double safetyOldVoxel; 2001 LocateGlobalPointWithinVolume(pGl 2002 safetyOldVoxel = 2003 fvoxelNav.ComputeSafety(localPo 2004 newSafety= safetyOldVoxel; 2005 #endif 2006 } 2007 else 2008 { 2009 newSafety=fnormalNav.ComputeSafet 2010 } 2011 break; 2012 case kParameterised: 2013 if( GetDaughtersRegularStructureId( 2014 { 2015 newSafety=fparamNav.ComputeSafety 2016 } 2017 else // Regular structure 2018 { 2019 newSafety=fregularNav.ComputeSafe 2020 } 2021 break; 2022 case kReplica: 2023 G4Exception("G4ITNavigator2::Comput 2024 FatalException, "Not ap 2025 break; 2026 case kExternal: 2027 G4Exception("G4ITNavigator2::Comput 2028 "GeomNav0001", FatalExc 2029 "Not applicable for ext 2030 break; 2031 } 2032 } 2033 else 2034 { 2035 newSafety = freplicaNav.ComputeSafety(p 2036 f 2037 } 2038 2039 if (keepState) 2040 { 2041 *fpNavigatorState = *savedState; 2042 delete savedState; 2043 // RestoreSavedState(); 2044 // This now overwrites the values of the 2045 } 2046 2047 // Remember last safety origin & value 2048 // 2049 // We overwrite the Safety 'sphere' - kee 2050 fPreviousSftOrigin = pGlobalpoint; 2051 fPreviousSafety = newSafety; 2052 } 2053 2054 #ifdef G4DEBUG_NAVIGATION 2055 if( fVerbose > 1 ) 2056 { 2057 G4cout << " ---- Exiting ComputeSafety 2058 if( fVerbose > 2 ) { PrintState(); } 2059 G4cout << " Returned value of Safety = 2060 } 2061 G4cout.precision(oldcoutPrec); 2062 #endif 2063 2064 return newSafety; 2065 } 2066 2067 2068 // ****************************************** 2069 // RecheckDistanceToCurrentBoundary 2070 // 2071 // Trial method for checking potential displa 2072 // Check position aDisplacedGlobalpoint, to s 2073 // current volume (mother). 2074 // If in mother, check distance to boundary a 2075 // If in entering daughter, check distance ba 2076 // NOTE: 2077 // Can be called only after ComputeStep() is 2078 // Deals only with current volume (and potent 2079 // Caveats 2080 // First VERSION: Does not consider daughter 2081 // neither for checking whether it has 2082 // nor for determining distance to exit 2083 // ****************************************** 2084 // 2085 G4bool G4ITNavigator2::RecheckDistanceToCurre 2086 const G4ThreeVector &aDi 2087 const G4ThreeVector &aNe 2088 const G4double ProposedM 2089 G4double *prDistance, 2090 G4double *prNewSafety) c 2091 { 2092 G4ThreeVector localPosition = ComputeLocal 2093 G4ThreeVector localDirection = ComputeLocal 2094 // G4double Step = kInfinity; 2095 2096 G4bool validExitNormal; 2097 G4ThreeVector exitNormal; 2098 // Check against mother solid 2099 G4VPhysicalVolume *motherPhysical = fHisto 2100 G4LogicalVolume *motherLogical = motherPhys 2101 2102 #ifdef CHECK_ORDER_OF_METHODS 2103 if( ! fLastTriedStepComputation ) 2104 { 2105 G4Exception("G4Navigator::RecheckDistanc 2106 "GeomNav0001", FatalExceptio 2107 "Method must be called after ComputeStep 2108 } 2109 #endif 2110 2111 EInside locatedDaug; // = kUndefined; 2112 G4double daughterStep= DBL_MAX; 2113 G4double daughterSafety= DBL_MAX; 2114 2115 if( fEnteredDaughter ) 2116 { 2117 if( motherLogical->CharacteriseDaughters 2118 2119 // Track arrived at boundary of a daught 2120 // the last call of ComputeStep(). 2121 // In case the proposed displaced point 2122 // it must backtrack at least to the e 2123 // NOTE: No check is made against other 2124 // assumed that the proposed displacem 2125 // this is not needed. 2126 2127 // Must check boundary of current daught 2128 // 2129 G4VPhysicalVolume *candPhysical= fBlocke 2130 G4LogicalVolume *candLogical= candPhysic 2131 G4VSolid *candSolid= candLogica 2132 2133 G4AffineTransform nextLevelTrf(candPhysi 2134 candPhysi 2135 2136 G4ThreeVector dgPosition= nextLevelTrf. 2137 G4ThreeVector dgDirection= nextLevelTrf. 2138 locatedDaug = candSolid->Inside(dgPositi 2139 2140 if( locatedDaug == kInside ) 2141 { 2142 // Reverse direction - and find first 2143 // Must backtrack 2144 G4double distanceBackOut = 2145 candSolid->DistanceToOut(dgPosition 2146 - dgDirect 2147 true, &val 2148 daughterStep= - distanceBackOut; 2149 // No check is made whether the par 2150 // at this location without encount 2151 // a different psurface of the curr 2152 if( prNewSafety != nullptr ) 2153 { 2154 daughterSafety= candSolid->Distanc 2155 } 2156 } 2157 else 2158 { 2159 if( locatedDaug == kOutside ) 2160 { 2161 // See whether it still intersect 2162 // 2163 daughterStep= candSolid->Distanc 2164 2165 if( prNewSafety != nullptr ) 2166 { 2167 daughterSafety= candSolid->Dist 2168 } 2169 } 2170 else 2171 { 2172 // The point remains on the surfac 2173 // 2174 daughterStep= 0.0; 2175 daughterSafety= 0.0; 2176 } 2177 } 2178 2179 // If trial point is in daughter (or on 2180 // answer, the rest is not relevant 2181 // 2182 if( locatedDaug != kOutside ) 2183 { 2184 *prDistance= daughterStep; 2185 if( prNewSafety != nullptr ) { *prNe 2186 return true; 2187 } 2188 // If ever extended, so that some type o 2189 // this would change 2190 } 2191 2192 G4VSolid *motherSolid= motherLogical->GetSo 2193 2194 G4double motherStep= DBL_MAX, motherSafety= 2195 2196 // Check distance to boundary of mother 2197 // 2198 if ( fHistory.GetTopVolumeType()!=kReplica 2199 { 2200 EInside locatedMoth = motherSolid->Insid 2201 2202 if( locatedMoth == kInside ) 2203 { 2204 motherSafety= motherSolid->DistanceTo 2205 if( ProposedMove >= motherSafety ) 2206 { 2207 motherStep= motherSolid->DistanceT 2208 2209 true, &validExit 2210 } 2211 else 2212 { 2213 motherStep= ProposedMove; 2214 } 2215 } 2216 else if( locatedMoth == kOutside) 2217 { 2218 motherSafety= motherSolid->DistanceTo 2219 if( ProposedMove >= motherSafety ) 2220 { 2221 motherStep= - motherSolid->Distan 2222 2223 } 2224 } 2225 else 2226 { 2227 motherSafety= 0.0; 2228 *prDistance= 0.0; // On surface - n 2229 if( prNewSafety != nullptr ) { *prNe 2230 return false; 2231 } 2232 } 2233 else 2234 { 2235 return false; 2236 } 2237 2238 *prDistance= std::min( motherStep, daughte 2239 if( prNewSafety != nullptr ) 2240 { 2241 *prNewSafety= std::min( motherSafety, da 2242 } 2243 return true; 2244 } 2245 2246 2247 // ****************************************** 2248 // CreateTouchableHistoryHandle 2249 // ****************************************** 2250 // 2251 G4TouchableHandle G4ITNavigator2::CreateTouch 2252 { 2253 CheckNavigatorStateIsValid(); 2254 return G4TouchableHandle( CreateTouchableHi 2255 } 2256 2257 // ****************************************** 2258 // PrintState 2259 // ****************************************** 2260 // 2261 void G4ITNavigator2::PrintState() const 2262 { 2263 CheckNavigatorStateIsValid(); 2264 G4long oldcoutPrec = G4cout.precision(4); 2265 if( fVerbose >= 4 ) 2266 { 2267 G4cout << "The current state of G4Navigat 2268 G4cout << " ValidExitNormal= " << fValid 2269 << " ExitNormal = " << fExitN 2270 << " Exiting = " << fExiti 2271 << " Entering = " << fEnter 2272 << " BlockedPhysicalVolume= " ; 2273 if (fBlockedPhysicalVolume==nullptr) 2274 G4cout << "None"; 2275 else 2276 G4cout << fBlockedPhysicalVolume->GetNa 2277 G4cout << G4endl 2278 << " BlockedReplicaNo = " << 2279 << " LastStepWasZero = " << 2280 << G4endl; 2281 } 2282 if( ( 1 < fVerbose) && (fVerbose < 4) ) 2283 { 2284 G4cout << G4endl; // Make sure to line up 2285 G4cout << std::setw(30) << " ExitNormal " 2286 << std::setw( 5) << " Valid " 2287 << std::setw( 9) << " Exiting " 2288 << std::setw( 9) << " Entering" 2289 << std::setw(15) << " Blocked:Volu 2290 << std::setw( 9) << " ReplicaNo" 2291 << std::setw( 8) << " LastStepZero 2292 << G4endl; 2293 G4cout << "( " << std::setw(7) << fExitNo 2294 << ", " << std::setw(7) << fExitNo 2295 << ", " << std::setw(7) << fExitNo 2296 << std::setw( 5) << fValidExitNor 2297 << std::setw( 9) << fExiting 2298 << std::setw( 9) << fEntering 2299 if ( fBlockedPhysicalVolume==nullptr ) 2300 { 2301 G4cout << std::setw(15) << "None"; 2302 } 2303 else 2304 { 2305 G4cout << std::setw(15)<< fBlockedPhysi 2306 } 2307 G4cout << std::setw( 9) << fBlockedRepli 2308 << std::setw( 8) << fLastStepWasZ 2309 << G4endl; 2310 } 2311 if( fVerbose > 2 ) 2312 { 2313 G4cout.precision(8); 2314 G4cout << " Current Localpoint = " << fLa 2315 G4cout << " PreviousSftOrigin = " << fPr 2316 G4cout << " PreviousSafety = " << fPr 2317 } 2318 G4cout.precision(oldcoutPrec); 2319 } 2320 2321 // ****************************************** 2322 // ComputeStepLog 2323 // ****************************************** 2324 // 2325 void G4ITNavigator2::ComputeStepLog(const G4T 2326 G4doub 2327 { 2328 CheckNavigatorStateIsValid(); 2329 // The following checks only make sense if 2330 // than the tolerance. 2331 2332 static const G4double fAccuracyForWarning 2333 fAccuracyForException 2334 2335 G4ThreeVector OriginalGlobalpoint = fHistor 2336 Transfo 2337 2338 G4double shiftOriginSafSq = (fPreviousSftOr 2339 2340 // Check that the starting point of this st 2341 // within the isotropic safety sphere of th 2342 // to a accuracy/precision given by fAccur 2343 // If so give warning. 2344 // If it fails by more than fAccuracyForE 2345 // 2346 if( shiftOriginSafSq >= sqr(fPreviousSafety 2347 { 2348 G4double shiftOrigin = std::sqrt(shiftOri 2349 G4double diffShiftSaf = shiftOrigin - fPr 2350 2351 if( diffShiftSaf > fAccuracyForWarning ) 2352 { 2353 G4long oldcoutPrec= G4cout.precision(8) 2354 G4long oldcerrPrec= G4cerr.precision(10 2355 std::ostringstream message, suggestion; 2356 message << "Accuracy error or slightly 2357 << G4endl 2358 << " The Step's starting po 2359 << std::sqrt(moveLenSq)/mm << " 2360 << " since the last call to 2361 << " This has resulted in m 2362 << shiftOrigin/mm << " mm " 2363 << " from the last point at whi 2364 << " was calculated " << G4 2365 << " which is more than the 2366 << fPreviousSafety/mm << " mm 2367 << " This difference is " 2368 << diffShiftSaf/mm << " mm." << 2369 << " The tolerated accuracy 2370 << fAccuracyForException/mm << 2371 2372 suggestion << " "; 2373 static G4ThreadLocal G4int warnNow = 0; 2374 if( ((++warnNow % 100) == 1) ) 2375 { 2376 message << G4endl 2377 << " This problem can be due 2378 << " - a process that has p 2379 << " larger than the current s 2380 << " - inaccuracy in the co 2381 suggestion << "We suggest that you " 2382 << " - find i) what part 2383 << " ii) through what part 2384 << " for example by r 2385 << G4endl 2386 << " /tracking/ver 2387 << " - check which proc 2388 << " this particle (and lo 2389 << G4endl 2390 << " - in case, create a 2391 << " of this event using:" 2392 << " /tracking/ver 2393 } 2394 G4Exception("G4ITNavigator2::ComputeSte 2395 "GeomNav1002", JustWarning, 2396 message, G4String(suggestio 2397 G4cout.precision(oldcoutPrec); 2398 G4cerr.precision(oldcerrPrec); 2399 } 2400 #ifdef G4DEBUG_NAVIGATION 2401 else 2402 { 2403 G4cerr << "WARNING - G4ITNavigator2::Co 2404 << " The Step's startin 2405 << std::sqrt(moveLenSq) << "," < 2406 << " which has taken it 2407 << " the current safety. " << G4 2408 } 2409 #endif 2410 } 2411 G4double safetyPlus = fPreviousSafety + fAc 2412 if ( shiftOriginSafSq > sqr(safetyPlus) ) 2413 { 2414 std::ostringstream message; 2415 message << "May lead to a crash or unreli 2416 << " Position has shifted 2417 << " notifying the navigator !" < 2418 << " Tolerated safety: " < 2419 << " Computed shift : " < 2420 G4Exception("G4ITNavigator2::ComputeStep( 2421 JustWarning, message); 2422 } 2423 } 2424 2425 // ****************************************** 2426 // Operator << 2427 // ****************************************** 2428 // 2429 std::ostream& operator << (std::ostream &os,c 2430 { 2431 // Old version did only the following: 2432 // os << "Current History: " << G4endl << n 2433 // Old behaviour is recovered for fVerbose 2434 2435 // Adapted from G4ITNavigator2::PrintState( 2436 2437 G4long oldcoutPrec = os.precision(4); 2438 if( n.fVerbose >= 4 ) 2439 { 2440 os << "The current state of G4ITNavigator 2441 os << " ValidExitNormal= " << n.fValidEx 2442 << " ExitNormal = " << n.fExitNormal 2443 << " Exiting = " << n.fExiting 2444 << " Entering = " << n.fEntering 2445 << " BlockedPhysicalVolume= " ; 2446 2447 if (n.fBlockedPhysicalVolume==nullptr) 2448 { 2449 os << "None"; 2450 } 2451 else 2452 { 2453 os << n.fBlockedPhysicalVolume->GetName 2454 } 2455 2456 os << G4endl 2457 << " BlockedReplicaNo = " << n.fBlo 2458 << " LastStepWasZero = " << n.fLa 2459 << G4endl; 2460 } 2461 if( ( 1 < n.fVerbose) && (n.fVerbose < 4) ) 2462 { 2463 os << G4endl; // Make sure to line up 2464 os << std::setw(30) << " ExitNormal " << 2465 << std::setw( 5) << " Valid " << " 2466 << std::setw( 9) << " Exiting " << " 2467 << std::setw( 9) << " Entering" << " 2468 << std::setw(15) << " Blocked:Volume " < 2469 << std::setw( 9) << " ReplicaNo" < 2470 << std::setw( 8) << " LastStepZero " < 2471 << G4endl; 2472 os << "( " << std::setw(7) << n.fExitNorm 2473 << ", " << std::setw(7) << n.fExitNormal. 2474 << ", " << std::setw(7) << n.fExitNormal. 2475 << std::setw( 5) << n.fValidExitNormal 2476 << std::setw( 9) << n.fExiting 2477 << std::setw( 9) << n.fEntering 2478 2479 if ( n.fBlockedPhysicalVolume==nullptr ) 2480 { os << std::setw(15) << "None"; } 2481 else 2482 { os << std::setw(15)<< n.fBlockedPhysica 2483 2484 os << std::setw( 9) << n.fBlockedReplica 2485 << std::setw( 8) << n.fLastStepWasZero 2486 << G4endl; 2487 } 2488 if( n.fVerbose > 2 ) 2489 { 2490 os.precision(8); 2491 os << " Current Localpoint = " << n.fLast 2492 os << " PreviousSftOrigin = " << n.fPrev 2493 os << " PreviousSafety = " << n.fPrev 2494 } 2495 if( n.fVerbose > 3 || n.fVerbose == 0 ) 2496 { 2497 os << "Current History: " << G4endl << n. 2498 } 2499 2500 os.precision(oldcoutPrec); 2501 return os; 2502 } 2503 2504 //------------------------------------------- 2505 2506 EInside 2507 G4ITNavigator2::InsideCurrentVolume(const G4T 2508 { 2509 const G4AffineTransform& transform = GetGlo 2510 G4ThreeVector localPoint(transform.Transfor 2511 2512 G4VSolid* solid = fHistory.GetTopVolume()-> 2513 2514 return solid->Inside(localPoint); 2515 } 2516 2517 //------------------------------------------- 2518 2519 void G4ITNavigator2::GetRandomInCurrentVolume 2520 { 2521 const G4AffineTransform& local2Global = Get 2522 G4VSolid* solid = fHistory.GetTopVolume()-> 2523 2524 G4VoxelLimits voxelLimits; // Defaults to 2525 G4double vmin, vmax; 2526 G4AffineTransform dummy; 2527 std::vector<std::vector<G4double> > fExtend 2528 2529 solid->CalculateExtent(kXAxis,voxelLimits,d 2530 fExtend[kXAxis][BoundingBox::kMin] = vmin; 2531 fExtend[kXAxis][BoundingBox::kMax] = vmax; 2532 2533 solid->CalculateExtent(kYAxis,voxelLimits,d 2534 fExtend[kYAxis][BoundingBox::kMin] = vmin; 2535 fExtend[kYAxis][BoundingBox::kMax] = vmax; 2536 2537 solid->CalculateExtent(kZAxis,voxelLimits,d 2538 fExtend[kZAxis][BoundingBox::kMin] = vmin; 2539 fExtend[kZAxis][BoundingBox::kMax] = vmax; 2540 2541 G4ThreeVector rndmPos; 2542 2543 while(true) 2544 { 2545 for(G4int i = 0 ; i < 3 ; ++i) 2546 { 2547 G4double min = fExtend[i][BoundingBox:: 2548 G4double max = fExtend[i][BoundingBox:: 2549 rndmPos[i] = G4UniformRand()*(max-min)+ 2550 } 2551 2552 if(solid->Inside(rndmPos) == kInside) bre 2553 } 2554 2555 _rndmPoint = local2Global.TransformPoint(rn 2556 } 2557 2558