Problem 983

The implementation of CalculateExtent() in the G4TessellatedSolid is incomplete.  The method checks if the affine transform is pure translation, and if so, it computes the extent of the solid correctly.  However, if there is a rotation component to the affine transform, it branches to an empty else block and always returns false.  This causes very confusing and erratic tracking behavior when geometry optimization is turned on (the default), since G4SmartVoxelHeader::BuildNodes() assumes a working CalculateExtent() method.  Rotated G4TessellatedSolid and G4ExtrudedSolid objects will disappear or be clipped in strange ways when viewed with the raytracer.

Here is one possible implementation for the method which is slower for pure translations, but is sufficiently general to handle all affine transformations:

G4bool G4TessellatedSolid::CalculateExtent(const EAxis pAxis,
         const G4VoxelLimits& pVoxelLimit, const G4AffineTransform& pTransform,
               G4double& pMin, G4double& pMax) const
{
    G4ThreeVectorList transVertexList(vertexList);

    // Put solid into transformed frame
    for (size_t i=0; i<vertexList.size(); i++)
      pTransform.ApplyPointTransform(transVertexList[i]);

    // Find min and max extent in each dimension
    G4ThreeVector minExtent(kInfinity, kInfinity, kInfinity);
    G4ThreeVector maxExtent(-kInfinity, -kInfinity, -kInfinity);
    for (size_t i=0; i<transVertexList.size(); i++)
    {
      for (int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++)
      {
        G4double coordinate = transVertexList[i][axis];
        if (coordinate < minExtent[axis])
          minExtent[axis] = coordinate;
        if (coordinate > maxExtent[axis])
          maxExtent[axis] = coordinate;
      }
    }
    
    // Check for containment and clamp to voxel boundaries
    for (int axis=G4ThreeVector::X; axis < G4ThreeVector::SIZE; axis++)
    {
      EAxis geomAxis; // G4 geom classes use different index type
      switch(axis)
      {
        case G4ThreeVector::X: geomAxis = kXAxis; break;
        case G4ThreeVector::Y: geomAxis = kYAxis; break;
        case G4ThreeVector::Z: geomAxis = kZAxis; break;
      }
      G4bool isLimited = pVoxelLimit.IsLimited(geomAxis);
      G4double voxelMinExtent = pVoxelLimit.GetMinExtent(geomAxis);
      G4double voxelMaxExtent = pVoxelLimit.GetMaxExtent(geomAxis);
      
      if (isLimited)
      {
        if ( minExtent[axis] > voxelMaxExtent+kCarTolerance ||
             maxExtent[axis] < voxelMinExtent-kCarTolerance    ) return false ;
        else
        {
          if (minExtent[axis] < voxelMinExtent)
          {
            minExtent[axis] = voxelMinExtent ;
          }
          if (maxExtent[axis] > voxelMaxExtent)
          {
            maxExtent[axis] = voxelMaxExtent;
          }
        }
      }
    }
    
    // Convert pAxis into G4ThreeVector index
    int clhepAxis;
    switch(pAxis)
    {
      case kXAxis: clhepAxis = G4ThreeVector::X; break;
      case kYAxis: clhepAxis = G4ThreeVector::Y; break;
      case kZAxis: clhepAxis = G4ThreeVector::Z; break;
      default: break;
    } 
    
    pMin = minExtent[clhepAxis] - kCarTolerance;
    pMax = maxExtent[clhepAxis] + kCarTolerance;
    
    return true;  
}