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