55 using namespace CLHEP;
64 :
G4VSolid( name ), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
79 :
G4VSolid(a), dx(0.), dy(0.), dz(0.), halfTol(0.),
80 fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
98 :
G4VSolid(rhs), dx(rhs.dx), dy(rhs.dy), dz(rhs.dz), halfTol(rhs.halfTol),
99 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
112 if (
this == &rhs) {
return *
this; }
121 halfTol = rhs.halfTol;
122 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
158 G4double rFudge = 1.0/std::cos(0.5*sigPhi);
170 cosPhi = std::cos(phi),
171 sinPhi = std::sin(phi);
173 v1( dxFudge*cosPhi, dyFudge*sinPhi, -
dz ),
180 if (numPhi == 1) phi = 0;
181 cosPhi = std::cos(phi),
182 sinPhi = std::sin(phi);
210 phiPoly.
SetNormal( (v1-v0).cross(w0-v0).unit() );
220 }
while( --numPhi > 0 );
299 if ( (distZ < halfTol ) && ( distR1 <= 1 ) )
304 if( (distR1 <= 1 ) && ( distR2 >= 1 ) )
310 if ( noSurfaces == 0 )
313 G4Exception(
"G4EllipticalTube::SurfaceNormal(p)",
"GeomSolids1002",
316 norm = ApproxSurfaceNormal(p);
318 else if ( noSurfaces == 1 ) { norm = sumnorm; }
319 else { norm = sumnorm.
unit(); }
329 G4EllipticalTube::ApproxSurfaceNormal(
const G4ThreeVector&
p )
const
342 if (distZ*distZ < distR2)
387 if (sigz < 0)
return kInfinity;
393 if (
CheckXY(p.
x(),p.
y(),-halfTol) <= 1.0)
return kInfinity;
406 yi = p.
y() + q*v.
y();
416 return (sigz < -halfTol) ? q : 0;
418 else if (xi*dy*dy*v.
x() + yi*
dx*
dx*v.
y() >= 0)
438 if (sigz > 0)
return kInfinity;
439 if (
CheckXY(p.
x(),p.
y(),-halfTol) <= 1.0)
return kInfinity;
445 yi = p.
y() + q*v.
y();
449 return (sigz > -halfTol) ? q : 0;
451 else if (xi*dy*dy*v.
x() + yi*
dx*
dx*v.
y() >= 0)
464 if (n==0)
return kInfinity;
469 if (std::fabs(p.
z()) <
dz+halfTol) {
470 if (
CheckXY( p.
x(), p.
y(), halfTol ) < 1.0)
475 if (p.
x()*dy*dy*v.
x() + p.
y()*
dx*
dx*v.
y() < 0)
return 0;
484 if (q[0] < 0)
return kInfinity;
494 if (zi < -
dz)
return kInfinity;
498 if (zi > +
dz)
return kInfinity;
525 if (
CheckXY( p.
x(), p.
y(), +halfTol ) < 1.0)
531 if (p.
z() < -
dz-halfTol)
533 else if (p.
z() >
dz+halfTol)
553 G4double tnorm = std::sqrt( tx*tx + ty*ty );
558 G4double distR = ( (p.
x()-xe)*ty - (p.
y()-ye)*tx )/tnorm;
568 return std::sqrt( (p.
z()+
dz)*(p.
z()+
dz) + distR*distR );
570 return std::sqrt( (p.
z()-
dz)*(p.
z()-
dz) + distR*distR );
593 if (calcNorm) { *validNorm =
true; }
607 sBest = -(p.
z()+
dz)/v.
z();
612 if (p.
z() < -
dz+halfTol)
614 if (calcNorm) { *norm = normHere; }
637 if (p.
z() > +
dz-halfTol)
639 if (calcNorm) { *norm = normHere; }
646 if (q < sBest) { sBest = q; nBest = normHere; }
657 if (sBest == kInfinity)
660 std::ostringstream message;
661 G4int oldprc = message.precision(16) ;
662 message <<
"Point p is outside !?" <<
G4endl
664 <<
" p.x() = " << p.
x()/
mm <<
" mm" <<
G4endl
665 <<
" p.y() = " << p.
y()/
mm <<
" mm" <<
G4endl
666 <<
" p.z() = " << p.
z()/
mm <<
" mm" <<
G4endl
668 <<
" v.x() = " << v.
x() <<
G4endl
669 <<
" v.y() = " << v.
y() <<
G4endl
670 <<
" v.z() = " << v.
z() <<
G4endl
671 <<
"Proposed distance :" <<
G4endl
672 <<
" snxt = " << sBest/
mm <<
" mm";
673 message.precision(oldprc) ;
674 G4Exception(
"G4EllipticalTube::DistanceToOut(p,v,...)",
677 if (calcNorm) { *norm = nBest; }
680 else if (q[n-1] > sBest)
682 if (calcNorm) { *norm = nBest; }
699 if (
CheckXY( p.
x(), p.
y(), -halfTol ) > 1.0)
704 if (p.
x()*dy*dy*v.
x() + p.
y()*
dx*
dx*v.
y() > 0)
return 0;
732 if (sBest < halfTol)
return 0;
739 if (radical < +
DBL_MIN)
return 0;
742 if (p.
x() < 0) xi = -xi;
747 radical = 1.0 - p.
x()*p.
x()/
dx/
dx;
748 if (radical < +
DBL_MIN)
return 0;
750 G4double yi = dy*std::sqrt( radical );
751 if (p.
y() < 0) yi = -yi;
760 G4double normi = std::sqrt( xdi*xdi + ydi*ydi );
761 if (normi < halfTol)
return 0;
765 G4double q = 0.5*(xdi*(p.
y()-yi) - ydi*(p.
x()-xi));
766 if (xi*yi < 0) q = -q;
768 if (q < sBest) sBest = q;
773 return sBest < halfTol ? 0 : sBest;
817 if (radical < -
DBL_MIN)
return 0;
828 radical = std::sqrt(radical);
830 G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
833 if (sa < sb) { ss[0] = sa; ss[1] = sb; }
else { ss[0] = sb; ss[1] = sa; }
843 return G4String(
"G4EllipticalTube");
861 if(fCubicVolume != 0.) {;}
871 if(fSurfaceArea != 0.) {;}
881 G4int oldprc = os.precision(16);
882 os <<
"-----------------------------------------------------------\n"
883 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
884 <<
" ===================================================\n"
885 <<
" Solid type: G4EllipticalTube\n"
887 <<
" length Z: " <<
dz/
mm <<
" mm \n"
888 <<
" surface equation in X and Y: \n"
889 <<
" (X / " <<
dx <<
")^2 + (Y / " << dy <<
")^2 = 1 \n"
890 <<
"-----------------------------------------------------------\n";
891 os.precision(oldprc);
905 G4double xRand, yRand, zRand, phi, cosphi, sinphi, zArea, cArea,
p, chose;
908 cosphi = std::cos(phi);
909 sinphi = std::sin(phi);
918 p = 2.*
pi*std::sqrt(0.5*(
dx*
dx+dy*dy));
929 if( (chose>=0) && (chose < cArea) )
933 else if( (chose >= cArea) && (chose < cArea + zArea) )
936 yRand = std::sqrt(1.-
sqr(xRand/
dx));
943 yRand = std::sqrt(1.-
sqr(xRand/
dx));
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
ThreeVector shoot(const G4int Ap, const G4int Af)
CLHEP::Hep3Vector G4ThreeVector
G4double GetSurfaceArea()
virtual ~G4EllipticalTube()
G4bool GetExtent(G4double &min, G4double &max) const
virtual G4double GetCubicVolume()
G4ThreeVector GetPointOnSurface() const
void SetNormal(const G4ThreeVector &newNormal)
virtual G4bool PartialClip(const G4VoxelLimits &voxelLimit, const EAxis IgnoreMe)
virtual void AddVertexInOrder(const G4ThreeVector vertex)
virtual void ClearAllVertices()
G4VisExtent GetExtent() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
HepPolyhedron & Transform(const G4Transform3D &t)
virtual void AddSolid(const G4Box &)=0
G4Polyhedron * GetPolyhedron() const
G4EllipticalTube & operator=(const G4EllipticalTube &rhs)
G4double CheckXY(const G4double x, const G4double y, const G4double toler) const
G4int IntersectXY(const G4ThreeVector &p, const G4ThreeVector &v, G4double s[2]) const
G4double GetCubicVolume()
void AddSurface(const G4ClippablePolygon &surface)
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
std::ostream & StreamInfo(std::ostream &os) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4Polyhedron * CreatePolyhedron() const
static G4int GetNumberOfRotationSteps()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4EllipticalTube(const G4String &name, G4double theDx, G4double theDy, G4double theDz)
G4GeometryType GetEntityType() const
virtual G4double GetSurfaceArea()
const G4int kMaxMeshSections
EInside Inside(const G4ThreeVector &p) const