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 // class G4VoxelNavigation Implementation 27 // 28 // Author: P.Kent, 1996 29 // 30 // -------------------------------------------------------------------- 31 #include "G4VoxelNavigation.hh" 32 #include "G4GeometryTolerance.hh" 33 #include "G4VoxelSafety.hh" 34 35 #include "G4AuxiliaryNavServices.hh" 36 37 #include <cassert> 38 #include <ostream> 39 40 // ******************************************************************** 41 // Constructor 42 // ******************************************************************** 43 // 44 G4VoxelNavigation::G4VoxelNavigation() 45 : fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis), 46 fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0), 47 fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.), 48 fVoxelNodeNoStack(kNavigatorVoxelStackMax,0), 49 fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)nullptr) 50 { 51 fLogger= new G4NavigationLogger("G4VoxelNavigation"); 52 fpVoxelSafety= new G4VoxelSafety(); 53 fHalfTolerance= 0.5*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 54 55 #ifdef G4DEBUG_NAVIGATION 56 SetVerboseLevel(5); // Reports most about daughter volumes 57 #endif 58 } 59 60 // ******************************************************************** 61 // Destructor 62 // ******************************************************************** 63 // 64 G4VoxelNavigation::~G4VoxelNavigation() 65 { 66 delete fpVoxelSafety; 67 delete fLogger; 68 } 69 70 // -------------------------------------------------------------------------- 71 // Input: 72 // exiting: : last step exited 73 // blockedPhysical : phys volume last exited (if exiting) 74 // blockedReplicaNo : copy/replica number of exited 75 // Output: 76 // entering : if true, found candidate volume to enter 77 // blockedPhysical : candidate phys volume to enter - if entering 78 // blockedReplicaNo : copy/replica number - if entering 79 // exiting: : will exit current (mother) volume 80 // In/Out 81 // -------------------------------------------------------------------------- 82 83 // ******************************************************************** 84 // ComputeStep 85 // ******************************************************************** 86 // 87 G4double 88 G4VoxelNavigation::ComputeStep( const G4ThreeVector& localPoint, 89 const G4ThreeVector& localDirection, 90 const G4double currentProposedStepLength, 91 G4double& newSafety, 92 /* const */ G4NavigationHistory& history, 93 G4bool& validExitNormal, 94 G4ThreeVector& exitNormal, 95 G4bool& exiting, 96 G4bool& entering, 97 G4VPhysicalVolume* (*pBlockedPhysical), 98 G4int& blockedReplicaNo ) 99 { 100 G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=nullptr; 101 G4LogicalVolume *motherLogical; 102 G4VSolid *motherSolid; 103 G4ThreeVector sampleDirection; 104 G4double ourStep=currentProposedStepLength, ourSafety; 105 G4double motherSafety, motherStep = DBL_MAX; 106 G4int localNoDaughters, sampleNo; 107 108 G4bool initialNode, noStep; 109 G4SmartVoxelNode *curVoxelNode; 110 G4long curNoVolumes, contentNo; 111 G4double voxelSafety; 112 113 motherPhysical = history.GetTopVolume(); 114 motherLogical = motherPhysical->GetLogicalVolume(); 115 motherSolid = motherLogical->GetSolid(); 116 117 // 118 // Compute mother safety 119 // 120 121 motherSafety = motherSolid->DistanceToOut(localPoint); 122 ourSafety = motherSafety; // Working isotropic safety 123 124 #ifdef G4VERBOSE 125 if ( fCheck ) 126 { 127 fLogger->PreComputeStepLog (motherPhysical, motherSafety, localPoint); 128 } 129 #endif 130 131 // 132 // Compute daughter safeties & intersections 133 // 134 135 // Exiting normal optimisation 136 // 137 if ( exiting && validExitNormal ) 138 { 139 if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine ) 140 { 141 // Block exited daughter volume 142 // 143 blockedExitedVol = *pBlockedPhysical; 144 ourSafety = 0; 145 } 146 } 147 exiting = false; 148 entering = false; 149 150 // For extra checking, get the distance to Mother early !! 151 G4bool motherValidExitNormal = false; 152 G4ThreeVector motherExitNormal(0.0, 0.0, 0.0); 153 154 #ifdef G4VERBOSE 155 if ( fCheck ) 156 { 157 // Compute early -- a) for validity 158 // b) to check against answer of daughters! 159 motherStep = motherSolid->DistanceToOut(localPoint, 160 localDirection, 161 true, 162 &motherValidExitNormal, 163 &motherExitNormal); 164 } 165 #endif 166 167 localNoDaughters = (G4int)motherLogical->GetNoDaughters(); 168 169 fBList.Enlarge(localNoDaughters); 170 fBList.Reset(); 171 172 initialNode = true; 173 noStep = true; 174 175 while (noStep) 176 { 177 curVoxelNode = fVoxelNode; 178 curNoVolumes = curVoxelNode->GetNoContained(); 179 for (contentNo=curNoVolumes-1; contentNo>=0; contentNo--) 180 { 181 sampleNo = curVoxelNode->GetVolume((G4int)contentNo); 182 if ( !fBList.IsBlocked(sampleNo) ) 183 { 184 fBList.BlockVolume(sampleNo); 185 samplePhysical = motherLogical->GetDaughter(sampleNo); 186 if ( samplePhysical!=blockedExitedVol ) 187 { 188 G4AffineTransform sampleTf(samplePhysical->GetRotation(), 189 samplePhysical->GetTranslation()); 190 sampleTf.Invert(); 191 const G4ThreeVector samplePoint = 192 sampleTf.TransformPoint(localPoint); 193 const G4VSolid *sampleSolid = 194 samplePhysical->GetLogicalVolume()->GetSolid(); 195 const G4double sampleSafety = 196 sampleSolid->DistanceToIn(samplePoint); 197 198 if ( sampleSafety<ourSafety ) 199 { 200 ourSafety = sampleSafety; 201 } 202 if ( sampleSafety<=ourStep ) 203 { 204 sampleDirection = sampleTf.TransformAxis(localDirection); 205 G4double sampleStep = 206 sampleSolid->DistanceToIn(samplePoint, sampleDirection); 207 #ifdef G4VERBOSE 208 if( fCheck ) 209 { 210 fLogger->PrintDaughterLog(sampleSolid, samplePoint, 211 sampleSafety, true, 212 sampleDirection, sampleStep); 213 } 214 #endif 215 if ( sampleStep<=ourStep ) 216 { 217 ourStep = sampleStep; 218 entering = true; 219 exiting = false; 220 *pBlockedPhysical = samplePhysical; 221 blockedReplicaNo = -1; 222 #ifdef G4VERBOSE 223 // Check to see that the resulting point is indeed in/on volume. 224 // This could be done only for successful candidate. 225 if ( fCheck ) 226 { 227 fLogger->AlongComputeStepLog (sampleSolid, samplePoint, 228 sampleDirection, localDirection, sampleSafety, sampleStep); 229 } 230 #endif 231 } 232 #ifdef G4VERBOSE 233 if ( fCheck && ( sampleStep < kInfinity ) 234 && ( sampleStep >= motherStep ) ) 235 { 236 // The intersection point with the daughter is after the exit 237 // point from the mother volume. Double check this !! 238 fLogger->CheckDaughterEntryPoint(sampleSolid, 239 samplePoint, sampleDirection, 240 motherSolid, 241 localPoint, localDirection, 242 motherStep, sampleStep); 243 } 244 #endif 245 } 246 #ifdef G4VERBOSE 247 else // ie if sampleSafety > outStep 248 { 249 if( fCheck ) 250 { 251 fLogger->PrintDaughterLog(sampleSolid, samplePoint, 252 sampleSafety, false, 253 G4ThreeVector(0.,0.,0.), -1.0 ); 254 } 255 } 256 #endif 257 } 258 } 259 } 260 if (initialNode) 261 { 262 initialNode = false; 263 voxelSafety = ComputeVoxelSafety(localPoint); 264 if ( voxelSafety<ourSafety ) 265 { 266 ourSafety = voxelSafety; 267 } 268 if ( currentProposedStepLength<ourSafety ) 269 { 270 // Guaranteed physics limited 271 // 272 noStep = false; 273 entering = false; 274 exiting = false; 275 *pBlockedPhysical = nullptr; 276 ourStep = kInfinity; 277 } 278 else 279 { 280 // 281 // Compute mother intersection if required 282 // 283 if ( motherSafety<=ourStep ) 284 { 285 // In case of check mode this is a duplicate call -- acceptable 286 motherStep = motherSolid->DistanceToOut(localPoint, localDirection, 287 true, &motherValidExitNormal, &motherExitNormal); 288 #ifdef G4VERBOSE 289 if ( fCheck ) 290 { 291 fLogger->PostComputeStepLog(motherSolid, localPoint, localDirection, 292 motherStep, motherSafety); 293 if( motherValidExitNormal ) 294 { 295 fLogger->CheckAndReportBadNormal(motherExitNormal, 296 localPoint, localDirection, 297 motherStep, motherSolid, 298 "From motherSolid::DistanceToOut" ); 299 } 300 } 301 #endif 302 if( (motherStep >= kInfinity) || (motherStep < 0.0) ) 303 { 304 #ifdef G4VERBOSE 305 if( fCheck ) // Error - indication of being outside solid !! 306 { 307 fLogger->ReportOutsideMother(localPoint, localDirection, 308 motherPhysical); 309 } 310 #endif 311 motherStep = 0.0; 312 ourStep = 0.0; 313 exiting = true; 314 entering = false; 315 316 // validExitNormal= motherValidExitNormal; 317 // exitNormal= motherExitNormal; 318 // Useful only if the point is very close to surface 319 // => but it would need to be rotated to grand-mother ref frame ! 320 validExitNormal= false; 321 322 *pBlockedPhysical = nullptr; // or motherPhysical ? 323 blockedReplicaNo = 0; // or motherReplicaNumber ? 324 325 newSafety = 0.0; 326 return ourStep; 327 } 328 329 if ( motherStep<=ourStep ) 330 { 331 ourStep = motherStep; 332 exiting = true; 333 entering = false; 334 335 // Exit normal: Natural location to set these;confirmed short step 336 // 337 validExitNormal = motherValidExitNormal; 338 exitNormal = motherExitNormal; 339 340 if ( validExitNormal ) 341 { 342 const G4RotationMatrix *rot = motherPhysical->GetRotation(); 343 if (rot != nullptr) 344 { 345 exitNormal *= rot->inverse(); 346 #ifdef G4VERBOSE 347 if( fCheck ) 348 { 349 fLogger->CheckAndReportBadNormal(exitNormal, // rotated 350 motherExitNormal, // original 351 *rot, 352 "From RotationMatrix" ); 353 } 354 #endif 355 } 356 } 357 } 358 else 359 { 360 validExitNormal = false; 361 } 362 } 363 } 364 newSafety = ourSafety; 365 } 366 if (noStep) 367 { 368 noStep = LocateNextVoxel(localPoint, localDirection, ourStep); 369 } 370 } // end -while (noStep)- loop 371 372 return ourStep; 373 } 374 375 // ******************************************************************** 376 // ComputeVoxelSafety 377 // 378 // Computes safety from specified point to voxel boundaries 379 // using already located point 380 // o collected boundaries for most derived level 381 // o adjacent boundaries for previous levels 382 // ******************************************************************** 383 // 384 G4double 385 G4VoxelNavigation::ComputeVoxelSafety(const G4ThreeVector& localPoint) const 386 { 387 G4SmartVoxelHeader *curHeader; 388 G4double voxelSafety, curNodeWidth; 389 G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta; 390 G4int minCurNodeNoDelta, maxCurNodeNoDelta; 391 G4int localVoxelDepth, curNodeNo; 392 EAxis curHeaderAxis; 393 394 localVoxelDepth = fVoxelDepth; 395 396 curHeader = fVoxelHeaderStack[localVoxelDepth]; 397 curHeaderAxis = fVoxelAxisStack[localVoxelDepth]; 398 curNodeNo = fVoxelNodeNoStack[localVoxelDepth]; 399 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth]; 400 401 // Compute linear intersection distance to boundaries of max/min 402 // to collected nodes at current level 403 // 404 curNodeOffset = curNodeNo*curNodeWidth; 405 maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo; 406 minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo(); 407 minCurCommonDelta = localPoint(curHeaderAxis) 408 - curHeader->GetMinExtent() - curNodeOffset; 409 maxCurCommonDelta = curNodeWidth-minCurCommonDelta; 410 411 if ( minCurNodeNoDelta<maxCurNodeNoDelta ) 412 { 413 voxelSafety = minCurNodeNoDelta*curNodeWidth; 414 voxelSafety += minCurCommonDelta; 415 } 416 else if (maxCurNodeNoDelta < minCurNodeNoDelta) 417 { 418 voxelSafety = maxCurNodeNoDelta*curNodeWidth; 419 voxelSafety += maxCurCommonDelta; 420 } 421 else // (maxCurNodeNoDelta == minCurNodeNoDelta) 422 { 423 voxelSafety = minCurNodeNoDelta*curNodeWidth; 424 voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta); 425 } 426 427 // Compute isotropic safety to boundaries of previous levels 428 // [NOT to collected boundaries] 429 430 // Loop checking, 07.10.2016, JA 431 while ( (localVoxelDepth>0) && (voxelSafety>0) ) 432 { 433 localVoxelDepth--; 434 curHeader = fVoxelHeaderStack[localVoxelDepth]; 435 curHeaderAxis = fVoxelAxisStack[localVoxelDepth]; 436 curNodeNo = fVoxelNodeNoStack[localVoxelDepth]; 437 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth]; 438 curNodeOffset = curNodeNo*curNodeWidth; 439 minCurCommonDelta = localPoint(curHeaderAxis) 440 - curHeader->GetMinExtent() - curNodeOffset; 441 maxCurCommonDelta = curNodeWidth-minCurCommonDelta; 442 443 if ( minCurCommonDelta<voxelSafety ) 444 { 445 voxelSafety = minCurCommonDelta; 446 } 447 if ( maxCurCommonDelta<voxelSafety ) 448 { 449 voxelSafety = maxCurCommonDelta; 450 } 451 } 452 if ( voxelSafety<0 ) 453 { 454 voxelSafety = 0; 455 } 456 457 return voxelSafety; 458 } 459 460 // ******************************************************************** 461 // LocateNextVoxel 462 // 463 // Finds the next voxel from the current voxel and point 464 // in the specified direction 465 // 466 // Returns false if all voxels considered 467 // [current Step ends inside same voxel or leaves all voxels] 468 // true otherwise 469 // [the information on the next voxel is put into the set of 470 // fVoxel* variables & "stacks"] 471 // ******************************************************************** 472 // 473 G4bool 474 G4VoxelNavigation::LocateNextVoxel(const G4ThreeVector& localPoint, 475 const G4ThreeVector& localDirection, 476 const G4double currentStep) 477 { 478 G4SmartVoxelHeader *workHeader=nullptr, *newHeader=nullptr; 479 G4SmartVoxelProxy *newProxy=nullptr; 480 G4SmartVoxelNode *newVoxelNode=nullptr; 481 G4ThreeVector targetPoint, voxelPoint; 482 G4double workNodeWidth, workMinExtent, workCoord; 483 G4double minVal, maxVal, newDistance=0.; 484 G4double newHeaderMin, newHeaderNodeWidth; 485 G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0; 486 EAxis workHeaderAxis, newHeaderAxis; 487 G4bool isNewVoxel = false; 488 489 G4double currentDistance = currentStep; 490 491 // Determine if end of Step within current voxel 492 // 493 for (depth=0; depth<fVoxelDepth; ++depth) 494 { 495 targetPoint = localPoint+localDirection*currentDistance; 496 newDistance = currentDistance; 497 workHeader = fVoxelHeaderStack[depth]; 498 workHeaderAxis = fVoxelAxisStack[depth]; 499 workNodeNo = fVoxelNodeNoStack[depth]; 500 workNodeWidth = fVoxelSliceWidthStack[depth]; 501 workMinExtent = workHeader->GetMinExtent(); 502 workCoord = targetPoint(workHeaderAxis); 503 minVal = workMinExtent+workNodeNo*workNodeWidth; 504 505 if ( minVal<=workCoord+fHalfTolerance ) 506 { 507 maxVal = minVal+workNodeWidth; 508 if ( maxVal<=workCoord-fHalfTolerance ) 509 { 510 // Must consider next voxel 511 // 512 newNodeNo = workNodeNo+1; 513 newHeader = workHeader; 514 newDistance = (maxVal-localPoint(workHeaderAxis)) 515 / localDirection(workHeaderAxis); 516 isNewVoxel = true; 517 newDepth = depth; 518 } 519 } 520 else 521 { 522 newNodeNo = workNodeNo-1; 523 newHeader = workHeader; 524 newDistance = (minVal-localPoint(workHeaderAxis)) 525 / localDirection(workHeaderAxis); 526 isNewVoxel = true; 527 newDepth = depth; 528 } 529 currentDistance = newDistance; 530 } 531 targetPoint = localPoint+localDirection*currentDistance; 532 533 // Check if end of Step within collected boundaries of current voxel 534 // 535 depth = fVoxelDepth; 536 { 537 workHeader = fVoxelHeaderStack[depth]; 538 workHeaderAxis = fVoxelAxisStack[depth]; 539 workNodeNo = fVoxelNodeNoStack[depth]; 540 workNodeWidth = fVoxelSliceWidthStack[depth]; 541 workMinExtent = workHeader->GetMinExtent(); 542 workCoord = targetPoint(workHeaderAxis); 543 minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth; 544 545 if ( minVal<=workCoord+fHalfTolerance ) 546 { 547 maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1) 548 *workNodeWidth; 549 if ( maxVal<=workCoord-fHalfTolerance ) 550 { 551 newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1; 552 newHeader = workHeader; 553 newDistance = (maxVal-localPoint(workHeaderAxis)) 554 / localDirection(workHeaderAxis); 555 isNewVoxel = true; 556 newDepth = depth; 557 } 558 } 559 else 560 { 561 newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1; 562 newHeader = workHeader; 563 newDistance = (minVal-localPoint(workHeaderAxis)) 564 / localDirection(workHeaderAxis); 565 isNewVoxel = true; 566 newDepth = depth; 567 } 568 currentDistance = newDistance; 569 } 570 if (isNewVoxel) 571 { 572 // Compute new voxel & adjust voxel stack 573 // 574 // newNodeNo=Candidate node no at 575 // newDepth =refinement depth of crossed voxel boundary 576 // newHeader=Header for crossed voxel 577 // newDistance=distance to crossed voxel boundary (along the track) 578 // 579 if ( (newNodeNo<0) || (newNodeNo>=G4int(newHeader->GetNoSlices()))) 580 { 581 // Leaving mother volume 582 // 583 isNewVoxel = false; 584 } 585 else 586 { 587 // Compute intersection point on the least refined 588 // voxel boundary that is hit 589 // 590 voxelPoint = localPoint+localDirection*newDistance; 591 fVoxelNodeNoStack[newDepth] = newNodeNo; 592 fVoxelDepth = newDepth; 593 newVoxelNode = nullptr; 594 while ( newVoxelNode == nullptr ) 595 { 596 newProxy = newHeader->GetSlice(newNodeNo); 597 if (newProxy->IsNode()) 598 { 599 newVoxelNode = newProxy->GetNode(); 600 } 601 else 602 { 603 ++fVoxelDepth; 604 newHeader = newProxy->GetHeader(); 605 newHeaderAxis = newHeader->GetAxis(); 606 newHeaderNoSlices = (G4int)newHeader->GetNoSlices(); 607 newHeaderMin = newHeader->GetMinExtent(); 608 newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin) 609 / newHeaderNoSlices; 610 newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin) 611 / newHeaderNodeWidth ); 612 // Rounding protection 613 // 614 if ( newNodeNo<0 ) 615 { 616 newNodeNo=0; 617 } 618 else if ( newNodeNo>=newHeaderNoSlices ) 619 { 620 newNodeNo = newHeaderNoSlices-1; 621 } 622 // Stack info for stepping 623 // 624 fVoxelAxisStack[fVoxelDepth] = newHeaderAxis; 625 fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices; 626 fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth; 627 fVoxelNodeNoStack[fVoxelDepth] = newNodeNo; 628 fVoxelHeaderStack[fVoxelDepth] = newHeader; 629 } 630 } 631 fVoxelNode = newVoxelNode; 632 } 633 } 634 return isNewVoxel; 635 } 636 637 // ******************************************************************** 638 // ComputeSafety 639 // 640 // Calculates the isotropic distance to the nearest boundary from the 641 // specified point in the local coordinate system. 642 // The localpoint utilised must be within the current volume. 643 // ******************************************************************** 644 // 645 G4double 646 G4VoxelNavigation::ComputeSafety(const G4ThreeVector& localPoint, 647 const G4NavigationHistory& history, 648 const G4double maxLength) 649 { 650 G4VPhysicalVolume *motherPhysical, *samplePhysical; 651 G4LogicalVolume *motherLogical; 652 G4VSolid *motherSolid; 653 G4double motherSafety, ourSafety; 654 G4int sampleNo; 655 G4SmartVoxelNode *curVoxelNode; 656 G4long curNoVolumes, contentNo; 657 G4double voxelSafety; 658 659 motherPhysical = history.GetTopVolume(); 660 motherLogical = motherPhysical->GetLogicalVolume(); 661 motherSolid = motherLogical->GetSolid(); 662 663 if( fBestSafety ) 664 { 665 return fpVoxelSafety->ComputeSafety( localPoint,*motherPhysical,maxLength ); 666 } 667 668 // 669 // Compute mother safety 670 // 671 672 motherSafety = motherSolid->DistanceToOut(localPoint); 673 ourSafety = motherSafety; // Working isotropic safety 674 675 if( motherSafety == 0.0 ) 676 { 677 #ifdef G4DEBUG_NAVIGATION 678 // Check that point is inside mother volume 679 EInside insideMother = motherSolid->Inside(localPoint); 680 681 if( insideMother == kOutside ) 682 { 683 G4ExceptionDescription message; 684 message << "Safety method called for location outside current Volume." << G4endl 685 << "Location for safety is Outside this volume. " << G4endl 686 << "The approximate distance to the solid " 687 << "(safety from outside) is: " 688 << motherSolid->DistanceToIn( localPoint ) << G4endl; 689 message << " Problem occurred with physical volume: " 690 << " Name: " << motherPhysical->GetName() 691 << " Copy No: " << motherPhysical->GetCopyNo() << G4endl 692 << " Local Point = " << localPoint << G4endl; 693 message << " Description of solid: " << G4endl 694 << *motherSolid << G4endl; 695 G4Exception("G4VoxelNavigation::ComputeSafety()", "GeomNav0003", 696 JustWarning, message); 697 } 698 699 // Following check is NOT for an issue - it is only for information 700 // It is allowed that a solid gives approximate safety - even zero. 701 // 702 if( insideMother == kInside ) // && fVerbose ) 703 { 704 G4ExceptionDescription messageIn; 705 706 messageIn << " Point is Inside, but safety is Zero ." << G4endl; 707 messageIn << " Inexact safety for volume " << motherPhysical->GetName() << G4endl 708 << " Solid: Name= " << motherSolid->GetName() 709 << " Type= " << motherSolid->GetEntityType() << G4endl; 710 messageIn << " Local point= " << localPoint << G4endl; 711 messageIn << " Solid parameters: " << G4endl << *motherSolid << G4endl; 712 G4Exception("G4VoxelNavigation::ComputeSafety()", "GeomNav0003", 713 JustWarning, messageIn); 714 } 715 #endif 716 // if( insideMother != kInside ) 717 return 0.0; 718 } 719 720 #ifdef G4VERBOSE 721 if( fCheck ) 722 { 723 fLogger->ComputeSafetyLog (motherSolid,localPoint,motherSafety,true,1); 724 } 725 #endif 726 // 727 // Compute daughter safeties 728 // 729 // Look only inside the current Voxel only (in the first version). 730 // 731 curVoxelNode = fVoxelNode; 732 curNoVolumes = curVoxelNode->GetNoContained(); 733 734 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- ) 735 { 736 sampleNo = curVoxelNode->GetVolume((G4int)contentNo); 737 samplePhysical = motherLogical->GetDaughter(sampleNo); 738 739 G4AffineTransform sampleTf(samplePhysical->GetRotation(), 740 samplePhysical->GetTranslation()); 741 sampleTf.Invert(); 742 const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint); 743 const G4VSolid* sampleSolid= samplePhysical->GetLogicalVolume()->GetSolid(); 744 G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint); 745 if ( sampleSafety<ourSafety ) 746 { 747 ourSafety = sampleSafety; 748 } 749 #ifdef G4VERBOSE 750 if( fCheck ) 751 { 752 fLogger->ComputeSafetyLog(sampleSolid, samplePoint, 753 sampleSafety, false, 0); 754 } 755 #endif 756 } 757 voxelSafety = ComputeVoxelSafety(localPoint); 758 if ( voxelSafety<ourSafety ) 759 { 760 ourSafety = voxelSafety; 761 } 762 return ourSafety; 763 } 764 765 void G4VoxelNavigation::RelocateWithinVolume( G4VPhysicalVolume* motherPhysical, 766 const G4ThreeVector& localPoint ) 767 { 768 auto motherLogical = motherPhysical->GetLogicalVolume(); 769 770 assert(motherLogical != nullptr); 771 772 if ( auto pVoxelHeader = motherLogical->GetVoxelHeader() ) 773 VoxelLocate( pVoxelHeader, localPoint ); 774 } 775 776 // ******************************************************************** 777 // SetVerboseLevel 778 // ******************************************************************** 779 // 780 void G4VoxelNavigation::SetVerboseLevel(G4int level) 781 { 782 if( fLogger != nullptr ) { fLogger->SetVerboseLevel(level); } 783 if( fpVoxelSafety != nullptr) { fpVoxelSafety->SetVerboseLevel(level); } 784 } 785