77 fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
78 fSide90(0), fSide180(0), fSide270(0),
79 fSurfaceArea(0.), fpPolyhedron(0)
99 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
100 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
101 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
102 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
106 if ( fDx1 != fDx2 && fDx3 != fDx4 )
108 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
111 std::ostringstream message;
112 message <<
"Not planar surface in untwisted Trapezoid: "
114 <<
"fDy2 is " << fDy2 <<
" but should be "
116 G4Exception(
"G4VTwistedFaceted::G4VTwistedFaceted()",
"GeomSolids0002",
122 if ( fDx1 == fDx2 && fDx3 == fDx4 )
129 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
131 std::ostringstream message;
132 message <<
"Not planar surface in untwisted Trapezoid: "
134 <<
"One endcap is rectangular, the other is a trapezoid." <<
G4endl
135 <<
"For planarity reasons they have to be rectangles or trapezoids "
137 G4Exception(
"G4VTwistedFaceted::G4VTwistedFaceted()",
"GeomSolids0002",
143 fPhiTwist = PhiTwist ;
148 fTAlph = std::tan(fAlph) ;
155 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
159 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
168 && ( std::fabs(fPhiTwist) > 2*kAngTolerance )
169 && ( std::fabs(fPhiTwist) <
pi/2 )
170 && ( std::fabs(fAlph) <
pi/2 )
171 && ( fTheta < pi/2 && fTheta >= 0 ) )
174 std::ostringstream message;
175 message <<
"Invalid dimensions. Too small, or twist angle too big: "
177 <<
"fDx 1-4 = " << fDx1/
cm <<
", " << fDx2/
cm <<
", "
178 << fDx3/
cm <<
", " << fDx4/
cm <<
" cm" <<
G4endl
179 <<
"fDy 1-2 = " << fDy1/
cm <<
", " << fDy2/
cm <<
", "
181 <<
"fDz = " << fDz/
cm <<
" cm" <<
G4endl
182 <<
" twistangle " << fPhiTwist/
deg <<
" deg" <<
G4endl
183 <<
" phi,theta = " << fPhi/
deg <<
", " << fTheta/
deg <<
" deg";
188 fCubicVolume = 2 * fDz * ( ( fDx1 + fDx2 ) * fDy1 + ( fDx3 + fDx4 ) * fDy2 );
196 :
G4VSolid(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
197 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fDx(0.), fDy(0.), fAlph(0.),
198 fTAlph(0.), fdeltaX(0.), fdeltaY(0.), fPhiTwist(0.),
199 fLowerEndcap(0), fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0),
200 fSide270(0), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
209 if (fLowerEndcap)
delete fLowerEndcap ;
210 if (fUpperEndcap)
delete fUpperEndcap ;
212 if (fSide0)
delete fSide0 ;
213 if (fSide90)
delete fSide90 ;
214 if (fSide180)
delete fSide180 ;
215 if (fSide270)
delete fSide270 ;
216 if (fpPolyhedron)
delete fpPolyhedron;
223 :
G4VSolid(rhs), fTheta(rhs.fTheta), fPhi(rhs.fPhi),
224 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fDy2(rhs.fDy2),
225 fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fDz), fDx(rhs.fDx), fDy(rhs.fDy),
226 fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdeltaX(rhs.fdeltaX),
227 fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTwist), fLowerEndcap(0),
228 fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0), fSide270(0),
229 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
231 fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
232 fLastDistanceToIn(rhs.fLastDistanceToIn),
233 fLastDistanceToOut(rhs.fLastDistanceToOut),
234 fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
235 fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
248 if (
this == &rhs) {
return *
this; }
256 fTheta = rhs.fTheta; fPhi = rhs.fPhi;
257 fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.fDx2; fDy2= rhs.fDy2;
258 fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fDz; fDx= rhs.fDx; fDy= rhs.fDy;
259 fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdeltaX= rhs.fdeltaX;
260 fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTwist; fLowerEndcap= 0;
261 fUpperEndcap= 0; fSide0= 0; fSide90= 0; fSide180= 0; fSide270= 0;
262 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
264 fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
265 fLastDistanceToIn= rhs.fLastDistanceToIn;
266 fLastDistanceToOut= rhs.fLastDistanceToOut;
267 fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
268 fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
282 G4Exception(
"G4VTwistedFaceted::ComputeDimensions()",
284 "G4VTwistedFaceted does not support Parameterisation.");
297 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
310 xMin = xoffset - maxRad ;
311 xMax = xoffset + maxRad ;
330 yMin = yoffset - maxRad ;
331 yMax = yoffset + maxRad ;
350 zMin = zoffset - fDz ;
351 zMax = zoffset + fDz ;
393 G4bool existsAfterClip = false ;
406 if (pVoxelLimit.
IsLimited(pAxis) ==
false)
408 if ( pMin != kInfinity || pMax != -kInfinity )
410 existsAfterClip = true ;
425 if ( pMin != kInfinity || pMax != -kInfinity )
427 existsAfterClip = true ;
463 existsAfterClip = true ;
469 return existsAfterClip;
483 vertices->reserve(8);
485 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
508 G4Exception(
"G4VTwistedFaceted::CreateRotatedVertices()",
510 "Error in allocation of vertices. Out of memory !");
523 if (fLastInside.p == p) {
524 return fLastInside.inside;
527 tmpin =
const_cast<EInside*
>(&(fLastInside.inside));
528 tmpp->
set(p.
x(), p.
y(), p.
z());
533 G4double phi = p.
z()/(2*fDz) * fPhiTwist ;
537 G4double px = p.
x() + fdeltaX * ( -phi/fPhiTwist) ;
538 G4double py = p.
y() + fdeltaY * ( -phi/fPhiTwist) ;
541 G4double posx = px * cphi - py * sphi ;
542 G4double posy = px * sphi + py * cphi ;
565 G4cout <<
"phi,theta = " << fPhi <<
" , " << fTheta <<
G4endl ;
608 return fLastInside.inside;
624 if (fLastNormal.p == p)
626 return fLastNormal.vec;
632 tmpp->
set(p.
x(), p.
y(), p.
z());
638 surfaces[0] = fSide0 ;
639 surfaces[1] = fSide90 ;
640 surfaces[2] = fSide180 ;
641 surfaces[3] = fSide270 ;
642 surfaces[4] = fLowerEndcap;
643 surfaces[5] = fUpperEndcap;
652 if (tmpdistance < distance)
654 distance = tmpdistance;
660 tmpsurface[0] = surfaces[besti];
661 *tmpnormal = tmpsurface[0]->
GetNormal(bestxx,
true);
663 return fLastNormal.vec;
686 if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v)
688 return fLastDistanceToIn.value;
692 tmpp =
const_cast<G4ThreeVector*
>(&(fLastDistanceToInWithV.p));
693 tmpv =
const_cast<G4ThreeVector*
>(&(fLastDistanceToInWithV.vec));
694 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToInWithV.value));
695 tmpp->
set(p.
x(), p.
y(), p.
z());
696 tmpv->
set(v.
x(), v.
y(), v.
z());
717 return fLastDistanceToInWithV.value;
731 surfaces[0] = fSide0;
732 surfaces[1] = fSide90 ;
733 surfaces[2] = fSide180 ;
734 surfaces[3] = fSide270 ;
735 surfaces[4] = fLowerEndcap;
736 surfaces[5] = fUpperEndcap;
741 for (i=0; i < 6 ; i++)
748 G4cout <<
"Solid DistanceToIn : distance = " << tmpdistance <<
G4endl ;
751 if (tmpdistance < distance)
753 distance = tmpdistance;
764 return fLastDistanceToInWithV.value;
781 if (fLastDistanceToIn.p == p)
783 return fLastDistanceToIn.value;
788 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToIn.value));
789 tmpp->
set(p.
x(), p.
y(), p.
z());
807 return fLastDistanceToIn.value;
820 surfaces[0] = fSide0;
821 surfaces[1] = fSide90 ;
822 surfaces[2] = fSide180 ;
823 surfaces[3] = fSide270 ;
824 surfaces[4] = fLowerEndcap;
825 surfaces[5] = fUpperEndcap;
833 if (tmpdistance < distance)
835 distance = tmpdistance;
840 return fLastDistanceToIn.value;
845 G4Exception(
"G4VTwistedFaceted::DistanceToIn(p)",
"GeomSolids0003",
877 if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v )
879 return fLastDistanceToOutWithV.value;
883 tmpp =
const_cast<G4ThreeVector*
>(&(fLastDistanceToOutWithV.p));
884 tmpv =
const_cast<G4ThreeVector*
>(&(fLastDistanceToOutWithV.vec));
885 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToOutWithV.value));
886 tmpp->
set(p.
x(), p.
y(), p.
z());
887 tmpv->
set(v.
x(), v.
y(), v.
z());
910 *norm = (blockedsurface->
GetNormal(p,
true));
915 return fLastDistanceToOutWithV.value;
927 surfaces[0] = fSide0;
928 surfaces[1] = fSide90 ;
929 surfaces[2] = fSide180 ;
930 surfaces[3] = fSide270 ;
931 surfaces[4] = fLowerEndcap;
932 surfaces[5] = fUpperEndcap;
938 for (i=0; i< 6 ; i++) {
940 if (tmpdistance < distance)
942 distance = tmpdistance;
952 *norm = (surfaces[besti]->
GetNormal(p,
true));
959 return fLastDistanceToOutWithV.value;
976 if (fLastDistanceToOut.p == p)
978 return fLastDistanceToOut.value;
983 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToOut.value));
984 tmpp->
set(p.
x(), p.
y(), p.
z());
1006 G4cout.precision(oldprc) ;
1007 G4Exception(
"G4VTwistedFaceted::DistanceToOut(p)",
"GeomSolids1002",
1015 retval = fLastDistanceToOut.value;
1029 surfaces[0] = fSide0;
1030 surfaces[1] = fSide90 ;
1031 surfaces[2] = fSide180 ;
1032 surfaces[3] = fSide270 ;
1033 surfaces[4] = fLowerEndcap;
1034 surfaces[5] = fUpperEndcap;
1039 for (i=0; i< 6; i++)
1042 if (tmpdistance < distance)
1044 distance = tmpdistance;
1048 *tmpdist = distance;
1050 retval = fLastDistanceToOut.value;
1056 G4Exception(
"G4VTwistedFaceted::DistanceToOut(p)",
"GeomSolids0003",
1074 G4int oldprc = os.precision(16);
1075 os <<
"-----------------------------------------------------------\n"
1076 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1077 <<
" ===================================================\n"
1078 <<
" Solid type: G4VTwistedFaceted\n"
1079 <<
" Parameters: \n"
1080 <<
" polar angle theta = " << fTheta/
degree <<
" deg" <<
G4endl
1081 <<
" azimuthal angle phi = " << fPhi/
degree <<
" deg" <<
G4endl
1082 <<
" tilt angle alpha = " << fAlph/
degree <<
" deg" <<
G4endl
1083 <<
" TWIST angle = " << fPhiTwist/
degree <<
" deg" <<
G4endl
1084 <<
" Half length along y (lower endcap) = " << fDy1/
cm <<
" cm"
1086 <<
" Half length along x (lower endcap, bottom) = " << fDx1/
cm <<
" cm"
1088 <<
" Half length along x (lower endcap, top) = " << fDx2/
cm <<
" cm"
1090 <<
" Half length along y (upper endcap) = " << fDy2/
cm <<
" cm"
1092 <<
" Half length along x (upper endcap, bottom) = " << fDx3/
cm <<
" cm"
1094 <<
" Half length along x (upper endcap, top) = " << fDx4/
cm <<
" cm"
1096 <<
"-----------------------------------------------------------\n";
1097 os.precision(oldprc);
1116 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
1127 void G4VTwistedFaceted::CreateSurfaces()
1132 if ( fDx1 == fDx2 && fDx3 == fDx4 )
1134 fSide0 =
new G4TwistBoxSide(
"0deg", fPhiTwist, fDz, fTheta, fPhi,
1135 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*
deg);
1137 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*
deg);
1142 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*
deg);
1144 fPhi+
pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*
deg);
1150 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*
deg);
1152 fPhi+
pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*
deg);
1157 fDz, fAlph, fPhi, fTheta, 1 );
1159 fDz, fAlph, fPhi, fTheta, -1 );
1163 fSide0->
SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
1164 fSide90->
SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
1165 fSide180->
SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
1166 fSide270->
SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
1167 fUpperEndcap->
SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1168 fLowerEndcap->
SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1178 return G4String(
"G4VTwistedFaceted");
1187 if (!fpPolyhedron ||
1191 delete fpPolyhedron;
1195 return fpPolyhedron;
1209 if ( z == fDz ) z -= 0.1*fDz ;
1210 if ( z == -fDz ) z += 0.1*fDz ;
1212 G4double phi = z/(2*fDz)*fPhiTwist ;
1214 return G4ThreeVector(fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z ) ;
1261 else if( (chose >= a1) && (chose < a1 + a2 ) )
1272 else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1282 else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1292 else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1325 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1326 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1330 typedef G4int G4int4[4];
1331 G4double3* xyz =
new G4double3[nnodes];
1332 G4int4* faces =
new G4int4[nfaces] ;
1334 fLowerEndcap->
GetFacets(k,k,xyz,faces,0) ;
1335 fUpperEndcap->
GetFacets(k,k,xyz,faces,1) ;
void set(double x, double y, double z)
void ClipCrossSection(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
ThreeVector shoot(const G4int Ap, const G4int Af)
virtual void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *)
G4ThreeVectorList * CreateRotatedVertices(const G4AffineTransform &pTransform) const
G4double GetMinYExtent() const
G4VTwistedFaceted(const G4String &pname, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph)
CLHEP::Hep3Vector G4ThreeVector
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
G4bool IsYLimited() const
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
virtual G4GeometryType GetEntityType() const
G4ThreeVector GetPointOnSurface() const
G4ThreeVector GetPointInSolid(G4double z) const
virtual G4Polyhedron * CreatePolyhedron() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool IsXLimited() const
G4VTwistedFaceted & operator=(const G4VTwistedFaceted &rhs)
virtual void AddSolid(const G4Box &)=0
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4bool IsValidNorm() const
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
G4double Xcoef(G4double u, G4double phi, G4double ftg) const
G4double GetMaxXExtent() const
virtual G4double GetSurfaceArea()=0
G4double GetMinZExtent() const
G4GLOB_DLL std::ostream G4cout
virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
void SetNeighbours(G4VTwistSurface *axis0min, G4VTwistSurface *axis1min, G4VTwistSurface *axis0max, G4VTwistSurface *axis1max)
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)=0
virtual ~G4VTwistedFaceted()
G4double GetValueB(G4double phi) const
std::vector< G4ThreeVector > G4ThreeVectorList
virtual G4double GetBoundaryMax(G4double)=0
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=false, G4bool *validnorm=0, G4ThreeVector *n=0) const
virtual G4double GetBoundaryMin(G4double)=0
G4double GetMinXExtent() const
G4double GetMaxZExtent() const
virtual G4Polyhedron * GetPolyhedron() const
static G4int GetNumberOfRotationSteps()
virtual G4VisExtent GetExtent() const
virtual EInside Inside(const G4ThreeVector &p) const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4double GetMaxYExtent() const
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
virtual std::ostream & StreamInfo(std::ostream &os) const
G4double GetMaxExtent(const EAxis pAxis) const
G4bool IsZLimited() const
G4double GetAngularTolerance() const
void ClipBetweenSections(G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
G4double GetMinExtent(const EAxis pAxis) const
static G4GeometryTolerance * GetInstance()