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