52 #if !defined(G4GEOM_USE_UTRAP)
66 using namespace CLHEP;
94 if ( pDz <= 0 || pDy1 <= 0 || pDx1 <= 0 ||
95 pDx2 <= 0 || pDy2 <= 0 || pDx3 <= 0 || pDx4 <= 0 )
97 std::ostringstream message;
98 message <<
"Invalid length parameters for Solid: " <<
GetName() <<
G4endl
100 << pDx1 <<
", " << pDx2 <<
", " << pDx3 <<
", " << pDx4 <<
G4endl
101 <<
" Y - " << pDy1 <<
", " << pDy2 <<
G4endl
140 && pt[0].
z() == pt[1].
z() && pt[0].
z() == pt[2].
z()
141 && pt[0].
z() == pt[3].
z()
143 && pt[4].
z() == pt[5].
z() && pt[4].
z() == pt[6].
z()
144 && pt[4].
z() == pt[7].
z()
146 && pt[0].y() == pt[1].y() && pt[2].y() == pt[3].y()
147 && pt[4].y() == pt[5].y() && pt[6].y() == pt[7].y()
148 && std::fabs( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) <
kCarTolerance
149 && std::fabs( pt[0].
x() + pt[1].
x() + pt[4].
x() + pt[5].
x() +
152 std::ostringstream message;
153 message <<
"Invalid vertice coordinates for Solid: " <<
GetName();
166 "Face at ~-Y not planar.");
175 std::ostringstream message;
176 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
187 std::ostringstream message;
188 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
198 std::ostringstream message;
199 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
205 fDy1 = ((pt[2]).y()-(pt[1]).y())*0.5;
206 fDx1 = ((pt[1]).
x()-(pt[0]).
x())*0.5;
207 fDx2 = ((pt[3]).
x()-(pt[2]).
x())*0.5;
210 fDy2 = ((pt[6]).y()-(pt[5]).y())*0.5;
211 fDx3 = ((pt[5]).
x()-(pt[4]).
x())*0.5;
212 fDx4 = ((pt[7]).
x()-(pt[6]).
x())*0.5;
231 if ( pZ<=0 || pY<=0 || pX<=0 || pLTX<=0 || pLTX>pX )
233 std::ostringstream message;
234 message <<
"Invalid length parameters for Solid: " <<
GetName();
277 std::ostringstream message;
278 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
288 std::ostringstream message;
289 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
299 std::ostringstream message;
300 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
310 std::ostringstream message;
311 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
329 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 )
331 std::ostringstream message;
332 message <<
"Invalid length parameters for Solid: " <<
GetName();
375 std::ostringstream message;
376 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
386 std::ostringstream message;
387 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
397 std::ostringstream message;
398 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
408 std::ostringstream message;
409 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
428 if ( pDz<=0 || pDy<=0 || pDx<=0 )
430 std::ostringstream message;
431 message <<
"Invalid length parameters for Solid: " <<
GetName();
474 std::ostringstream message;
475 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
485 std::ostringstream message;
486 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
496 std::ostringstream message;
497 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
507 std::ostringstream message;
508 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
521 :
G4CSGSolid (pName), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
522 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
523 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
534 :
G4CSGSolid(a), fDz(1.), fTthetaCphi(0.), fTthetaSphi(0.),
535 fDy1(1.), fDx1(1.), fDx2(1.), fTalpha1(0.),
536 fDy2(1.), fDx3(1.), fDx4(1.), fTalpha2(0.)
555 fTthetaCphi(rhs.fTthetaCphi), fTthetaSphi(rhs.fTthetaSphi),
556 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fTalpha1(rhs.fTalpha1),
557 fDy2(rhs.fDy2), fDx3(rhs.fDx3), fDx4(rhs.fDx4), fTalpha2(rhs.fTalpha2)
559 for (
size_t i=0; i<4; ++i)
576 if (
this == &rhs) {
return *
this; }
588 for (
size_t i=0; i<4; ++i)
616 if ( pDz<=0 || pDy1<=0 || pDx1<=0 || pDx2<=0 || pDy2<=0 || pDx3<=0 || pDx4<=0 )
618 std::ostringstream message;
619 message <<
"Invalid Length Parameters for Solid: " <<
GetName() <<
G4endl
621 << pDx1 <<
", " << pDx2 <<
", " << pDx3 <<
", " << pDx4 <<
G4endl
622 <<
" Y - " << pDy1 <<
", " << pDy2 <<
G4endl
624 G4Exception(
"G4Trap::SetAllParameters()",
"GeomSolids0002",
679 std::ostringstream message;
680 message <<
"Face at ~-Y not planar for Solid: " <<
GetName();
681 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
690 std::ostringstream message;
691 message <<
"Face at ~+Y not planar for Solid: " <<
GetName();
692 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
701 std::ostringstream message;
702 message <<
"Face at ~-X not planar for Solid: " <<
GetName();
703 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
712 std::ostringstream message;
713 message <<
"Face at ~+X not planar for Solid: " <<
GetName();
714 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
744 Vcross = v12.cross(v13);
766 a = +(p4.y() - p2.y())*(p3.z() - p1.z())
767 - (p3.y() - p1.y())*(p4.z() - p2.z());
769 b = -(p4.x() - p2.x())*(p3.z() - p1.z())
770 + (p3.x() - p1.x())*(p4.z() - p2.z());
772 c = +(p4.x() - p2.x())*(p3.y() - p1.y())
773 - (p3.x() - p1.x())*(p4.y() - p2.y());
775 sd = std::sqrt( a*a + b*b + c*c );
785 std::ostringstream message;
786 message <<
"Invalid parameters: norm.mod() <= 0, for Solid: "
788 G4Exception(
"G4Trap::MakePlanes()",
"GeomSolids0002",
793 plane.
d = -( plane.
a*p1.x() + plane.
b*p1.y() + plane.
c*p1.z() );
822 G4double xMin, xMax, yMin, yMax, zMin, zMax;
880 temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMin-pt[0].
z())
881 /(pt[4].
z()-pt[0].z()) ;
882 temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())*(zMax-pt[0].
z())
883 /(pt[4].
z()-pt[0].z()) ;
884 temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMin-pt[2].
z())
885 /(pt[6].
z()-pt[2].z()) ;
886 temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())*(zMax-pt[2].
z())
887 /(pt[6].
z()-pt[2].z()) ;
892 for( i = 0 ; i < 4 ; i++ )
894 if( temp[i] > yMax ) yMax = temp[i] ;
895 if( temp[i] < yMin ) yMin = temp[i] ;
916 temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
917 *(zMin-pt[0].
z())/(pt[4].
z()-pt[0].z()) ;
918 temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
919 *(zMax-pt[0].
z())/(pt[4].
z()-pt[0].z()) ;
920 temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
921 *(zMin-pt[2].
z())/(pt[6].
z()-pt[2].z()) ;
922 temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
923 *(zMax-pt[2].
z())/(pt[6].
z()-pt[2].z()) ;
924 temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
925 *(zMin-pt[3].
z())/(pt[7].
z()-pt[3].z()) ;
926 temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
927 *(zMax-pt[3].
z())/(pt[7].
z()-pt[3].z()) ;
928 temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
929 *(zMin-pt[1].
z())/(pt[5].
z()-pt[1].z()) ;
930 temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
931 *(zMax-pt[1].
z())/(pt[5].
z()-pt[1].z()) ;
936 for( i = 0 ; i < 8 ; i++ )
938 if( temp[i] > xMax) xMax = temp[i] ;
939 if( temp[i] < xMin) xMin = temp[i] ;
987 G4bool existsAfterClip = false ;
999 for(
G4int nv = 0 ; nv < 8 ; nv++ )
1001 if( (*vertices)[nv].
x() > xMax ) xMax = (*vertices)[nv].x();
1002 if( (*vertices)[nv].y() > yMax ) yMax = (*vertices)[nv].y();
1003 if( (*vertices)[nv].
z() > zMax ) zMax = (*vertices)[nv].z();
1005 if( (*vertices)[nv].
x() < xMin ) xMin = (*vertices)[nv].x();
1006 if( (*vertices)[nv].y() < yMin ) yMin = (*vertices)[nv].y();
1007 if( (*vertices)[nv].
z() < zMin ) zMin = (*vertices)[nv].z();
1091 existsAfterClip=
true;
1099 flag = existsAfterClip ;
1118 for ( i = 0;i < 4;i++ )
1132 for ( i = 0; i < 4; i++ )
1152 G4int i, noSurfaces = 0;
1157 for (i = 0; i < 4; i++)
1166 distz = std::fabs( std::fabs( p.z() ) -
fDz );
1191 if (distmx <= delta)
1201 if (distmy <= delta)
1209 if ( p.z() >= 0.) sumnorm += nZ;
1212 if ( noSurfaces == 0 )
1215 G4Exception(
"G4Trap::SurfaceNormal(p)",
"GeomSolids1002",
1220 else if ( noSurfaces == 1 ) norm = sumnorm;
1221 else norm = sumnorm.unit();
1244 safez=std::fabs(std::fabs(p.z())-
fDz);
1288 smin = (-
fDz-p.z())/v.z();
1295 else if (v.z() < 0 )
1297 max = -
fDz - p.z() ;
1301 smin=(
fDz-p.z())/v.z();
1397 safe=std::fabs(p.z())-
fDz;
1402 if (Dist > safe) safe=Dist;
1679 std::ostringstream message;
1680 G4int oldprc = message.precision(16);
1681 message <<
"Undefined side for valid surface normal to solid."
1683 <<
"Position:" << G4endl << G4endl
1684 <<
"p.x() = " << p.x()/
mm <<
" mm" << G4endl
1685 <<
"p.y() = " << p.y()/
mm <<
" mm" << G4endl
1686 <<
"p.z() = " << p.z()/
mm <<
" mm" << G4endl << G4endl
1687 <<
"Direction:" << G4endl << G4endl
1688 <<
"v.x() = " << v.x() << G4endl
1689 <<
"v.y() = " << v.y() << G4endl
1690 <<
"v.z() = " << v.z() << G4endl << G4endl
1691 <<
"Proposed distance :" << G4endl << G4endl
1692 <<
"snxt = " << snxt/
mm <<
" mm" <<
G4endl;
1693 message.precision(oldprc);
1694 G4Exception(
"G4Trap::DistanceToOut(p,v,..)",
"GeomSolids1002",
1722 G4cout.precision(oldprc) ;
1724 "GeomSolids1002",
JustWarning,
"Point p is outside !?" );
1728 safe=
fDz-std::fabs(p.z());
1736 if (Dist<safe) safe=Dist;
1759 vertices->reserve(8);
1791 "Error in allocation of vertices. Out of memory !");
1811 return new G4Trap(*
this);
1820 G4int oldprc = os.precision(16);
1821 os <<
"-----------------------------------------------------------\n"
1822 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1823 <<
" ===================================================\n"
1824 <<
" Solid type: G4Trap\n"
1825 <<
" Parameters: \n"
1826 <<
" half length Z: " <<
fDz/
mm <<
" mm \n"
1827 <<
" half length Y of face -fDz: " <<
fDy1/
mm <<
" mm \n"
1828 <<
" half length X of side -fDy1, face -fDz: " <<
fDx1/
mm <<
" mm \n"
1829 <<
" half length X of side +fDy1, face -fDz: " <<
fDx2/
mm <<
" mm \n"
1830 <<
" half length Y of face +fDz: " <<
fDy2/
mm <<
" mm \n"
1831 <<
" half length X of side -fDy2, face +fDz: " <<
fDx3/
mm <<
" mm \n"
1832 <<
" half length X of side +fDy2, face +fDz: " <<
fDx4/
mm <<
" mm \n"
1835 <<
" std::tan(alpha), -fDz: " <<
fTalpha1/
degree <<
" degrees \n"
1836 <<
" std::tan(alpha), +fDz: " <<
fTalpha2/
degree <<
" degrees \n"
1837 <<
" trap side plane equations:\n"
1846 <<
"-----------------------------------------------------------\n";
1847 os.precision(oldprc);
1862 G4double lambda1, lambda2, chose, aOne, aTwo;
1871 w.z()*v.x() - w.x()*v.z(),
1872 w.x()*v.y() - w.y()*v.x());
1874 aOne = 0.5*Area.mag();
1877 t.z()*u.x() - t.x()*u.z(),
1878 t.x()*u.y() - t.y()*u.x());
1880 aTwo = 0.5*Area.mag();
1886 if( (chose>=0.) && (chose < aOne) )
1890 return (p2+lambda1*v+lambda2*w);
1898 return (p0+lambda1*t+lambda2*u);
1907 G4double aOne, aTwo, aThree, aFour, aFive, aSix, chose;
1938 if( (chose>=0.) && (chose<aOne) )
1940 else if( (chose>=aOne) && (chose<aOne+aTwo) )
1942 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThree) )
1944 else if( (chose>=aOne+aTwo+aThree) && (chose<aOne+aTwo+aThree+aFour) )
1946 else if( (chose>=aOne+aTwo+aThree+aFour)
1947 && (chose<aOne+aTwo+aThree+aFour+aFive) )
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 const G4double kInfinity
G4double GetMinYExtent() const
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)
G4bool IsYLimited() const
G4Trap & operator=(const G4Trap &rhs)
const G4double w[NPOINTSGL]
G4bool IsXLimited() const
virtual void AddSolid(const G4Box &)=0
static double normal(HepRandomEngine *eptr)
G4double GetMaxXExtent() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
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
std::vector< G4ThreeVector > G4ThreeVectorList
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
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
G4double GetMinXExtent() const
std::ostream & StreamInfo(std::ostream &os) const
G4double GetMaxZExtent() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double x[NPOINTSGL]
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
static const double degree
G4ThreeVector GetPointOnSurface() const
G4double GetMaxYExtent() const
const G4double kCoplanar_Tolerance
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)
EInside Inside(const G4ThreeVector &p) const
G4CSGSolid & operator=(const G4CSGSolid &rhs)
G4bool IsZLimited() const
void DescribeYourselfTo(G4VGraphicsScene &scene) const