55 #if !defined(G4GEOM_USE_UTRAP)
70 using namespace CLHEP;
98 if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
99 pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
101 std::ostringstream message;
102 message <<
"Invalid length parameters for Solid: " <<
GetName() <<
G4endl
104 << pDx1 <<
", " << pDx2 <<
", " << pDx3 <<
", " << pDx4 <<
G4endl
105 <<
" Y - " << pDy1 <<
", " << pDy2 <<
G4endl
144 && pt[0].z() == pt[1].z() && pt[0].z() == pt[2].z()
145 && pt[0].z() == pt[3].z()
147 && pt[4].z() == pt[5].z() && pt[4].z() == pt[6].z()
148 && pt[4].z() == pt[7].z()
150 && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
151 && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
152 && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) <
kCarTolerance
153 && std::fabs( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() +
154 pt[2].x() + pt[3].x() + pt[6].x() + pt[7].x() ) <
kCarTolerance ) )
156 std::ostringstream message;
157 message <<
"Invalid vertice coordinates for Solid: " <<
GetName();
170 "Face at ~-Y not planar.");
179 std::ostringstream message;
180 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
191 std::ostringstream message;
192 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
202 std::ostringstream message;
203 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
209 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
210 fDx1 = ((pt[1]).x()-(pt[0]).x())*0.5;
211 fDx2 = ((pt[3]).x()-(pt[2]).x())*0.5;
212 fTalpha1 = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/
fDy1;
214 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
215 fDx3 = ((pt[5]).x()-(pt[4]).x())*0.5;
216 fDx4 = ((pt[7]).x()-(pt[6]).x())*0.5;
217 fTalpha2 = ((pt[6]).x()+(pt[7]).x()-(pt[5]).x()-(pt[4]).x())*0.25/
fDy2;
235 if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
237 std::ostringstream message;
238 message <<
"Invalid length parameters for Solid: " <<
GetName();
281 std::ostringstream message;
282 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
292 std::ostringstream message;
293 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
303 std::ostringstream message;
304 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
314 std::ostringstream message;
315 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
333 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
335 std::ostringstream message;
336 message <<
"Invalid length parameters for Solid: " <<
GetName();
379 std::ostringstream message;
380 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
390 std::ostringstream message;
391 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
401 std::ostringstream message;
402 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
412 std::ostringstream message;
413 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
432 if ( pDz<=0 || pDy<=0 || pDx<=0 )
434 std::ostringstream message;
435 message <<
"Invalid length parameters for Solid: " <<
GetName();
478 std::ostringstream message;
479 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
489 std::ostringstream message;
490 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
500 std::ostringstream message;
501 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
511 std::ostringstream message;
512 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
525 :
G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
526 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
527 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
538 :
G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
539 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
540 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
559 fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
560 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
561 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
563 for (
size_t i=0; i<4; ++i)
580 if (
this == &rhs) {
return *
this; }
592 for (
size_t i=0; i<4; ++i)
620 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
622 std::ostringstream message;
623 message <<
"Invalid Length Parameters for Solid: " <<
GetName() <<
G4endl
625 << pDx1 <<
", " << pDx2 <<
", " << pDx3 <<
", " << pDx4 <<
G4endl
626 <<
" Y - " << pDy1 <<
", " << pDy2 <<
G4endl
628 G4Exception(
"G4Trap::SetAllParameters()",
"GeomSolids0002",
683 std::ostringstream message;
684 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
685 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
694 std::ostringstream message;
695 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
696 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
705 std::ostringstream message;
706 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
707 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
716 std::ostringstream message;
717 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
718 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
748 Vcross = v12.cross(v13);
770 a = +(p4.y() - p2.y())*(p3.z() - p1.z())
771 - (p3.y() - p1.y())*(p4.z() - p2.z());
773 b = -(p4.x() - p2.x())*(p3.z() - p1.z())
774 + (p3.x() - p1.x())*(p4.z() - p2.z());
776 c = +(p4.x() - p2.x())*(p3.y() - p1.y())
777 - (p3.x() - p1.x())*(p4.y() - p2.y());
779 sd = std::sqrt( a*a + b*b + c*c );
789 std::ostringstream message;
790 message <<
"Invalid parameters: norm.mod() <= 0, for Solid: "
792 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
797 plane.
d = -( plane.
a*p1.x() + plane.
b*p1.y() + plane.
c*p1.z() );
836 std::min(-x0-x1-dx1,-x0+x1-dx2),x0-x2-dx3),x0+x2-dx4);
840 std::max(-x0-x1+dx1,-x0+x1+dx2),x0-x2+dx3),x0+x2+dx4);
846 pMin.set(xmin,ymin,-dz);
847 pMax.set(xmax,ymax, dz);
851 if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
853 std::ostringstream message;
854 message <<
"Bad bounding box (min >= max) for solid: "
856 <<
"\npMin = " << pMin
857 <<
"\npMax = " << pMax;
880 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
884 return exist = (pMin < pMax) ?
true :
false;
903 baseA[0].set(-x0-x1-dx1,-y0-dy1,-dz);
904 baseA[1].set(-x0-x1+dx1,-y0-dy1,-dz);
905 baseA[2].set(-x0+x1+dx2,-y0+dy1,-dz);
906 baseA[3].set(-x0+x1-dx2,-y0+dy1,-dz);
908 baseB[0].set( x0-x2-dx3, y0-dy2, dz);
909 baseB[1].set( x0-x2+dx3, y0-dy2, dz);
910 baseB[2].set( x0+x2+dx4, y0+dy2, dz);
911 baseB[3].set( x0+x2-dx4, y0+dy2, dz);
913 std::vector<const G4ThreeVectorList *> polygons(2);
914 polygons[0] = &baseA;
915 polygons[1] = &baseB;
935 for ( i = 0;i < 4;i++ )
949 for ( i = 0; i < 4; i++ )
969 G4int i, noSurfaces = 0;
974 for (i = 0; i < 4; i++)
983 distz = std::fabs( std::fabs( p.z() ) -
fDz );
1008 if (distmx <= delta)
1018 if (distmy <= delta)
1026 if ( p.z() >= 0.) sumnorm += nZ;
1029 if ( noSurfaces == 0 )
1032 G4Exception(
"G4Trap::SurfaceNormal(p)",
"GeomSolids1002",
1037 else if ( noSurfaces == 1 ) norm = sumnorm;
1038 else norm = sumnorm.unit();
1061 safez=std::fabs(std::fabs(p.z())-
fDz);
1105 smin = (-
fDz-p.z())/v.z();
1112 else if (v.z() < 0 )
1114 max = -
fDz - p.z() ;
1118 smin=(
fDz-p.z())/v.z();
1214 safe=std::fabs(p.z())-
fDz;
1219 if (Dist > safe) safe=Dist;
1496 std::ostringstream message;
1497 G4int oldprc = message.precision(16);
1498 message <<
"Undefined side for valid surface normal to solid."
1500 <<
"Position:" << G4endl << G4endl
1501 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1502 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1503 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1504 <<
"Direction:" << G4endl << G4endl
1505 <<
"v.x() = " << v.x() << G4endl
1506 <<
"v.y() = " << v.y() << G4endl
1507 <<
"v.z() = " << v.z() << G4endl << G4endl
1508 <<
"Proposed distance :" << G4endl << G4endl
1509 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl;
1510 message.precision(oldprc);
1511 G4Exception(
"G4Trap::DistanceToOut(p,v,..)",
"GeomSolids1002",
1539 G4cout.precision(oldprc) ;
1541 "GeomSolids1002",
JustWarning,
"Point p is outside !?" );
1545 safe=
fDz-std::fabs(p.z());
1553 if (Dist<safe) safe=Dist;
1575 return new G4Trap(*
this);
1584 G4int oldprc = os.precision(16);
1585 os <<
"-----------------------------------------------------------\n"
1586 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1587 <<
" ===================================================\n"
1588 <<
" Solid type: G4Trap\n"
1589 <<
" Parameters: \n"
1590 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1591 <<
" half length Y of face -fDz: " <<
fDy1/
mm <<
" mm \n"
1592 <<
" half length X of side -fDy1, face -fDz: " <<
fDx1/
mm <<
" mm \n"
1593 <<
" half length X of side +fDy1, face -fDz: " <<
fDx2/
mm <<
" mm \n"
1594 <<
" half length Y of face +fDz: " <<
fDy2/
mm <<
" mm \n"
1595 <<
" half length X of side -fDy2, face +fDz: " <<
fDx3/
mm <<
" mm \n"
1596 <<
" half length X of side +fDy2, face +fDz: " <<
fDx4/
mm <<
" mm \n"
1599 <<
" std::tan(alpha), -fDz: " <<
fTalpha1/
degree <<
" degrees \n"
1600 <<
" std::tan(alpha), +fDz: " <<
fTalpha2/
degree <<
" degrees \n"
1601 <<
" trap side plane equations:\n"
1610 <<
"-----------------------------------------------------------\n";
1611 os.precision(oldprc);
1626 G4double lambda1, lambda2, chose, aOne, aTwo;
1635 w.z()*v.x() - w.x()*v.z(),
1636 w.x()*v.y() - w.y()*v.x());
1638 aOne = 0.5*Area.mag();
1641 t.z()*u.x() - t.x()*u.z(),
1642 t.x()*u.y() - t.y()*u.x());
1644 aTwo = 0.5*Area.mag();
1650 if( (chose>=0.) && (chose < aOne) )
1654 return (p2+lambda1*v+lambda2*w);
1662 return (p0+lambda1*t+lambda2*u);
1671 G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1702 if( (chose>=0.) && (chose<aOne) )
1704 else if( (chose>=aOne) && (chose<aOne+aTwo) )
1706 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
1708 else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
1710 else if( (chose>=aOne+aTwo+aThree+aFour)
1711 && (chose<aOne+aTwo+aThree+aFour+aFive) )
G4double GetXHalfLength4() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4double &area) const
G4GeometryType GetEntityType() const
static constexpr double mm
G4double GetYHalfLength2() const
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4bool fRebuildPolyhedron
G4Trap(const G4String &pName, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
std::vector< ExP01TrackerHit * > a
G4double GetZHalfLength() const
G4Trap & operator=(const G4Trap &rhs)
virtual void AddSolid(const G4Box &)=0
G4double GetXHalfLength2() const
G4double GetTanAlpha2() const
static double normal(HepRandomEngine *eptr)
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double GetXHalfLength1() const
G4GLOB_DLL std::ostream G4cout
G4double GetXHalfLength3() const
static constexpr double degree
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
void SetAllParameters(G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
std::vector< G4ThreeVector > G4ThreeVectorList
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
std::ostream & StreamInfo(std::ostream &os) const
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
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4ThreeVector GetPointOnSurface() const
const G4double kCoplanar_Tolerance
G4double GetYHalfLength1() const
G4Polyhedron * CreatePolyhedron() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool MakePlane(const G4ThreeVector &p1, const G4ThreeVector &p2, const G4ThreeVector &p3, const G4ThreeVector &p4, TrapSidePlane &plane)
G4double GetTanAlpha1() const
EInside Inside(const G4ThreeVector &p) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
void DescribeYourselfTo(G4VGraphicsScene &scene) const