Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4VoxelNavigation.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// class G4VoxelNavigation Implementation
27//
28// Author: P.Kent, 1996
29//
30// --------------------------------------------------------------------
31#include "G4VoxelNavigation.hh"
33#include "G4VoxelSafety.hh"
34
36
37#include <cassert>
38#include <ostream>
39
40// ********************************************************************
41// Constructor
42// ********************************************************************
43//
45 : fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
46 fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
47 fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
48 fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
49 fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)nullptr)
50{
51 fLogger= new G4NavigationLogger("G4VoxelNavigation");
54
55#ifdef G4DEBUG_NAVIGATION
56 SetVerboseLevel(5); // Reports most about daughter volumes
57#endif
58}
59
60// ********************************************************************
61// Destructor
62// ********************************************************************
63//
69
70// --------------------------------------------------------------------------
71// Input:
72// exiting: : last step exited
73// blockedPhysical : phys volume last exited (if exiting)
74// blockedReplicaNo : copy/replica number of exited
75// Output:
76// entering : if true, found candidate volume to enter
77// blockedPhysical : candidate phys volume to enter - if entering
78// blockedReplicaNo : copy/replica number - if entering
79// exiting: : will exit current (mother) volume
80// In/Out
81// --------------------------------------------------------------------------
82
83// ********************************************************************
84// ComputeStep
85// ********************************************************************
86//
89 const G4ThreeVector& localDirection,
90 const G4double currentProposedStepLength,
91 G4double& newSafety,
92 /* const */ G4NavigationHistory& history,
93 G4bool& validExitNormal,
94 G4ThreeVector& exitNormal,
95 G4bool& exiting,
96 G4bool& entering,
97 G4VPhysicalVolume* (*pBlockedPhysical),
98 G4int& blockedReplicaNo )
99{
100 G4VPhysicalVolume *motherPhysical, *samplePhysical, *blockedExitedVol=nullptr;
101 G4LogicalVolume *motherLogical;
102 G4VSolid *motherSolid;
103 G4ThreeVector sampleDirection;
104 G4double ourStep=currentProposedStepLength, ourSafety;
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->GetLogicalVolume();
115 motherSolid = motherLogical->GetSolid();
116
117 //
118 // Compute mother safety
119 //
120
121 motherSafety = motherSolid->DistanceToOut(localPoint);
122 ourSafety = motherSafety; // Working isotropic safety
123
124#ifdef G4VERBOSE
125 if ( fCheck )
126 {
127 fLogger->PreComputeStepLog (motherPhysical, motherSafety, localPoint);
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)>=kMinExitingNormalCosine )
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 Mother early !!
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 answer of daughters!
159 motherStep = motherSolid->DistanceToOut(localPoint,
160 localDirection,
161 true,
162 &motherValidExitNormal,
163 &motherExitNormal);
164 }
165#endif
166
167 localNoDaughters = (G4int)motherLogical->GetNoDaughters();
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->GetNoContained();
179 for (contentNo=curNoVolumes-1; contentNo>=0; contentNo--)
180 {
181 sampleNo = curVoxelNode->GetVolume((G4int)contentNo);
182 if ( !fBList.IsBlocked(sampleNo) )
183 {
184 fBList.BlockVolume(sampleNo);
185 samplePhysical = motherLogical->GetDaughter(sampleNo);
186 if ( samplePhysical!=blockedExitedVol )
187 {
188 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
189 samplePhysical->GetTranslation());
190 sampleTf.Invert();
191 const G4ThreeVector samplePoint =
192 sampleTf.TransformPoint(localPoint);
193 const G4VSolid *sampleSolid =
194 samplePhysical->GetLogicalVolume()->GetSolid();
195 const G4double sampleSafety =
196 sampleSolid->DistanceToIn(samplePoint);
197
198 if ( sampleSafety<ourSafety )
199 {
200 ourSafety = sampleSafety;
201 }
202 if ( sampleSafety<=ourStep )
203 {
204 sampleDirection = sampleTf.TransformAxis(localDirection);
205 G4double sampleStep =
206 sampleSolid->DistanceToIn(samplePoint, sampleDirection);
207#ifdef G4VERBOSE
208 if( fCheck )
209 {
210 fLogger->PrintDaughterLog(sampleSolid, samplePoint,
211 sampleSafety, true,
212 sampleDirection, sampleStep);
213 }
214#endif
215 if ( sampleStep<=ourStep )
216 {
217 ourStep = sampleStep;
218 entering = true;
219 exiting = false;
220 *pBlockedPhysical = samplePhysical;
221 blockedReplicaNo = -1;
222#ifdef G4VERBOSE
223 // Check to see that the resulting point is indeed in/on volume.
224 // This could be done only for successful candidate.
225 if ( fCheck )
226 {
227 fLogger->AlongComputeStepLog (sampleSolid, samplePoint,
228 sampleDirection, localDirection, sampleSafety, sampleStep);
229 }
230#endif
231 }
232#ifdef G4VERBOSE
233 if ( fCheck && ( sampleStep < kInfinity )
234 && ( sampleStep >= motherStep ) )
235 {
236 // The intersection point with the daughter is after the exit
237 // point from the mother volume. Double check this !!
238 fLogger->CheckDaughterEntryPoint(sampleSolid,
239 samplePoint, sampleDirection,
240 motherSolid,
241 localPoint, localDirection,
242 motherStep, sampleStep);
243 }
244#endif
245 }
246#ifdef G4VERBOSE
247 else // ie if sampleSafety > outStep
248 {
249 if( fCheck )
250 {
251 fLogger->PrintDaughterLog(sampleSolid, samplePoint,
252 sampleSafety, false,
253 G4ThreeVector(0.,0.,0.), -1.0 );
254 }
255 }
256#endif
257 }
258 }
259 }
260 if (initialNode)
261 {
262 initialNode = false;
263 voxelSafety = ComputeVoxelSafety(localPoint);
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 required
282 //
283 if ( motherSafety<=ourStep )
284 {
285 // In case of check mode this is a duplicate call -- acceptable
286 motherStep = motherSolid->DistanceToOut(localPoint, localDirection,
287 true, &motherValidExitNormal, &motherExitNormal);
288#ifdef G4VERBOSE
289 if ( fCheck )
290 {
291 fLogger->PostComputeStepLog(motherSolid, localPoint, localDirection,
292 motherStep, motherSafety);
293 if( motherValidExitNormal )
294 {
295 fLogger->CheckAndReportBadNormal(motherExitNormal,
296 localPoint, localDirection,
297 motherStep, motherSolid,
298 "From motherSolid::DistanceToOut" );
299 }
300 }
301#endif
302 if( (motherStep >= kInfinity) || (motherStep < 0.0) )
303 {
304#ifdef G4VERBOSE
305 if( fCheck ) // Error - indication of being outside solid !!
306 {
307 fLogger->ReportOutsideMother(localPoint, localDirection,
308 motherPhysical);
309 }
310#endif
311 motherStep = 0.0;
312 ourStep = 0.0;
313 exiting = true;
314 entering = false;
315
316 // validExitNormal= motherValidExitNormal;
317 // exitNormal= motherExitNormal;
318 // Useful only if the point is very close to surface
319 // => but it would need to be rotated to grand-mother ref frame !
320 validExitNormal= false;
321
322 *pBlockedPhysical = nullptr; // or motherPhysical ?
323 blockedReplicaNo = 0; // or motherReplicaNumber ?
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 to set these;confirmed short step
336 //
337 validExitNormal = motherValidExitNormal;
338 exitNormal = motherExitNormal;
339
340 if ( validExitNormal )
341 {
342 const G4RotationMatrix *rot = motherPhysical->GetRotation();
343 if (rot != nullptr)
344 {
345 exitNormal *= rot->inverse();
346#ifdef G4VERBOSE
347 if( fCheck )
348 {
349 fLogger->CheckAndReportBadNormal(exitNormal, // rotated
350 motherExitNormal, // original
351 *rot,
352 "From RotationMatrix" );
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, localDirection, ourStep);
369 }
370 } // end -while (noStep)- loop
371
372 return ourStep;
373}
374
375// ********************************************************************
376// ComputeVoxelSafety
377//
378// Computes safety from specified point to voxel boundaries
379// using already located point
380// o collected boundaries for most derived level
381// o adjacent boundaries for previous levels
382// ********************************************************************
383//
386{
387 G4SmartVoxelHeader *curHeader;
388 G4double voxelSafety, curNodeWidth;
389 G4double curNodeOffset, minCurCommonDelta, maxCurCommonDelta;
390 G4int minCurNodeNoDelta, maxCurNodeNoDelta;
391 G4int localVoxelDepth, curNodeNo;
392 EAxis curHeaderAxis;
393
394 localVoxelDepth = fVoxelDepth;
395
396 curHeader = fVoxelHeaderStack[localVoxelDepth];
397 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
398 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
399 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
400
401 // Compute linear intersection distance to boundaries of max/min
402 // to collected nodes at current level
403 //
404 curNodeOffset = curNodeNo*curNodeWidth;
405 maxCurNodeNoDelta = fVoxelNode->GetMaxEquivalentSliceNo()-curNodeNo;
406 minCurNodeNoDelta = curNodeNo-fVoxelNode->GetMinEquivalentSliceNo();
407 minCurCommonDelta = localPoint(curHeaderAxis)
408 - curHeader->GetMinExtent() - curNodeOffset;
409 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
410
411 if ( minCurNodeNoDelta<maxCurNodeNoDelta )
412 {
413 voxelSafety = minCurNodeNoDelta*curNodeWidth;
414 voxelSafety += minCurCommonDelta;
415 }
416 else if (maxCurNodeNoDelta < minCurNodeNoDelta)
417 {
418 voxelSafety = maxCurNodeNoDelta*curNodeWidth;
419 voxelSafety += maxCurCommonDelta;
420 }
421 else // (maxCurNodeNoDelta == minCurNodeNoDelta)
422 {
423 voxelSafety = minCurNodeNoDelta*curNodeWidth;
424 voxelSafety += std::min(minCurCommonDelta,maxCurCommonDelta);
425 }
426
427 // Compute isotropic safety to boundaries of previous levels
428 // [NOT to collected boundaries]
429
430 // Loop checking, 07.10.2016, JA
431 while ( (localVoxelDepth>0) && (voxelSafety>0) )
432 {
433 localVoxelDepth--;
434 curHeader = fVoxelHeaderStack[localVoxelDepth];
435 curHeaderAxis = fVoxelAxisStack[localVoxelDepth];
436 curNodeNo = fVoxelNodeNoStack[localVoxelDepth];
437 curNodeWidth = fVoxelSliceWidthStack[localVoxelDepth];
438 curNodeOffset = curNodeNo*curNodeWidth;
439 minCurCommonDelta = localPoint(curHeaderAxis)
440 - curHeader->GetMinExtent() - curNodeOffset;
441 maxCurCommonDelta = curNodeWidth-minCurCommonDelta;
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 and point
464// in the specified direction
465//
466// Returns false if all voxels considered
467// [current Step ends inside same voxel or leaves all voxels]
468// true otherwise
469// [the information on the next voxel is put into the set of
470// fVoxel* variables & "stacks"]
471// ********************************************************************
472//
473G4bool
475 const G4ThreeVector& localDirection,
476 const G4double currentStep)
477{
478 G4SmartVoxelHeader *workHeader=nullptr, *newHeader=nullptr;
479 G4SmartVoxelProxy *newProxy=nullptr;
480 G4SmartVoxelNode *newVoxelNode=nullptr;
481 G4ThreeVector targetPoint, voxelPoint;
482 G4double workNodeWidth, workMinExtent, workCoord;
483 G4double minVal, maxVal, newDistance=0.;
484 G4double newHeaderMin, newHeaderNodeWidth;
485 G4int depth=0, newDepth=0, workNodeNo=0, newNodeNo=0, newHeaderNoSlices=0;
486 EAxis workHeaderAxis, newHeaderAxis;
487 G4bool isNewVoxel = false;
488
489 G4double currentDistance = currentStep;
490
491 // Determine if end of Step within current voxel
492 //
493 for (depth=0; depth<fVoxelDepth; ++depth)
494 {
495 targetPoint = localPoint+localDirection*currentDistance;
496 newDistance = currentDistance;
497 workHeader = fVoxelHeaderStack[depth];
498 workHeaderAxis = fVoxelAxisStack[depth];
499 workNodeNo = fVoxelNodeNoStack[depth];
500 workNodeWidth = fVoxelSliceWidthStack[depth];
501 workMinExtent = workHeader->GetMinExtent();
502 workCoord = targetPoint(workHeaderAxis);
503 minVal = workMinExtent+workNodeNo*workNodeWidth;
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(workHeaderAxis))
515 / localDirection(workHeaderAxis);
516 isNewVoxel = true;
517 newDepth = depth;
518 }
519 }
520 else
521 {
522 newNodeNo = workNodeNo-1;
523 newHeader = workHeader;
524 newDistance = (minVal-localPoint(workHeaderAxis))
525 / localDirection(workHeaderAxis);
526 isNewVoxel = true;
527 newDepth = depth;
528 }
529 currentDistance = newDistance;
530 }
531 targetPoint = localPoint+localDirection*currentDistance;
532
533 // Check if end of Step within collected boundaries of current voxel
534 //
535 depth = fVoxelDepth;
536 {
537 workHeader = fVoxelHeaderStack[depth];
538 workHeaderAxis = fVoxelAxisStack[depth];
539 workNodeNo = fVoxelNodeNoStack[depth];
540 workNodeWidth = fVoxelSliceWidthStack[depth];
541 workMinExtent = workHeader->GetMinExtent();
542 workCoord = targetPoint(workHeaderAxis);
543 minVal = workMinExtent+fVoxelNode->GetMinEquivalentSliceNo()*workNodeWidth;
544
545 if ( minVal<=workCoord+fHalfTolerance )
546 {
547 maxVal = workMinExtent+(fVoxelNode->GetMaxEquivalentSliceNo()+1)
548 *workNodeWidth;
549 if ( maxVal<=workCoord-fHalfTolerance )
550 {
551 newNodeNo = fVoxelNode->GetMaxEquivalentSliceNo()+1;
552 newHeader = workHeader;
553 newDistance = (maxVal-localPoint(workHeaderAxis))
554 / localDirection(workHeaderAxis);
555 isNewVoxel = true;
556 newDepth = depth;
557 }
558 }
559 else
560 {
561 newNodeNo = fVoxelNode->GetMinEquivalentSliceNo()-1;
562 newHeader = workHeader;
563 newDistance = (minVal-localPoint(workHeaderAxis))
564 / localDirection(workHeaderAxis);
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 voxel boundary
576 // newHeader=Header for crossed voxel
577 // newDistance=distance to crossed voxel boundary (along the track)
578 //
579 if ( (newNodeNo<0) || (newNodeNo>=G4int(newHeader->GetNoSlices())))
580 {
581 // Leaving mother volume
582 //
583 isNewVoxel = false;
584 }
585 else
586 {
587 // Compute intersection point on the least refined
588 // voxel boundary that is hit
589 //
590 voxelPoint = localPoint+localDirection*newDistance;
591 fVoxelNodeNoStack[newDepth] = newNodeNo;
592 fVoxelDepth = newDepth;
593 newVoxelNode = nullptr;
594 while ( newVoxelNode == nullptr )
595 {
596 newProxy = newHeader->GetSlice(newNodeNo);
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->GetNoSlices();
607 newHeaderMin = newHeader->GetMinExtent();
608 newHeaderNodeWidth = (newHeader->GetMaxExtent()-newHeaderMin)
609 / newHeaderNoSlices;
610 newNodeNo = G4int( (voxelPoint(newHeaderAxis)-newHeaderMin)
611 / newHeaderNodeWidth );
612 // Rounding protection
613 //
614 if ( newNodeNo<0 )
615 {
616 newNodeNo=0;
617 }
618 else if ( newNodeNo>=newHeaderNoSlices )
619 {
620 newNodeNo = newHeaderNoSlices-1;
621 }
622 // Stack info for stepping
623 //
624 fVoxelAxisStack[fVoxelDepth] = newHeaderAxis;
625 fVoxelNoSlicesStack[fVoxelDepth] = newHeaderNoSlices;
626 fVoxelSliceWidthStack[fVoxelDepth] = newHeaderNodeWidth;
627 fVoxelNodeNoStack[fVoxelDepth] = newNodeNo;
628 fVoxelHeaderStack[fVoxelDepth] = newHeader;
629 }
630 }
631 fVoxelNode = newVoxelNode;
632 }
633 }
634 return isNewVoxel;
635}
636
637// ********************************************************************
638// ComputeSafety
639//
640// Calculates the isotropic distance to the nearest boundary from the
641// specified point in the local coordinate system.
642// The localpoint utilised must be within the current volume.
643// ********************************************************************
644//
647 const G4NavigationHistory& history,
648 const G4double maxLength)
649{
650 G4VPhysicalVolume *motherPhysical, *samplePhysical;
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->GetLogicalVolume();
661 motherSolid = motherLogical->GetSolid();
662
663 if( fBestSafety )
664 {
665 return fpVoxelSafety->ComputeSafety( localPoint,*motherPhysical,maxLength );
666 }
667
668 //
669 // Compute mother safety
670 //
671
672 motherSafety = motherSolid->DistanceToOut(localPoint);
673 ourSafety = motherSafety; // Working isotropic safety
674
675 if( motherSafety == 0.0 )
676 {
677#ifdef G4DEBUG_NAVIGATION
678 // Check that point is inside mother volume
679 EInside insideMother = motherSolid->Inside(localPoint);
680
681 if( insideMother == kOutside )
682 {
684 message << "Safety method called for location outside current Volume." << G4endl
685 << "Location for safety is Outside this volume. " << G4endl
686 << "The approximate distance to the solid "
687 << "(safety from outside) is: "
688 << motherSolid->DistanceToIn( localPoint ) << G4endl;
689 message << " Problem occurred with physical volume: "
690 << " Name: " << motherPhysical->GetName()
691 << " Copy No: " << motherPhysical->GetCopyNo() << G4endl
692 << " Local Point = " << localPoint << G4endl;
693 message << " Description of solid: " << G4endl
694 << *motherSolid << G4endl;
695 G4Exception("G4VoxelNavigation::ComputeSafety()", "GeomNav0003",
696 JustWarning, message);
697 }
698
699 // Following check is NOT for an issue - it is only for information
700 // It is allowed that a solid gives approximate safety - even zero.
701 //
702 if( insideMother == kInside ) // && fVerbose )
703 {
704 G4ExceptionDescription messageIn;
705
706 messageIn << " Point is Inside, but safety is Zero ." << G4endl;
707 messageIn << " Inexact safety for volume " << motherPhysical->GetName() << G4endl
708 << " Solid: Name= " << motherSolid->GetName()
709 << " Type= " << motherSolid->GetEntityType() << G4endl;
710 messageIn << " Local point= " << localPoint << G4endl;
711 messageIn << " Solid parameters: " << G4endl << *motherSolid << G4endl;
712 G4Exception("G4VoxelNavigation::ComputeSafety()", "GeomNav0003",
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,localPoint,motherSafety,true,1);
724 }
725#endif
726 //
727 // Compute daughter safeties
728 //
729 // Look only inside the current Voxel only (in the first version).
730 //
731 curVoxelNode = fVoxelNode;
732 curNoVolumes = curVoxelNode->GetNoContained();
733
734 for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
735 {
736 sampleNo = curVoxelNode->GetVolume((G4int)contentNo);
737 samplePhysical = motherLogical->GetDaughter(sampleNo);
738
739 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
740 samplePhysical->GetTranslation());
741 sampleTf.Invert();
742 const G4ThreeVector samplePoint = sampleTf.TransformPoint(localPoint);
743 const G4VSolid* sampleSolid= samplePhysical->GetLogicalVolume()->GetSolid();
744 G4double sampleSafety = sampleSolid->DistanceToIn(samplePoint);
745 if ( sampleSafety<ourSafety )
746 {
747 ourSafety = sampleSafety;
748 }
749#ifdef G4VERBOSE
750 if( fCheck )
751 {
752 fLogger->ComputeSafetyLog(sampleSolid, samplePoint,
753 sampleSafety, false, 0);
754 }
755#endif
756 }
757 voxelSafety = ComputeVoxelSafety(localPoint);
758 if ( voxelSafety<ourSafety )
759 {
760 ourSafety = voxelSafety;
761 }
762 return ourSafety;
763}
764
766 const G4ThreeVector& localPoint )
767{
768 auto motherLogical = motherPhysical->GetLogicalVolume();
769
770 assert(motherLogical != nullptr);
771
772 if ( auto pVoxelHeader = motherLogical->GetVoxelHeader() )
773 VoxelLocate( pVoxelHeader, localPoint );
774}
775
776// ********************************************************************
777// SetVerboseLevel
778// ********************************************************************
779//
781{
782 if( fLogger != nullptr ) { fLogger->SetVerboseLevel(level); }
783 if( fpVoxelSafety != nullptr) { fpVoxelSafety->SetVerboseLevel(level); }
784}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
double G4double
Definition G4Types.hh:83
long G4long
Definition G4Types.hh:87
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
#define G4endl
Definition G4ios.hh:67
double dot(const Hep3Vector &) const
HepRotation inverse() const
G4AffineTransform & Invert()
G4ThreeVector TransformPoint(const G4ThreeVector &vec) const
G4ThreeVector TransformAxis(const G4ThreeVector &axis) const
void BlockVolume(const G4int v)
void Enlarge(const G4int nv)
G4bool IsBlocked(const G4int v) const
G4double GetSurfaceTolerance() const
static G4GeometryTolerance * GetInstance()
G4VSolid * GetSolid() const
std::size_t GetNoDaughters() const
G4VPhysicalVolume * GetDaughter(const std::size_t i) const
G4VPhysicalVolume * GetTopVolume() const
void SetVerboseLevel(G4int level)
void PreComputeStepLog(const G4VPhysicalVolume *motherPhysical, G4double motherSafety, const G4ThreeVector &localPoint) const
void CheckDaughterEntryPoint(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, const G4ThreeVector &sampleDirection, const G4VSolid *motherSolid, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double motherStep, G4double sampleStep) const
void PrintDaughterLog(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, G4double sampleSafety, G4bool onlySafety, const G4ThreeVector &sampleDirection, G4double sampleStep) const
G4bool CheckAndReportBadNormal(const G4ThreeVector &unitNormal, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double step, const G4VSolid *solid, const char *msg) const
void ComputeSafetyLog(const G4VSolid *solid, const G4ThreeVector &point, G4double safety, G4bool isMotherVolume, G4int banner=-1) const
void PostComputeStepLog(const G4VSolid *motherSolid, const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, G4double motherStep, G4double motherSafety) const
void ReportOutsideMother(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4VPhysicalVolume *motherPV, G4double tDist=30.0 *CLHEP::cm) const
void AlongComputeStepLog(const G4VSolid *sampleSolid, const G4ThreeVector &samplePoint, const G4ThreeVector &sampleDirection, const G4ThreeVector &localDirection, G4double sampleSafety, G4double sampleStep) const
G4double GetMinExtent() const
EAxis GetAxis() const
G4int GetMaxEquivalentSliceNo() const
G4int GetVolume(G4int pVolumeNo) const
std::size_t GetNoContained() const
G4int GetMinEquivalentSliceNo() const
G4bool IsNode() const
G4SmartVoxelNode * GetNode() const
G4SmartVoxelHeader * GetHeader() const
const G4RotationMatrix * GetRotation() const
const G4ThreeVector GetTranslation() const
G4LogicalVolume * GetLogicalVolume() const
virtual G4int GetCopyNo() const =0
const G4String & GetName() const
G4String GetName() const
virtual EInside Inside(const G4ThreeVector &p) const =0
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const =0
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
virtual G4GeometryType GetEntityType() const =0
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4NavigationHistory &history, const G4double pMaxLength=DBL_MAX) override
virtual void SetVerboseLevel(G4int level) override
G4NavigationLogger * fLogger
std::vector< G4double > fVoxelSliceWidthStack
std::vector< EAxis > fVoxelAxisStack
G4VoxelSafety * fpVoxelSafety
G4SmartVoxelNode * fVoxelNode
virtual G4double ComputeStep(const G4ThreeVector &globalPoint, const G4ThreeVector &globalDirection, const G4double currentProposedStepLength, G4double &newSafety, G4NavigationHistory &history, G4bool &validExitNormal, G4ThreeVector &exitNormal, G4bool &exiting, G4bool &entering, G4VPhysicalVolume *(*pBlockedPhysical), G4int &blockedReplicaNo) override
G4bool LocateNextVoxel(const G4ThreeVector &localPoint, const G4ThreeVector &localDirection, const G4double currentStep)
G4SmartVoxelNode * VoxelLocate(G4SmartVoxelHeader *pHead, const G4ThreeVector &localPoint)
virtual void RelocateWithinVolume(G4VPhysicalVolume *motherPhysical, const G4ThreeVector &localPoint) override
std::vector< G4int > fVoxelNodeNoStack
std::vector< G4SmartVoxelHeader * > fVoxelHeaderStack
std::vector< G4int > fVoxelNoSlicesStack
G4double ComputeVoxelSafety(const G4ThreeVector &localPoint) const
void SetVerboseLevel(G4int level)
G4double ComputeSafety(const G4ThreeVector &localPoint, const G4VPhysicalVolume &currentPhysical, G4double maxLength=DBL_MAX)
EAxis
Definition geomdefs.hh:54
@ kXAxis
Definition geomdefs.hh:55
EInside
Definition geomdefs.hh:67
@ kInside
Definition geomdefs.hh:70
@ kOutside
Definition geomdefs.hh:68
#define DBL_MAX
Definition templates.hh:62