84 :
G4VSolid(pname), fRebuildPolyhedron(false), fpPolyhedron(0),
85 fLowerEndcap(0), fUpperEndcap(0), fSide0(0),
86 fSide90(0), fSide180(0), fSide270(0),
107 fDxDown = ( fDx1 > fDx2 ? fDx1 : fDx2 ) ;
108 fDxUp = ( fDx3 > fDx4 ? fDx3 : fDx4 ) ;
109 fDx = ( fDxUp > fDxDown ? fDxUp : fDxDown ) ;
110 fDy = ( fDy1 > fDy2 ? fDy1 : fDy2 ) ;
114 if ( fDx1 != fDx2 && fDx3 != fDx4 )
116 pDytmp = fDy1 * ( fDx3 - fDx4 ) / ( fDx1 - fDx2 ) ;
119 std::ostringstream message;
120 message <<
"Not planar surface in untwisted Trapezoid: "
122 <<
"fDy2 is " << fDy2 <<
" but should be "
124 G4Exception(
"G4VTwistedFaceted::G4VTwistedFaceted()",
"GeomSolids0002",
130 if ( fDx1 == fDx2 && fDx3 == fDx4 )
137 if ( ( fDx1 == fDx2 && fDx3 != fDx4 ) || ( fDx1 != fDx2 && fDx3 == fDx4 ) )
139 std::ostringstream message;
140 message <<
"Not planar surface in untwisted Trapezoid: "
142 <<
"One endcap is rectangular, the other is a trapezoid." <<
G4endl
143 <<
"For planarity reasons they have to be rectangles or trapezoids "
145 G4Exception(
"G4VTwistedFaceted::G4VTwistedFaceted()",
"GeomSolids0002",
151 fPhiTwist = PhiTwist ;
156 fTAlph = std::tan(fAlph) ;
163 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
167 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
176 && ( std::fabs(fPhiTwist) > 2*kAngTolerance )
177 && ( std::fabs(fPhiTwist) <
pi/2 )
178 && ( std::fabs(fAlph) <
pi/2 )
179 && ( fTheta < pi/2 && fTheta >= 0 ) )
182 std::ostringstream message;
183 message <<
"Invalid dimensions. Too small, or twist angle too big: "
185 <<
"fDx 1-4 = " << fDx1/
cm <<
", " << fDx2/
cm <<
", "
186 << fDx3/
cm <<
", " << fDx4/
cm <<
" cm" <<
G4endl
187 <<
"fDy 1-2 = " << fDy1/
cm <<
", " << fDy2/
cm <<
", "
189 <<
"fDz = " << fDz/
cm <<
" cm" <<
G4endl
190 <<
" twistangle " << fPhiTwist/
deg <<
" deg" <<
G4endl
191 <<
" phi,theta = " << fPhi/
deg <<
", " << fTheta/
deg <<
" deg";
196 fCubicVolume = 2 * fDz * ( ( fDx1 + fDx2 ) * fDy1 + ( fDx3 + fDx4 ) * fDy2 );
204 :
G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
205 fTheta(0.), fPhi(0.), fDy1(0.),
206 fDx1(0.), fDx2(0.), fDy2(0.), fDx3(0.), fDx4(0.),
207 fDz(0.), fDx(0.), fDy(0.), fAlph(0.),
208 fTAlph(0.), fdeltaX(0.), fdeltaY(0.), fPhiTwist(0.),
209 fLowerEndcap(0), fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0),
210 fSide270(0), fCubicVolume(0.), fSurfaceArea(0.)
220 if (fLowerEndcap) {
delete fLowerEndcap ; }
221 if (fUpperEndcap) {
delete fUpperEndcap ; }
223 if (fSide0) {
delete fSide0 ; }
224 if (fSide90) {
delete fSide90 ; }
225 if (fSide180) {
delete fSide180 ; }
226 if (fSide270) {
delete fSide270 ; }
235 :
G4VSolid(rhs), fRebuildPolyhedron(false), fpPolyhedron(0),
236 fTheta(rhs.fTheta), fPhi(rhs.fPhi),
237 fDy1(rhs.fDy1), fDx1(rhs.fDx1), fDx2(rhs.fDx2), fDy2(rhs.fDy2),
238 fDx3(rhs.fDx3), fDx4(rhs.fDx4), fDz(rhs.fDz), fDx(rhs.fDx), fDy(rhs.fDy),
239 fAlph(rhs.fAlph), fTAlph(rhs.fTAlph), fdeltaX(rhs.fdeltaX),
240 fdeltaY(rhs.fdeltaY), fPhiTwist(rhs.fPhiTwist), fLowerEndcap(0),
241 fUpperEndcap(0), fSide0(0), fSide90(0), fSide180(0), fSide270(0),
242 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
243 fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
244 fLastDistanceToIn(rhs.fLastDistanceToIn),
245 fLastDistanceToOut(rhs.fLastDistanceToOut),
246 fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
247 fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
260 if (
this == &rhs) {
return *
this; }
268 fTheta = rhs.fTheta; fPhi = rhs.fPhi;
269 fDy1= rhs.fDy1; fDx1= rhs.fDx1; fDx2= rhs.fDx2; fDy2= rhs.fDy2;
270 fDx3= rhs.fDx3; fDx4= rhs.fDx4; fDz= rhs.fDz; fDx= rhs.fDx; fDy= rhs.fDy;
271 fAlph= rhs.fAlph; fTAlph= rhs.fTAlph; fdeltaX= rhs.fdeltaX;
272 fdeltaY= rhs.fdeltaY; fPhiTwist= rhs.fPhiTwist; fLowerEndcap= 0;
273 fUpperEndcap= 0; fSide0= 0; fSide90= 0; fSide180= 0; fSide270= 0;
274 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
277 fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
278 fLastDistanceToIn= rhs.fLastDistanceToIn;
279 fLastDistanceToOut= rhs.fLastDistanceToOut;
280 fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
281 fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
296 G4Exception(
"G4VTwistedFaceted::ComputeDimensions()",
298 "G4VTwistedFaceted does not support Parameterisation.");
308 G4double maxRad = std::sqrt(fDx*fDx + fDy*fDy);
309 pMin.
set(-maxRad,-maxRad,-fDz);
310 pMax.
set( maxRad, maxRad, fDz);
343 if (fLastInside.p == p) {
344 return fLastInside.inside;
347 tmpin =
const_cast<EInside*
>(&(fLastInside.inside));
348 tmpp->
set(p.
x(), p.
y(), p.
z());
353 G4double phi = p.
z()/(2*fDz) * fPhiTwist ;
357 G4double px = p.
x() + fdeltaX * ( -phi/fPhiTwist) ;
358 G4double py = p.
y() + fdeltaY * ( -phi/fPhiTwist) ;
361 G4double posx = px * cphi - py * sphi ;
362 G4double posy = px * sphi + py * cphi ;
385 G4cout <<
"phi,theta = " << fPhi <<
" , " << fTheta <<
G4endl ;
428 return fLastInside.inside;
445 if (fLastNormal.p == p)
447 return fLastNormal.vec;
453 tmpp->
set(p.
x(), p.
y(), p.
z());
459 surfaces[0] = fSide0 ;
460 surfaces[1] = fSide90 ;
461 surfaces[2] = fSide180 ;
462 surfaces[3] = fSide270 ;
463 surfaces[4] = fLowerEndcap;
464 surfaces[5] = fUpperEndcap;
473 if (tmpdistance < distance)
475 distance = tmpdistance;
481 tmpsurface[0] = surfaces[besti];
482 *tmpnormal = tmpsurface[0]->
GetNormal(bestxx,
true);
484 return fLastNormal.vec;
508 if (fLastDistanceToInWithV.p == p && fLastDistanceToInWithV.vec == v)
510 return fLastDistanceToIn.value;
514 tmpp =
const_cast<G4ThreeVector*
>(&(fLastDistanceToInWithV.p));
515 tmpv =
const_cast<G4ThreeVector*
>(&(fLastDistanceToInWithV.vec));
516 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToInWithV.value));
517 tmpp->
set(p.
x(), p.
y(), p.
z());
518 tmpv->
set(v.
x(), v.
y(), v.
z());
539 return fLastDistanceToInWithV.value;
553 surfaces[0] = fSide0;
554 surfaces[1] = fSide90 ;
555 surfaces[2] = fSide180 ;
556 surfaces[3] = fSide270 ;
557 surfaces[4] = fLowerEndcap;
558 surfaces[5] = fUpperEndcap;
563 for (i=0; i < 6 ; i++)
570 G4cout <<
"Solid DistanceToIn : distance = " << tmpdistance <<
G4endl ;
573 if (tmpdistance < distance)
575 distance = tmpdistance;
586 return fLastDistanceToInWithV.value;
606 if (fLastDistanceToIn.p == p)
608 return fLastDistanceToIn.value;
613 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToIn.value));
614 tmpp->
set(p.
x(), p.
y(), p.
z());
632 return fLastDistanceToIn.value;
645 surfaces[0] = fSide0;
646 surfaces[1] = fSide90 ;
647 surfaces[2] = fSide180 ;
648 surfaces[3] = fSide270 ;
649 surfaces[4] = fLowerEndcap;
650 surfaces[5] = fUpperEndcap;
658 if (tmpdistance < distance)
660 distance = tmpdistance;
665 return fLastDistanceToIn.value;
670 G4Exception(
"G4VTwistedFaceted::DistanceToIn(p)",
"GeomSolids0003",
702 if (fLastDistanceToOutWithV.p == p && fLastDistanceToOutWithV.vec == v )
704 return fLastDistanceToOutWithV.value;
708 tmpp =
const_cast<G4ThreeVector*
>(&(fLastDistanceToOutWithV.p));
709 tmpv =
const_cast<G4ThreeVector*
>(&(fLastDistanceToOutWithV.vec));
710 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToOutWithV.value));
711 tmpp->
set(p.
x(), p.
y(), p.
z());
712 tmpv->
set(v.
x(), v.
y(), v.
z());
735 *norm = (blockedsurface->
GetNormal(p,
true));
740 return fLastDistanceToOutWithV.value;
752 surfaces[0] = fSide0;
753 surfaces[1] = fSide90 ;
754 surfaces[2] = fSide180 ;
755 surfaces[3] = fSide270 ;
756 surfaces[4] = fLowerEndcap;
757 surfaces[5] = fUpperEndcap;
763 for (i=0; i< 6 ; i++) {
765 if (tmpdistance < distance)
767 distance = tmpdistance;
777 *norm = (surfaces[besti]->
GetNormal(p,
true));
784 return fLastDistanceToOutWithV.value;
804 if (fLastDistanceToOut.p == p)
806 return fLastDistanceToOut.value;
811 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToOut.value));
812 tmpp->
set(p.
x(), p.
y(), p.
z());
834 G4cout.precision(oldprc) ;
835 G4Exception(
"G4VTwistedFaceted::DistanceToOut(p)",
"GeomSolids1002",
843 retval = fLastDistanceToOut.value;
857 surfaces[0] = fSide0;
858 surfaces[1] = fSide90 ;
859 surfaces[2] = fSide180 ;
860 surfaces[3] = fSide270 ;
861 surfaces[4] = fLowerEndcap;
862 surfaces[5] = fUpperEndcap;
870 if (tmpdistance < distance)
872 distance = tmpdistance;
878 retval = fLastDistanceToOut.value;
884 G4Exception(
"G4VTwistedFaceted::DistanceToOut(p)",
"GeomSolids0003",
902 G4int oldprc = os.precision(16);
903 os <<
"-----------------------------------------------------------\n"
904 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
905 <<
" ===================================================\n"
906 <<
" Solid type: G4VTwistedFaceted\n"
908 <<
" polar angle theta = " << fTheta/
degree <<
" deg" <<
G4endl
909 <<
" azimuthal angle phi = " << fPhi/
degree <<
" deg" <<
G4endl
910 <<
" tilt angle alpha = " << fAlph/
degree <<
" deg" <<
G4endl
911 <<
" TWIST angle = " << fPhiTwist/
degree <<
" deg" <<
G4endl
912 <<
" Half length along y (lower endcap) = " << fDy1/
cm <<
" cm"
914 <<
" Half length along x (lower endcap, bottom) = " << fDx1/
cm <<
" cm"
916 <<
" Half length along x (lower endcap, top) = " << fDx2/
cm <<
" cm"
918 <<
" Half length along y (upper endcap) = " << fDy2/
cm <<
" cm"
920 <<
" Half length along x (upper endcap, bottom) = " << fDx3/
cm <<
" cm"
922 <<
" Half length along x (upper endcap, top) = " << fDx4/
cm <<
" cm"
924 <<
"-----------------------------------------------------------\n";
925 os.precision(oldprc);
945 G4double maxRad = std::sqrt( fDx*fDx + fDy*fDy);
956 void G4VTwistedFaceted::CreateSurfaces()
961 if ( fDx1 == fDx2 && fDx3 == fDx4 )
964 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 0.*
deg);
966 fDy1, fDx1, fDx1, fDy2, fDx3, fDx3, fAlph, 180.*
deg);
971 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*
deg);
973 fPhi+
pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*
deg);
979 fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*
deg);
981 fPhi+
pi, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*
deg);
986 fDz, fAlph, fPhi, fTheta, 1 );
988 fDz, fAlph, fPhi, fTheta, -1 );
992 fSide0->
SetNeighbours( fSide270 , fLowerEndcap , fSide90 , fUpperEndcap );
993 fSide90->
SetNeighbours( fSide0 , fLowerEndcap , fSide180 , fUpperEndcap );
994 fSide180->
SetNeighbours(fSide90 , fLowerEndcap , fSide270 , fUpperEndcap );
995 fSide270->
SetNeighbours(fSide180 , fLowerEndcap , fSide0 , fUpperEndcap );
996 fUpperEndcap->
SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
997 fLowerEndcap->
SetNeighbours( fSide180, fSide270 , fSide0 , fSide90 );
1007 return G4String(
"G4VTwistedFaceted");
1042 if ( z == fDz ) z -= 0.1*fDz ;
1043 if ( z == -fDz ) z += 0.1*fDz ;
1045 G4double phi = z/(2*fDz)*fPhiTwist ;
1047 return G4ThreeVector(fdeltaX * phi/fPhiTwist, fdeltaY * phi/fPhiTwist, z ) ;
1094 else if( (chose >= a1) && (chose < a1 + a2 ) )
1105 else if( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1115 else if( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1125 else if( (chose >= a1 + a2 + a3 + a4 ) && (chose < a1 + a2 + a3 + a4 + a5 ) )
1158 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
1159 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
1163 typedef G4int G4int4[4];
1164 G4double3* xyz =
new G4double3[nnodes];
1165 G4int4* faces =
new G4int4[nfaces] ;
1167 fLowerEndcap->
GetFacets(k,k,xyz,faces,0) ;
1168 fUpperEndcap->
GetFacets(k,k,xyz,faces,1) ;
void set(double x, double y, double z)
ThreeVector shoot(const G4int Ap, const G4int Af)
virtual void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *)
static constexpr double mm
static const G4double kInfinity
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 fRebuildPolyhedron
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
virtual void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4VTwistedFaceted & operator=(const G4VTwistedFaceted &rhs)
virtual void AddSolid(const G4Box &)=0
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
#define G4MUTEX_INITIALIZER
G4bool IsValidNorm() const
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
static constexpr double twopi
G4double Xcoef(G4double u, G4double phi, G4double ftg) const
static double normal(HepRandomEngine *eptr)
virtual G4double GetSurfaceArea()=0
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)
static constexpr double degree
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)=0
static constexpr double cm
virtual ~G4VTwistedFaceted()
G4double GetValueB(G4double phi) const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
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
virtual G4Polyhedron * GetPolyhedron() const
static G4int GetNumberOfRotationSteps()
virtual G4VisExtent GetExtent() const
virtual EInside Inside(const G4ThreeVector &p) const
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
G4VSolid & operator=(const G4VSolid &rhs)
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
static constexpr double pi
static constexpr double deg
virtual std::ostream & StreamInfo(std::ostream &os) const
G4Polyhedron * fpPolyhedron
G4double GetAngularTolerance() const
static G4GeometryTolerance * GetInstance()