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