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