| 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 // $Id: G4VoxelSafety.cc,v $
26 // 27 //
27 // Author: John Apostolakis 28 // Author: John Apostolakis
28 // First version: 31 May 2010 29 // First version: 31 May 2010
29 // 30 //
30 // ------------------------------------------- 31 // --------------------------------------------------------------------
31 #include "G4VoxelSafety.hh" 32 #include "G4VoxelSafety.hh"
32 33
33 #include "G4GeometryTolerance.hh" 34 #include "G4GeometryTolerance.hh"
34 35
35 #include "G4SmartVoxelProxy.hh" 36 #include "G4SmartVoxelProxy.hh"
36 #include "G4SmartVoxelNode.hh" 37 #include "G4SmartVoxelNode.hh"
37 #include "G4SmartVoxelHeader.hh" 38 #include "G4SmartVoxelHeader.hh"
38 39
39 // ******************************************* 40 // ********************************************************************
40 // Constructor 41 // Constructor
41 // - copied from G4VoxelNavigation (1st v 42 // - copied from G4VoxelNavigation (1st version)
42 // ******************************************* 43 // ********************************************************************
43 // 44 //
44 G4VoxelSafety::G4VoxelSafety() 45 G4VoxelSafety::G4VoxelSafety()
45 : fVoxelAxisStack(kNavigatorVoxelStackMax,kX << 46 : fBlockList(),
>> 47 fpMotherLogical(0),
>> 48 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),
>> 55 fCheck(false),
>> 56 fVerbose(0)
50 { 57 {
51 kCarTolerance = G4GeometryTolerance::GetInst 58 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
52 } 59 }
53 60
54 // ******************************************* 61 // ********************************************************************
55 // Destructor 62 // Destructor
56 // ******************************************* 63 // ********************************************************************
57 // 64 //
58 G4VoxelSafety::~G4VoxelSafety() = default; << 65 G4VoxelSafety::~G4VoxelSafety()
>> 66 {
>> 67 }
59 68
60 // ******************************************* 69 // ********************************************************************
61 // ComputeSafety 70 // ComputeSafety
62 // 71 //
63 // Calculates the isotropic distance to the ne 72 // Calculates the isotropic distance to the nearest boundary from the
64 // specified point in the local coordinate sys 73 // specified point in the local coordinate system.
65 // The localpoint utilised must be within the 74 // The localpoint utilised must be within the current volume.
66 // ******************************************* 75 // ********************************************************************
67 // 76 //
68 G4double 77 G4double
69 G4VoxelSafety::ComputeSafety(const G4ThreeVect << 78 G4VoxelSafety::ComputeSafety(const G4ThreeVector& localPoint,
70 const G4VPhysical 79 const G4VPhysicalVolume& currentPhysical,
71 G4double ma << 80 G4double maxLength)
72 { 81 {
73 G4LogicalVolume *motherLogical; 82 G4LogicalVolume *motherLogical;
74 G4VSolid *motherSolid; 83 G4VSolid *motherSolid;
75 G4SmartVoxelHeader *motherVoxelHeader; 84 G4SmartVoxelHeader *motherVoxelHeader;
76 G4double motherSafety, ourSafety; 85 G4double motherSafety, ourSafety;
77 G4int localNoDaughters; 86 G4int localNoDaughters;
78 G4double daughterSafety; 87 G4double daughterSafety;
79 88
80 motherLogical = currentPhysical.GetLogicalVo 89 motherLogical = currentPhysical.GetLogicalVolume();
81 fpMotherLogical= motherLogical; // For use 90 fpMotherLogical= motherLogical; // For use by the other methods
82 motherSolid = motherLogical->GetSolid(); 91 motherSolid = motherLogical->GetSolid();
83 motherVoxelHeader = motherLogical->GetVoxelH << 92 motherVoxelHeader= motherLogical->GetVoxelHeader();
84 93
85 #ifdef G4VERBOSE 94 #ifdef G4VERBOSE
86 if( fVerbose > 0 ) 95 if( fVerbose > 0 )
87 { 96 {
88 G4cout << "*** G4VoxelSafety::ComputeSafet 97 G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl;
89 } 98 }
90 #endif 99 #endif
91 100
92 // Check that point is inside mother volume 101 // Check that point is inside mother volume
93 // 102 //
94 EInside insideMother = motherSolid->Inside( << 103 EInside insideMother= motherSolid->Inside(localPoint);
95 if( insideMother != kInside ) 104 if( insideMother != kInside )
96 { 105 {
97 #ifdef G4DEBUG_NAVIGATION 106 #ifdef G4DEBUG_NAVIGATION
98 if( insideMother == kOutside ) 107 if( insideMother == kOutside )
99 { 108 {
100 std::ostringstream message; 109 std::ostringstream message;
101 message << "Safety method called for loc 110 message << "Safety method called for location outside current Volume."
102 << G4endl 111 << G4endl
103 << "Location for safety is Outsi 112 << "Location for safety is Outside this volume. " << G4endl
104 << "The approximate distance to 113 << "The approximate distance to the solid "
105 << "(safety from outside) is: " 114 << "(safety from outside) is: "
106 << motherSolid->DistanceToIn( lo 115 << motherSolid->DistanceToIn( localPoint ) << G4endl;
107 message << " Problem occurred with phys 116 message << " Problem occurred with physical volume: "
108 << " Name: " << currentPhysical. 117 << " Name: " << currentPhysical.GetName()
109 << " Copy No: " << currentPhysic 118 << " Copy No: " << currentPhysical.GetCopyNo() << G4endl
110 << " Local Point = " << local 119 << " Local Point = " << localPoint << G4endl;
111 message << " Description of solid: " << 120 message << " Description of solid: " << G4endl
112 << *motherSolid << G4endl; 121 << *motherSolid << G4endl;
113 G4Exception("G4VoxelSafety::ComputeSafet 122 G4Exception("G4VoxelSafety::ComputeSafety()", "GeomNav0003",
114 FatalException, message); 123 FatalException, message);
115 } 124 }
116 #endif 125 #endif
117 return 0.0; 126 return 0.0;
118 } 127 }
119 128
120 // First limit: mother safety - distance t 129 // First limit: mother safety - distance to outer boundaries
121 // 130 //
122 motherSafety = motherSolid->DistanceToOut(lo 131 motherSafety = motherSolid->DistanceToOut(localPoint);
123 ourSafety = motherSafety; // 132 ourSafety = motherSafety; // Working isotropic safety
124 133
125 #ifdef G4VERBOSE 134 #ifdef G4VERBOSE
126 if(( fCheck ) ) // && ( fVerbose == 1 )) 135 if(( fCheck ) ) // && ( fVerbose == 1 ))
127 { 136 {
128 G4cout << " Invoked DistanceToOut(p) fo 137 G4cout << " Invoked DistanceToOut(p) for mother solid: "
129 << motherSolid->GetName() 138 << motherSolid->GetName()
130 << ". Solid replied: " << motherSaf 139 << ". Solid replied: " << motherSafety << G4endl
131 << " For local point p: " << loc 140 << " For local point p: " << localPoint
132 << ", to be considered as 'mother s 141 << ", to be considered as 'mother safety'." << G4endl;
133 } 142 }
134 #endif 143 #endif
135 localNoDaughters = (G4int)motherLogical->Get << 144 localNoDaughters = motherLogical->GetNoDaughters();
136 145
137 fBlockList.Enlarge(localNoDaughters); 146 fBlockList.Enlarge(localNoDaughters);
138 fBlockList.Reset(); 147 fBlockList.Reset();
139 148
140 fVoxelDepth = -1; // Resets the depth -- mu 149 fVoxelDepth = -1; // Resets the depth -- must be done for next method
141 daughterSafety= SafetyForVoxelHeader(motherV 150 daughterSafety= SafetyForVoxelHeader(motherVoxelHeader, localPoint, maxLength,
142 current << 151 currentPhysical, 0, ourSafety);
143 ourSafety= std::min( motherSafety, daughterS 152 ourSafety= std::min( motherSafety, daughterSafety );
144 153
145 return ourSafety; 154 return ourSafety;
146 } 155 }
147 156
148 // ******************************************* 157 // ********************************************************************
149 // SafetyForVoxelNode 158 // SafetyForVoxelNode
150 // 159 //
151 // Calculate the safety for volumes included i 160 // Calculate the safety for volumes included in current Voxel Node
152 // ******************************************* 161 // ********************************************************************
153 // 162 //
154 G4double 163 G4double
155 G4VoxelSafety::SafetyForVoxelNode( const G4Sma << 164 G4VoxelSafety::SafetyForVoxelNode( const G4SmartVoxelNode *curVoxelNode,
156 const G4Thr << 165 const G4ThreeVector& localPoint )
157 { 166 {
158 G4double ourSafety = DBL_MAX; << 167 G4double ourSafety= DBL_MAX;
159 168
160 G4long curNoVolumes, contentNo; << 169 G4int curNoVolumes, contentNo, sampleNo;
161 G4int sampleNo; << 170 G4VPhysicalVolume *samplePhysical;
162 G4VPhysicalVolume* samplePhysical; <<
163 171
164 G4double sampleSafety = 0.0; << 172 G4double sampleSafety=0.0;
165 G4ThreeVector samplePoint; 173 G4ThreeVector samplePoint;
166 G4VSolid* ptrSolid = nullptr; << 174 G4VSolid* ptrSolid=0;
167 175
168 curNoVolumes = curVoxelNode->GetNoContained 176 curNoVolumes = curVoxelNode->GetNoContained();
169 177
170 for ( contentNo=curNoVolumes-1; contentNo>= 178 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
171 { 179 {
172 sampleNo = curVoxelNode->GetVolume((G4in << 180 sampleNo = curVoxelNode->GetVolume(contentNo);
173 if ( !fBlockList.IsBlocked(sampleNo) ) 181 if ( !fBlockList.IsBlocked(sampleNo) )
174 { 182 {
175 fBlockList.BlockVolume(sampleNo); 183 fBlockList.BlockVolume(sampleNo);
176 184
177 samplePhysical = fpMotherLogical->GetD 185 samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
178 G4AffineTransform sampleTf(samplePhysi 186 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
179 samplePhysi 187 samplePhysical->GetTranslation());
180 sampleTf.Invert(); 188 sampleTf.Invert();
181 samplePoint = sampleTf.TransformPoint( << 189 samplePoint = sampleTf.TransformPoint(localPoint);
182 ptrSolid = samplePhysical->GetLogicalV << 190 ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid();
183 191
184 sampleSafety = ptrSolid->DistanceToIn( 192 sampleSafety = ptrSolid->DistanceToIn(samplePoint);
185 ourSafety = std::min( sampleSafety, ou << 193 ourSafety = std::min( sampleSafety, ourSafety );
186 #ifdef G4VERBOSE 194 #ifdef G4VERBOSE
187 if(( fCheck ) && ( fVerbose == 1 )) 195 if(( fCheck ) && ( fVerbose == 1 ))
188 { 196 {
>> 197 // ReportSolidSafetyToIn( MethodName, solid, value, point );
189 G4cout << "*** G4VoxelSafety::Safety 198 G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
190 << " Invoked DistanceToIn( 199 << " Invoked DistanceToIn(p) for daughter solid: "
191 << ptrSolid->GetName() 200 << ptrSolid->GetName()
192 << ". Solid replied: " << sam 201 << ". Solid replied: " << sampleSafety << G4endl
193 << " For local point p: " 202 << " For local point p: " << samplePoint
194 << ", to be considered as 'da 203 << ", to be considered as 'daughter safety'." << G4endl;
195 } 204 }
196 #endif 205 #endif
197 } 206 }
198 } // end for contentNo 207 } // end for contentNo
199 208
200 return ourSafety; 209 return ourSafety;
201 } 210 }
202 211
203 // ******************************************* 212 // ********************************************************************
204 // SafetyForVoxelHeader 213 // SafetyForVoxelHeader
205 // 214 //
206 // Cycles through levels of headers to process 215 // Cycles through levels of headers to process each node level
207 // Obtained by modifying VoxelLocate (to cycle 216 // Obtained by modifying VoxelLocate (to cycle through Node Headers)
208 // ******************************************* 217 // *********************************************************************
209 // 218 //
210 G4double 219 G4double
211 G4VoxelSafety::SafetyForVoxelHeader( const G4S 220 G4VoxelSafety::SafetyForVoxelHeader( const G4SmartVoxelHeader* pHeader,
212 const G4T << 221 const G4ThreeVector& localPoint,
213 G4d << 222 G4double maxLength,
214 const G4V 223 const G4VPhysicalVolume& currentPhysical, //Debug
215 G4d << 224 G4double distUpperDepth_Sq,
216 G4d << 225 G4double previousMinSafety
217 ) 226 )
218 { 227 {
219 const G4SmartVoxelHeader* const targetVoxelH << 228 const G4SmartVoxelHeader * const targetVoxelHeader=pHeader;
220 G4SmartVoxelNode* targetVoxelNode = nullptr; << 229 G4SmartVoxelNode *targetVoxelNode=0;
221 230
222 const G4SmartVoxelProxy* sampleProxy; << 231 const G4SmartVoxelProxy *sampleProxy;
223 EAxis targetHeaderAxis; 232 EAxis targetHeaderAxis;
224 G4double targetHeaderMin, targetHeaderMax, t 233 G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
225 G4int targetHeaderNoSlices; 234 G4int targetHeaderNoSlices;
226 G4int targetNodeNo; 235 G4int targetNodeNo;
227 236
228 G4double minSafety = previousMinSafety; << 237 G4double minSafety= previousMinSafety;
229 G4double ourSafety = DBL_MAX; << 238 G4double ourSafety= DBL_MAX;
230 unsigned int checkedNum= 0; 239 unsigned int checkedNum= 0;
231 240
232 ++fVoxelDepth; << 241 fVoxelDepth++;
233 // fVoxelDepth set by ComputeSafety or prev 242 // fVoxelDepth set by ComputeSafety or previous level call
234 243
235 targetHeaderAxis = targetVoxelHeader->G 244 targetHeaderAxis = targetVoxelHeader->GetAxis();
236 targetHeaderNoSlices = (G4int)targetVoxelHe << 245 targetHeaderNoSlices = targetVoxelHeader->GetNoSlices();
237 targetHeaderMin = targetVoxelHeader->G 246 targetHeaderMin = targetVoxelHeader->GetMinExtent();
238 targetHeaderMax = targetVoxelHeader->G 247 targetHeaderMax = targetVoxelHeader->GetMaxExtent();
239 248
240 targetHeaderNodeWidth = (targetHeaderMax-tar 249 targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
241 / targetHeaderNoSlic 250 / targetHeaderNoSlices;
242 251
243 G4double localCrd = localPoint(targetHeaderA << 252 G4double localCrd= localPoint(targetHeaderAxis);
244 253
245 const auto candNodeNo = G4int( (localCrd-ta << 254 const G4int candNodeNo= G4int( (localCrd-targetHeaderMin)
246 / targetHeade << 255 / targetHeaderNodeWidth );
247 // Ensure that it is between 0 and targetHea 256 // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
248 257
249 #ifdef G4DEBUG_VOXELISATION 258 #ifdef G4DEBUG_VOXELISATION
250 if( candNodeNo < 0 || candNodeNo > targetHea 259 if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
251 { 260 {
252 G4ExceptionDescription ed; 261 G4ExceptionDescription ed;
253 ed << " Potential ERROR." 262 ed << " Potential ERROR."
254 << " Point is outside range of Voxel i 263 << " Point is outside range of Voxel in current coordinate" << G4endl;
255 ed << " Node number of point " << localP 264 ed << " Node number of point " << localPoint
256 << "is outside the range. " << G4endl; 265 << "is outside the range. " << G4endl;
257 ed << " Voxel node Num= " << candNodeN 266 ed << " Voxel node Num= " << candNodeNo << " versus minimum= " << 0
258 << " and maximum= " << targetHeaderNoS 267 << " and maximum= " << targetHeaderNoSlices-1 << G4endl;
259 ed << " Axis = " << targetHeaderAxis 268 ed << " Axis = " << targetHeaderAxis
260 << " No of slices = " << targetHeader 269 << " No of slices = " << targetHeaderNoSlices << G4endl;
261 ed << " Local coord = " << localCrd 270 ed << " Local coord = " << localCrd
262 << " Voxel Min = " << targetHeaderMin 271 << " Voxel Min = " << targetHeaderMin
263 << " Max = " << targetHeade 272 << " Max = " << targetHeaderMax << G4endl;
264 G4LogicalVolume *pLogical= currentPhysica 273 G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume();
265 ed << " Current volume (physical) = " << 274 ed << " Current volume (physical) = " << currentPhysical.GetName()
266 << " (logical) = " << pLogical->GetN 275 << " (logical) = " << pLogical->GetName() << G4endl;
267 G4VSolid* pSolid= pLogical->GetSolid(); 276 G4VSolid* pSolid= pLogical->GetSolid();
268 ed << " Solid type = " << pSolid->GetEnt 277 ed << " Solid type = " << pSolid->GetEntityType() << G4endl;
269 ed << *pSolid << G4endl; 278 ed << *pSolid << G4endl;
270 G4Exception("G4VoxelSafety::SafetyForVoxe 279 G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003",
271 JustWarning, ed, 280 JustWarning, ed,
272 "Point is outside range of Vo 281 "Point is outside range of Voxel in current coordinate");
273 } 282 }
274 #endif 283 #endif
275 284
276 const G4int pointNodeNo = 285 const G4int pointNodeNo =
277 std::max( 0, std::min( candNodeN 286 std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) );
278 287
279 #ifdef G4VERBOSE 288 #ifdef G4VERBOSE
280 if( fVerbose > 2 ) 289 if( fVerbose > 2 )
281 { 290 {
282 G4cout << G4endl; 291 G4cout << G4endl;
283 G4cout << "**** G4VoxelSafety::SafetyForVo 292 G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader " << G4endl;
284 G4cout << " Called at Depth = " << fV 293 G4cout << " Called at Depth = " << fVoxelDepth;
285 G4cout << " Calculated pointNodeNo= " << p 294 G4cout << " Calculated pointNodeNo= " << pointNodeNo
286 << " from position= " << localPoi 295 << " from position= " << localPoint(targetHeaderAxis)
287 << " min= " << targetHeaderMin 296 << " min= " << targetHeaderMin
288 << " max= " << targetHeaderMax 297 << " max= " << targetHeaderMax
289 << " width= " << targetHeaderNode 298 << " width= " << targetHeaderNodeWidth
290 << " no-slices= " << targetHeaderN 299 << " no-slices= " << targetHeaderNoSlices
291 << " axis= " << targetHeaderAxis 300 << " axis= " << targetHeaderAxis << G4endl;
292 } 301 }
293 else if (fVerbose == 1) 302 else if (fVerbose == 1)
294 { 303 {
295 G4cout << " VoxelSafety: Depth = " << fVo 304 G4cout << " VoxelSafety: Depth = " << fVoxelDepth
296 << " Number of Slices = " << target 305 << " Number of Slices = " << targetHeaderNoSlices
297 << " Header (address) = " << target 306 << " Header (address) = " << targetVoxelHeader << G4endl;
298 } 307 }
299 #endif 308 #endif
300 309
301 // Stack info for stepping 310 // Stack info for stepping
302 // 311 //
303 fVoxelAxisStack[fVoxelDepth] = targetHeaderA 312 fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
304 fVoxelNoSlicesStack[fVoxelDepth] = targetHea 313 fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
305 fVoxelSliceWidthStack[fVoxelDepth] = targetH 314 fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
306 315
307 fVoxelHeaderStack[fVoxelDepth] = pHeader; 316 fVoxelHeaderStack[fVoxelDepth] = pHeader;
308 317
309 G4int trialUp = -1, trialDown = -1; << 318 G4int trialUp= -1, trialDown= -1;
310 G4double distUp = DBL_MAX, distDown = DBL_MA << 319 G4double distUp= DBL_MAX, distDown= DBL_MAX;
311 320
312 // Using Equivalent voxels - this is pre-ini 321 // Using Equivalent voxels - this is pre-initialisation only
313 // 322 //
314 G4int nextUp = pointNodeNo+1; << 323 G4int nextUp= pointNodeNo+1;
315 G4int nextDown = pointNodeNo-1; << 324 G4int nextDown= pointNodeNo-1;
316 325
317 G4int nextNodeNo = pointNodeNo; << 326 G4int nextNodeNo= pointNodeNo;
318 G4double distAxis; // Distance in current A 327 G4double distAxis; // Distance in current Axis
319 distAxis = 0.0; // Starting in node contain << 328 distAxis= 0.0; // Starting in node containing local Coordinate
320 329
321 G4bool nextIsInside = false; << 330 G4bool nextIsInside= false;
322 331
323 G4double distMaxInterest= std::min( previous 332 G4double distMaxInterest= std::min( previousMinSafety, maxLength);
324 // We will not look beyond this distance. 333 // We will not look beyond this distance.
325 // This distance will be updated to reflec 334 // This distance will be updated to reflect the
326 // max ( minSafety, maxLength ) at each st 335 // max ( minSafety, maxLength ) at each step
327 336
328 targetNodeNo = pointNodeNo; << 337 targetNodeNo= pointNodeNo;
329 do 338 do
330 { 339 {
331 G4double nodeSafety = DBL_MAX, headerSafe << 340 G4double nodeSafety= DBL_MAX, headerSafety= DBL_MAX;
332 fVoxelNodeNoStack[fVoxelDepth] = targetNo 341 fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
333 342
334 ++checkedNum; << 343 checkedNum++;
335 344
336 sampleProxy = targetVoxelHeader->GetSlice 345 sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
337 346
338 #ifdef G4DEBUG_NAVIGATION 347 #ifdef G4DEBUG_NAVIGATION
339 assert( sampleProxy != 0); 348 assert( sampleProxy != 0);
340 if( fVerbose > 2 ) 349 if( fVerbose > 2 )
341 { 350 {
342 G4cout << " -Checking node " << targetN 351 G4cout << " -Checking node " << targetNodeNo
343 << " is proxy with address " << 352 << " is proxy with address " << sampleProxy << G4endl;
344 } 353 }
345 #endif 354 #endif
346 355
347 if ( sampleProxy == nullptr ) << 356 if ( sampleProxy == 0 )
348 { 357 {
349 G4ExceptionDescription ed; 358 G4ExceptionDescription ed;
350 ed << " Problem for node number= " << t 359 ed << " Problem for node number= " << targetNodeNo
351 << " Number of slides = " << targ 360 << " Number of slides = " << targetHeaderNoSlices
352 << G4endl; 361 << G4endl;
353 G4Exception( "G4VoxelSafety::SafetyForV 362 G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003",
354 FatalException, ed, 363 FatalException, ed,
355 "Problem sampleProxy is Ze 364 "Problem sampleProxy is Zero. Failure in loop.");
356 } 365 }
357 else if ( sampleProxy->IsNode() ) 366 else if ( sampleProxy->IsNode() )
358 { 367 {
359 targetVoxelNode = sampleProxy->GetNode( 368 targetVoxelNode = sampleProxy->GetNode();
360 369
361 // Deal with the node here [ i.e. the l 370 // Deal with the node here [ i.e. the last level ]
362 // 371 //
363 nodeSafety= SafetyForVoxelNode( targetV 372 nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
364 #ifdef G4DEBUG_NAVIGATION 373 #ifdef G4DEBUG_NAVIGATION
365 if( fVerbose > 2 ) 374 if( fVerbose > 2 )
366 { 375 {
367 G4cout << " -- It is a Node "; 376 G4cout << " -- It is a Node ";
368 G4cout << " its safety= " << node 377 G4cout << " its safety= " << nodeSafety
369 << " our level Saf = " << our 378 << " our level Saf = " << ourSafety
370 << " MinSaf= " << minSafety 379 << " MinSaf= " << minSafety << G4endl;
371 } 380 }
372 #endif 381 #endif
373 ourSafety= std::min( ourSafety, nodeSa 382 ourSafety= std::min( ourSafety, nodeSafety );
374 383
375 trialUp = targetVoxelNode->GetMaxEquiv << 384 trialUp = targetVoxelNode->GetMaxEquivalentSliceNo()+1;
376 trialDown = targetVoxelNode->GetMinEqu 385 trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1;
377 } 386 }
378 else 387 else
379 { 388 {
380 const G4SmartVoxelHeader* pNewVoxelHea << 389 const G4SmartVoxelHeader *pNewVoxelHeader = sampleProxy->GetHeader();
381 390
382 G4double distCombined_Sq; 391 G4double distCombined_Sq;
383 distCombined_Sq = distUpperDepth_Sq + 392 distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
384 393
385 #ifdef G4DEBUG_NAVIGATION 394 #ifdef G4DEBUG_NAVIGATION
386 if( fVerbose > 2 ) 395 if( fVerbose > 2 )
387 { 396 {
388 G4double distCombined= std::sqrt( di 397 G4double distCombined= std::sqrt( distCombined_Sq );
389 G4double distUpperDepth= std::sqrt ( 398 G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
390 G4cout << " -- It is a Header " << G 399 G4cout << " -- It is a Header " << G4endl;
391 G4cout << " Recurse to deal with ne 400 G4cout << " Recurse to deal with next level, fVoxelDepth= "
392 << fVoxelDepth+1 << G4endl; 401 << fVoxelDepth+1 << G4endl;
393 G4cout << " Distances: Upper= " << 402 G4cout << " Distances: Upper= " << distUpperDepth
394 << " Axis= " << distAxis 403 << " Axis= " << distAxis
395 << " Combined= " << distCombi 404 << " Combined= " << distCombined << G4endl;
396 } 405 }
397 #endif 406 #endif
398 407
399 // Recurse to deal with lower levels 408 // Recurse to deal with lower levels
400 // 409 //
401 headerSafety= SafetyForVoxelHeader( pN 410 headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint,
402 ma 411 maxLength, currentPhysical,
403 di 412 distCombined_Sq, minSafety);
404 ourSafety = std::min( ourSafety, heade << 413 ourSafety= std::min( ourSafety, headerSafety );
405 414
406 #ifdef G4DEBUG_NAVIGATION 415 #ifdef G4DEBUG_NAVIGATION
407 if( fVerbose > 2 ) 416 if( fVerbose > 2 )
408 { 417 {
409 G4cout << " >> Header safety = 418 G4cout << " >> Header safety = " << headerSafety
410 << " our level Saf = " << our 419 << " our level Saf = " << ourSafety << G4endl;
411 } 420 }
412 #endif 421 #endif
413 trialUp = pNewVoxelHeader->GetMaxEqu 422 trialUp = pNewVoxelHeader->GetMaxEquivalentSliceNo()+1;
414 trialDown = pNewVoxelHeader->GetMinEqu 423 trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1;
415 } 424 }
416 minSafety = std::min( minSafety, ourSafet << 425 minSafety= std::min( minSafety, ourSafety );
417 426
418 // Find next closest Voxel 427 // Find next closest Voxel
419 // - first try: by simple subtraction 428 // - first try: by simple subtraction
420 // - later: using distance (TODO - t 429 // - later: using distance (TODO - tbc)
421 // 430 //
422 if( targetNodeNo >= pointNodeNo ) 431 if( targetNodeNo >= pointNodeNo )
423 { 432 {
424 nextUp = trialUp; << 433 nextUp = trialUp;
425 // distUp = std::max( targetHeaderMa 434 // distUp = std::max( targetHeaderMax-localCrd, 0.0 );
426 G4double lowerEdgeOfNext = targetHeade 435 G4double lowerEdgeOfNext = targetHeaderMin
427 + nextUp * ta 436 + nextUp * targetHeaderNodeWidth;
428 distUp = lowerEdgeOfNext-localCrd ; 437 distUp = lowerEdgeOfNext-localCrd ;
429 if( distUp < 0.0 ) 438 if( distUp < 0.0 )
430 { 439 {
431 distUp = DBL_MAX; // On the wrong 440 distUp = DBL_MAX; // On the wrong side - must not be considered
432 } 441 }
433 #ifdef G4DEBUG_NAVIGATION 442 #ifdef G4DEBUG_NAVIGATION
434 if( fVerbose > 2 ) 443 if( fVerbose > 2 )
435 { 444 {
436 G4cout << " > Updated nextUp= " << ne 445 G4cout << " > Updated nextUp= " << nextUp << G4endl;
437 } 446 }
438 #endif 447 #endif
439 } 448 }
440 449
441 if( targetNodeNo <= pointNodeNo ) 450 if( targetNodeNo <= pointNodeNo )
442 { 451 {
443 nextDown = trialDown; 452 nextDown = trialDown;
444 // distDown = std::max( localCrd-targ 453 // distDown = std::max( localCrd-targetHeaderMin, 0.0);
445 G4double upperEdgeOfNext = targetHead 454 G4double upperEdgeOfNext = targetHeaderMin
446 + (nextDown+ 455 + (nextDown+1) * targetHeaderNodeWidth;
447 distDown = localCrd-upperEdgeOfNext; 456 distDown = localCrd-upperEdgeOfNext;
448 if( distDown < 0.0 ) 457 if( distDown < 0.0 )
449 { 458 {
450 distDown= DBL_MAX; // On the wrong 459 distDown= DBL_MAX; // On the wrong side - must not be considered
451 } 460 }
452 #ifdef G4DEBUG_NAVIGATION 461 #ifdef G4DEBUG_NAVIGATION
453 if( fVerbose > 2 ) 462 if( fVerbose > 2 )
454 { 463 {
455 G4cout << " > Updated nextDown= " << 464 G4cout << " > Updated nextDown= " << nextDown << G4endl;
456 } 465 }
457 #endif 466 #endif
458 } 467 }
459 468
460 #ifdef G4DEBUG_NAVIGATION 469 #ifdef G4DEBUG_NAVIGATION
461 if( fVerbose > 2 ) 470 if( fVerbose > 2 )
462 { 471 {
463 G4cout << " Node= " << pointNodeNo 472 G4cout << " Node= " << pointNodeNo
464 << " Up: next= " << nextUp << 473 << " Up: next= " << nextUp << " d# "
465 << nextUp - pointNodeNo 474 << nextUp - pointNodeNo
466 << " trialUp= " << trialUp << " 475 << " trialUp= " << trialUp << " d# "
467 << trialUp - pointNodeNo 476 << trialUp - pointNodeNo
468 << " Down: next= " << nextDown < 477 << " Down: next= " << nextDown << " d# "
469 << targetNodeNo - nextDown 478 << targetNodeNo - nextDown
470 << " trialDwn= " << trialDown << 479 << " trialDwn= " << trialDown << " d# "
471 << targetNodeNo - trialDown 480 << targetNodeNo - trialDown
472 << " condition (next is Inside)= 481 << " condition (next is Inside)= " << nextIsInside
473 << G4endl; 482 << G4endl;
474 } 483 }
475 #endif 484 #endif
476 485
477 G4bool UpIsClosest; 486 G4bool UpIsClosest;
478 UpIsClosest = distUp < distDown; << 487 UpIsClosest= distUp < distDown;
479 488
480 if( (nextUp < targetHeaderNoSlices) 489 if( (nextUp < targetHeaderNoSlices)
481 && (UpIsClosest || (nextDown < 0)) ) 490 && (UpIsClosest || (nextDown < 0)) )
482 { 491 {
483 nextNodeNo = nextUp; << 492 nextNodeNo=nextUp;
484 distAxis = distUp; 493 distAxis = distUp;
485 ++nextUp; // Default 494 ++nextUp; // Default
486 #ifdef G4VERBOSE 495 #ifdef G4VERBOSE
487 if( fVerbose > 2 ) 496 if( fVerbose > 2 )
488 { 497 {
489 G4cout << " > Chose Up. Depth= " < 498 G4cout << " > Chose Up. Depth= " << fVoxelDepth
490 << " Nodes: next= " << next 499 << " Nodes: next= " << nextNodeNo
491 << " new nextUp= " << nextU 500 << " new nextUp= " << nextUp
492 << " Dist = " << distAxis << 501 << " Dist = " << distAxis << G4endl;
493 } 502 }
494 #endif 503 #endif
495 } 504 }
496 else 505 else
497 { 506 {
498 nextNodeNo = nextDown; << 507 nextNodeNo=nextDown;
499 distAxis = distDown; 508 distAxis = distDown;
500 --nextDown; // A default value 509 --nextDown; // A default value
501 #ifdef G4VERBOSE 510 #ifdef G4VERBOSE
502 if( fVerbose > 2 ) 511 if( fVerbose > 2 )
503 { 512 {
504 G4cout << " > Chose Down. Depth= " < 513 G4cout << " > Chose Down. Depth= " << fVoxelDepth
505 << " Nodes: next= " << nextNo 514 << " Nodes: next= " << nextNodeNo
506 << " new nextDown= " << nextD 515 << " new nextDown= " << nextDown
507 << " Dist = " << distAxis << G 516 << " Dist = " << distAxis << G4endl;
508 } 517 }
509 #endif 518 #endif
510 } 519 }
511 520
512 nextIsInside = (nextNodeNo >= 0) && (next 521 nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
513 if( nextIsInside ) 522 if( nextIsInside )
514 { 523 {
515 targetNodeNo= nextNodeNo; 524 targetNodeNo= nextNodeNo;
516 525
517 #ifdef G4DEBUG_NAVIGATION 526 #ifdef G4DEBUG_NAVIGATION
518 assert( targetVoxelHeader->GetSlice(nex 527 assert( targetVoxelHeader->GetSlice(nextNodeNo) != 0 );
519 G4bool bContinue = (distAxis<minSafety) << 528 G4bool bContinue= (distAxis<minSafety);
520 if( !bContinue ) 529 if( !bContinue )
521 { 530 {
522 if( fVerbose > 2 ) 531 if( fVerbose > 2 )
523 { 532 {
524 G4cout << " Can skip remaining at d 533 G4cout << " Can skip remaining at depth " << targetHeaderAxis
525 << " >> distAxis= " << dist 534 << " >> distAxis= " << distAxis
526 << " minSafety= " << minSafe 535 << " minSafety= " << minSafety << G4endl;
527 } 536 }
528 } 537 }
529 #endif 538 #endif
530 } 539 }
531 else 540 else
532 { 541 {
533 #ifdef G4DEBUG_NAVIGATION 542 #ifdef G4DEBUG_NAVIGATION
534 if( fVerbose > 2) 543 if( fVerbose > 2)
535 { 544 {
536 G4cout << " VoxSaf> depth= " << fVoxe 545 G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl;
537 G4cout << " VoxSaf> No more candidate 546 G4cout << " VoxSaf> No more candidates: nodeDown= " << nextDown
538 << " nodeUp= " << nextUp 547 << " nodeUp= " << nextUp
539 << " lastSlice= " << targetHea 548 << " lastSlice= " << targetHeaderNoSlices << G4endl;
540 } 549 }
541 #endif 550 #endif
542 } 551 }
543 552
544 // This calculation can be 'hauled'-up to 553 // This calculation can be 'hauled'-up to where minSafety is calculated
545 // 554 //
546 distMaxInterest = std::min( minSafety, ma 555 distMaxInterest = std::min( minSafety, maxLength );
547 556
548 } while ( nextIsInside && ( distAxis*distAxi 557 } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
549 < distMaxInterest* 558 < distMaxInterest*distMaxInterest ) );
550 559
551 #ifdef G4VERBOSE 560 #ifdef G4VERBOSE
552 if( fVerbose > 0 ) 561 if( fVerbose > 0 )
553 { 562 {
554 G4cout << " Ended for targetNodeNo -- chec 563 G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of "
555 << targetHeaderNoSlices << " slices 564 << targetHeaderNoSlices << " slices." << G4endl;
556 G4cout << " ===== Returning from SafetyFor 565 G4cout << " ===== Returning from SafetyForVoxelHeader "
557 << " Depth= " << fVoxelDepth << G4 566 << " Depth= " << fVoxelDepth << G4endl
558 << G4endl; 567 << G4endl;
559 } 568 }
560 #endif 569 #endif
561 570
562 // Go back one level 571 // Go back one level
563 fVoxelDepth--; 572 fVoxelDepth--;
564 573
565 return ourSafety; 574 return ourSafety;
566 } 575 }
567 576