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 // Author: John Apostolakis 28 // First version: 31 May 2010 29 // 30 // -------------------------------------------------------------------- 31 #include "G4VoxelSafety.hh" 32 33 #include "G4GeometryTolerance.hh" 34 35 #include "G4SmartVoxelProxy.hh" 36 #include "G4SmartVoxelNode.hh" 37 #include "G4SmartVoxelHeader.hh" 38 39 // ******************************************************************** 40 // Constructor 41 // - copied from G4VoxelNavigation (1st version) 42 // ******************************************************************** 43 // 44 G4VoxelSafety::G4VoxelSafety() 45 : fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis), 46 fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0), 47 fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.), 48 fVoxelNodeNoStack(kNavigatorVoxelStackMax,0), 49 fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)nullptr) 50 { 51 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 52 } 53 54 // ******************************************************************** 55 // Destructor 56 // ******************************************************************** 57 // 58 G4VoxelSafety::~G4VoxelSafety() = default; 59 60 // ******************************************************************** 61 // ComputeSafety 62 // 63 // Calculates the isotropic distance to the nearest boundary from the 64 // specified point in the local coordinate system. 65 // The localpoint utilised must be within the current volume. 66 // ******************************************************************** 67 // 68 G4double 69 G4VoxelSafety::ComputeSafety(const G4ThreeVector& localPoint, 70 const G4VPhysicalVolume& currentPhysical, 71 G4double maxLength) 72 { 73 G4LogicalVolume *motherLogical; 74 G4VSolid *motherSolid; 75 G4SmartVoxelHeader *motherVoxelHeader; 76 G4double motherSafety, ourSafety; 77 G4int localNoDaughters; 78 G4double daughterSafety; 79 80 motherLogical = currentPhysical.GetLogicalVolume(); 81 fpMotherLogical= motherLogical; // For use by the other methods 82 motherSolid = motherLogical->GetSolid(); 83 motherVoxelHeader = motherLogical->GetVoxelHeader(); 84 85 #ifdef G4VERBOSE 86 if( fVerbose > 0 ) 87 { 88 G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl; 89 } 90 #endif 91 92 // Check that point is inside mother volume 93 // 94 EInside insideMother = motherSolid->Inside(localPoint); 95 if( insideMother != kInside ) 96 { 97 #ifdef G4DEBUG_NAVIGATION 98 if( insideMother == kOutside ) 99 { 100 std::ostringstream message; 101 message << "Safety method called for location outside current Volume." 102 << G4endl 103 << "Location for safety is Outside this volume. " << G4endl 104 << "The approximate distance to the solid " 105 << "(safety from outside) is: " 106 << motherSolid->DistanceToIn( localPoint ) << G4endl; 107 message << " Problem occurred with physical volume: " 108 << " Name: " << currentPhysical.GetName() 109 << " Copy No: " << currentPhysical.GetCopyNo() << G4endl 110 << " Local Point = " << localPoint << G4endl; 111 message << " Description of solid: " << G4endl 112 << *motherSolid << G4endl; 113 G4Exception("G4VoxelSafety::ComputeSafety()", "GeomNav0003", 114 FatalException, message); 115 } 116 #endif 117 return 0.0; 118 } 119 120 // First limit: mother safety - distance to outer boundaries 121 // 122 motherSafety = motherSolid->DistanceToOut(localPoint); 123 ourSafety = motherSafety; // Working isotropic safety 124 125 #ifdef G4VERBOSE 126 if(( fCheck ) ) // && ( fVerbose == 1 )) 127 { 128 G4cout << " Invoked DistanceToOut(p) for mother solid: " 129 << motherSolid->GetName() 130 << ". Solid replied: " << motherSafety << G4endl 131 << " For local point p: " << localPoint 132 << ", to be considered as 'mother safety'." << G4endl; 133 } 134 #endif 135 localNoDaughters = (G4int)motherLogical->GetNoDaughters(); 136 137 fBlockList.Enlarge(localNoDaughters); 138 fBlockList.Reset(); 139 140 fVoxelDepth = -1; // Resets the depth -- must be done for next method 141 daughterSafety= SafetyForVoxelHeader(motherVoxelHeader, localPoint, maxLength, 142 currentPhysical, 0.0, ourSafety); 143 ourSafety= std::min( motherSafety, daughterSafety ); 144 145 return ourSafety; 146 } 147 148 // ******************************************************************** 149 // SafetyForVoxelNode 150 // 151 // Calculate the safety for volumes included in current Voxel Node 152 // ******************************************************************** 153 // 154 G4double 155 G4VoxelSafety::SafetyForVoxelNode( const G4SmartVoxelNode* curVoxelNode, 156 const G4ThreeVector& localPoint ) 157 { 158 G4double ourSafety = DBL_MAX; 159 160 G4long curNoVolumes, contentNo; 161 G4int sampleNo; 162 G4VPhysicalVolume* samplePhysical; 163 164 G4double sampleSafety = 0.0; 165 G4ThreeVector samplePoint; 166 G4VSolid* ptrSolid = nullptr; 167 168 curNoVolumes = curVoxelNode->GetNoContained(); 169 170 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- ) 171 { 172 sampleNo = curVoxelNode->GetVolume((G4int)contentNo); 173 if ( !fBlockList.IsBlocked(sampleNo) ) 174 { 175 fBlockList.BlockVolume(sampleNo); 176 177 samplePhysical = fpMotherLogical->GetDaughter(sampleNo); 178 G4AffineTransform sampleTf(samplePhysical->GetRotation(), 179 samplePhysical->GetTranslation()); 180 sampleTf.Invert(); 181 samplePoint = sampleTf.TransformPoint(localPoint); 182 ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid(); 183 184 sampleSafety = ptrSolid->DistanceToIn(samplePoint); 185 ourSafety = std::min( sampleSafety, ourSafety ); 186 #ifdef G4VERBOSE 187 if(( fCheck ) && ( fVerbose == 1 )) 188 { 189 G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl 190 << " Invoked DistanceToIn(p) for daughter solid: " 191 << ptrSolid->GetName() 192 << ". Solid replied: " << sampleSafety << G4endl 193 << " For local point p: " << samplePoint 194 << ", to be considered as 'daughter safety'." << G4endl; 195 } 196 #endif 197 } 198 } // end for contentNo 199 200 return ourSafety; 201 } 202 203 // ******************************************************************** 204 // SafetyForVoxelHeader 205 // 206 // Cycles through levels of headers to process each node level 207 // Obtained by modifying VoxelLocate (to cycle through Node Headers) 208 // ********************************************************************* 209 // 210 G4double 211 G4VoxelSafety::SafetyForVoxelHeader( const G4SmartVoxelHeader* pHeader, 212 const G4ThreeVector& localPoint, 213 G4double maxLength, 214 const G4VPhysicalVolume& currentPhysical, //Debug 215 G4double distUpperDepth_Sq, 216 G4double previousMinSafety 217 ) 218 { 219 const G4SmartVoxelHeader* const targetVoxelHeader = pHeader; 220 G4SmartVoxelNode* targetVoxelNode = nullptr; 221 222 const G4SmartVoxelProxy* sampleProxy; 223 EAxis targetHeaderAxis; 224 G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth; 225 G4int targetHeaderNoSlices; 226 G4int targetNodeNo; 227 228 G4double minSafety = previousMinSafety; 229 G4double ourSafety = DBL_MAX; 230 unsigned int checkedNum= 0; 231 232 ++fVoxelDepth; 233 // fVoxelDepth set by ComputeSafety or previous level call 234 235 targetHeaderAxis = targetVoxelHeader->GetAxis(); 236 targetHeaderNoSlices = (G4int)targetVoxelHeader->GetNoSlices(); 237 targetHeaderMin = targetVoxelHeader->GetMinExtent(); 238 targetHeaderMax = targetVoxelHeader->GetMaxExtent(); 239 240 targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin) 241 / targetHeaderNoSlices; 242 243 G4double localCrd = localPoint(targetHeaderAxis); 244 245 const auto candNodeNo = G4int( (localCrd-targetHeaderMin) 246 / targetHeaderNodeWidth ); 247 // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1 248 249 #ifdef G4DEBUG_VOXELISATION 250 if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 ) 251 { 252 G4ExceptionDescription ed; 253 ed << " Potential ERROR." 254 << " Point is outside range of Voxel in current coordinate" << G4endl; 255 ed << " Node number of point " << localPoint 256 << "is outside the range. " << G4endl; 257 ed << " Voxel node Num= " << candNodeNo << " versus minimum= " << 0 258 << " and maximum= " << targetHeaderNoSlices-1 << G4endl; 259 ed << " Axis = " << targetHeaderAxis 260 << " No of slices = " << targetHeaderNoSlices << G4endl; 261 ed << " Local coord = " << localCrd 262 << " Voxel Min = " << targetHeaderMin 263 << " Max = " << targetHeaderMax << G4endl; 264 G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume(); 265 ed << " Current volume (physical) = " << currentPhysical.GetName() 266 << " (logical) = " << pLogical->GetName() << G4endl; 267 G4VSolid* pSolid= pLogical->GetSolid(); 268 ed << " Solid type = " << pSolid->GetEntityType() << G4endl; 269 ed << *pSolid << G4endl; 270 G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003", 271 JustWarning, ed, 272 "Point is outside range of Voxel in current coordinate"); 273 } 274 #endif 275 276 const G4int pointNodeNo = 277 std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) ); 278 279 #ifdef G4VERBOSE 280 if( fVerbose > 2 ) 281 { 282 G4cout << G4endl; 283 G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader " << G4endl; 284 G4cout << " Called at Depth = " << fVoxelDepth; 285 G4cout << " Calculated pointNodeNo= " << pointNodeNo 286 << " from position= " << localPoint(targetHeaderAxis) 287 << " min= " << targetHeaderMin 288 << " max= " << targetHeaderMax 289 << " width= " << targetHeaderNodeWidth 290 << " no-slices= " << targetHeaderNoSlices 291 << " axis= " << targetHeaderAxis << G4endl; 292 } 293 else if (fVerbose == 1) 294 { 295 G4cout << " VoxelSafety: Depth = " << fVoxelDepth 296 << " Number of Slices = " << targetHeaderNoSlices 297 << " Header (address) = " << targetVoxelHeader << G4endl; 298 } 299 #endif 300 301 // Stack info for stepping 302 // 303 fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis; 304 fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices; 305 fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth; 306 307 fVoxelHeaderStack[fVoxelDepth] = pHeader; 308 309 G4int trialUp = -1, trialDown = -1; 310 G4double distUp = DBL_MAX, distDown = DBL_MAX; 311 312 // Using Equivalent voxels - this is pre-initialisation only 313 // 314 G4int nextUp = pointNodeNo+1; 315 G4int nextDown = pointNodeNo-1; 316 317 G4int nextNodeNo = pointNodeNo; 318 G4double distAxis; // Distance in current Axis 319 distAxis = 0.0; // Starting in node containing local Coordinate 320 321 G4bool nextIsInside = false; 322 323 G4double distMaxInterest= std::min( previousMinSafety, maxLength); 324 // We will not look beyond this distance. 325 // This distance will be updated to reflect the 326 // max ( minSafety, maxLength ) at each step 327 328 targetNodeNo = pointNodeNo; 329 do 330 { 331 G4double nodeSafety = DBL_MAX, headerSafety = DBL_MAX; 332 fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo; 333 334 ++checkedNum; 335 336 sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo); 337 338 #ifdef G4DEBUG_NAVIGATION 339 assert( sampleProxy != 0); 340 if( fVerbose > 2 ) 341 { 342 G4cout << " -Checking node " << targetNodeNo 343 << " is proxy with address " << sampleProxy << G4endl; 344 } 345 #endif 346 347 if ( sampleProxy == nullptr ) 348 { 349 G4ExceptionDescription ed; 350 ed << " Problem for node number= " << targetNodeNo 351 << " Number of slides = " << targetHeaderNoSlices 352 << G4endl; 353 G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003", 354 FatalException, ed, 355 "Problem sampleProxy is Zero. Failure in loop."); 356 } 357 else if ( sampleProxy->IsNode() ) 358 { 359 targetVoxelNode = sampleProxy->GetNode(); 360 361 // Deal with the node here [ i.e. the last level ] 362 // 363 nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint); 364 #ifdef G4DEBUG_NAVIGATION 365 if( fVerbose > 2 ) 366 { 367 G4cout << " -- It is a Node "; 368 G4cout << " its safety= " << nodeSafety 369 << " our level Saf = " << ourSafety 370 << " MinSaf= " << minSafety << G4endl; 371 } 372 #endif 373 ourSafety= std::min( ourSafety, nodeSafety ); 374 375 trialUp = targetVoxelNode->GetMaxEquivalentSliceNo()+1; 376 trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1; 377 } 378 else 379 { 380 const G4SmartVoxelHeader* pNewVoxelHeader = sampleProxy->GetHeader(); 381 382 G4double distCombined_Sq; 383 distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis; 384 385 #ifdef G4DEBUG_NAVIGATION 386 if( fVerbose > 2 ) 387 { 388 G4double distCombined= std::sqrt( distCombined_Sq ); 389 G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq ); 390 G4cout << " -- It is a Header " << G4endl; 391 G4cout << " Recurse to deal with next level, fVoxelDepth= " 392 << fVoxelDepth+1 << G4endl; 393 G4cout << " Distances: Upper= " << distUpperDepth 394 << " Axis= " << distAxis 395 << " Combined= " << distCombined << G4endl; 396 } 397 #endif 398 399 // Recurse to deal with lower levels 400 // 401 headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint, 402 maxLength, currentPhysical, 403 distCombined_Sq, minSafety); 404 ourSafety = std::min( ourSafety, headerSafety ); 405 406 #ifdef G4DEBUG_NAVIGATION 407 if( fVerbose > 2 ) 408 { 409 G4cout << " >> Header safety = " << headerSafety 410 << " our level Saf = " << ourSafety << G4endl; 411 } 412 #endif 413 trialUp = pNewVoxelHeader->GetMaxEquivalentSliceNo()+1; 414 trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1; 415 } 416 minSafety = std::min( minSafety, ourSafety ); 417 418 // Find next closest Voxel 419 // - first try: by simple subtraction 420 // - later: using distance (TODO - tbc) 421 // 422 if( targetNodeNo >= pointNodeNo ) 423 { 424 nextUp = trialUp; 425 // distUp = std::max( targetHeaderMax-localCrd, 0.0 ); 426 G4double lowerEdgeOfNext = targetHeaderMin 427 + nextUp * targetHeaderNodeWidth; 428 distUp = lowerEdgeOfNext-localCrd ; 429 if( distUp < 0.0 ) 430 { 431 distUp = DBL_MAX; // On the wrong side - must not be considered 432 } 433 #ifdef G4DEBUG_NAVIGATION 434 if( fVerbose > 2 ) 435 { 436 G4cout << " > Updated nextUp= " << nextUp << G4endl; 437 } 438 #endif 439 } 440 441 if( targetNodeNo <= pointNodeNo ) 442 { 443 nextDown = trialDown; 444 // distDown = std::max( localCrd-targetHeaderMin, 0.0); 445 G4double upperEdgeOfNext = targetHeaderMin 446 + (nextDown+1) * targetHeaderNodeWidth; 447 distDown = localCrd-upperEdgeOfNext; 448 if( distDown < 0.0 ) 449 { 450 distDown= DBL_MAX; // On the wrong side - must not be considered 451 } 452 #ifdef G4DEBUG_NAVIGATION 453 if( fVerbose > 2 ) 454 { 455 G4cout << " > Updated nextDown= " << nextDown << G4endl; 456 } 457 #endif 458 } 459 460 #ifdef G4DEBUG_NAVIGATION 461 if( fVerbose > 2 ) 462 { 463 G4cout << " Node= " << pointNodeNo 464 << " Up: next= " << nextUp << " d# " 465 << nextUp - pointNodeNo 466 << " trialUp= " << trialUp << " d# " 467 << trialUp - pointNodeNo 468 << " Down: next= " << nextDown << " d# " 469 << targetNodeNo - nextDown 470 << " trialDwn= " << trialDown << " d# " 471 << targetNodeNo - trialDown 472 << " condition (next is Inside)= " << nextIsInside 473 << G4endl; 474 } 475 #endif 476 477 G4bool UpIsClosest; 478 UpIsClosest = distUp < distDown; 479 480 if( (nextUp < targetHeaderNoSlices) 481 && (UpIsClosest || (nextDown < 0)) ) 482 { 483 nextNodeNo = nextUp; 484 distAxis = distUp; 485 ++nextUp; // Default 486 #ifdef G4VERBOSE 487 if( fVerbose > 2 ) 488 { 489 G4cout << " > Chose Up. Depth= " << fVoxelDepth 490 << " Nodes: next= " << nextNodeNo 491 << " new nextUp= " << nextUp 492 << " Dist = " << distAxis << G4endl; 493 } 494 #endif 495 } 496 else 497 { 498 nextNodeNo = nextDown; 499 distAxis = distDown; 500 --nextDown; // A default value 501 #ifdef G4VERBOSE 502 if( fVerbose > 2 ) 503 { 504 G4cout << " > Chose Down. Depth= " << fVoxelDepth 505 << " Nodes: next= " << nextNodeNo 506 << " new nextDown= " << nextDown 507 << " Dist = " << distAxis << G4endl; 508 } 509 #endif 510 } 511 512 nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices); 513 if( nextIsInside ) 514 { 515 targetNodeNo= nextNodeNo; 516 517 #ifdef G4DEBUG_NAVIGATION 518 assert( targetVoxelHeader->GetSlice(nextNodeNo) != 0 ); 519 G4bool bContinue = (distAxis<minSafety); 520 if( !bContinue ) 521 { 522 if( fVerbose > 2 ) 523 { 524 G4cout << " Can skip remaining at depth " << targetHeaderAxis 525 << " >> distAxis= " << distAxis 526 << " minSafety= " << minSafety << G4endl; 527 } 528 } 529 #endif 530 } 531 else 532 { 533 #ifdef G4DEBUG_NAVIGATION 534 if( fVerbose > 2) 535 { 536 G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl; 537 G4cout << " VoxSaf> No more candidates: nodeDown= " << nextDown 538 << " nodeUp= " << nextUp 539 << " lastSlice= " << targetHeaderNoSlices << G4endl; 540 } 541 #endif 542 } 543 544 // This calculation can be 'hauled'-up to where minSafety is calculated 545 // 546 distMaxInterest = std::min( minSafety, maxLength ); 547 548 } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq 549 < distMaxInterest*distMaxInterest ) ); 550 551 #ifdef G4VERBOSE 552 if( fVerbose > 0 ) 553 { 554 G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of " 555 << targetHeaderNoSlices << " slices." << G4endl; 556 G4cout << " ===== Returning from SafetyForVoxelHeader " 557 << " Depth= " << fVoxelDepth << G4endl 558 << G4endl; 559 } 560 #endif 561 562 // Go back one level 563 fVoxelDepth--; 564 565 return ourSafety; 566 } 567