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