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