59 #if !defined(G4GEOM_USE_UPOLYHEDRA)
73 using namespace CLHEP;
92 std::ostringstream message;
93 message <<
"Solid must have at least one side - " <<
GetName() <<
G4endl
94 <<
" No sides specified !";
95 G4Exception(
"G4Polyhedra::G4Polyhedra()",
"GeomSolids0002",
104 { phiTotal =
twopi; }
105 G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
121 for (i=0; i<numZPlanes; i++)
123 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
125 if( (rInner[i] > rOuter[i+1])
126 ||(rInner[i+1] > rOuter[i]) )
129 std::ostringstream message;
130 message <<
"Cannot create a Polyhedra with no contiguous segments."
132 <<
" Segments are not contiguous !" <<
G4endl
133 <<
" rMin[" << i <<
"] = " << rInner[i]
134 <<
" -- rMax[" << i+1 <<
"] = " << rOuter[i+1] <<
G4endl
135 <<
" rMin[" << i+1 <<
"] = " << rInner[i+1]
136 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
137 G4Exception(
"G4Polyhedra::G4Polyhedra()",
"GeomSolids0002",
152 rz->
ScaleA( 1/convertRad );
157 Create( phiStart, phiTotal, theNumSide, rz );
177 Create( phiStart, phiTotal, theNumSide, rz );
201 if (rz->
Amin() < 0.0)
203 std::ostringstream message;
204 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
205 <<
" All R values must be >= 0 !";
206 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
216 std::ostringstream message;
217 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
218 <<
" R/Z cross section is zero or near zero: " << rzArea;
219 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
226 std::ostringstream message;
227 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
228 <<
" Too few unique R/Z values !";
229 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
235 std::ostringstream message;
236 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
237 <<
" R/Z segments cross !";
238 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
263 endPhi = phiStart+phiTotal;
286 next->
r = iterRZ.
GetA();
287 next->
z = iterRZ.
GetB();
288 }
while( ++next, iterRZ.
Next() );
315 if (corner->
r < 1/kInfinity && next->
r < 1/kInfinity)
continue;
340 }
while( prev=corner, corner=next, corner >
corners );
370 phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
371 original_parameters(0), enclosingCylinder(0)
403 if (
this == &source)
return *
this;
471 std::ostringstream message;
472 message <<
"Solid " <<
GetName() <<
" built using generic construct."
473 <<
G4endl <<
"Not applicable to the generic construct !";
474 G4Exception(
"G4Polyhedra::Reset()",
"GeomSolids1001",
588 G4int oldprc = os.precision(16);
589 os <<
"-----------------------------------------------------------\n"
590 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
591 <<
" ===================================================\n"
592 <<
" Solid type: G4Polyhedra\n"
595 <<
" ending phi angle : " <<
endPhi/
degree <<
" degrees \n";
600 os <<
" number of Z planes: " << numPlanes <<
"\n"
602 for (i=0; i<numPlanes; i++)
604 os <<
" Z plane " << i <<
": "
607 os <<
" Tangent distances to inner surface (Rmin): \n";
608 for (i=0; i<numPlanes; i++)
610 os <<
" Z plane " << i <<
": "
613 os <<
" Tangent distances to outer surface (Rmax): \n";
614 for (i=0; i<numPlanes; i++)
616 os <<
" Z plane " << i <<
": "
620 os <<
" number of RZ points: " <<
numCorner <<
"\n"
621 <<
" RZ values (corners): \n";
627 os <<
"-----------------------------------------------------------\n";
628 os.precision(oldprc);
642 G4double lambda1, lambda2, chose,aOne,aTwo;
653 if( (chose>=0.) && (chose < aOne) )
657 return (p2+lambda1*v+lambda2*w);
662 return (p0+lambda1*t+lambda2*u);
681 return (p2 + lambda1*w + lambda2*v);
693 G4double chose, totArea=0., Achose1, Achose2,
694 rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
696 area, aTop=0., aBottom=0., zVal=0.;
699 std::vector<G4double> aVector1;
700 std::vector<G4double> aVector2;
701 std::vector<G4double> aVector3;
709 for(j=0; j<numPlanes-1; j++)
715 area = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+
b)*cosksi;
716 aVector1.push_back(area);
719 for(j=0; j<numPlanes-1; j++)
725 area = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+
b)*cosksi;
726 aVector2.push_back(area);
729 for(j=0; j<numPlanes-1; j++)
740 else { aVector3.push_back(0.); }
743 for(j=0; j<numPlanes-1; j++)
745 totArea +=
numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
755 aTop = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+
b)*cosksi;
763 aBottom = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+
b)*cosksi;
767 Achose2 =
numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
770 if( (chose >= 0.) && (chose < aTop + aBottom) )
773 rang = std::floor((chose-
startPhi)/ksi-0.01);
774 if(rang<0) { rang=0; }
775 rang = std::fabs(rang);
776 sinphi1 = std::sin(
startPhi+rang*ksi);
777 sinphi2 = std::sin(
startPhi+(rang+1)*ksi);
778 cosphi1 = std::cos(
startPhi+rang*ksi);
779 cosphi2 = std::cos(
startPhi+(rang+1)*ksi);
781 if(chose>=0. && chose<aTop)
801 for (j=0; j<numPlanes-1; j++)
803 if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
807 Achose1 +=
numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
808 Achose2 = Achose1 +
numSide*(aVector1[j+1]+aVector2[j+1])
818 totArea =
numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
821 if( (chose>=0.) && (chose<
numSide*aVector1[j]) )
824 rang = std::floor((chose-
startPhi)/ksi-0.01);
825 if(rang<0) { rang=0; }
826 rang = std::fabs(rang);
829 sinphi1 = std::sin(
startPhi+rang*ksi);
830 sinphi2 = std::sin(
startPhi+(rang+1)*ksi);
831 cosphi1 = std::cos(
startPhi+rang*ksi);
832 cosphi2 = std::cos(
startPhi+(rang+1)*ksi);
844 else if ( (chose >=
numSide*aVector1[j])
845 && (chose <=
numSide*(aVector1[j]+aVector2[j])) )
848 rang = std::floor((chose-
startPhi)/ksi-0.01);
849 if(rang<0) { rang=0; }
850 rang = std::fabs(rang);
853 sinphi1 = std::sin(
startPhi+rang*ksi);
854 sinphi2 = std::sin(
startPhi+(rang+1)*ksi);
855 cosphi1 = std::cos(
startPhi+rang*ksi);
856 cosphi2 = std::cos(
startPhi+(rang+1)*ksi);
870 if( (chose>=0.) && (chose < 1.) )
945 typedef G4int int4[4];
952 std::vector<G4bool> chopped(
numCorner,
false);
953 std::vector<G4int*> triQuads;
956 while (remaining >= 3)
960 G4int A = -1, B = -1, C = -1;
961 G4int iStepper = iStarter;
964 if (A < 0) { A = iStepper; }
965 else if (B < 0) { B = iStepper; }
966 else if (C < 0) { C = iStepper; }
969 if (++iStepper >=
numCorner) iStepper = 0;
971 while (chopped[iStepper]);
973 while (C < 0 && iStepper != iStarter);
988 triQuads.push_back(tq);
996 if (++iStarter >=
numCorner) { iStarter = 0; }
998 while (chopped[iStarter]);
1006 faces_vec =
new int4[nFaces];
1010 for (
G4int iEnd = 0; iEnd < 2; ++iEnd)
1012 for (
size_t i = 0; i < triQuads.size(); ++i)
1025 a = triQuads[i][0] + addition;
1026 b = triQuads[i][2] + addition;
1027 c = triQuads[i][1] + addition;
1029 G4int ab = std::abs(b - a);
1030 G4int bc = std::abs(c - b);
1031 G4int ca = std::abs(a - c);
1032 faces_vec[iface][0] = (ab == 1 || ab ==
d)? a: -a;
1033 faces_vec[iface][1] = (bc == 1 || bc ==
d)? b: -b;
1034 faces_vec[iface][2] = (ca == 1 || ca ==
d)? c: -c;
1035 faces_vec[iface][3] = 0;
1042 xyz =
new double3[nNodes];
1050 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1051 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1053 if (iCorner < numCorner - 1)
1055 faces_vec[iface][0] = ixyz + 1;
1056 faces_vec[iface][1] = ixyz + numCorner + 1;
1057 faces_vec[iface][2] = ixyz + numCorner + 2;
1058 faces_vec[iface][3] = ixyz + 2;
1062 faces_vec[iface][0] = ixyz + 1;
1063 faces_vec[iface][1] = ixyz + numCorner + 1;
1064 faces_vec[iface][2] = ixyz + 2;
1065 faces_vec[iface][3] = ixyz - numCorner + 2;
1077 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1078 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1087 xyz =
new double3[nNodes];
1088 faces_vec =
new int4[nFaces];
1092 G4int ixyz = 0, iface = 0;
1097 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1098 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1100 if (iSide < numSide - 1)
1102 if (iCorner < numCorner - 1)
1104 faces_vec[iface][0] = ixyz + 1;
1105 faces_vec[iface][1] = ixyz + numCorner + 1;
1106 faces_vec[iface][2] = ixyz + numCorner + 2;
1107 faces_vec[iface][3] = ixyz + 2;
1111 faces_vec[iface][0] = ixyz + 1;
1112 faces_vec[iface][1] = ixyz + numCorner + 1;
1113 faces_vec[iface][2] = ixyz + 2;
1114 faces_vec[iface][3] = ixyz - numCorner + 2;
1119 if (iCorner < numCorner - 1)
1121 faces_vec[iface][0] = ixyz + 1;
1122 faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
1123 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1124 faces_vec[iface][3] = ixyz + 2;
1128 faces_vec[iface][0] = ixyz + 1;
1129 faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
1130 faces_vec[iface][2] = ixyz - nFaces + 2;
1131 faces_vec[iface][3] = ixyz - numCorner + 2;
1142 delete [] faces_vec;
1146 std::ostringstream message;
1147 message <<
"Problem creating G4Polyhedron for: " <<
GetName();
1148 G4Exception(
"G4Polyhedra::CreatePolyhedron()",
"GeomSolids1002",
1164 G4bool isConvertible=
true;
1170 std::vector<G4double>
Z;
1171 std::vector<G4double> Rmin;
1172 std::vector<G4double> Rmax;
1174 G4int countPlanes=1;
1185 Rmax.push_back (
corners[1].r);icurr=1;
1187 else if (Zprev ==
corners[numPlanes-1].z)
1189 Rmin.push_back(
corners[numPlanes-1].
r);
1196 Rmax.push_back (
corners[0].r);
1201 G4int inextr=0, inextl=0;
1202 for (
G4int i=0; i < numPlanes-2; i++)
1205 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1222 Rmin.push_back(
corners[icurl].r);
1223 Rmax.push_back(Rmax[countPlanes-2]);
1224 Rmax[countPlanes-2]=
corners[icurl].
r;
1228 Rmin.push_back(
corners[inextl].r);
1229 Rmax.push_back(
corners[icurl].r);
1235 Rmax.push_back (
corners[icurr].r + (Zleft-
corners[icurr].z)/difZr
1240 isConvertible=
false;
break;
1242 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1250 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1253 Rmax.push_back (
corners[inextr].r);
1257 Z.push_back(Zright);
1266 Rmin.push_back(
corners[icurr].r);
1267 Rmax.push_back(
corners[inextr].r);
1271 Rmin.push_back(
corners[inextr].r);
1272 Rmax.push_back(
corners[icurr].r);
1273 Rmax[countPlanes-2]=
corners[inextr].
r;
1282 Rmin.push_back (
corners[icurr].r);
1294 isConvertible=
false;
break;
1304 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1309 Rmin.push_back(
corners[inextl].r);
1314 Rmin.push_back(
corners[inextl].r);
1327 for(
G4int j=0; j < countPlanes; j++)
1341 std::ostringstream message;
1343 <<
"cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1344 G4Exception(
"G4Polyhedra::SetOriginalParameters()",
1353 for(
G4int j=0; j < numPlanes; j++)
void CopyStuff(const G4Polyhedra &source)
G4bool CrossesItself(G4double tolerance)
ThreeVector shoot(const G4int Ap, const G4int Af)
CLHEP::Hep3Vector G4ThreeVector
G4GeometryType GetEntityType() const
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
G4PolyhedraHistorical * original_parameters
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4PolyhedraSideRZ * corners
const G4Polyhedra & operator=(const G4Polyhedra &source)
G4bool MustBeOutside(const G4ThreeVector &p) const
void Create(G4double phiStart, G4double phiTotal, G4int numSide, G4ReduciblePolygon *rz)
G4ThreeVector GetPointOnSurface() const
G4bool RemoveDuplicateVertices(G4double tolerance)
G4EnclosingCylinder * enclosingCylinder
G4bool RemoveRedundantVertices(G4double tolerance)
G4ThreeVector GetPointOnSurfaceGeneric() const
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4int NumVertices() const
void ScaleA(G4double scale)
const G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
EInside Inside(const G4ThreeVector &p) const
G4Polyhedron * CreatePolyhedron() const
void SetOriginalParameters(G4PolyhedraHistorical *pars)
G4ThreeVector GetPointOnTriangle(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2) const
virtual EInside Inside(const G4ThreeVector &p) const
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
std::ostream & StreamInfo(std::ostream &os) const
G4Polyhedra(const G4String &name, G4double phiStart, G4double phiTotal, G4int numSide, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])