63 using namespace CLHEP;
72 :
G4VSolid( name ), fCubicVolume(0.), fSurfaceArea(0.),
73 fRebuildPolyhedron(false), fpPolyhedron(0)
88 :
G4VSolid(a), dx(0.), dy(0.), dz(0.), halfTol(0.),
89 fCubicVolume(0.), fSurfaceArea(0.),
90 fRebuildPolyhedron(false), fpPolyhedron(0)
108 :
G4VSolid(rhs), dx(rhs.dx), dy(rhs.dy), dz(rhs.dz), halfTol(rhs.halfTol),
109 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
110 fRebuildPolyhedron(false), fpPolyhedron(0)
122 if (
this == &rhs) {
return *
this; }
168 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
172 return exist = (pMin < pMax) ?
true :
false;
177 const G4int NSTEPS = 48;
180 G4double sinHalf = std::sin(0.5*ang);
181 G4double cosHalf = std::cos(0.5*ang);
182 G4double sinStep = 2.*sinHalf*cosHalf;
183 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
190 for (
G4int k=0; k<NSTEPS; ++k)
192 baseA[k].set(sx*cosCur,sy*sinCur,-
dz);
193 baseB[k].set(sx*cosCur,sy*sinCur,
dz);
196 sinCur = sinCur*cosStep + cosCur*sinStep;
197 cosCur = cosCur*cosStep - sinTmp*sinStep;
200 std::vector<const G4ThreeVectorList *> polygons(2);
201 polygons[0] = &baseA;
202 polygons[1] = &baseB;
261 if ( (distZ <
halfTol ) && ( distR1 <= 1 ) )
266 if( (distR1 <= 1 ) && ( distR2 >= 1 ) )
272 if ( noSurfaces == 0 )
275 G4Exception(
"G4EllipticalTube::SurfaceNormal(p)",
"GeomSolids1002",
280 else if ( noSurfaces == 1 ) { norm = sumnorm; }
281 else { norm = sumnorm.
unit(); }
304 if (distZ*distZ < distR2)
368 yi = p.
y() + q*v.
y();
378 return (sigz < -
halfTol) ? q : 0;
380 else if (xi*
dy*
dy*v.
x() + yi*
dx*
dx*v.
y() >= 0)
407 yi = p.
y() + q*v.
y();
411 return (sigz > -
halfTol) ? q : 0;
413 else if (xi*
dy*
dy*v.
x() + yi*
dx*
dx*v.
y() >= 0)
437 if (p.
x()*
dy*
dy*v.
x() + p.
y()*
dx*
dx*v.
y() < 0)
return 0;
515 G4double tnorm = std::sqrt( tx*tx + ty*ty );
520 G4double distR = ( (p.
x()-xe)*ty - (p.
y()-ye)*tx )/tnorm;
530 return std::sqrt( (p.
z()+
dz)*(p.
z()+
dz) + distR*distR );
532 return std::sqrt( (p.
z()-
dz)*(p.
z()-
dz) + distR*distR );
555 if (calcNorm) { *validNorm =
true; }
569 sBest = -(p.
z()+
dz)/v.
z();
576 if (calcNorm) { *norm = normHere; }
601 if (calcNorm) { *norm = normHere; }
608 if (q < sBest) { sBest = q; nBest = normHere; }
623 G4int oldprc = message.precision(16) ;
624 message <<
"Point p is outside !?" <<
G4endl
626 <<
" p.x() = " << p.
x()/
mm <<
" mm" <<
G4endl
627 <<
" p.y() = " << p.
y()/
mm <<
" mm" <<
G4endl
628 <<
" p.z() = " << p.
z()/
mm <<
" mm" <<
G4endl
630 <<
" v.x() = " << v.
x() <<
G4endl
631 <<
" v.y() = " << v.
y() <<
G4endl
632 <<
" v.z() = " << v.
z() <<
G4endl
633 <<
"Proposed distance :" <<
G4endl
634 <<
" snxt = " << sBest/
mm <<
" mm";
635 message.precision(oldprc) ;
636 G4Exception(
"G4EllipticalTube::DistanceToOut(p,v,...)",
639 if (calcNorm) { *norm = nBest; }
642 else if (q[n-1] > sBest)
644 if (calcNorm) { *norm = nBest; }
666 if (p.
x()*
dy*
dy*v.
x() + p.
y()*
dx*
dx*v.
y() > 0)
return 0;
701 if (radical < +
DBL_MIN)
return 0;
704 if (p.
x() < 0) xi = -xi;
709 radical = 1.0 - p.
x()*p.
x()/
dx/
dx;
710 if (radical < +
DBL_MIN)
return 0;
713 if (p.
y() < 0) yi = -yi;
722 G4double normi = std::sqrt( xdi*xdi + ydi*ydi );
727 G4double q = 0.5*(xdi*(p.
y()-yi) - ydi*(p.
x()-xi));
728 if (xi*yi < 0) q = -q;
730 if (q < sBest) sBest = q;
735 return sBest <
halfTol ? 0 : sBest;
779 if (radical < -
DBL_MIN)
return 0;
790 radical = std::sqrt(radical);
792 G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
795 if (sa < sb) { ss[0] = sa; ss[1] = sb; }
else { ss[0] = sb; ss[1] = sa; }
805 return G4String(
"G4EllipticalTube");
842 G4int oldprc = os.precision(16);
843 os <<
"-----------------------------------------------------------\n"
844 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
845 <<
" ===================================================\n"
846 <<
" Solid type: G4EllipticalTube\n"
848 <<
" length Z: " <<
dz/
mm <<
" mm \n"
849 <<
" surface equation in X and Y: \n"
850 <<
" (X / " <<
dx <<
")^2 + (Y / " <<
dy <<
")^2 = 1 \n"
851 <<
"-----------------------------------------------------------\n";
852 os.precision(oldprc);
866 G4double xRand, yRand, zRand, phi, cosphi, sinphi, zArea, cArea,
p, chose;
869 cosphi = std::cos(phi);
870 sinphi = std::sin(phi);
890 if( (chose>=0) && (chose < cArea) )
894 else if( (chose >= cArea) && (chose < cArea + zArea) )
897 yRand = std::sqrt(1.-
sqr(xRand/
dx));
904 yRand = std::sqrt(1.-
sqr(xRand/
dx));
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void set(double x, double y, double z)
virtual ~G4EllipticalTube()
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
virtual void AddSolid(const G4Box &)=0
static const G4double kInfinity
static constexpr double mm
G4EllipticalTube & operator=(const G4EllipticalTube &rhs)
EInside Inside(const G4ThreeVector &p) const
void message(RunManager *runmanager)
G4bool fRebuildPolyhedron
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4Polyhedron * fpPolyhedron
G4VSolid & operator=(const G4VSolid &rhs)
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4ThreeVector GetPointOnSurface() const
std::vector< G4ThreeVector > G4ThreeVectorList
#define G4MUTEX_INITIALIZER
G4double GetCubicVolume()
G4VisExtent GetExtent() const
HepPolyhedron & Transform(const G4Transform3D &t)
G4EllipticalTube(const G4String &name, G4double theDx, G4double theDy, G4double theDz)
static constexpr double twopi
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4Polyhedron * CreatePolyhedron() const
G4double CheckXY(const G4double x, const G4double y, const G4double toler) const
ThreeVector shoot(const G4int Ap, const G4int Af)
static G4int GetNumberOfRotationSteps()
G4double GetSurfaceArea()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
static constexpr double pi
G4int IntersectXY(const G4ThreeVector &p, const G4ThreeVector &v, G4double s[2]) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4GeometryType GetEntityType() const
std::ostream & StreamInfo(std::ostream &os) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4Polyhedron * GetPolyhedron() const