G4SmartVoxelHeader Class Reference

#include <G4SmartVoxelHeader.hh>

Public Member Functions
	G4SmartVoxelHeader (G4LogicalVolume *pVolume, G4int pSlice=0)
	~G4SmartVoxelHeader ()
G4int	GetMaxEquivalentSliceNo () const
void	SetMaxEquivalentSliceNo (G4int pMax)
G4int	GetMinEquivalentSliceNo () const
void	SetMinEquivalentSliceNo (G4int pMin)
EAxis	GetAxis () const
EAxis	GetParamAxis () const
G4double	GetMaxExtent () const
G4double	GetMinExtent () const
G4int	GetNoSlices () const
G4SmartVoxelProxy *	GetSlice (G4int n) const
G4bool	AllSlicesEqual () const
G4bool	operator== (const G4SmartVoxelHeader &pHead) const
	G4SmartVoxelHeader (G4LogicalVolume *pVolume, const G4VoxelLimits &pLimits, const G4VolumeNosVector *pCandidates, G4int pSlice=0)

Protected Member Functions
void	BuildVoxels (G4LogicalVolume *pVolume)
void	BuildReplicaVoxels (G4LogicalVolume *pVolume)
void	BuildConsumedNodes (G4int nReplicas)
void	BuildVoxelsWithinLimits (G4LogicalVolume *pVolume, G4VoxelLimits pLimits, const G4VolumeNosVector *pCandidates)
void	BuildEquivalentSliceNos ()
void	CollectEquivalentNodes ()
void	CollectEquivalentHeaders ()
G4ProxyVector *	BuildNodes (G4LogicalVolume *pVolume, G4VoxelLimits pLimits, const G4VolumeNosVector *pCandidates, EAxis pAxis)
G4double	CalculateQuality (G4ProxyVector *pSlice)
void	RefineNodes (G4LogicalVolume *pVolume, G4VoxelLimits pLimits)

Protected Attributes
G4int	fminEquivalent
G4int	fmaxEquivalent
EAxis	faxis
EAxis	fparamAxis
G4double	fmaxExtent
G4double	fminExtent
G4ProxyVector	fslices

Friends
std::ostream &	operator<< (std::ostream &s, const G4SmartVoxelHeader &h)

Detailed Description

Definition at line 78 of file G4SmartVoxelHeader.hh.

Constructor & Destructor Documentation

G4SmartVoxelHeader::G4SmartVoxelHeader	(	G4LogicalVolume *	pVolume,
		G4int	pSlice = 0
	)

Definition at line 66 of file G4SmartVoxelHeader.cc.

  : fminEquivalent(pSlice),
    fmaxEquivalent(pSlice),
    fparamAxis(kUndefined)
{
  G4int nDaughters = pVolume->GetNoDaughters();
  G4VoxelLimits limits;   // Create `unlimited' limits object

  // Determine whether daughter is replicated
  //
  if ((nDaughters!=1) || (!pVolume->GetDaughter(0)->IsReplicated()))
  {
    // Daughter not replicated => conventional voxel Build
    // where each daughters extents are computed
    //
    BuildVoxels(pVolume);
  }
  else
  {
    // Single replicated daughter
    //
    BuildReplicaVoxels(pVolume);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4SmartVoxelHeader::~G4SmartVoxelHeader

(

)

Definition at line 127 of file G4SmartVoxelHeader.cc.

{
  // Manually destroy underlying nodes/headers
  // Delete collected headers and nodes once only
  //
  G4int node, proxy, maxNode=fslices.size();
  G4SmartVoxelProxy *lastProxy=0;
  G4SmartVoxelNode *dyingNode, *lastNode=0;
  G4SmartVoxelHeader *dyingHeader, *lastHeader=0;

  for (node=0; node<maxNode; node++)
  {
    if (fslices[node]->IsHeader())
    {
      dyingHeader = fslices[node]->GetHeader();
      if (lastHeader!=dyingHeader)
      {
        lastHeader = dyingHeader;
        lastNode = 0;
        delete dyingHeader;
      }
    }
    else
    {
      dyingNode = fslices[node]->GetNode();
      if (dyingNode!=lastNode)
      {
        lastNode=dyingNode;
        lastHeader=0;
        delete dyingNode;
      }
    }
  }
  // Delete proxies
  //
  for (proxy=0; proxy<maxNode; proxy++)
  {
    if (fslices[proxy]!=lastProxy)
    {
      lastProxy = fslices[proxy];
      delete lastProxy;
    }
  }
  // Don't need to clear slices
  // fslices.clear();
}

G4SmartVoxelHeader::G4SmartVoxelHeader	(	G4LogicalVolume *	pVolume,
		const G4VoxelLimits &	pLimits,
		const G4VolumeNosVector *	pCandidates,
		G4int	pSlice = 0
	)

Definition at line 100 of file G4SmartVoxelHeader.cc.

  : fminEquivalent(pSlice),
    fmaxEquivalent(pSlice),
    fparamAxis(kUndefined)
{
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "**** G4SmartVoxelHeader::G4SmartVoxelHeader" << G4endl
         << "     Limits " << pLimits << G4endl
         << "     Candidate #s = " ;
  for (size_t i=0;i<pCandidates->size();i++)
  {
    G4cout << (*pCandidates)[i] << " ";
  }
  G4cout << G4endl;
#endif   

  BuildVoxelsWithinLimits(pVolume,pLimits,pCandidates);
}

Here is the call graph for this function:

Member Function Documentation

G4bool G4SmartVoxelHeader::AllSlicesEqual

(

)

const

Definition at line 1281 of file G4SmartVoxelHeader.cc.

{
  G4int noSlices = fslices.size();
  G4SmartVoxelProxy* refProxy;

  if (noSlices>1)
  {
    refProxy=fslices[0];
    for (G4int i=1; i<noSlices; i++)
    {
      if (refProxy!=fslices[i])
      {
        return false;
      }
    }
  }
  return true;
}

void G4SmartVoxelHeader::BuildConsumedNodes ( G4int  nReplicas )

protected

Definition at line 395 of file G4SmartVoxelHeader.cc.

{
  G4int nNode, nVol;
  G4SmartVoxelNode *pNode;
  G4SmartVoxelProxy *pProxyNode;

  // Create and fill nodes in temporary G4NodeVector (on stack)
  //
  G4NodeVector nodeList;
  nodeList.reserve(nReplicas);
  for (nNode=0; nNode<nReplicas; nNode++)
  {
    pNode=new G4SmartVoxelNode(nNode);
    if (!pNode)
    {
      G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "GeomMgt0003",
                  FatalException, "Node allocation error.");
    }
    nodeList.push_back(pNode);
  }
  for (nVol=0; nVol<nReplicas; nVol++)
  {
    nodeList[nVol]->Insert(nVol);   // Insert replication of number
  }                                 // identical to voxel number

  // Create & fill proxy List `in place' by modifying instance data fslices
  //
  fslices.clear();
  for (nNode=0; nNode<nReplicas; nNode++)
  {
    pProxyNode = new G4SmartVoxelProxy(nodeList[nNode]);
    if (!pProxyNode)
    {
      G4Exception("G4SmartVoxelHeader::BuildConsumedNodes()", "GeomMgt0003",
                  FatalException, "Proxy node allocation error.");
    }
    fslices.push_back(pProxyNode);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4SmartVoxelHeader::BuildEquivalentSliceNos ( )

protected

Definition at line 583 of file G4SmartVoxelHeader.cc.

{
  G4int sliceNo, minNo, maxNo, equivNo;
  G4int maxNode = fslices.size();
  G4SmartVoxelNode *startNode, *sampleNode;
  for (sliceNo=0; sliceNo<maxNode; sliceNo++)
  {
    minNo = sliceNo;

    // Get first node (see preconditions - will throw exception if a header)
    //
    startNode = fslices[minNo]->GetNode();

    // Find max equivalent
    //
    for (equivNo=minNo+1; equivNo<maxNode; equivNo++)
    {
      sampleNode = fslices[equivNo]->GetNode();
      if (!((*startNode) == (*sampleNode))) { break; }
    }
    maxNo = equivNo-1;
    if (maxNo != minNo)
    {
      // Set min and max nos
      //
      for (equivNo=minNo; equivNo<=maxNo; equivNo++)
      {
        sampleNode = fslices[equivNo]->GetNode();
        sampleNode->SetMinEquivalentSliceNo(minNo);
        sampleNode->SetMaxEquivalentSliceNo(maxNo);
      }
      // Advance outer loop to end of equivalent group
      //
      sliceNo = maxNo;
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4ProxyVector * G4SmartVoxelHeader::BuildNodes ( G4LogicalVolume *  pVolume, G4VoxelLimits  pLimits, const G4VolumeNosVector *  pCandidates, EAxis  pAxis  )

protected

Definition at line 750 of file G4SmartVoxelHeader.cc.

{
  G4double motherMinExtent= kInfinity, motherMaxExtent= -kInfinity,
           targetMinExtent= kInfinity, targetMaxExtent= -kInfinity;
  G4VPhysicalVolume *pDaughter=0;
  G4VPVParameterisation *pParam=0;
  G4VSolid *targetSolid;
  G4AffineTransform targetTransform;
  G4bool replicated;
  G4int nCandidates = pCandidates->size();
  G4int nVol, nNode, targetVolNo;
  G4VoxelLimits noLimits;
    
#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
         << "     Limits = " << pLimits << G4endl
         << "       Axis = " << pAxis << G4endl
         << " Candidates = " << nCandidates << G4endl;
#endif

  // Compute extent of logical volume's solid along this axis
  // NOTE: results stored locally and not preserved/reused
  //
  G4VSolid* outerSolid = pVolume->GetSolid();
  const G4AffineTransform origin;
  if( !outerSolid->CalculateExtent(pAxis, pLimits, origin,
                                   motherMinExtent, motherMaxExtent) )
  {
    outerSolid->CalculateExtent(pAxis, noLimits, origin,
                                motherMinExtent, motherMaxExtent);
  }
  G4VolumeExtentVector minExtents(nCandidates,0.);
  G4VolumeExtentVector maxExtents(nCandidates,0.);

  if ( (pVolume->GetNoDaughters()==1)
    && (pVolume->GetDaughter(0)->IsReplicated()==true) )
  {
    // Replication data not required: only parameterisation object 
    // and volume no. List used
    //
    pDaughter = pVolume->GetDaughter(0);
    pParam = pDaughter->GetParameterisation();
    if (!pParam)
    {
      std::ostringstream message;
      message << "PANIC! - Missing parameterisation." << G4endl
              << "         Replicated volume with no parameterisation object !";
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                  FatalException, message);
      return 0;
    }

    // Setup daughter's transformations
    //
    targetTransform = G4AffineTransform(pDaughter->GetRotation(),
                                        pDaughter->GetTranslation());
    replicated = true;
  }
    else
  {
    replicated = false;
  }
    
  // Compute extents
  //
  for (nVol=0; nVol<nCandidates; nVol++)
  {
    targetVolNo=(*pCandidates)[nVol];
    if (replicated == false)
    {
      pDaughter=pVolume->GetDaughter(targetVolNo);

      // Setup daughter's transformations
      //
      targetTransform = G4AffineTransform(pDaughter->GetRotation(),
                                          pDaughter->GetTranslation());
      // Get underlying (and setup) solid
      //
      targetSolid = pDaughter->GetLogicalVolume()->GetSolid();
    }
    else
    {
      // Find  solid
      //
      targetSolid = pParam->ComputeSolid(targetVolNo,pDaughter);

      // Setup solid
      //
      targetSolid->ComputeDimensions(pParam,targetVolNo,pDaughter);

      // Setup transform
      //
      pParam->ComputeTransformation(targetVolNo,pDaughter);
      targetTransform = G4AffineTransform(pDaughter->GetRotation(),
                                          pDaughter->GetTranslation());
    }
    // Calculate extents
    //
    if(!targetSolid->CalculateExtent(pAxis, pLimits, targetTransform,
                                     targetMinExtent, targetMaxExtent))
    {
      targetSolid->CalculateExtent(pAxis, noLimits, targetTransform,
                                   targetMinExtent,targetMaxExtent);
    }
    minExtents[nVol] = targetMinExtent;
    maxExtents[nVol] = targetMaxExtent;

#ifdef G4GEOMETRY_VOXELDEBUG
   G4cout << "---------------------------------------------------" << G4endl
          << "     Volume = " << pDaughter->GetName() << G4endl
          << " Min Extent = " << targetMinExtent << G4endl
          << " Max Extent = " << targetMaxExtent << G4endl
          << "---------------------------------------------------" << G4endl;
#endif

    // Check not entirely outside mother when processing toplevel nodes
    //
    if ( (!pLimits.IsLimited()) && ((targetMaxExtent<=motherMinExtent)
                                  ||(targetMinExtent>=motherMaxExtent)) )
    {
      std::ostringstream message;
      message << "PANIC! - Overlapping daughter with mother volume." << G4endl
              << "         Daughter physical volume "
              << pDaughter->GetName() << G4endl
              << "         is entirely outside mother logical volume "
              << pVolume->GetName() << " !!";
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0002",
                  FatalException, message);
    }

#ifdef G4GEOMETRY_VOXELDEBUG
    // Check for straddling volumes when debugging.
    // If a volume is >kStraddlePercent percent over the mother
    // boundary, print a warning.
    //
    if (!pLimits.IsLimited())
    {
      G4double width;
      G4int kStraddlePercent=5;
      width = maxExtents[nVol]-minExtents[nVol];
      if ( (((motherMinExtent-minExtents[nVol])*100/width) > kStraddlePercent)
         ||(((maxExtents[nVol]-motherMaxExtent)*100/width) > kStraddlePercent) )
      {
        G4cout << "**** G4SmartVoxelHeader::BuildNodes" << G4endl
               << "     WARNING : Daughter # " << nVol
               << " name = " << pDaughter->GetName() << G4endl
               << "     Crosses mother boundary of logical volume, name = " 
               << pVolume->GetName() << G4endl
               << "     by more than " << kStraddlePercent 
               << "%" << G4endl;
      }
    }
#endif
  }

  // Extents of all daughters known

  // Calculate minimum slice width, only including volumes inside the limits
  //
  G4double minWidth = kInfinity;
  G4double currentWidth;
  for (nVol=0; nVol<nCandidates; nVol++)
  {
    // currentWidth should -always- be a positive value. Inaccurate computed extent
    // from the solid or situations of malformed geometries (overlaps) may lead to
    // negative values and therefore unpredictable crashes !
    //
    currentWidth = std::abs(maxExtents[nVol]-minExtents[nVol]);
    if ( (currentWidth<minWidth)
      && (maxExtents[nVol]>=pLimits.GetMinExtent(pAxis))
      && (minExtents[nVol]<=pLimits.GetMaxExtent(pAxis)) )
    {
      minWidth = currentWidth;
    }
  }

  // No. of Nodes formula - nearest integer to
  // mother width/half min daughter width +1
  //
  G4double noNodesExactD = ((motherMaxExtent-motherMinExtent)*2.0/minWidth)+1.0;

  // Compare with "smartless quality", i.e. the average number of slices
  // used per contained volume.
  //
  G4double smartlessComputed = noNodesExactD / nCandidates;
  G4double smartlessUser = pVolume->GetSmartless();
  G4double smartless = (smartlessComputed <= smartlessUser)
                       ? smartlessComputed : smartlessUser;
  G4double noNodesSmart = smartless*nCandidates;
  G4int    noNodesExactI = G4int(noNodesSmart);
  G4int    noNodes = ((noNodesSmart-noNodesExactI)>=0.5)
                     ? noNodesExactI+1 : noNodesExactI;
  if( noNodes == 0 ) { noNodes=1; }

#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "     Smartless computed = " << smartlessComputed << G4endl
         << "     Smartless volume = " << smartlessUser
         << " => # Smartless = " << smartless << G4endl;
  G4cout << "     Min width = " << minWidth
         << " => # Nodes = " << noNodes << G4endl;
#endif

  if (noNodes>kMaxVoxelNodes)
  {
    noNodes=kMaxVoxelNodes;
#ifdef G4GEOMETRY_VOXELDEBUG
    G4cout << "     Nodes Clipped to = " << kMaxVoxelNodes << G4endl;
#endif   
  }
  G4double nodeWidth = (motherMaxExtent-motherMinExtent)/noNodes;

  // Create G4VoxelNodes. Will Add proxies before setting fslices
  //
  G4NodeVector* nodeList = new G4NodeVector();
  if (!nodeList)
  {
    G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                FatalException, "NodeList allocation error.");
    return 0;
  }
  nodeList->reserve(noNodes);

  for (nNode=0; nNode<noNodes; nNode++)
  {
    G4SmartVoxelNode *pNode;
    pNode = new G4SmartVoxelNode(nNode);
    if (!pNode)
    {
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                  FatalException, "Node allocation error.");
      return 0;
    }
    nodeList->push_back(pNode);
  }

  // All nodes created (empty)

  // Fill nodes: Step through extent lists
  //
  for (nVol=0; nVol<nCandidates; nVol++)
  {
    G4int nodeNo, minContainingNode, maxContainingNode;
    minContainingNode = G4int((minExtents[nVol]-motherMinExtent)/nodeWidth);
    maxContainingNode = G4int((maxExtents[nVol]-motherMinExtent)/nodeWidth);

    // Only add nodes that are inside the limits of the axis
    //
    if ( (maxContainingNode>=0) && (minContainingNode<noNodes) )
    {
      // If max extent is on max boundary => maxContainingNode=noNodes:
      // should be one less as nodeList has noNodes entries
      //
      if (maxContainingNode>=noNodes)
      {
        maxContainingNode = noNodes-1;
      }
      //
      // Protection against protruding volumes
      //
      if (minContainingNode<0)
      {
        minContainingNode=0;
      }
      for (nodeNo=minContainingNode; nodeNo<=maxContainingNode; nodeNo++)
      {
        (*nodeList)[nodeNo]->Insert((*pCandidates)[nVol]);
      }
    }
  }

  // All nodes filled

  // Create proxy List : caller has deletion responsibility
  // (but we must delete nodeList *itself* - not the contents)
  //
  G4ProxyVector* proxyList = new G4ProxyVector();
  if (!proxyList)
  {
    G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                FatalException, "Proxy list allocation error.");
    return 0;
  }
  proxyList->reserve(noNodes);

  //
  // Fill proxy List
  //
  for (nNode=0; nNode<noNodes; nNode++)
  {
    // Get rid of possible excess capacity in the internal node vector
    //
    ((*nodeList)[nNode])->Shrink();
    G4SmartVoxelProxy* pProxyNode = new G4SmartVoxelProxy((*nodeList)[nNode]);
    if (!pProxyNode)
    {
      G4Exception("G4SmartVoxelHeader::BuildNodes()", "GeomMgt0003",
                  FatalException, "Proxy node allocation failed.");
      return 0;
    }
    proxyList->push_back(pProxyNode);
  }
  delete nodeList;
  return proxyList;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4SmartVoxelHeader::BuildReplicaVoxels ( G4LogicalVolume *  pVolume )

protected

Definition at line 267 of file G4SmartVoxelHeader.cc.

{
  G4VPhysicalVolume *pDaughter=0;

  // Replication data
  //
  EAxis axis;
  G4int nReplicas;
  G4double width,offset;
  G4bool consuming;

  // Consistency check: pVolume should contain single replicated volume
  //
  if ( (pVolume->GetNoDaughters()==1)
    && (pVolume->GetDaughter(0)->IsReplicated()==true) )
  {
    // Obtain replication data
    //
    pDaughter=pVolume->GetDaughter(0);
    pDaughter->GetReplicationData(axis,nReplicas,width,offset,consuming);
    fparamAxis = axis;
    if ( consuming==false )
    {
      G4VoxelLimits limits;   // Create `unlimited' limits object
      G4VolumeNosVector targetList;
      targetList.reserve(nReplicas);
      for (G4int i=0; i<nReplicas; i++)
      {
        targetList.push_back(i);
      }
      if (axis != kUndefined)
      {
        // Apply voxelisation along the specified axis only

        G4ProxyVector* pSlices=BuildNodes(pVolume,limits,&targetList,axis);
        faxis = axis;
        fslices = *pSlices;
        delete pSlices;

        // Calculate and set min and max extents given our axis
        //
        const G4AffineTransform origin;
        pVolume->GetSolid()->CalculateExtent(faxis, limits, origin,
                                             fminExtent, fmaxExtent);
        // Calculate equivalent nos
        //
        BuildEquivalentSliceNos();
        CollectEquivalentNodes();   // Collect common nodes
      }
      else
      {
        // Build voxels similarly as for normal placements considering
        // all three cartesian axes.

        BuildVoxelsWithinLimits(pVolume, limits, &targetList);
      }
    }
    else
    {
      // Replication is consuming -> Build voxels directly
      //
      // o Cartesian axes - range is -width*nREplicas/2 to +width*nREplicas/2
      //                    nReplicas replications result
      // o Radial axis (rho) = range is 0 to width*nReplicas
      //                    nReplicas replications result
      // o Phi axi       - range is offset to offset+width*nReplicas radians
      //
      // Equivalent slices no computation & collection not required - all
      // slices are different
      //
      switch (axis)
      {
        case kXAxis:
        case kYAxis:
        case kZAxis:
          fminExtent = -width*nReplicas*0.5;
          fmaxExtent =  width*nReplicas*0.5;
          break;
        case kRho:
          fminExtent = offset;
          fmaxExtent = width*nReplicas+offset;
          break;
        case kPhi:
          fminExtent = offset;
          fmaxExtent = offset+width*nReplicas;
          break;
        default:
          G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels()",
                      "GeomMgt0002", FatalException, "Illegal axis.");
          break;
      }  
      faxis = axis;   // Set axis
      BuildConsumedNodes(nReplicas);
      if ( (axis==kXAxis) || (axis==kYAxis) || (axis==kZAxis) )
      {
        // Sanity check on extent
        //
        G4double emin = kInfinity, emax = -kInfinity;
        G4VoxelLimits limits;
        G4AffineTransform origin;
        pVolume->GetSolid()->CalculateExtent(axis, limits, origin, emin, emax);
        if ( (std::fabs((emin-fminExtent)/fminExtent) +
              std::fabs((emax-fmaxExtent)/fmaxExtent)) > 0.05)
        {
          std::ostringstream message;
          message << "Sanity check: wrong solid extent." << G4endl
                  << "        Replicated geometry, logical volume: "
                  << pVolume->GetName();
          G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels",
                      "GeomMgt0002", FatalException, message);
        }
      }
    }
  }
  else
  {
    G4Exception("G4SmartVoxelHeader::BuildReplicaVoxels", "GeomMgt0002",
                FatalException, "Only one replicated daughter is allowed !");
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4SmartVoxelHeader::BuildVoxels ( G4LogicalVolume *  pVolume )

protected

Definition at line 247 of file G4SmartVoxelHeader.cc.

{
  G4VoxelLimits limits;   // Create `unlimited' limits object
  G4int nDaughters = pVolume->GetNoDaughters();

  G4VolumeNosVector targetList;
  targetList.reserve(nDaughters);
  for (G4int i=0; i<nDaughters; i++)
  {
    targetList.push_back(i);
  }
  BuildVoxelsWithinLimits(pVolume, limits, &targetList);
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4SmartVoxelHeader::BuildVoxelsWithinLimits ( G4LogicalVolume *  pVolume, G4VoxelLimits  pLimits, const G4VolumeNosVector *  pCandidates  )

protected

Definition at line 443 of file G4SmartVoxelHeader.cc.

{
  // Choose best axis for slicing by:
  // 1. Trying all unlimited cartesian axes
  // 2. Select axis which gives greatest no slices

  G4ProxyVector *pGoodSlices=0, *pTestSlices, *tmpSlices;
  G4double goodSliceScore=kInfinity, testSliceScore;
  EAxis goodSliceAxis = kXAxis;
  EAxis testAxis      = kXAxis;
  G4int node, maxNode, iaxis;
  G4VoxelLimits noLimits;

  // Try all non-limited cartesian axes
  //
  for (iaxis=0; iaxis<3; iaxis++)
  {
    switch(iaxis)
    {
      case 0:
        testAxis = kXAxis;
        break;
      case 1:
        testAxis = kYAxis;
        break;
      case 2:
        testAxis = kZAxis;
        break;
    }
    if (!pLimits.IsLimited(testAxis))
    {
      pTestSlices = BuildNodes(pVolume,pLimits,pCandidates,testAxis);
      testSliceScore = CalculateQuality(pTestSlices);
      if ( (!pGoodSlices) || (testSliceScore<goodSliceScore) )
      {
        goodSliceAxis  = testAxis;
        goodSliceScore = testSliceScore;
        tmpSlices      = pGoodSlices;
        pGoodSlices    = pTestSlices;
        pTestSlices    = tmpSlices;
      }
      if (pTestSlices)
      {
        // Destroy pTestSlices and all its contents
        //
        maxNode=pTestSlices->size();
        for (node=0; node<maxNode; node++)
        {
          delete (*pTestSlices)[node]->GetNode();
        }
        G4SmartVoxelProxy* tmpProx;
        while (pTestSlices->size()>0)  // Loop checking, 06.08.2015, G.Cosmo
        {
          tmpProx = pTestSlices->back();
          pTestSlices->pop_back();
          for (G4ProxyVector::iterator i=pTestSlices->begin();
                                       i!=pTestSlices->end(); )
          {
            if (*i==tmpProx)
            {
              i = pTestSlices->erase(i);
            }
            else
            {
              ++i;
            }
          }
          if ( tmpProx ) { delete tmpProx; }
        }
        delete pTestSlices;
      }
    }
  }
  // Check for error case.. when limits already 3d,
  // so cannot select a new axis
  //
  if (!pGoodSlices)
  {
    G4Exception("G4SmartVoxelHeader::BuildVoxelsWithinLimits()",
                "GeomMgt0002", FatalException,
                "Cannot select more than 3 axis for optimisation.");
    return;
  }

  // 
  // We have selected pGoodSlices, with a score testSliceScore
  //

  // Store chosen axis, slice ptr
  //
  fslices=*pGoodSlices; // Set slice information, copy ptrs in collection
  delete pGoodSlices;   // Destroy slices vector, but not contained
                        // proxies or nodes
  faxis=goodSliceAxis;

#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << G4endl << "     Volume = " << pVolume->GetName()
         << G4endl << "     Selected axis = " << faxis << G4endl;
  for (size_t islice=0; islice<fslices.size(); islice++)
  {
    G4cout << "     Node #" << islice << " = {";
    for (G4int j=0; j<fslices[islice]->GetNode()->GetNoContained(); j++)
    {
      G4cout << " " << fslices[islice]->GetNode()->GetVolume(j);
    }
    G4cout << " }" << G4endl;
  }
  G4cout << G4endl;
#endif

  // Calculate and set min and max extents given our axis
  //
  G4VSolid* outerSolid = pVolume->GetSolid();
  const G4AffineTransform origin;
  if(!outerSolid->CalculateExtent(faxis,pLimits,origin,fminExtent,fmaxExtent))
  {
    outerSolid->CalculateExtent(faxis,noLimits,origin,fminExtent,fmaxExtent);
  }

  // Calculate equivalent nos
  //
  BuildEquivalentSliceNos();
  CollectEquivalentNodes();     // Collect common nodes
  RefineNodes(pVolume,pLimits); // Refine nodes creating headers

  // No common headers can exist because collapsed by construction
}

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4SmartVoxelHeader::CalculateQuality ( G4ProxyVector *  pSlice )

protected

Definition at line 1073 of file G4SmartVoxelHeader.cc.

{
  G4double quality;
  G4int nNodes = pSlice->size();
  G4int noContained, maxContained=0, sumContained=0, sumNonEmptyNodes=0;
  G4SmartVoxelNode *node;

  for (G4int i=0; i<nNodes; i++)
  {
    if ((*pSlice)[i]->IsNode())
    {
      // Definitely a node. Add info to running totals
      //
      node = (*pSlice)[i]->GetNode();
      noContained = node->GetNoContained();
      if (noContained)
      {
        sumNonEmptyNodes++;
        sumContained += noContained;
        //
        // Calc maxContained for statistics
        //
        if (noContained>maxContained)
        {
          maxContained = noContained;
        }
      }
    }
    else
    {
      G4Exception("G4SmartVoxelHeader::CalculateQuality()", "GeomMgt0001",
                  FatalException, "Not applicable to replicated volumes.");
    }
  }

  // Calculate quality with protection against no non-empty nodes
  //
  if (sumNonEmptyNodes)
  {
    quality = sumContained/sumNonEmptyNodes;
  }
  else
  {
    quality = kInfinity;
  }

#ifdef G4GEOMETRY_VOXELDEBUG
  G4cout << "**** G4SmartVoxelHeader::CalculateQuality" << G4endl
         << "     Quality = " << quality << G4endl
         << "     Nodes = " << nNodes 
         << " of which " << sumNonEmptyNodes << " non empty" << G4endl
         << "     Max Contained = " << maxContained << G4endl;
#endif

  return quality;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4SmartVoxelHeader::CollectEquivalentHeaders ( )

protected

Definition at line 677 of file G4SmartVoxelHeader.cc.

{
  G4int sliceNo, maxNo, equivNo;
  G4int maxNode = fslices.size();
  G4SmartVoxelHeader *equivHeader, *sampleHeader;
  G4SmartVoxelProxy *equivProxy;

  for (sliceNo=0; sliceNo<maxNode; sliceNo++)
  {
    equivProxy = fslices[sliceNo];
    if (equivProxy->IsHeader())
    {
      equivHeader = equivProxy->GetHeader();
      maxNo = equivHeader->GetMaxEquivalentSliceNo();
      if (maxNo != sliceNo)
      {
        // Attempt collection between sliceNo and maxNo inclusive:
        // look for common headers. All slices between sliceNo and maxNo
        // are guaranteed to be headers but may not have equal contents
        //
#ifdef G4GEOMETRY_VOXELDEBUG
        G4cout << "**** G4SmartVoxelHeader::CollectEquivalentHeaders" << G4endl
               << "     Collecting Headers =";
#endif
        for (equivNo=sliceNo+1; equivNo<=maxNo; equivNo++)
        {
          sampleHeader = fslices[equivNo]->GetHeader();
          if ( (*sampleHeader) == (*equivHeader) )
          {
#ifdef G4GEOMETRY_VOXELDEBUG
            G4cout << " " << equivNo;
#endif
            // Delete sampleHeader + proxy and replace with equivHeader/Proxy
            //
            delete sampleHeader;
            delete fslices[equivNo];
            fslices[equivNo] = equivProxy;
          }
          else
          {
            // Not equal. Set this header to be
            // the current header for comparisons
            //
            equivProxy  = fslices[equivNo];
            equivHeader = equivProxy->GetHeader();
          }

        }
#ifdef G4GEOMETRY_VOXELDEBUG
        G4cout << G4endl;
#endif
        // Skip past examined slices
        //
        sliceNo = maxNo;
      }
    }
  }
}

Here is the call graph for this function:

void G4SmartVoxelHeader::CollectEquivalentNodes ( )

protected

Definition at line 630 of file G4SmartVoxelHeader.cc.

{
  G4int sliceNo, maxNo, equivNo;
  G4int maxNode=fslices.size();
  G4SmartVoxelNode *equivNode;
  G4SmartVoxelProxy *equivProxy;

  for (sliceNo=0; sliceNo<maxNode; sliceNo++)
  {
    equivProxy=fslices[sliceNo];

    // Assumption (see preconditions): all slices are nodes
    //
    equivNode = equivProxy->GetNode();
    maxNo = equivNode->GetMaxEquivalentSliceNo();
    if (maxNo != sliceNo)
    {
#ifdef G4GEOMETRY_VOXELDEBUG
      G4cout << "**** G4SmartVoxelHeader::CollectEquivalentNodes" << G4endl
             << "     Collecting Nodes = " 
             << sliceNo << " - " << maxNo << G4endl;
#endif
      // Do collection between sliceNo and maxNo inclusive
      //
      for (equivNo=sliceNo+1; equivNo<=maxNo; equivNo++)
      {
        delete fslices[equivNo]->GetNode();
        delete fslices[equivNo];
        fslices[equivNo] = equivProxy;
      }
      sliceNo = maxNo;
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

EAxis G4SmartVoxelHeader::GetAxis

(

)

const

Here is the caller graph for this function:

G4int G4SmartVoxelHeader::GetMaxEquivalentSliceNo

(

)

const

Here is the caller graph for this function:

G4double G4SmartVoxelHeader::GetMaxExtent

(

)

const

Here is the caller graph for this function:

G4int G4SmartVoxelHeader::GetMinEquivalentSliceNo

(

)

const

Here is the caller graph for this function:

G4double G4SmartVoxelHeader::GetMinExtent

(

)

const

Here is the caller graph for this function:

G4int G4SmartVoxelHeader::GetNoSlices

(

)

const

Here is the caller graph for this function:

EAxis G4SmartVoxelHeader::GetParamAxis

(

)

const

G4SmartVoxelProxy* G4SmartVoxelHeader::GetSlice

(

G4int

n

)

const

Here is the caller graph for this function:

G4bool G4SmartVoxelHeader::operator==

(

const G4SmartVoxelHeader &

pHead

)

const

Definition at line 183 of file G4SmartVoxelHeader.cc.

{
  if ( (GetAxis()      == pHead.GetAxis())
    && (GetNoSlices()  == pHead.GetNoSlices())
    && (GetMinExtent() == pHead.GetMinExtent())
    && (GetMaxExtent() == pHead.GetMaxExtent()) )
  {
    G4int node, maxNode;
    G4SmartVoxelProxy *leftProxy, *rightProxy;
    G4SmartVoxelHeader *leftHeader, *rightHeader;
    G4SmartVoxelNode *leftNode, *rightNode;

    maxNode=GetNoSlices();
    for (node=0; node<maxNode; node++)
    {
      leftProxy  = GetSlice(node);
      rightProxy = pHead.GetSlice(node);
      if (leftProxy->IsHeader())
      {
        if (rightProxy->IsNode())
        {
          return false;
        }
        else
        {
          leftHeader  = leftProxy->GetHeader();
          rightHeader = rightProxy->GetHeader();
          if (!(*leftHeader==*rightHeader))
          {
            return false;
          }
        }
      }
      else
      {
        if (rightProxy->IsHeader())
        {
          return false;
        }
        else
        {
          leftNode  = leftProxy->GetNode();
          rightNode = rightProxy->GetNode();
          if (!(*leftNode==*rightNode))
          {
            return false;
          }
        }
      }
    }
    return true;
  }
  else
  {
    return false;
  }
}

Here is the call graph for this function:

void G4SmartVoxelHeader::RefineNodes ( G4LogicalVolume *  pVolume, G4VoxelLimits  pLimits  )

protected

Definition at line 1139 of file G4SmartVoxelHeader.cc.

{
  G4int refinedDepth=0, minVolumes;
  G4int maxNode = fslices.size();

  if (pLimits.IsXLimited()) 
  {
    refinedDepth++;
  }
  if (pLimits.IsYLimited()) 
  {
    refinedDepth++;
  }
  if (pLimits.IsZLimited()) 
  {
    refinedDepth++;
  }

  // Calculate minimum number of volumes necessary to refine
  //
  switch (refinedDepth)
  {
    case 0:
      minVolumes=kMinVoxelVolumesLevel2;
      break;
    case 1:
      minVolumes=kMinVoxelVolumesLevel3;
      break;
    default:
      minVolumes=10000;   // catch refinedDepth=3 and errors
      break;
  }

  if (refinedDepth<2)
  {
    G4int targetNo, noContainedDaughters, minNo, maxNo, replaceNo, i;
    G4double sliceWidth = (fmaxExtent-fminExtent)/maxNode;
    G4VoxelLimits newLimits;
    G4SmartVoxelNode* targetNode;
    G4SmartVoxelProxy* targetNodeProxy;
    G4SmartVoxelHeader* replaceHeader;
    G4SmartVoxelProxy* replaceHeaderProxy;
    G4VolumeNosVector* targetList;
    G4SmartVoxelProxy* lastProxy;
      
    for (targetNo=0; targetNo<maxNode; targetNo++)
    {
      // Assume all slices are nodes (see preconditions)
      //
      targetNodeProxy = fslices[targetNo];
      targetNode = targetNodeProxy->GetNode();

      if (targetNode->GetNoContained() >= minVolumes)
      {
        noContainedDaughters = targetNode->GetNoContained();
        targetList = new G4VolumeNosVector();
        if (!targetList)
        {
          G4Exception("G4SmartVoxelHeader::RefineNodes()",
                      "GeomMgt0003", FatalException,
                      "Target volume node list allocation error.");
          return;
        }
        targetList->reserve(noContainedDaughters);
        for (i=0; i<noContainedDaughters; i++)
        {
          targetList->push_back(targetNode->GetVolume(i));
        }
        minNo = targetNode->GetMinEquivalentSliceNo();
        maxNo = targetNode->GetMaxEquivalentSliceNo();

#ifdef G4GEOMETRY_VOXELDEBUG
        G4cout << "**** G4SmartVoxelHeader::RefineNodes" << G4endl
               << "     Refining nodes " << minNo 
               << " - " << maxNo << " inclusive" << G4endl;
#endif
        if (minNo > maxNo)    // Delete node and list to be replaced
        {                     // and avoid further action ...
          delete targetNode;
          delete targetList;
          return;
        }

        // Delete node proxies at start of collected sets of nodes/headers
        //
        lastProxy=0;
        for (replaceNo=minNo; replaceNo<=maxNo; replaceNo++)
        {
          if (lastProxy != fslices[replaceNo])
          {
            lastProxy=fslices[replaceNo];
            delete lastProxy;
          }
        }
        // Delete node to be replaced
        //
        delete targetNode;

        // Create new headers + proxies and replace in fslices
        //
        newLimits = pLimits;
        newLimits.AddLimit(faxis,fminExtent+sliceWidth*minNo,
                           fminExtent+sliceWidth*(maxNo+1));
        replaceHeader = new G4SmartVoxelHeader(pVolume,newLimits,
                                               targetList,replaceNo);
        if (!replaceHeader)
        {
          G4Exception("G4SmartVoxelHeader::RefineNodes()", "GeomMgt0003",
                      FatalException, "Refined VoxelHeader allocation error.");
          return;
        }
        replaceHeader->SetMinEquivalentSliceNo(minNo);
        replaceHeader->SetMaxEquivalentSliceNo(maxNo);
        replaceHeaderProxy = new G4SmartVoxelProxy(replaceHeader);
        if (!replaceHeaderProxy)
        {
          G4Exception("G4SmartVoxelHeader::RefineNodes()", "GeomMgt0003",
                      FatalException, "Refined VoxelProxy allocation error.");
          return;
        }
        for (replaceNo=minNo; replaceNo<=maxNo; replaceNo++)
        {
          fslices[replaceNo] = replaceHeaderProxy;
        }
        // Finished replacing current `equivalent' group
        //
        delete targetList;
        targetNo=maxNo;
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void G4SmartVoxelHeader::SetMaxEquivalentSliceNo

(

G4int

pMax

)

Here is the caller graph for this function:

void G4SmartVoxelHeader::SetMinEquivalentSliceNo

(

G4int

pMin

)

Here is the caller graph for this function:

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  s, const G4SmartVoxelHeader &  h  )

friend

Definition at line 1304 of file G4SmartVoxelHeader.cc.

{
  os << "Axis = " << G4int(h.faxis) << G4endl;
  G4SmartVoxelProxy *collectNode=0, *collectHead=0;
  G4int collectNodeNo=0;
  G4int collectHeadNo=0;
  size_t i, j;
  G4bool haveHeaders=false;

  for (i=0; i<h.fslices.size(); i++)
  {
    os << "Slice #" << i << " = ";
    if (h.fslices[i]->IsNode())
    {
      if (h.fslices[i]!=collectNode)
      {
        os << "{";
        for (G4int k=0; k<h.fslices[i]->GetNode()->GetNoContained(); k++)
        {
          os << " " << h.fslices[i]->GetNode()->GetVolume(k);
        }
        os << " }" << G4endl;
        collectNode = h.fslices[i];
        collectNodeNo = i;
      }
      else
      {
        os << "As slice #" << collectNodeNo << G4endl;
      }
    }
    else
    {
      haveHeaders=true;
      if (h.fslices[i] != collectHead)
      {
        os << "Header" << G4endl;
        collectHead = h.fslices[i];
        collectHeadNo = i;
      }
      else
      {
        os << "As slice #" << collectHeadNo << G4endl;
      }
    }
  }

  if (haveHeaders)
  {
    collectHead=0;
    for (j=0; j<h.fslices.size(); j++)
    {
      if (h.fslices[j]->IsHeader())
      {
        os << "Header at Slice #" << j << " = ";
        if (h.fslices[j] != collectHead)
        {
          os << G4endl 
             << (*(h.fslices[j]->GetHeader()));
          collectHead = h.fslices[j];
          collectHeadNo = j;
        }
        else
        {
          os << "As slice #" << collectHeadNo << G4endl;
        }
      }
    }
  }
  return os;
}

Member Data Documentation

EAxis G4SmartVoxelHeader::faxis

protected

Definition at line 191 of file G4SmartVoxelHeader.hh.

G4int G4SmartVoxelHeader::fmaxEquivalent

protected

Definition at line 188 of file G4SmartVoxelHeader.hh.

G4double G4SmartVoxelHeader::fmaxExtent

protected

Definition at line 194 of file G4SmartVoxelHeader.hh.

G4int G4SmartVoxelHeader::fminEquivalent

protected

Definition at line 187 of file G4SmartVoxelHeader.hh.

G4double G4SmartVoxelHeader::fminExtent

protected

Definition at line 195 of file G4SmartVoxelHeader.hh.

EAxis G4SmartVoxelHeader::fparamAxis

protected

Definition at line 191 of file G4SmartVoxelHeader.hh.

G4ProxyVector G4SmartVoxelHeader::fslices

protected

Definition at line 198 of file G4SmartVoxelHeader.hh.

---

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

• source/geant4.10.03.p03/source/geometry/management/include/G4SmartVoxelHeader.hh
• source/geant4.10.03.p03/source/geometry/management/src/G4SmartVoxelHeader.cc