G4VoxelSafety is an utility class for the handling isotropic safety in volumes containing only G4PVPlacement daughter volumes for which voxels have been constructed. More...

#include <G4VoxelSafety.hh>

Public Member Functions
	G4VoxelSafety ()
	~G4VoxelSafety ()
G4double	ComputeSafety (const G4ThreeVector &localPoint, const G4VPhysicalVolume &currentPhysical, G4double maxLength=DBL_MAX)
G4int	GetVerboseLevel () const
void	SetVerboseLevel (G4int level)

Protected Member Functions
G4double	SafetyForVoxelHeader (const G4SmartVoxelHeader *pHead, const G4ThreeVector &localPoint, G4double maxLength, const G4VPhysicalVolume &currentPhysical, G4double distUpperDepth=0.0, G4double previousMinSafety=DBL_MAX)
G4double	SafetyForVoxelNode (const G4SmartVoxelNode *curVoxelNode, const G4ThreeVector &localPoint)

Detailed Description

G4VoxelSafety is an utility class for the handling isotropic safety in volumes containing only G4PVPlacement daughter volumes for which voxels have been constructed.

Definition at line 60 of file G4VoxelSafety.hh.

Constructor & Destructor Documentation

◆ G4VoxelSafety()

G4VoxelSafety::G4VoxelSafety ( )

Constructor and Destructor.

Definition at line 44 of file G4VoxelSafety.cc.

  : fVoxelAxisStack(kNavigatorVoxelStackMax,kXAxis),
    fVoxelNoSlicesStack(kNavigatorVoxelStackMax,0),
    fVoxelSliceWidthStack(kNavigatorVoxelStackMax,0.),
    fVoxelNodeNoStack(kNavigatorVoxelStackMax,0),
    fVoxelHeaderStack(kNavigatorVoxelStackMax,(G4SmartVoxelHeader*)nullptr)
{
  kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
}

◆ ~G4VoxelSafety()

G4VoxelSafety::~G4VoxelSafety ( )

default

Member Function Documentation

◆ ComputeSafety()

G4double G4VoxelSafety::ComputeSafety	(	const G4ThreeVector &	localPoint,
		const G4VPhysicalVolume &	currentPhysical,
		G4double	maxLength = DBL_MAX )

Calculates the isotropic distance to the nearest boundary from the specified point in the local coordinate system. The localpoint utilised must be within the current volume.

Parameters

[in]	localPoint	Local point.
[in]	currentPhysical	Current physical volume.
[in]	maxLength	Maximum length beyond which volumes are not checked.

Returns: Isotropic distance of given point to closest surface.

Definition at line 69 of file G4VoxelSafety.cc.

{
  G4LogicalVolume *motherLogical;
  G4VSolid *motherSolid;
  G4SmartVoxelHeader *motherVoxelHeader;
  G4double motherSafety, ourSafety;
  G4int localNoDaughters;
  G4double daughterSafety;
 
  motherLogical = currentPhysical.GetLogicalVolume();
  fpMotherLogical= motherLogical;   // For use by the other methods
  motherSolid = motherLogical->GetSolid();
  motherVoxelHeader = motherLogical->GetVoxelHeader();
 
#ifdef G4VERBOSE  
  if( fVerbose > 0 )
  { 
    G4cout << "*** G4VoxelSafety::ComputeSafety(): ***" << G4endl; 
  }
#endif
 
  // Check that point is inside mother volume
  //
  EInside  insideMother = motherSolid->Inside(localPoint); 
  if( insideMother != kInside  )
  { 
#ifdef G4DEBUG_NAVIGATION
    if( insideMother == kOutside )
    {
      std::ostringstream message;
      message << "Safety method called for location outside current Volume."
              << G4endl
              << "Location for safety is Outside this volume. " << G4endl
              << "The approximate distance to the solid "
              << "(safety from outside) is: " 
              << motherSolid->DistanceToIn( localPoint ) << G4endl;
      message << "  Problem occurred with physical volume: "
              << " Name: " << currentPhysical.GetName()
              << " Copy No: " << currentPhysical.GetCopyNo() << G4endl
              << "    Local Point = " << localPoint << G4endl;
      message << "  Description of solid: " << G4endl
              << *motherSolid << G4endl;
      G4Exception("G4VoxelSafety::ComputeSafety()", "GeomNav0003",
                  FatalException, message);
    }
#endif
    return 0.0;
  }   
 
  //  First limit:  mother safety - distance to outer boundaries
  //
  motherSafety = motherSolid->DistanceToOut(localPoint);
  ourSafety = motherSafety;                 // Working isotropic safety
   
#ifdef G4VERBOSE
  if(( fCheck ) ) // && ( fVerbose == 1 ))
  {
    G4cout << "    Invoked DistanceToOut(p) for mother solid: "
           << motherSolid->GetName()
           << ". Solid replied: " << motherSafety << G4endl
           << "    For local point p: " << localPoint
           << ", to be considered as 'mother safety'." << G4endl;
  }
#endif
  localNoDaughters = (G4int)motherLogical->GetNoDaughters();
 
  fBlockList.Enlarge(localNoDaughters);
  fBlockList.Reset();
 
  fVoxelDepth = -1;  // Resets the depth -- must be done for next method
  daughterSafety= SafetyForVoxelHeader(motherVoxelHeader, localPoint, maxLength,
                                       currentPhysical, 0.0, ourSafety);
  ourSafety= std::min( motherSafety, daughterSafety ); 
 
  return ourSafety;
}

◆ GetVerboseLevel()

G4int G4VoxelSafety::GetVerboseLevel ( ) const

inline

Verbosity control.

Note: If level>0 && G4VERBOSE, printout can occur.

Definition at line 87 of file G4VoxelSafety.hh.

87{ return fVerbose; }

◆ SafetyForVoxelHeader()

G4double G4VoxelSafety::SafetyForVoxelHeader	(	const G4SmartVoxelHeader *	pHead,
		const G4ThreeVector &	localPoint,
		G4double	maxLength,
		const G4VPhysicalVolume &	currentPhysical,
		G4double	distUpperDepth = 0.0,
		G4double	previousMinSafety = DBL_MAX )

protected

Cycles through levels of headers to process each node level.

Parameters

[in]	pHead	Voxel header.
[in]	localPoint	Local point.
[in]	maxLength	Maximum length beyond which volumes are not checked.
[in]	currentPhysical	Current volume (used for debug printout).
[in]	distUpperDepth	Upper square distance from voxel.
[in]	previousMinSafety	Minimum distance beyond which not to look.

Returns: Isotropic distance of the point to closest volume in all nodes.

Definition at line 211 of file G4VoxelSafety.cc.

{
  const G4SmartVoxelHeader* const targetVoxelHeader = pHeader;
  G4SmartVoxelNode* targetVoxelNode = nullptr;
 
  const G4SmartVoxelProxy* sampleProxy;
  EAxis    targetHeaderAxis;
  G4double targetHeaderMin, targetHeaderMax, targetHeaderNodeWidth;
  G4int    targetHeaderNoSlices;
  G4int    targetNodeNo;
 
  G4double minSafety = previousMinSafety;
  G4double ourSafety = DBL_MAX;
  unsigned int checkedNum= 0;
  
  ++fVoxelDepth;
  // fVoxelDepth  set by ComputeSafety or previous level call
 
  targetHeaderAxis =      targetVoxelHeader->GetAxis();
  targetHeaderNoSlices =  (G4int)targetVoxelHeader->GetNoSlices();
  targetHeaderMin =       targetVoxelHeader->GetMinExtent();
  targetHeaderMax =       targetVoxelHeader->GetMaxExtent();
 
  targetHeaderNodeWidth = (targetHeaderMax-targetHeaderMin)
                          / targetHeaderNoSlices;
 
  G4double localCrd = localPoint(targetHeaderAxis);
 
  const auto  candNodeNo = G4int( (localCrd-targetHeaderMin)
                                 / targetHeaderNodeWidth );
  // Ensure that it is between 0 and targetHeader->GetMaxExtent() - 1
 
#ifdef G4DEBUG_VOXELISATION  
  if( candNodeNo < 0 || candNodeNo > targetHeaderNoSlices-1 )
  {
     G4ExceptionDescription ed;
     ed << " Potential ERROR."
        << " Point is outside range of Voxel in current coordinate" << G4endl;
     ed << "  Node number of point " << localPoint
        << "is outside the range. " << G4endl;
     ed << "    Voxel node Num= " << candNodeNo << " versus  minimum= " << 0
        << " and maximum= " << targetHeaderNoSlices-1 << G4endl;
     ed << "    Axis = " << targetHeaderAxis
        << "  No of slices = " << targetHeaderNoSlices << G4endl;
     ed << "    Local coord = " << localCrd
        << "  Voxel Min = " << targetHeaderMin
        <<            " Max = " << targetHeaderMax << G4endl;
     G4LogicalVolume *pLogical= currentPhysical.GetLogicalVolume();
     ed << "  Current volume (physical) = " << currentPhysical.GetName()
        << "   (logical) = " << pLogical->GetName()     << G4endl;
     G4VSolid* pSolid= pLogical->GetSolid();
     ed << "  Solid type = " << pSolid->GetEntityType() << G4endl;
     ed << *pSolid << G4endl;
     G4Exception("G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav1003",
                 JustWarning, ed,
                 "Point is outside range of Voxel in current coordinate");
  }
#endif
  
  const G4int pointNodeNo =
              std::max( 0, std::min( candNodeNo, targetHeaderNoSlices-1 ) );
  
#ifdef G4VERBOSE  
  if( fVerbose > 2 )
  { 
    G4cout << G4endl;
    G4cout << "**** G4VoxelSafety::SafetyForVoxelHeader  " << G4endl;
    G4cout << "  Called at Depth     = " << fVoxelDepth;
    G4cout << " Calculated pointNodeNo= " << pointNodeNo
           << "  from position= " <<  localPoint(targetHeaderAxis)
           << "  min= "    << targetHeaderMin
           << "  max= "    << targetHeaderMax
           << "  width= "  << targetHeaderNodeWidth 
           << "  no-slices= " << targetHeaderNoSlices
           << "  axis=  "  << targetHeaderAxis   << G4endl;
  }
  else if (fVerbose == 1)
  {
    G4cout << " VoxelSafety: Depth  = " << fVoxelDepth 
           << " Number of Slices = " << targetHeaderNoSlices
           << " Header (address) = " << targetVoxelHeader  << G4endl;
  }
#endif
 
  // Stack info for stepping
  //
  fVoxelAxisStack[fVoxelDepth] = targetHeaderAxis;
  fVoxelNoSlicesStack[fVoxelDepth] = targetHeaderNoSlices;
  fVoxelSliceWidthStack[fVoxelDepth] = targetHeaderNodeWidth;
 
  fVoxelHeaderStack[fVoxelDepth] = pHeader;
 
  G4int   trialUp = -1,     trialDown = -1;
  G4double distUp = DBL_MAX, distDown = DBL_MAX;
 
  // Using Equivalent voxels - this is pre-initialisation only
  //
  G4int nextUp =   pointNodeNo+1;
  G4int nextDown = pointNodeNo-1;
 
  G4int    nextNodeNo = pointNodeNo;
  G4double distAxis;  // Distance in current Axis
  distAxis = 0.0;  // Starting in node containing local Coordinate
 
  G4bool nextIsInside = false;
 
  G4double distMaxInterest= std::min( previousMinSafety, maxLength);
    // We will not look beyond this distance.
    // This distance will be updated to reflect the
    // max ( minSafety, maxLength ) at each step
 
  targetNodeNo = pointNodeNo;
  do
  {
     G4double nodeSafety = DBL_MAX, headerSafety = DBL_MAX;
     fVoxelNodeNoStack[fVoxelDepth] = targetNodeNo;
    
     ++checkedNum; 
    
     sampleProxy = targetVoxelHeader->GetSlice(targetNodeNo);
 
#ifdef G4DEBUG_NAVIGATION
     assert( sampleProxy != 0);
     if( fVerbose > 2 ) 
     {
       G4cout << " -Checking node " << targetNodeNo
              << " is proxy with address " << sampleProxy << G4endl;
     }
#endif 
 
     if ( sampleProxy == nullptr )
     {
       G4ExceptionDescription ed;
       ed << " Problem for node number= " << targetNodeNo
          << "    Number of slides = " << targetHeaderNoSlices
          << G4endl;
       G4Exception( "G4VoxelSafety::SafetyForVoxelHeader()", "GeomNav0003",
                    FatalException, ed,
                    "Problem sampleProxy is Zero. Failure in loop.");
     }
     else if ( sampleProxy->IsNode() )
     {
       targetVoxelNode = sampleProxy->GetNode();
 
       // Deal with the node here [ i.e. the last level ]
       //
       nodeSafety= SafetyForVoxelNode( targetVoxelNode, localPoint);
#ifdef G4DEBUG_NAVIGATION
       if( fVerbose > 2 )
       {
          G4cout << " -- It is a Node ";
          G4cout << "    its safety= " << nodeSafety
                 << " our level Saf = " << ourSafety
                 << "  MinSaf= " << minSafety << G4endl;
       }
#endif
        ourSafety= std::min( ourSafety, nodeSafety );
       
        trialUp = targetVoxelNode->GetMaxEquivalentSliceNo()+1;
        trialDown = targetVoxelNode->GetMinEquivalentSliceNo()-1;
     }
     else  
     {
        const G4SmartVoxelHeader* pNewVoxelHeader = sampleProxy->GetHeader();
 
        G4double distCombined_Sq;
        distCombined_Sq = distUpperDepth_Sq + distAxis*distAxis;
 
#ifdef G4DEBUG_NAVIGATION
        if( fVerbose > 2 )
        {
          G4double distCombined= std::sqrt( distCombined_Sq );
          G4double distUpperDepth= std::sqrt ( distUpperDepth_Sq );
          G4cout << " -- It is a Header " << G4endl;
          G4cout << "  Recurse to deal with next level, fVoxelDepth= " 
                 << fVoxelDepth+1 << G4endl;
          G4cout << "  Distances:  Upper= " << distUpperDepth
                 << " Axis= " << distAxis
                 << " Combined= " << distCombined << G4endl;
        }
#endif
        
        // Recurse to deal with lower levels
        //
        headerSafety= SafetyForVoxelHeader( pNewVoxelHeader, localPoint,
                                            maxLength, currentPhysical,
                                            distCombined_Sq, minSafety);
        ourSafety = std::min( ourSafety, headerSafety );
       
#ifdef G4DEBUG_NAVIGATION
        if( fVerbose > 2 )
        { 
          G4cout << "      >> Header safety = " << headerSafety
                 << " our level Saf = " << ourSafety << G4endl;
        }
#endif
        trialUp =   pNewVoxelHeader->GetMaxEquivalentSliceNo()+1;
        trialDown = pNewVoxelHeader->GetMinEquivalentSliceNo()-1;
     }
     minSafety = std::min( minSafety, ourSafety );
 
     // Find next closest Voxel
     //    - first try: by simple subtraction
     //    - later:  using distance  (TODO - tbc)
     //
     if( targetNodeNo >= pointNodeNo )
     {
        nextUp = trialUp;
        // distUp   = std::max( targetHeaderMax-localCrd,  0.0 );
        G4double lowerEdgeOfNext = targetHeaderMin
                                 + nextUp * targetHeaderNodeWidth;
        distUp = lowerEdgeOfNext-localCrd ;
        if( distUp < 0.0 )
        {
           distUp = DBL_MAX;  // On the wrong side - must not be considered
        }
#ifdef G4DEBUG_NAVIGATION
       if( fVerbose > 2 )
       {
         G4cout << " > Updated nextUp= " << nextUp << G4endl;
       }
#endif
     }
     
     if( targetNodeNo <= pointNodeNo ) 
     {
         nextDown = trialDown;
         // distDown = std::max( localCrd-targetHeaderMin,  0.0);
         G4double upperEdgeOfNext = targetHeaderMin
                                  + (nextDown+1) * targetHeaderNodeWidth;
         distDown = localCrd-upperEdgeOfNext;
         if( distDown < 0.0 )
         {
            distDown= DBL_MAX; // On the wrong side - must not be considered 
         }
#ifdef G4DEBUG_NAVIGATION
       if( fVerbose > 2 )
       {
         G4cout << " > Updated nextDown= " << nextDown << G4endl;
       }
#endif
     }
 
#ifdef G4DEBUG_NAVIGATION
     if( fVerbose > 2 )
     {
       G4cout << " Node= " << pointNodeNo
              << " Up:   next= " << nextUp  << " d# "
              << nextUp - pointNodeNo
              << " trialUp=  " << trialUp << " d# "
              << trialUp - pointNodeNo
              << " Down: next= " << nextDown << " d# "
              << targetNodeNo - nextDown
              << " trialDwn= " << trialDown << " d# "
              << targetNodeNo - trialDown
              << " condition (next is Inside)= " << nextIsInside
              << G4endl;
     }
#endif     
     
     G4bool UpIsClosest;
     UpIsClosest = distUp < distDown;
     
     if( (nextUp < targetHeaderNoSlices)
         && (UpIsClosest || (nextDown < 0)) )
     {
       nextNodeNo = nextUp;
       distAxis = distUp;
       ++nextUp; // Default
#ifdef G4VERBOSE
       if( fVerbose > 2 )
       {
          G4cout << " > Chose Up.   Depth= " << fVoxelDepth
                 << "   Nodes: next= " << nextNodeNo
                 << " new nextUp=   " << nextUp
                 << " Dist = " << distAxis << G4endl;
       }
#endif
     }
     else
     {
       nextNodeNo = nextDown;
       distAxis = distDown;
       --nextDown; // A default value
#ifdef G4VERBOSE
       if( fVerbose > 2 )
       {
         G4cout << " > Chose Down.  Depth= " << fVoxelDepth
                << "  Nodes: next= " << nextNodeNo
                << " new nextDown=  " << nextDown
                << " Dist = " << distAxis << G4endl;
       }
#endif
     }
 
     nextIsInside = (nextNodeNo >= 0) && (nextNodeNo < targetHeaderNoSlices);
     if( nextIsInside )
     {
       targetNodeNo= nextNodeNo;
 
#ifdef G4DEBUG_NAVIGATION
       assert( targetVoxelHeader->GetSlice(nextNodeNo) != 0 );
       G4bool bContinue = (distAxis<minSafety);
       if( !bContinue )
       {
         if( fVerbose > 2 )
         {
           G4cout << " Can skip remaining at depth " << targetHeaderAxis
                  << " >>  distAxis= " << distAxis
                  << " minSafety= " << minSafety << G4endl;
         }
       }
#endif
     }
     else
     {
#ifdef G4DEBUG_NAVIGATION
       if( fVerbose > 2)
       {
         G4cout << " VoxSaf> depth= " << fVoxelDepth << G4endl;
         G4cout << " VoxSaf> No more candidates:  nodeDown= " << nextDown
                << " nodeUp= " << nextUp
                << " lastSlice= " << targetHeaderNoSlices << G4endl;
       }
#endif
     }
 
     // This calculation can be 'hauled'-up to where minSafety is calculated
     //
     distMaxInterest = std::min( minSafety, maxLength ); 
 
  } while ( nextIsInside && ( distAxis*distAxis + distUpperDepth_Sq
                            < distMaxInterest*distMaxInterest ) ); 
 
#ifdef G4VERBOSE
  if( fVerbose > 0 )
  { 
    G4cout << " Ended for targetNodeNo -- checked " << checkedNum << " out of "
           << targetHeaderNoSlices << " slices." << G4endl;
    G4cout << " ===== Returning from SafetyForVoxelHeader " 
           << "  Depth= " << fVoxelDepth << G4endl
           << G4endl;  
  }
#endif
 
  // Go back one level
  fVoxelDepth--; 
  
  return ourSafety;
}

Referenced by ComputeSafety(), and SafetyForVoxelHeader().

◆ SafetyForVoxelNode()

G4double G4VoxelSafety::SafetyForVoxelNode	(	const G4SmartVoxelNode *	curVoxelNode,
		const G4ThreeVector &	localPoint )

protected

Calculates the safety for volumes included in current Voxel Node.

Parameters

[in]	curVoxelNode	Voxel node.
[in]	localPoint	Local point.

Returns: Isotropic distance of given point to closest volume in node.

Definition at line 155 of file G4VoxelSafety.cc.

{
   G4double ourSafety = DBL_MAX;
 
   G4long curNoVolumes, contentNo;
   G4int sampleNo;
   G4VPhysicalVolume* samplePhysical;
 
   G4double      sampleSafety = 0.0; 
   G4ThreeVector samplePoint;
   G4VSolid*     ptrSolid = nullptr;
 
   curNoVolumes = curVoxelNode->GetNoContained();
 
   for ( contentNo=curNoVolumes-1; contentNo>=0; contentNo-- )
   {
      sampleNo = curVoxelNode->GetVolume((G4int)contentNo);
      if ( !fBlockList.IsBlocked(sampleNo) ) 
      { 
        fBlockList.BlockVolume(sampleNo);
 
        samplePhysical = fpMotherLogical->GetDaughter(sampleNo);
        G4AffineTransform sampleTf(samplePhysical->GetRotation(),
                                   samplePhysical->GetTranslation());
        sampleTf.Invert();
        samplePoint = sampleTf.TransformPoint(localPoint);
        ptrSolid = samplePhysical->GetLogicalVolume()->GetSolid();
 
        sampleSafety = ptrSolid->DistanceToIn(samplePoint);
        ourSafety = std::min( sampleSafety, ourSafety ); 
#ifdef G4VERBOSE
        if(( fCheck ) && ( fVerbose == 1 ))
        {
          G4cout << "*** G4VoxelSafety::SafetyForVoxelNode(): ***" << G4endl
                 << "    Invoked DistanceToIn(p) for daughter solid: "
                 << ptrSolid->GetName()
                 << ". Solid replied: " << sampleSafety << G4endl
                 << "    For local point p: " << samplePoint
                 << ", to be considered as 'daughter safety'." << G4endl;
        }
#endif
      }
   }  // end for contentNo
 
   return ourSafety; 
}

Referenced by SafetyForVoxelHeader().

◆ SetVerboseLevel()

void G4VoxelSafety::SetVerboseLevel ( G4int level )

inline

Definition at line 88 of file G4VoxelSafety.hh.

88{ fVerbose = level; }

The documentation for this class was generated from the following files:

Public Member Functions

Protected Member Functions

Detailed Description

Constructor & Destructor Documentation

◆ G4VoxelSafety()

◆ ~G4VoxelSafety()

Member Function Documentation

◆ ComputeSafety()

◆ GetVerboseLevel()

◆ SafetyForVoxelHeader()

◆ SafetyForVoxelNode()

◆ SetVerboseLevel()