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