42 #if !defined(G4GEOM_USE_UPOLYCONE)
58 using namespace CLHEP;
85 for (i=0; i<numZPlanes; i++)
87 if(rInner[i]>rOuter[i])
91 message <<
"Cannot create a Polycone with rInner > rOuter for the same Z"
93 <<
" rInner > rOuter for the same Z !" <<
G4endl
94 <<
" rMin[" << i <<
"] = " << rInner[i]
95 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
96 G4Exception(
"G4Polycone::G4Polycone()",
"GeomSolids0002",
99 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
101 if( (rInner[i] > rOuter[i+1])
102 ||(rInner[i+1] > rOuter[i]) )
106 message <<
"Cannot create a Polycone with no contiguous segments."
108 <<
" Segments are not contiguous !" <<
G4endl
109 <<
" rMin[" << i <<
"] = " << rInner[i]
110 <<
" -- rMax[" << i+1 <<
"] = " << rOuter[i+1] <<
G4endl
111 <<
" rMin[" << i+1 <<
"] = " << rInner[i+1]
112 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
113 G4Exception(
"G4Polycone::G4Polycone()",
"GeomSolids0002",
131 Create( phiStart, phiTotal, rz );
151 Create( phiStart, phiTotal, rz );
161 message <<
"Polycone " <<
GetName() <<
"cannot be converted" <<
G4endl
162 <<
"to Polycone with (Rmin,Rmaz,Z) parameters!";
163 G4Exception(
"G4Polycone::G4Polycone()",
"GeomSolids0002",
169 <<
"to optimized polycone with (Rmin,Rmaz,Z) parameters !"
188 if (rz->
Amin() < 0.0)
191 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
192 <<
" All R values must be >= 0 !";
193 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
205 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
206 <<
" R/Z cross section is zero or near zero: " << rzArea;
207 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
215 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
216 <<
" Too few unique R/Z values !";
217 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
224 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
225 <<
" R/Z segments cross !";
226 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
236 if (phiTotal <= 0 || phiTotal >
twopi-1
E-10)
253 endPhi = phiStart+phiTotal;
272 next->
r = iterRZ.
GetA();
273 next->
z = iterRZ.
GetB();
274 }
while( ++next, iterRZ.
Next() );
306 if (corner->
z > next->
z)
323 }
while( prev=corner, corner=next, corner >
corners );
352 :
G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
353 numCorner(0), corners(0),original_parameters(0),
385 if (
this == &source)
return *
this;
535 if (corner.
r < rmin) rmin = corner.
r;
536 if (corner.
r > rmax) rmax = corner.
r;
537 if (corner.
z < zmin) zmin = corner.
z;
538 if (corner.
z > zmax) zmax = corner.
z;
548 pMin.
set(vmin.
x(),vmin.
y(),zmin);
549 pMax.
set(vmax.
x(),vmax.
y(),zmax);
553 pMin.
set(-rmax,-rmax, zmin);
554 pMax.
set( rmax, rmax, zmax);
559 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
562 message <<
"Bad bounding box (min >= max) for solid: "
564 <<
"\npMin = " << pMin
565 <<
"\npMax = " << pMax;
566 G4Exception(
"G4Polycone::BoundingLimits()",
"GeomMgt0001",
589 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
593 return exist = (pMin < pMax) ?
true :
false;
602 std::vector<G4int> iout;
614 if (area < 0.)
std::reverse(contourRZ.begin(),contourRZ.end());
620 message <<
"Triangulation of RZ contour has failed for solid: "
622 <<
"\nExtent has been calculated using boundary box";
629 const G4int NSTEPS = 24;
635 G4int ksteps = (dphi <= astep) ? 1 : (
G4int)((dphi-
deg)/astep) + 1;
638 G4double sinHalf = std::sin(0.5*ang);
639 G4double cosHalf = std::cos(0.5*ang);
640 G4double sinStep = 2.*sinHalf*cosHalf;
641 G4double cosStep = 1. - 2.*sinHalf*sinHalf;
649 std::vector<const G4ThreeVectorList *> polygons;
650 polygons.resize(ksteps+2);
652 for (
G4int k=0; k<ksteps+2; ++k) pols[k].resize(6);
653 for (
G4int k=0; k<ksteps+2; ++k) polygons[k] = &pols[k];
660 G4int ntria = triangles.size()/3;
661 for (
G4int i=0; i<ntria; ++i)
664 for (
G4int k=0; k<3; ++k)
666 G4int e0 = i3+k, e1 = (k<2) ? e0+1 : i3;
669 r0[k2+0] = triangles[e0].x(); z0[k2+0] = triangles[e0].y();
670 r0[k2+1] = triangles[e1].x(); z0[k2+1] = triangles[e1].y();
674 if (z0[k2+1] - z0[k2+0] <= 0)
continue;
680 G4double sinCur = sinStart*cosHalf + cosStart*sinHalf;
681 G4double cosCur = cosStart*cosHalf - sinStart*sinHalf;
682 for (
G4int j=0; j<6; ++j) pols[0][j].set(r0[j]*cosStart,r0[j]*sinStart,z0[j]);
683 for (
G4int k=1; k<ksteps+1; ++k)
685 for (
G4int j=0; j<6; ++j) pols[k][j].set(r1[j]*cosCur,r1[j]*sinCur,z0[j]);
687 sinCur = sinCur*cosStep + cosCur*sinStep;
688 cosCur = cosCur*cosStep - sinTmp*sinStep;
690 for (
G4int j=0; j<6; ++j) pols[ksteps+1][j].set(r0[j]*cosEnd,r0[j]*sinEnd,z0[j]);
695 if (!benv.
CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax))
continue;
696 if (emin < pMin) pMin = emin;
697 if (emax > pMax) pMax =
emax;
698 if (eminlim > pMin && emaxlim < pMax)
return true;
700 return (pMin < pMax);
734 G4int oldprc = os.precision(16);
735 os <<
"-----------------------------------------------------------\n"
736 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
737 <<
" ===================================================\n"
738 <<
" Solid type: G4Polycone\n"
741 <<
" ending phi angle : " <<
endPhi/
degree <<
" degrees \n";
745 os <<
" number of Z planes: " << numPlanes <<
"\n"
747 for (i=0; i<numPlanes; i++)
749 os <<
" Z plane " << i <<
": "
752 os <<
" Tangent distances to inner surface (Rmin): \n";
753 for (i=0; i<numPlanes; i++)
755 os <<
" Z plane " << i <<
": "
758 os <<
" Tangent distances to outer surface (Rmax): \n";
759 for (i=0; i<numPlanes; i++)
761 os <<
" Z plane " << i <<
": "
765 os <<
" number of RZ points: " <<
numCorner <<
"\n"
766 <<
" RZ values (corners): \n";
772 os <<
"-----------------------------------------------------------\n";
773 os.precision(oldprc);
791 G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
792 G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
793 G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
796 rone = (fRmax1-fRmax2)/(2.*fDz);
797 rtwo = (fRmin1-fRmin2)/(2.*fDz);
798 if(fRmax1==fRmax2){qone=0.;}
801 qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
803 if(fRmin1==fRmin2){qtwo=0.;}
806 qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
808 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(
sqr(fRmin1-fRmin2)+
sqr(zTwo-zOne));
809 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(
sqr(fRmax1-fRmax2)+
sqr(zTwo-zOne));
810 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
811 totArea = Aone+Atwo+2.*Afive;
814 cosu = std::cos(phi);
815 sinu = std::sin(phi);
818 if( (startPhi == 0) && (
endPhi ==
twopi) ) { Afive = 0; }
820 if( (chose >= 0) && (chose < Aone) )
826 rone*sinu*(qone-zRand), zRand);
835 else if(chose >= Aone && chose < Aone + Atwo)
841 rtwo*sinu*(qtwo-zRand), zRand);
850 else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
853 rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
854 rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
857 rRand1*std::sin(startPhi), zRand);
862 rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
863 rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
866 rRand1*std::sin(
endPhi), zRand);
883 G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
884 aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
885 fDz = std::fabs(0.5*(zTwo-zOne));
889 aOne = 2.*fDz*fDPhi*fRMax;
890 aTwo = 2.*fDz*fDPhi*fRMin;
891 aFou = 2.*fDz*(fRMax-fRMin);
892 totArea = aOne+aTwo+2.*aFou;
894 cosphi = std::cos(phi);
895 sinphi = std::sin(phi);
902 if( (chose >= 0) && (chose < aOne) )
904 xRand = fRMax*cosphi;
905 yRand = fRMax*sinphi;
909 else if( (chose >= aOne) && (chose < aOne + aTwo) )
911 xRand = fRMin*cosphi;
912 yRand = fRMin*sinphi;
916 else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
918 xRand = rRand*std::cos(fSPhi+fDPhi);
919 yRand = rRand*std::sin(fSPhi+fDPhi);
926 xRand = rRand*std::cos(fSPhi+fDPhi);
927 yRand = rRand*std::sin(fSPhi+fDPhi);
942 G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
944 cosphi = std::cos(phi);
945 sinphi = std::sin(phi);
949 rRand1 = fRMin1; A1=0.;
954 A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
958 rRand2=fRMax1; Atot=A1;
963 Atot = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
967 if(rCh>A1) { rRand1=rRand2; }
969 xRand = rRand1*cosphi;
970 yRand = rRand1*sinphi;
989 if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
993 return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
1002 G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
1007 cosphi = std::cos(phi);
1008 sinphi = std::sin(phi);
1014 std::vector<G4double> areas;
1015 std::vector<G4ThreeVector> points;
1020 for(i=0; i<numPlanes-1; i++)
1045 areas.push_back(Area);
1052 totArea += (areas[0]+areas[numPlanes]);
1055 if( (chose>=0.) && (chose<areas[0]) )
1061 for (i=0; i<numPlanes-1; i++)
1063 Achose1 += areas[i];
1064 Achose2 = (Achose1+areas[i+1]);
1065 if(chose>=Achose1 && chose<Achose2)
1105 G4bool isConvertible=
true;
1111 std::vector<G4double>
Z;
1112 std::vector<G4double> Rmin;
1113 std::vector<G4double> Rmax;
1115 G4int countPlanes=1;
1126 Rmax.push_back (
corners[1].r);icurr=1;
1128 else if (Zprev ==
corners[numPlanes-1].z)
1130 Rmin.push_back(
corners[numPlanes-1].r);
1131 Rmax.push_back (
corners[0].r);
1137 Rmax.push_back (
corners[0].r);
1142 G4int inextr=0, inextl=0;
1143 for (
G4int i=0; i < numPlanes-2; i++)
1146 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1163 Rmin.push_back(
corners[inextl].r);
1164 Rmax.push_back(
corners[icurr].r);
1168 Rmin.push_back(
corners[inextl].r);
1177 Rmin.push_back(
corners[icurl].r);
1178 Rmax.push_back(
corners[icurr].r);
1182 Rmin.push_back(
corners[icurl].r);
1189 isConvertible=
false;
break;
1191 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1199 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1201 Rmin.push_back(
corners[inextl].r);
1202 Rmax.push_back(
corners[inextr].r);
1206 Z.push_back(Zright);
1215 Rmax.push_back(
corners[inextr].r);
1216 Rmin.push_back(
corners[icurr].r);
1220 Rmin.push_back(
corners[icurl].r + (Zright-
corners[icurl].z)/difZl
1222 Rmax.push_back(
corners[inextr].r);
1230 Rmax.push_back(
corners[inextr].r);
1231 Rmin.push_back (
corners[icurr].r);
1235 Rmax.push_back(
corners[inextr].r);
1243 isConvertible=
false;
break;
1253 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1255 Rmax.push_back(
corners[inextr].r);
1256 Rmin.push_back(
corners[inextl].r);
1267 for(
G4int j=0; j < countPlanes; j++)
1283 <<
"cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1284 G4Exception(
"G4Polycone::SetOriginalParameters()",
"GeomSolids0002",
1292 for(
G4int j=0; j < numPlanes; j++)
1302 return isConvertible;
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4PolyconeSideRZ GetCorner(G4int index) const
void set(double x, double y, double z)
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
CLHEP::Hep3Vector G4ThreeVector
G4double GetMaxExtent(const EAxis pAxis) const
std::vector< ExP01TrackerHit * > a
void SetOriginalParameters(G4PolyconeHistorical *pars)
static const G4double kInfinity
G4Polycone(const G4String &name, G4double phiStart, G4double phiTotal, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
CLHEP::Hep2Vector G4TwoVector
void message(RunManager *runmanager)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
static const G4double emax
void CopyStuff(const G4Polycone &source)
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double GetSinStartPhi() const
G4bool MustBeOutside(const G4ThreeVector &p) const
virtual EInside Inside(const G4ThreeVector &p) const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4EnclosingCylinder * enclosingCylinder
G4ThreeVector GetPointOnCut(G4double fRMin1, G4double fRMax1, G4double fRMin2, G4double fRMax2, G4double zOne, G4double zTwo, G4double &totArea) const
std::vector< G4ThreeVector > G4ThreeVectorList
G4double GetMinExtent(const EAxis pAxis) const
G4double GetCosStartPhi() const
static constexpr double deg
G4int GetNumRZCorner() const
G4bool RemoveDuplicateVertices(G4double tolerance)
G4ThreeVector GetPointOnRing(G4double fRMin, G4double fRMax, G4double fRMin2, G4double fRMax2, G4double zOne) const
G4Polyhedron * CreatePolyhedron() const
static constexpr double twopi
void Create(G4double phiStart, G4double phiTotal, G4ReduciblePolygon *rz)
G4ThreeVector GetPointOnSurface() const
std::ostream & StreamInfo(std::ostream &os) const
G4bool CrossesItself(G4double tolerance)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double GetEndPhi() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4int NumVertices() const
EInside Inside(const G4ThreeVector &p) const
G4double GetStartPhi() const
G4bool RemoveRedundantVertices(G4double tolerance)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
G4Polycone & operator=(const G4Polycone &source)
G4PolyconeHistorical * original_parameters
static constexpr double degree
G4GLOB_DLL std::ostream G4cout
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4ThreeVector GetPointOnCone(G4double fRmin1, G4double fRmax1, G4double fRmin2, G4double fRmax2, G4double zOne, G4double zTwo, G4double &totArea) const
static constexpr double pi
G4ThreeVector GetPointOnTubs(G4double fRMin, G4double fRMax, G4double zOne, G4double zTwo, G4double &totArea) const
G4GeometryType GetEntityType() const
G4PolyconeSideRZ * corners
G4double GetCosEndPhi() const
G4bool fRebuildPolyhedron
G4double GetSinEndPhi() const
G4Polyhedron * fpPolyhedron