69 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
70 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
74 G4Exception(
"G4TwistedTubs::G4TwistedTubs()",
"GeomSolids0002",
78 G4double sinhalftwist = std::sin(0.5 * twistedangle);
80 G4double endinnerradX = endinnerrad * sinhalftwist;
81 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
82 - endinnerradX * endinnerradX );
84 G4double endouterradX = endouterrad * sinhalftwist;
85 G4double outerrad = std::sqrt( endouterrad * endouterrad
86 - endouterradX * endouterradX );
89 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
101 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
102 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
107 std::ostringstream message;
108 message <<
"Invalid number of segments." <<
G4endl
109 <<
" nseg = " << nseg;
110 G4Exception(
"G4TwistedTubs::G4TwistedTubs()",
"GeomSolids0002",
115 G4Exception(
"G4TwistedTubs::G4TwistedTubs()",
"GeomSolids0002",
119 G4double sinhalftwist = std::sin(0.5 * twistedangle);
121 G4double endinnerradX = endinnerrad * sinhalftwist;
122 G4double innerrad = std::sqrt( endinnerrad * endinnerrad
123 - endinnerradX * endinnerradX );
125 G4double endouterradX = endouterrad * sinhalftwist;
126 G4double outerrad = std::sqrt( endouterrad * endouterrad
127 - endouterradX * endouterradX );
130 fDPhi = totphi / nseg;
131 SetFields(twistedangle, innerrad, outerrad, -halfzlen, halfzlen);
143 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
144 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
148 G4Exception(
"G4TwistedTubs::G4TwistedTubs()",
"GeomSolids0002",
152 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
165 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr),
166 fFormerTwisted(nullptr), fInnerHype(nullptr), fOuterHype(nullptr)
170 std::ostringstream message;
171 message <<
"Invalid number of segments." <<
G4endl
172 <<
" nseg = " << nseg;
173 G4Exception(
"G4TwistedTubs::G4TwistedTubs()",
"GeomSolids0002",
178 G4Exception(
"G4TwistedTubs::G4TwistedTubs()",
"GeomSolids0002",
182 fDPhi = totphi / nseg;
183 SetFields(twistedangle, innerrad, outerrad, negativeEndz, positiveEndz);
192 fLowerEndcap(nullptr), fUpperEndcap(nullptr),
193 fLatterTwisted(nullptr), fFormerTwisted(nullptr),
194 fInnerHype(nullptr), fOuterHype(nullptr)
205 delete fLatterTwisted;
206 delete fFormerTwisted;
209 delete fpPolyhedron; fpPolyhedron =
nullptr;
216 :
G4VSolid(rhs), fPhiTwist(rhs.fPhiTwist),
217 fInnerRadius(rhs.fInnerRadius), fOuterRadius(rhs.fOuterRadius),
218 fDPhi(rhs.fDPhi), fZHalfLength(rhs.fZHalfLength),
219 fInnerStereo(rhs.fInnerStereo), fOuterStereo(rhs.fOuterStereo),
220 fTanInnerStereo(rhs.fTanInnerStereo), fTanOuterStereo(rhs.fTanOuterStereo),
221 fKappa(rhs.fKappa), fInnerRadius2(rhs.fInnerRadius2),
222 fOuterRadius2(rhs.fOuterRadius2), fTanInnerStereo2(rhs.fTanInnerStereo2),
223 fTanOuterStereo2(rhs.fTanOuterStereo2),
224 fLowerEndcap(nullptr), fUpperEndcap(nullptr), fLatterTwisted(nullptr), fFormerTwisted(nullptr),
225 fInnerHype(nullptr), fOuterHype(nullptr),
226 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
227 fLastInside(rhs.fLastInside), fLastNormal(rhs.fLastNormal),
228 fLastDistanceToIn(rhs.fLastDistanceToIn),
229 fLastDistanceToOut(rhs.fLastDistanceToOut),
230 fLastDistanceToInWithV(rhs.fLastDistanceToInWithV),
231 fLastDistanceToOutWithV(rhs.fLastDistanceToOutWithV)
233 for (
auto i=0; i<2; ++i)
235 fEndZ[i] = rhs.fEndZ[i];
236 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
237 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
238 fEndPhi[i] = rhs.fEndPhi[i];
239 fEndZ2[i] = rhs.fEndZ2[i];
252 if (
this == &rhs) {
return *
this; }
260 fPhiTwist= rhs.fPhiTwist;
261 fInnerRadius= rhs.fInnerRadius; fOuterRadius= rhs.fOuterRadius;
262 fDPhi= rhs.fDPhi; fZHalfLength= rhs.fZHalfLength;
263 fInnerStereo= rhs.fInnerStereo; fOuterStereo= rhs.fOuterStereo;
264 fTanInnerStereo= rhs.fTanInnerStereo; fTanOuterStereo= rhs.fTanOuterStereo;
265 fKappa= rhs.fKappa; fInnerRadius2= rhs.fInnerRadius2;
266 fOuterRadius2= rhs.fOuterRadius2; fTanInnerStereo2= rhs.fTanInnerStereo2;
267 fTanOuterStereo2= rhs.fTanOuterStereo2;
268 fLowerEndcap= fUpperEndcap= fLatterTwisted= fFormerTwisted=
nullptr;
269 fInnerHype= fOuterHype=
nullptr;
270 fCubicVolume= rhs.fCubicVolume; fSurfaceArea= rhs.fSurfaceArea;
271 fLastInside= rhs.fLastInside; fLastNormal= rhs.fLastNormal;
272 fLastDistanceToIn= rhs.fLastDistanceToIn;
273 fLastDistanceToOut= rhs.fLastDistanceToOut;
274 fLastDistanceToInWithV= rhs.fLastDistanceToInWithV;
275 fLastDistanceToOutWithV= rhs.fLastDistanceToOutWithV;
277 for (
auto i=0; i<2; ++i)
279 fEndZ[i] = rhs.fEndZ[i];
280 fEndInnerRadius[i] = rhs.fEndInnerRadius[i];
281 fEndOuterRadius[i] = rhs.fEndOuterRadius[i];
282 fEndPhi[i] = rhs.fEndPhi[i];
283 fEndZ2[i] = rhs.fEndZ2[i];
287 fRebuildPolyhedron =
false;
288 delete fpPolyhedron; fpPolyhedron =
nullptr;
302 "G4TwistedTubs does not support Parameterisation.");
324 if (dphi <= 0 || totalphi >= CLHEP::twopi)
326 pMin.
set(-rmax,-rmax, zmin);
327 pMax.
set( rmax, rmax, zmax);
333 pMin.
set(vmin.
x(), vmin.
y(), zmin);
334 pMax.
set(vmax.
x(), vmax.
y(), zmax);
339 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
341 std::ostringstream message;
342 message <<
"Bad bounding box (min >= max) for solid: "
344 <<
"\npMin = " << pMin
345 <<
"\npMax = " << pMax;
346 G4Exception(
"G4TwistedTubs::BoundingLimits()",
"GeomMgt0001",
386 if (fLastInside.p == p)
388 return fLastInside.inside;
393 tmpinside =
const_cast<EInside*
>(&(fLastInside.inside));
394 tmpp->
set(p.
x(), p.
y(), p.
z());
401 if ((outerhypearea ==
kOutside) || (distanceToOut < -halftol))
411 if (distanceToOut <= halftol)
421 return fLastInside.inside;
436 if (fLastNormal.p == p)
438 return fLastNormal.vec;
441 auto tmpnormal =
const_cast<G4ThreeVector*
>(&(fLastNormal.vec));
443 tmpp->set(p.
x(), p.
y(), p.
z());
448 surfaces[0] = fLatterTwisted;
449 surfaces[1] = fFormerTwisted;
450 surfaces[2] = fInnerHype;
451 surfaces[3] = fOuterHype;
452 surfaces[4] = fLowerEndcap;
453 surfaces[5] = fUpperEndcap;
458 for (
auto i=0; i<6; ++i)
461 if (tmpdistance < distance)
463 distance = tmpdistance;
469 tmpsurface[0] = surfaces[besti];
470 *tmpnormal = tmpsurface[0]->
GetNormal(bestxx,
true);
472 return fLastNormal.vec;
495 if ((fLastDistanceToInWithV.p == p) && (fLastDistanceToInWithV.vec == v))
497 return fLastDistanceToIn.value;
501 tmpp =
const_cast<G4ThreeVector*
>(&(fLastDistanceToInWithV.p));
502 tmpv =
const_cast<G4ThreeVector*
>(&(fLastDistanceToInWithV.vec));
503 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToInWithV.value));
504 tmpp->
set(p.
x(), p.
y(), p.
z());
505 tmpv->
set(v.
x(), v.
y(), v.
z());
528 return fLastDistanceToInWithV.value;
542 surfaces[0] = fLowerEndcap;
543 surfaces[1] = fUpperEndcap;
544 surfaces[2] = fLatterTwisted;
545 surfaces[3] = fFormerTwisted;
546 surfaces[4] = fInnerHype;
547 surfaces[5] = fOuterHype;
551 for (
const auto & surface : surfaces)
553 G4double tmpdistance = surface->DistanceToIn(p, v, xx);
554 if (tmpdistance < distance)
556 distance = tmpdistance;
562 return fLastDistanceToInWithV.value;
580 if (fLastDistanceToIn.p == p)
582 return fLastDistanceToIn.value;
587 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToIn.value));
588 tmpp->
set(p.
x(), p.
y(), p.
z());
604 return fLastDistanceToIn.value;
613 surfaces[0] = fLowerEndcap;
614 surfaces[1] = fUpperEndcap;
615 surfaces[2] = fLatterTwisted;
616 surfaces[3] = fFormerTwisted;
617 surfaces[4] = fInnerHype;
618 surfaces[5] = fOuterHype;
622 for (
const auto & surface : surfaces)
624 G4double tmpdistance = surface->DistanceTo(p, xx);
625 if (tmpdistance < distance)
627 distance = tmpdistance;
632 return fLastDistanceToIn.value;
636 G4Exception(
"G4TwistedTubs::DistanceToIn(p)",
"GeomSolids0003",
666 if ((fLastDistanceToOutWithV.p == p) && (fLastDistanceToOutWithV.vec == v) )
668 return fLastDistanceToOutWithV.value;
672 tmpp =
const_cast<G4ThreeVector*
>(&(fLastDistanceToOutWithV.p));
673 tmpv =
const_cast<G4ThreeVector*
>(&(fLastDistanceToOutWithV.vec));
674 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToOutWithV.value));
675 tmpp->
set(p.
x(), p.
y(), p.
z());
676 tmpv->
set(v.
x(), v.
y(), v.
z());
701 *norm = (blockedsurface->
GetNormal(p,
true));
705 return fLastDistanceToOutWithV.value;
719 surfaces[0] = fLatterTwisted;
720 surfaces[1] = fFormerTwisted;
721 surfaces[2] = fInnerHype;
722 surfaces[3] = fOuterHype;
723 surfaces[4] = fLowerEndcap;
724 surfaces[5] = fUpperEndcap;
729 for (
auto i=0; i<6; ++i)
732 if (tmpdistance < distance)
734 distance = tmpdistance;
744 *norm = (surfaces[besti]->
GetNormal(p,
true));
751 return fLastDistanceToOutWithV.value;
770 if (fLastDistanceToOut.p == p)
772 return fLastDistanceToOut.value;
777 tmpdist =
const_cast<G4double*
>(&(fLastDistanceToOut.value));
778 tmpp->
set(p.
x(), p.
y(), p.
z());
795 return fLastDistanceToOut.value;
804 surfaces[0] = fLatterTwisted;
805 surfaces[1] = fFormerTwisted;
806 surfaces[2] = fInnerHype;
807 surfaces[3] = fOuterHype;
808 surfaces[4] = fLowerEndcap;
809 surfaces[5] = fUpperEndcap;
813 for (
const auto & surface : surfaces)
815 G4double tmpdistance = surface->DistanceTo(p, xx);
816 if (tmpdistance < distance)
818 distance = tmpdistance;
824 return fLastDistanceToOut.value;
828 G4Exception(
"G4TwistedTubs::DistanceToOut(p)",
"GeomSolids0003",
844 G4long oldprc = os.precision(16);
845 os <<
"-----------------------------------------------------------\n"
846 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
847 <<
" ===================================================\n"
848 <<
" Solid type: G4TwistedTubs\n"
850 <<
" -ve end Z : " << fEndZ[0]/mm <<
" mm \n"
851 <<
" +ve end Z : " << fEndZ[1]/mm <<
" mm \n"
852 <<
" inner end radius(-ve z): " << fEndInnerRadius[0]/mm <<
" mm \n"
853 <<
" inner end radius(+ve z): " << fEndInnerRadius[1]/mm <<
" mm \n"
854 <<
" outer end radius(-ve z): " << fEndOuterRadius[0]/mm <<
" mm \n"
855 <<
" outer end radius(+ve z): " << fEndOuterRadius[1]/mm <<
" mm \n"
856 <<
" inner radius (z=0) : " << fInnerRadius/mm <<
" mm \n"
857 <<
" outer radius (z=0) : " << fOuterRadius/mm <<
" mm \n"
858 <<
" twisted angle : " << fPhiTwist/degree <<
" degrees \n"
859 <<
" inner stereo angle : " << fInnerStereo/degree <<
" degrees \n"
860 <<
" outer stereo angle : " << fOuterStereo/degree <<
" degrees \n"
861 <<
" phi-width of a piece : " << fDPhi/degree <<
" degrees \n"
862 <<
"-----------------------------------------------------------\n";
863 os.precision(oldprc);
886 return { pmin.
x(),pmax.
x(),
898 G4double absPhiTwist = std::abs(fPhiTwist);
899 G4double dA = std::max(fDPhi,absPhiTwist);
905 const G4int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
906 const G4int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
910 typedef G4int G4int4[4];
911 auto xyz =
new G4double3[nnodes];
912 auto faces =
new G4int4[nfaces] ;
913 fLowerEndcap->
GetFacets(k,k,xyz,faces,0) ;
914 fUpperEndcap->
GetFacets(k,k,xyz,faces,1) ;
916 fFormerTwisted->
GetFacets(k,n,xyz,faces,3) ;
918 fLatterTwisted->
GetFacets(k,n,xyz,faces,5) ;
920 ph->createPolyhedron(nnodes,nfaces,xyz,faces);
933 if (fpPolyhedron ==
nullptr ||
934 fRebuildPolyhedron ||
941 fRebuildPolyhedron =
false;
950void G4TwistedTubs::CreateSurfaces()
958 fEndInnerRadius, fEndOuterRadius,
959 fDPhi, fEndPhi, fEndZ, -1) ;
962 fEndInnerRadius, fEndOuterRadius,
963 fDPhi, fEndPhi, fEndZ, 1) ;
966 rotHalfDPhi.
rotateZ(0.5*fDPhi);
969 fEndInnerRadius, fEndOuterRadius,
970 fDPhi, fEndPhi, fEndZ,
971 fInnerRadius, fOuterRadius, fKappa,
974 fEndInnerRadius, fEndOuterRadius,
975 fDPhi, fEndPhi, fEndZ,
976 fInnerRadius, fOuterRadius, fKappa,
980 fEndInnerRadius, fEndOuterRadius,
981 fDPhi, fEndPhi, fEndZ,
982 fInnerRadius, fOuterRadius,fKappa,
983 fTanInnerStereo, fTanOuterStereo, -1) ;
985 fEndInnerRadius, fEndOuterRadius,
986 fDPhi, fEndPhi, fEndZ,
987 fInnerRadius, fOuterRadius,fKappa,
988 fTanInnerStereo, fTanOuterStereo, 1) ;
994 fOuterHype, fFormerTwisted);
996 fOuterHype, fFormerTwisted);
998 fOuterHype, fUpperEndcap);
1000 fOuterHype, fUpperEndcap);
1002 fFormerTwisted, fUpperEndcap);
1004 fFormerTwisted, fUpperEndcap);
1013 return {
"G4TwistedTubs"};
1029 if (fCubicVolume == 0.)
1042 fCubicVolume = (2.*(Z1 - Z0)*(Aout + Ain)*(Aout - Ain)
1043 + Z1*(R1out + R1in)*(R1out - R1in)
1044 - Z0*(R0out + R0in)*(R0out - R0in))*DPhi/6.;
1046 return fCubicVolume;
1055 if (z == 0)
return 0.;
1064 G4double k = std::sqrt(aa + cc)/cc;
1066 area = 0.5*a*(h*std::sqrt(1. + kh*kh) + std::asinh(kh)/k);
1077 if (
GetDPhi() >= CLHEP::twopi || r <= 0 || z == 0)
return 0.;
1088 G4double sqroot = std::sqrt(pp + qq + 1);
1090 0.5*p*(
pp + 3.)*std::atanh(q/sqroot) +
1091 0.5*q*(qq + 3.)*std::atanh(p/sqroot) +
1092 std::atan(sqroot/(pq)) - CLHEP::halfpi)*h*a/(3.*pq);
1102 if (fSurfaceArea == 0.)
1114 G4double base0 = 0.5*dphi*(Rout0*Rout0 - Rinn0*Rinn0);
1115 G4double inner0 = GetLateralArea(Ainn, Rinn0, z0);
1116 G4double outer0 = GetLateralArea(Aout, Rout0, z0);
1118 GetPhiCutArea(Aout, Rout0, z0) - GetPhiCutArea(Ainn, Rinn0, z0);
1124 if (std::abs(z0) != std::abs(z1))
1126 base1 = 0.5*dphi*(Rout1*Rout1 - Rinn1*Rinn1);
1127 inner1 = GetLateralArea(Ainn, Rinn1, z1);
1128 outer1 = GetLateralArea(Aout, Rout1, z1);
1130 GetPhiCutArea(Aout, Rout1, z1) - GetPhiCutArea(Ainn, Rinn1, z1);
1132 fSurfaceArea = base0 + base1 +
1134 (inner0 + inner1 + outer0 + outer1 + 2.*(cut0 + cut1)) :
1135 std::abs(inner0 - inner1 + outer0 - outer1 + 2.*(cut0 - cut1)));
1137 return fSurfaceArea;
1146 G4double z = G4RandFlat::shoot(fEndZ[0],fEndZ[1]);
1158 G4double chose = G4RandFlat::shoot(0.,a1 + a2 + a3 + a4 + a5 + a6) ;
1165 phi = G4RandFlat::shoot(phimin,phimax) ;
1170 else if ( (chose >= a1) && (chose < a1 + a2 ) )
1175 phi = G4RandFlat::shoot(phimin,phimax) ;
1180 else if ( (chose >= a1 + a2 ) && (chose < a1 + a2 + a3 ) )
1185 x = G4RandFlat::shoot(xmin,xmax) ;
1190 else if ( (chose >= a1 + a2 + a3 ) && (chose < a1 + a2 + a3 + a4 ) )
1195 x = G4RandFlat::shoot(xmin,xmax) ;
1199 else if( (chose >= a1 + a2 + a3 + a4 )&&(chose < a1 + a2 + a3 + a4 + a5 ) )
1203 r = std::sqrt(G4RandFlat::shoot()*(
sqr(rmax)-
sqr(rmin))+
sqr(rmin));
1207 phi = G4RandFlat::shoot(phimin,phimax) ;
1215 r = rmin + (rmax-rmin)*std::sqrt(G4RandFlat::shoot());
1219 phi = G4RandFlat::shoot(phimin,phimax) ;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
#define G4MUTEX_INITIALIZER
void set(double x, double y, double z)
HepRotation & rotateZ(double delta)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double GetRadialTolerance() const
static G4GeometryTolerance * GetInstance()
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcnorm=false, G4bool *validnorm=nullptr, G4ThreeVector *n=nullptr) const override
G4double GetOuterRadius() const
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const override
G4double GetZHalfLength() const
G4double GetPhiTwist() const
G4Polyhedron * CreatePolyhedron() const override
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const override
G4TwistedTubs & operator=(const G4TwistedTubs &rhs)
G4double GetEndInnerRadius() const
G4double GetEndOuterRadius() const
G4VSolid * Clone() const override
G4Polyhedron * GetPolyhedron() const override
~G4TwistedTubs() override
G4GeometryType GetEntityType() const override
G4double GetEndPhi(G4int i) const
G4TwistedTubs(const G4String &pname, G4double twistedangle, G4double endinnerrad, G4double endouterrad, G4double halfzlen, G4double dphi)
std::ostream & StreamInfo(std::ostream &os) const override
G4double GetSurfaceArea() override
G4double GetCubicVolume() override
void ComputeDimensions(G4VPVParameterisation *, const G4int, const G4VPhysicalVolume *) override
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const override
G4VisExtent GetExtent() const override
G4ThreeVector GetPointOnSurface() const override
void DescribeYourselfTo(G4VGraphicsScene &scene) const override
G4double GetEndZ(G4int i) const
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const override
G4double GetInnerRadius() const
EInside Inside(const G4ThreeVector &p) const override
virtual void AddSolid(const G4Box &)=0
G4VSolid & operator=(const G4VSolid &rhs)
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)=0
void SetNeighbours(G4VTwistSurface *ax0min, G4VTwistSurface *ax1min, G4VTwistSurface *ax0max, G4VTwistSurface *ax1max)
virtual G4double GetBoundaryMin(G4double)=0
G4bool IsValidNorm() const
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 G4double GetBoundaryMax(G4double)=0
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
virtual G4double GetSurfaceArea()=0
static G4int GetNumberOfRotationSteps()