UOrb Class Reference

#include <UOrb.hh>

Inheritance diagram for UOrb:

Public Member Functions
	UOrb ()
	UOrb (const std::string &name, double pRmax)
	~UOrb ()
	UOrb (const UOrb &rhs)
UOrb &	operator= (const UOrb &rhs)
double	GetRadius () const
void	SetRadius (double newRmax)
EnumInside	Inside (const UVector3 &aPo6int) const
double	SafetyFromInside (const UVector3 &aPoint, bool aAccurate=false) const
double	SafetyFromOutside (const UVector3 &aPoint, bool aAccurate=false) const
double	DistanceToIn (const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
double	DistanceToOut (const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
bool	Normal (const UVector3 &aPoint, UVector3 &aNormal) const
void	Extent (UVector3 &aMin, UVector3 &aMax) const
double	Capacity ()
double	SurfaceArea ()
UGeometryType	GetEntityType () const
void	ComputeBBox (UBBox *, bool)
void	GetParametersList (int, double *) const
VUSolid *	Clone () const
double	GetRadialTolerance ()
UVector3	GetPointOnSurface () const
std::ostream &	StreamInfo (std::ostream &os) const
Public Member Functions inherited from VUSolid
	VUSolid ()
	VUSolid (const std::string &name)
virtual	~VUSolid ()
double	GetCarTolerance () const
double	GetRadTolerance () const
double	GetAngTolerance () const
void	SetCarTolerance (double eps)
void	SetRadTolerance (double eps)
void	SetAngTolerance (double eps)
virtual void	ExtentAxis (EAxisType aAxis, double &aMin, double &aMax) const
const std::string &	GetName () const
void	SetName (const std::string &aName)
virtual void	SamplePointsInside (int, UVector3 *) const
virtual void	SamplePointsOnSurface (int, UVector3 *) const
virtual void	SamplePointsOnEdge (int, UVector3 *) const
double	EstimateCubicVolume (int nStat, double epsilon) const
double	EstimateSurfaceArea (int nStat, double ell) const

Additional Inherited Members
Public Types inherited from VUSolid
enum	EnumInside { eInside =0, eSurface =1, eOutside =2 }
enum	EAxisType { eXaxis =0, eYaxis =1, eZaxis =2 }
Static Public Member Functions inherited from VUSolid
static double	Tolerance ()
Static Protected Attributes inherited from VUSolid
static double	fgTolerance = 1.0E-9
static double	frTolerance = 1.0E-9
static double	faTolerance = 1.0E-9

Detailed Description

Definition at line 28 of file UOrb.hh.

Constructor & Destructor Documentation

UOrb::UOrb ( )

inline

Definition at line 32 of file UOrb.hh.

Referenced by Clone().

: VUSolid(), fR(0), fRTolerance(0) {}

UOrb::UOrb	(	const std::string &	name,
		double	pRmax
	)

Definition at line 26 of file UOrb.cc.

References UUtils::Exception(), FatalErrorInArguments, VUSolid::fgTolerance, VUSolid::frTolerance, and G4INCL::Math::max().

  : VUSolid(name), fR(r), fCubicVolume(0), fSurfaceArea(0)
{
  const double epsilon = 2.e-11;  // relative tolerance of fR

  // Check radius
  //
  if (r < 10 * VUSolid::fgTolerance) // cartesian tolerance
  {
    UUtils::Exception("G4Orb::G4Orb()", "InvalidSetup", FatalErrorInArguments, 1, "Invalid radius > 10*kCarTolerance.");
  }
  // VUSolid::fRTolerance is radial tolerance (note: half of G4 tolerance)
  fRTolerance =  max(VUSolid::frTolerance, epsilon * r);
}

UOrb::~UOrb ( )

inline

Definition at line 34 of file UOrb.hh.

{}

UOrb::UOrb

(

const UOrb &

rhs

)

Definition at line 496 of file UOrb.cc.

  : VUSolid(rhs), fR(rhs.fR), fRTolerance(rhs.fRTolerance), fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea)
{
}

Member Function Documentation

double UOrb::Capacity ( )

inlinevirtual

Implements VUSolid.

Definition at line 102 of file UOrb.hh.

{
  if (fCubicVolume != 0.)
  {
    ;
  }
  else
  {
    fCubicVolume = (4 * UUtils::kPi / 3) * fR * fR * fR;
  }
  return fCubicVolume;
}

VUSolid * UOrb::Clone ( ) const

virtual

Implements VUSolid.

Definition at line 489 of file UOrb.cc.

References UOrb().

{
  return new UOrb(*this);
}

void UOrb::ComputeBBox ( UBBox *  , bool    )

inlinevirtual

Implements VUSolid.

Definition at line 67 of file UOrb.hh.

{}

double UOrb::DistanceToIn ( const UVector3 &  aPoint, const UVector3 &  aDirection, double  aPstep = UUtils::kInfinity  ) const

virtual

Implements VUSolid.

Definition at line 84 of file UOrb.cc.

References test::c, UUtils::Exception(), python.hepunit::rad, Warning, UVector3::x, UVector3::y, and UVector3::z.

{
  double snxt = UUtils::kInfinity;      // snxt = default return value

  double rad, pDotV3d; // , tolORMax2, tolIRMax2;
  double c, d2, s = UUtils::kInfinity;

  const double dRmax = 100.*fR;

  // General Precalcs

  rad    = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
  pDotV3d = p.x * v.x + p.y * v.y + p.z * v.z;

  // Radial Precalcs

  // tolORMax2 = (fR+fRTolerance*0.5)*(fR+fRTolerance*0.5);
  // tolIRMax2 = (fR-fRTolerance*0.5)*(fR-fRTolerance*0.5);

  // Outer spherical shell intersection
  // - Only if outside tolerant fR
  // - Check for if inside and outer G4Orb heading through solid (-> 0)
  // - No intersect -> no intersection with G4Orb
  //
  // Shell eqn: x^2+y^2+z^2 = RSPH^2
  //
  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
  //
  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
  // =>      rad2        +2s(pDotV3d)       +s^2                =R^2
  //
  // => s=-pDotV3d+-sqrt(pDotV3d^2-(rad2-R^2))

  c = (rad - fR) * (rad + fR); // c = (rad2-R^2))

  if (rad > fR - fRTolerance * 0.5) // not inside in terms of Inside(p)
  {
    if (c > fRTolerance * fR)
    {
      // If outside tolerant boundary of outer G4Orb in terms of c
      // [ should be sqrt(rad2) - fR > fRTolerance*0.5 ]

      d2 = pDotV3d * pDotV3d - c;

      if (d2 >= 0)
      {
        s = -pDotV3d - sqrt(d2); // ok! = [ ( -2 p dot v) +- sqrt [(-2p dot v)2 - 4*(rad - fR)*(rad + fR)] ] /  2
        // pDotV3d must be positive always, if not use alternative http://en.wikipedia.org/wiki/Quadratic_equation#Alternative_quadratic_formula
        if (s >= 0)
        {
          if (s > dRmax)   // Avoid rounding errors due to precision issues seen on
          {
            // 64 bits systems. Split long distances and recompute
            double fTerm = s - fmod(s, dRmax);
            s = fTerm + DistanceToIn(p + fTerm * v, v);
          }
          return snxt = s;
        }
      }
      else    // No intersection with UOrb
      {
        return snxt = UUtils::kInfinity;
      }
    }
    else // not outside in terms of c
    {
      if (c > -fRTolerance * fR)  // on surface
      {
        d2 = pDotV3d * pDotV3d - c;
        if ((d2 < fRTolerance * fR) || (pDotV3d >= 0)) // pDotV3d = cos si >= 0
        {
          return snxt = UUtils::kInfinity;
        }
        else
        {
          return snxt = 0.;
        }
      }
    }
  }
#ifdef UDEBUG
  else // inside ???
  {
    UUtils::Exception("UOrb::DistanceToIn(p,v)", "Notification", Warning, 1, "Point p is inside !?");
  }
#endif

  return snxt;
}

double UOrb::DistanceToOut ( const UVector3 &  aPoint, const UVector3 &  aDirection, UVector3 &  aNormalVector, bool &  aConvex, double  aPstep = UUtils::kInfinity  ) const

virtual

Implements VUSolid.

Definition at line 203 of file UOrb.cc.

References test::c, UUtils::Exception(), python.hepunit::rad, Warning, UVector3::x, UVector3::y, and UVector3::z.

{
  double snxt = 0;     // snxt: distance to next surface, is default return value
  bool notOutside = false;
  convex = true; // orb is always convex, if we leave surface of Orb, we will neber bump on the orb again ...

  double rad2, pDotV3d;
  double xi, yi, zi;    // Intersection point
  double c, d2;

  rad2    = p.x * p.x + p.y * p.y + p.z * p.z;
  pDotV3d = p.x * v.x + p.y * v.y + p.z * v.z;

  // Radial Intersection from UOrb::DistanceToIn
  //
  // Outer spherical shell intersection
  // - Only if outside tolerant fR
  // - Check for if inside and outer UOrb heading through solid (-> 0)
  // - No intersect -> no intersection with UOrb
  //
  // Shell eqn: x^2+y^2+z^2=RSPH^2
  //
  // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2
  //
  // => (px^2+py^2+pz^2) +2s(pxvx+pyvy+pzvz)+s^2(vx^2+vy^2+vz^2)=R^2
  // =>      rad2        +2s(pDotV3d)       +s^2                =R^2
  //
  // => s=-pDotV3d+-sqrt(pDotV3d^2-(rad2-R^2))

  const double rPlus = fR + fRTolerance;
  double rad = sqrt(rad2);

  if (rad <= rPlus)
  {
    c = (rad - fR) * (rad + fR); // rad2 - fR2

    if (c < fRTolerance * fR)
    {
      // Within tolerant Outer radius
      //
      // The test is
      //     rad  - fR < 0.5*fRTolerance
      // =>  rad  < fR + 0.5*kRadTol
      // =>  rad2 < (fR + 0.5*kRadTol)^2
      // =>  rad2 < fR^2 + 2.*0.5*fR*kRadTol + 0.25*kRadTol*kRadTol
      // =>  rad2 - fR^2    <~    fR*kRadTol

      d2 = pDotV3d * pDotV3d - c;

      if ((c > -2 * fRTolerance * fR) &&      // => point is on tolerant surface (i.e. within +- tolerance)
          ((pDotV3d >= 0)   || (d2 < 0)))         // if (pDotV3d >= 0 ) => leaving outside from Rmax; i.e. from surface
        // not re-entering
        // if (d2 < 0) => it means the point is already outside
      {
        // if(calcNorm) // NOTE: we do not have this variable, calcNorm is true always
        {
          //          *validNorm = true; // NOTE: we do not have this variable, probably always true
          n = UVector3(p.x / fR, p.y / fR, p.z / fR);
        }
        return snxt = 0;
      }
      else
      {
        // we are inside, with + version of quadratic eq. solution we calculate solution for distance
        snxt = -pDotV3d + sqrt(d2);    // second root since inside Rmax
        // the solution is safe because pDotV3d is negative
        // c alternative formula, see http://en.wikipedia.org/wiki/Quadratic_equation#Alternative_quadratic_formula
        // is not neccessary in this case
        notOutside = true;
      }
    }
  }
  else // p is outside ???
  {
    // Rule 2: DistanceToOut
    // Surface points = 0
    // Outside points: If pointing outwards (dot product with normal positive) return 0, otherwise ignore first surface, another surface should always be on the direction line, use instead distance to this one.

    //    double res = DistanceToOutForOutsidePoints(p, v, n);

    cout.precision(16);
    cout << endl;
    //    DumpInfo();
    cout << "Position:"  << endl << endl;
    cout << "p.x() = "   << p.x << endl;
    cout << "p.y() = "   << p.y << endl;
    cout << "p.z() = "   << p.z << endl << endl;
    cout << "Rp = " << sqrt(p.x * p.x + p.y * p.y + p.z * p.z) << endl << endl;
    cout << "Direction:" << endl << endl;
    cout << "v.x() = "   << v.x << endl;
    cout << "v.y() = "   << v.y << endl;
    cout << "v.z() = "   << v.z << endl << endl;
    cout << "Proposed distance :" << endl << endl;
    cout << "snxt = "    << snxt << endl << endl;



    cout.precision(6);
    UUtils::Exception("UOrb::DistanceToOut(p,v,..)", "Notification",
                      Warning, 1, "Logic error: snxt = kInfinity ???");

  }

  // if (calcNorm)    // Output switch operator
  {
    if (notOutside)
    {
      xi = p.x + snxt * v.x; // we move to the point on surface, then return normal at that point which for orb, see method bool UOrb::Normal( const UVector3& p, UVector3 &n)
      yi = p.y + snxt * v.y;
      zi = p.z + snxt * v.z;
      n = UVector3(xi / fR, yi / fR, zi / fR); // we return normalized vector
    }
    else
    {

      cout.precision(16);
      cout << endl;
      //        DumpInfo();
      cout << "Position:"  << endl << endl;
      cout << "p.x() = "   << p.x << " mm" << endl;
      cout << "p.y() = "   << p.y << " mm" << endl;
      cout << "p.z() = "   << p.z << " mm" << endl << endl;
      cout << "Direction:" << endl << endl;
      cout << "v.x() = "   << v.x << endl;
      cout << "v.y() = "   << v.y << endl;
      cout << "v.z() = "   << v.z << endl << endl;
      cout << "Proposed distance :" << endl << endl;
      cout << "snxt = "    << snxt << " mm" << endl << endl;
      cout.precision(6);


      UUtils::Exception("UOrb::DistanceToOut(p,v,..)", "Notification", Warning, 1, "Undefined side for valid surface normal to solid.");
    }
  }
  return snxt;
}

void UOrb::Extent ( UVector3 &  aMin, UVector3 &  aMax  ) const

virtual

Returns extent of the solid along a given cartesian axis OK

Implements VUSolid.

Definition at line 443 of file UOrb.cc.

References UVector3::Set().

{
  aMin.Set(-fR);
  aMax.Set(fR);
}

UGeometryType UOrb::GetEntityType ( ) const

virtual

Implements VUSolid.

Definition at line 538 of file UOrb.cc.

{
   return "Orb";
}

void UOrb::GetParametersList ( int  , double *  aArray  ) const

virtual

Implements VUSolid.

Definition at line 530 of file UOrb.cc.

References GetRadius().

{
  aArray[0] = GetRadius();
}

UVector3 UOrb::GetPointOnSurface ( ) const

virtual

Implements VUSolid.

Definition at line 474 of file UOrb.cc.

References UUtils::Random(), and UUtils::sqr().

{
  //  generate a random number from zero to 2UUtils::kPi...
  //
  double phi      = UUtils::Random(0., 2.*UUtils::kPi);
  double cosphi  = std::cos(phi);
  double sinphi  = std::sin(phi);

  // generate a random point uniform in area
  double costheta = UUtils::Random(-1., 1.);
  double sintheta = std::sqrt(1. - UUtils::sqr(costheta));

  return UVector3(fR * sintheta * cosphi, fR * sintheta * sinphi, fR * costheta);
}

double UOrb::GetRadialTolerance ( )

inline

Definition at line 74 of file UOrb.hh.

    {
      return fRTolerance;
    }

double UOrb::GetRadius ( ) const

inline

Definition at line 93 of file UOrb.hh.

Referenced by GetParametersList(), and G4UOrb::GetRadius().

{
  return fR;
}

VUSolid::EnumInside UOrb::Inside ( const UVector3 &  p ) const

virtual

Return whether point inside/outside/on surface Split into radius checks

Classify point location with respect to solid: o eInside - inside the solid o eSurface - close to surface within tolerance o eOutside - outside the solid

Implements VUSolid.

Definition at line 53 of file UOrb.cc.

References VUSolid::eInside, VUSolid::eOutside, VUSolid::eSurface, UVector3::x, UVector3::y, and UVector3::z.

Referenced by SafetyFromInside().

{
  double rad2 = p.x * p.x + p.y * p.y + p.z * p.z;
  //  if (false) double rad = sqrt(rad2);

  double tolRMax = fR - fRTolerance * 0.5;

  // Check radial surface
  double tolRMax2 = tolRMax * tolRMax;
  if (rad2 <= tolRMax2)
    return eInside;
  else
  {
    tolRMax = fR + fRTolerance * 0.5;
    tolRMax2 = tolRMax * tolRMax;
    if (rad2 <= tolRMax2)
      return eSurface;
    else
      return eOutside;
  }
}

bool UOrb::Normal ( const UVector3 &  p, UVector3 &  n  ) const

virtual

Return unit normal of surface closest to p

From http://lists.trolltech.com/qt-interest/2002-09/thread01124-0.html :

does anybody here have an algorithm to calculate the normal vector in a given point in space (x, y, z) in a sphere? I know that it's not about qt but i'll like very mutch the help.

It's simply the connecting vector from the centre of the sphere to the point (other way around for inward normals) obtained through vector subtraction, normalized to unity.

You really should get an algebra book though, as you are bound to encounter more of these problems in a 3d application.

Implements VUSolid.

Definition at line 421 of file UOrb.cc.

References python.hepunit::rad, UVector3::x, UVector3::y, and UVector3::z.

{
  double rad2 = p.x * p.x + p.y * p.y + p.z * p.z;
  double rad = sqrt(rad2);

  n = UVector3(p.x / rad, p.y / rad, p.z / rad);

  double tolRMaxP = fR + fRTolerance;
  double tolRMaxM = fR - fRTolerance;

  // Check radial surface
  bool result = ((rad2 <= tolRMaxP * tolRMaxP) && (rad2 >= tolRMaxM * tolRMaxM)); // means we are on surface
  return result;
}

UOrb & UOrb::operator=

(

const UOrb &

rhs

)

Definition at line 505 of file UOrb.cc.

{
  // Check assignment to self
  //
  if (this == &rhs)
  {
    return *this;
  }

  // Copy base class data
  //
  VUSolid::operator=(rhs);

  // Copy data
  //
  fR = rhs.fR;
  fRTolerance = rhs.fRTolerance;
  fCubicVolume = rhs.fCubicVolume;
  fSurfaceArea = rhs.fSurfaceArea;
  return *this;
}

double UOrb::SafetyFromInside ( const UVector3 &  aPoint, bool  aAccurate = false  ) const

virtual

Implements VUSolid.

Definition at line 350 of file UOrb.cc.

References UUtils::Exception(), Inside(), kOutside, python.hepunit::rad, Warning, UVector3::x, UVector3::y, and UVector3::z.

{

  /////////////////////////////////////////////////////////////////////////
  //
  // Calculate distance (<=actual) to closest surface of shape from inside

  double safe = 0.0, rad = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);

#ifdef UDEBUG
  if (Inside(p) == kOutside)
  {
    //     int oldprc = cout.precision(16);
    cout << endl;
    //     DumpInfo();
    cout << "Position:"  << endl << endl;
    cout << "p.x = "   << p.x << endl;
    cout << "p.y = "   << p.y << endl;
    cout << "p.z = "   << p.z << endl << endl;
    //     cout.precision(oldprc);
    UUtils::Exception("UOrb::DistanceToOut(p)", "Notification", Warning, 1,
                      "Point p is outside !?");
  }
#endif

  safe = fR - rad;
  if (safe < 0.) safe = 0.;
  return safe;
}

double UOrb::SafetyFromOutside ( const UVector3 &  aPoint, bool  aAccurate = false  ) const

virtual

Implements VUSolid.

Definition at line 393 of file UOrb.cc.

References python.hepunit::rad, UVector3::x, UVector3::y, and UVector3::z.

{
  double safe = 0.0;
  double rad  = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
  safe = rad - fR;
  if (safe < 0)
  {
    safe = 0.;
  }
  return safe;
}

void UOrb::SetRadius ( double  newRmax )

inline

Definition at line 97 of file UOrb.hh.

Referenced by G4UOrb::SetRadius().

{
  fR = newRmax;
}

std::ostream & UOrb::StreamInfo ( std::ostream &  os ) const

virtual

Implements VUSolid.

Definition at line 454 of file UOrb.cc.

References VUSolid::GetName().

{
  int oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "       *** Dump for solid - " << GetName() << " ***\n"
     << "       ===================================================\n"
     << " Solid type: UOrb\n"
     << " Parameters: \n"

     << "       outer radius: " << fR << " mm \n"
     << "-----------------------------------------------------------\n";
  os.precision(oldprc);

  return os;
}

double UOrb::SurfaceArea ( )

inlinevirtual

Implements VUSolid.

Definition at line 115 of file UOrb.hh.

{
  if (fSurfaceArea != 0.)
  {
    ;
  }
  else
  {
    fSurfaceArea = (4 * UUtils::kPi) * fR * fR;
  }
  return fSurfaceArea;
}

---

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

• UOrb.hh
• UOrb.cc