G4Box.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// Implementation for G4Box class
//
//  30.06.95 - P.Kent: First version
//  20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
//  18.04.17 - E.Tcherniaev: complete revision, speed-up
// --------------------------------------------------------------------
 
#include "G4Box.hh"
 
#if !defined(G4GEOM_USE_UBOX)
 
#include "G4SystemOfUnits.hh"
#include "G4VoxelLimits.hh"
#include "G4AffineTransform.hh"
#include "G4BoundingEnvelope.hh"
#include "G4QuickRand.hh"
 
#include "G4VPVParameterisation.hh"
 
#include "G4VGraphicsScene.hh"
#include "G4VisExtent.hh"
#include "G4AutoLock.hh"
 
namespace
{
  G4Mutex boxMutex = G4MUTEX_INITIALIZER;
}
 
////////////////////////////////////////////////////////////////////////
//
// Constructor - check & set half widths
 
G4Box::G4Box(const G4String& pName,
                   G4double pX,
                   G4double pY,
                   G4double pZ)
  : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
{
  delta = 0.5*kCarTolerance;
  if (pX < 2*kCarTolerance ||
      pY < 2*kCarTolerance ||
      pZ < 2*kCarTolerance)  // limit to thickness of surfaces
  {
    std::ostringstream message;
    message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
            << "     hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
    G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
  }
}
 
//////////////////////////////////////////////////////////////////////////
//
// Fake default constructor - sets only member data and allocates memory
//                            for usage restricted to object persistency
 
G4Box::G4Box( __void__& a )
  : G4CSGSolid(a)
{
}
 
//////////////////////////////////////////////////////////////////////////
//
// Assignment operator
 
G4Box& G4Box::operator = (const G4Box& rhs)
{
   // Check assignment to self
   //
   if (this == &rhs)  { return *this; }
 
   // Copy base class data
   //
   G4CSGSolid::operator=(rhs);
 
   // Copy data
   //
   fDx = rhs.fDx;
   fDy = rhs.fDy;
   fDz = rhs.fDz;
   delta = rhs.delta;
 
   return *this;
}
 
//////////////////////////////////////////////////////////////////////////
//
//  Set X dimension
 
void G4Box::SetXHalfLength(G4double dx)
{
  if(dx > 2*kCarTolerance)  // limit to thickness of surfaces
  {
    fDx = dx;
  }
  else
  {
    std::ostringstream message;
    message << "Dimension X too small for solid: " << GetName() << "!"
            << G4endl
            << "       hX = " << dx;
    G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
                FatalException, message);
  }
  fCubicVolume = 0.;
  fSurfaceArea = 0.;
  fRebuildPolyhedron = true;
}
 
//////////////////////////////////////////////////////////////////////////
//
//  Set Y dimension
 
void G4Box::SetYHalfLength(G4double dy)
{
  if(dy > 2*kCarTolerance)  // limit to thickness of surfaces
  {
    fDy = dy;
  }
  else
  {
    std::ostringstream message;
    message << "Dimension Y too small for solid: " << GetName() << "!\n"
            << "       hY = " << dy;
    G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
                FatalException, message);
  }
  fCubicVolume = 0.;
  fSurfaceArea = 0.;
  fRebuildPolyhedron = true;
}
 
//////////////////////////////////////////////////////////////////////////
//
//  Set Z dimension
 
void G4Box::SetZHalfLength(G4double dz)
{
  if(dz > 2*kCarTolerance)  // limit to thickness of surfaces
  {
    fDz = dz;
  }
  else
  {
    std::ostringstream message;
    message << "Dimension Z too small for solid: " << GetName() << "!\n"
            << "       hZ = " << dz;
    G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
                FatalException, message);
  }
  fCubicVolume = 0.;
  fSurfaceArea = 0.;
  fRebuildPolyhedron = true;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Dispatch to parameterisation for replication mechanism dimension
// computation & modification.
 
void G4Box::ComputeDimensions(G4VPVParameterisation* p,
                              const G4int n,
                              const G4VPhysicalVolume* pRep)
{
  p->ComputeDimensions(*this,n,pRep);
}
 
//////////////////////////////////////////////////////////////////////////
//
// Get bounding box
 
void G4Box::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
{
  pMin.set(-fDx,-fDy,-fDz);
  pMax.set( fDx, fDy, fDz);
 
  // Check correctness of the bounding box
  //
  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
  {
    std::ostringstream message;
    message << "Bad bounding box (min >= max) for solid: "
            << GetName() << " !"
            << "\npMin = " << pMin
            << "\npMax = " << pMax;
    G4Exception("G4Box::BoundingLimits()", "GeomMgt0001", JustWarning, message);
    DumpInfo();
  }
}
 
//////////////////////////////////////////////////////////////////////////
//
// Calculate extent under transform and specified limit
 
G4bool G4Box::CalculateExtent(const EAxis pAxis,
                              const G4VoxelLimits& pVoxelLimit,
                              const G4AffineTransform& pTransform,
                                    G4double& pMin, G4double& pMax) const
{
  G4ThreeVector bmin, bmax;
 
  // Get bounding box
  BoundingLimits(bmin,bmax);
 
  // Find extent
  G4BoundingEnvelope bbox(bmin,bmax);
  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
}
 
//////////////////////////////////////////////////////////////////////////
//
// Return whether point inside/outside/on surface, using tolerance
 
EInside G4Box::Inside(const G4ThreeVector& p) const
{
  G4double dist = std::max(std::max(
                  std::abs(p.x())-fDx,
                  std::abs(p.y())-fDy),
                  std::abs(p.z())-fDz);
  return (dist > delta) ? kOutside :
    ((dist > -delta) ? kSurface : kInside);
}
 
//////////////////////////////////////////////////////////////////////////
//
// Detect the side(s) and return corresponding normal
 
G4ThreeVector G4Box::SurfaceNormal(const G4ThreeVector& p) const
{
  G4double px = p.x(), py = p.y(), pz = p.z();
  G4ThreeVector norm(0.,0.,0.);
  if (std::abs(std::abs(px)-fDx) <= delta) { norm.setX(std::copysign(1.,px)); }
  if (std::abs(std::abs(py)-fDy) <= delta) { norm.setY(std::copysign(1.,py)); }
  if (std::abs(std::abs(pz)-fDz) <= delta) { norm.setZ(std::copysign(1.,pz)); }
 
  G4double nside = norm.mag2(); // number of sides = magnitude squared
  if (nside == 1)
  {
    return norm;
  }
  if (nside > 1)
  {
    return norm.unit(); // edge or corner
  }
 
  // Point is not on the surface
  //
#ifdef G4CSGDEBUG
  std::ostringstream message;
  G4int oldprc = message.precision(16);
  message << "Point p is not on surface (!?) of solid: "
          << GetName() << G4endl;
  message << "Position:\n";
  message << "   p.x() = " << p.x()/mm << " mm\n";
  message << "   p.y() = " << p.y()/mm << " mm\n";
  message << "   p.z() = " << p.z()/mm << " mm";
  G4cout.precision(oldprc);
  G4Exception("G4Box::SurfaceNormal(p)", "GeomSolids1002",
              JustWarning, message );
  DumpInfo();
#endif
 
  return ApproxSurfaceNormal(p);
}
 
//////////////////////////////////////////////////////////////////////////
//
// Algorithm for SurfaceNormal() following the original specification
// for points not on the surface
 
G4ThreeVector G4Box::ApproxSurfaceNormal(const G4ThreeVector& p) const
{
  G4double distx = std::abs(p.x()) - fDx;
  G4double disty = std::abs(p.y()) - fDy;
  G4double distz = std::abs(p.z()) - fDz;
 
  if (distx >= disty && distx >= distz)
  {
    return {std::copysign(1.,p.x()), 0., 0.};
  }
  if (disty >= distx && disty >= distz)
  {
    return {0., std::copysign(1.,p.y()), 0.};
  }
 
  return {0., 0., std::copysign(1.,p.z())};
}
 
//////////////////////////////////////////////////////////////////////////
//
// Calculate distance to box from an outside point
// - return kInfinity if no intersection
//
 
G4double G4Box::DistanceToIn(const G4ThreeVector& p,
                             const G4ThreeVector& v) const
{
  // Check if point is on the surface and traveling away
  //
  if ((std::abs(p.x())-fDx) >= -delta && p.x()*v.x() >= 0) { return kInfinity; }
  if ((std::abs(p.y())-fDy) >= -delta && p.y()*v.y() >= 0) { return kInfinity; }
  if ((std::abs(p.z())-fDz) >= -delta && p.z()*v.z() >= 0) { return kInfinity; }
 
  // Find intersection
  //
  G4double invx = (v.x() == 0) ? DBL_MAX : -1./v.x();
  G4double dx = std::copysign(fDx,invx);
  G4double txmin = (p.x() - dx)*invx;
  G4double txmax = (p.x() + dx)*invx;
 
  G4double invy = (v.y() == 0) ? DBL_MAX : -1./v.y();
  G4double dy = std::copysign(fDy,invy);
  G4double tymin = std::max(txmin,(p.y() - dy)*invy);
  G4double tymax = std::min(txmax,(p.y() + dy)*invy);
 
  G4double invz = (v.z() == 0) ? DBL_MAX : -1./v.z();
  G4double dz = std::copysign(fDz,invz);
  G4double tmin = std::max(tymin,(p.z() - dz)*invz);
  G4double tmax = std::min(tymax,(p.z() + dz)*invz);
 
  if (tmax <= tmin + delta) { return kInfinity; } // touch or no hit
 
  return (tmin < delta) ? 0. : tmin;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Appoximate distance to box.
// Returns largest perpendicular distance to the closest x/y/z sides of
// the box, which is the most fast estimation of the shortest distance to box
// - If inside return 0
 
G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
{
  G4double dist = std::max(std::max(
                  std::abs(p.x())-fDx,
                  std::abs(p.y())-fDy),
                  std::abs(p.z())-fDz);
  return (dist > 0) ? dist : 0.;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Calculate distance to surface of the box from inside and
// find normal at exit point, if required
// - when leaving the surface, return 0
 
G4double G4Box::DistanceToOut(const G4ThreeVector& p,
                              const G4ThreeVector& v,
                              const G4bool calcNorm,
                              G4bool* validNorm, G4ThreeVector* n) const
{
  G4double px = p.x(), vx = v.x();
  G4double py = p.y(), vy = v.y();
  G4double pz = p.z(), vz = v.z();
 
  if (!calcNorm) // calculation of normal is not needed
  {
    if ((std::abs(px) - fDx) >= -delta && px*vx > 0) { return 0.; }
    if ((std::abs(py) - fDy) >= -delta && py*vy > 0) { return 0.; }
    if ((std::abs(pz) - fDz) >= -delta && pz*vz > 0) { return 0.; }
    G4double tx = (vx == 0) ? DBL_MAX : (std::copysign(fDx,vx) - px)/vx;
    G4double ty = (vy == 0) ? DBL_MAX : (std::copysign(fDy,vy) - py)/vy;
    G4double tz = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - pz)/vz;
    G4double tmax = std::min(std::min(tx,ty),tz);
    return tmax;
  }
 
  *validNorm = true;
  // Check if point is on the surface and traveling away
  if ((std::abs(px) - fDx) >= -delta && px*vx > 0)
  {
    n->set(std::copysign(1.,px), 0., 0.);
    return 0.;
  }
  if ((std::abs(py) - fDy) >= -delta && py*vy > 0)
  {
    n->set(0., std::copysign(1.,py), 0.);
    return 0.;
  }
  if ((std::abs(pz) - fDz) >= -delta && pz*vz > 0)
  {
    n->set(0., 0., std::copysign(1.,pz));
    return 0.;
  }
 
  // Find intersection
  G4double tx = (vx == 0) ? DBL_MAX : (std::copysign(fDx,vx) - px)/vx;
  G4double ty = (vy == 0) ? DBL_MAX : (std::copysign(fDy,vy) - py)/vy;
  G4double tz = (vz == 0) ? DBL_MAX : (std::copysign(fDz,vz) - pz)/vz;
  G4double tmax = std::min(std::min(tx, ty), tz);
 
  // Find normal
  G4bool pickZ = (tmax == tz);
  G4bool pickX = (!pickZ) && (tmax == tx);
  G4bool pickY = (!pickZ) && (!pickX);
  G4double nz = std::copysign((G4double)pickZ, vz);
  G4double nx = std::copysign((G4double)pickX, vx);
  G4double ny = std::copysign((G4double)pickY, vy);
  n->set(nx, ny, nz);
  return tmax;
}
 
////////////////////////////////////////////////////////////////////////////
//
// Calculate exact shortest distance to any boundary from inside
// - if outside return 0
 
G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
{
#ifdef G4CSGDEBUG
  if( Inside(p) == kOutside )
  {
    std::ostringstream message;
    G4int oldprc = message.precision(16);
    message << "Point p is outside (!?) of solid: " << GetName() << G4endl;
    message << "Position:\n";
    message << "   p.x() = " << p.x()/mm << " mm\n";
    message << "   p.y() = " << p.y()/mm << " mm\n";
    message << "   p.z() = " << p.z()/mm << " mm";
    G4cout.precision(oldprc);
    G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
                JustWarning, message );
    DumpInfo();
  }
#endif
  G4double dist = std::min(std::min(
                  fDx-std::abs(p.x()),
                  fDy-std::abs(p.y())),
                  fDz-std::abs(p.z()));
  return (dist > 0) ? dist : 0.;
}
 
//////////////////////////////////////////////////////////////////////////
//
// GetEntityType
 
G4GeometryType G4Box::GetEntityType() const
{
  return {"G4Box"};
}
 
//////////////////////////////////////////////////////////////////////////
//
// IsFaceted
 
G4bool G4Box::IsFaceted() const
{
  return true;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Stream object contents to an output stream
 
std::ostream& G4Box::StreamInfo(std::ostream& os) const
{
  G4long oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << "Solid type: G4Box\n"
     << "Parameters: \n"
     << "   half length X: " << fDx/mm << " mm \n"
     << "   half length Y: " << fDy/mm << " mm \n"
     << "   half length Z: " << fDz/mm << " mm \n"
     << "-----------------------------------------------------------\n";
  os.precision(oldprc);
  return os;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Return a point randomly and uniformly selected on the surface
 
G4ThreeVector G4Box::GetPointOnSurface() const
{
  G4double sxy = fDx*fDy, sxz = fDx*fDz, syz = fDy*fDz;
  G4double select = (sxy + sxz + syz)*G4QuickRand();
  G4double u = 2.*G4QuickRand() - 1.;
  G4double v = 2.*G4QuickRand() - 1.;
 
  G4double x, y, z;
  if (select < sxy)
  {
    x = u*fDx;
    y = v*fDy;
    z = (select < 0.5*sxy) ? -fDz : fDz;
  }
  else if (select < sxy + sxz)
  {
    x = u*fDx;
    y = (select < sxy + 0.5*sxz) ? -fDy : fDy;
    z = v*fDz;
  }
  else
  {
    x = (select < sxy + sxz + 0.5*syz) ? -fDx : fDx;
    y = u*fDy;
    z = v*fDz;
  }
  return { x, y, z };
}
 
//////////////////////////////////////////////////////////////////////////
//
// Computes/returns volume capacity
//
G4double G4Box::GetCubicVolume()
{
  if (fCubicVolume == 0)
  {
    G4AutoLock l(&boxMutex);
    fCubicVolume = 8*fDx*fDy*fDz;
    l.unlock();
  }
  return fCubicVolume;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Computes/returns surface area
//
G4double G4Box::GetSurfaceArea()
{
  if (fSurfaceArea == 0)
  {
    G4AutoLock l(&boxMutex);
    fSurfaceArea = 8*(fDx*fDy+fDx*fDz+fDy*fDz);
    l.unlock();
  }
  return fSurfaceArea;
}
 
//////////////////////////////////////////////////////////////////////////
//
// Make a clone of the object
//
G4VSolid* G4Box::Clone() const
{
  return new G4Box(*this);
}
 
//////////////////////////////////////////////////////////////////////////
//
// Methods for visualisation
 
void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const
{
  scene.AddSolid (*this);
}
 
G4VisExtent G4Box::GetExtent() const
{
  return { -fDx, fDx, -fDy, fDy, -fDz, fDz };
}
 
G4Polyhedron* G4Box::CreatePolyhedron () const
{
  return new G4PolyhedronBox (fDx, fDy, fDz);
}
#endif