59 #if !defined(G4GEOM_USE_UPOLYHEDRA)
77 using namespace CLHEP;
97 message <<
"Solid must have at least one side - " <<
GetName() <<
G4endl
98 <<
" No sides specified !";
99 G4Exception(
"G4Polyhedra::G4Polyhedra()",
"GeomSolids0002",
108 { phiTotal =
twopi; }
109 G4double convertRad = std::cos(0.5*phiTotal/theNumSide);
125 for (i=0; i<numZPlanes; i++)
127 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
129 if( (rInner[i] > rOuter[i+1])
130 ||(rInner[i+1] > rOuter[i]) )
134 message <<
"Cannot create a Polyhedra with no contiguous segments."
136 <<
" Segments are not contiguous !" <<
G4endl
137 <<
" rMin[" << i <<
"] = " << rInner[i]
138 <<
" -- rMax[" << i+1 <<
"] = " << rOuter[i+1] <<
G4endl
139 <<
" rMin[" << i+1 <<
"] = " << rInner[i+1]
140 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
141 G4Exception(
"G4Polyhedra::G4Polyhedra()",
"GeomSolids0002",
156 rz->
ScaleA( 1/convertRad );
161 Create( phiStart, phiTotal, theNumSide, rz );
182 message <<
"Solid must have at least one side - " <<
GetName() <<
G4endl
183 <<
" No sides specified !";
184 G4Exception(
"G4Polyhedra::G4Polyhedra()",
"GeomSolids0002",
190 Create( phiStart, phiTotal, theNumSide, rz );
214 if (rz->
Amin() < 0.0)
217 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
218 <<
" All R values must be >= 0 !";
219 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
231 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
232 <<
" R/Z cross section is zero or near zero: " << rzArea;
233 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
241 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
242 <<
" Too few unique R/Z values !";
243 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
250 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
251 <<
" R/Z segments cross !";
252 G4Exception(
"G4Polyhedra::Create()",
"GeomSolids0002",
278 endPhi = phiStart+phiTotal;
302 next->
r = iterRZ.
GetA();
303 next->
z = iterRZ.
GetB();
304 }
while( ++next, iterRZ.
Next() );
356 }
while( prev=corner, corner=next, corner >
corners );
386 phiIsOpen(false), genericPgon(false), numCorner(0), corners(0),
387 original_parameters(0), enclosingCylinder(0)
419 if (
this == &source)
return *
this;
491 message <<
"Solid " <<
GetName() <<
" built using generic construct."
492 <<
G4endl <<
"Not applicable to the generic construct !";
493 G4Exception(
"G4Polyhedra::Reset()",
"GeomSolids1001",
584 if (corner.
r < rmin) rmin = corner.
r;
585 if (corner.
r > rmax) rmax = corner.
r;
586 if (corner.
z < zmin) zmin = corner.
z;
587 if (corner.
z > zmax) zmax = corner.
z;
601 G4double xmin = rmin*cosCur, xmax = xmin;
603 for (
G4int k=0; k<ksteps+1; ++k)
606 if (x < xmin) xmin =
x;
607 if (x > xmax) xmax =
x;
609 if (y < ymin) ymin =
y;
610 if (y > ymax) ymax =
y;
614 if (xx < xmin) xmin =
xx;
615 if (xx > xmax) xmax =
xx;
617 if (yy < ymin) ymin = yy;
618 if (yy > ymax) ymax = yy;
621 sinCur = sinCur*cosStep + cosCur*sinStep;
622 cosCur = cosCur*cosStep - sinTmp*sinStep;
624 pMin.
set(xmin,ymin,zmin);
625 pMax.
set(xmax,ymax,zmax);
629 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
632 message <<
"Bad bounding box (min >= max) for solid: "
634 <<
"\npMin = " << pMin
635 <<
"\npMax = " << pMax;
636 G4Exception(
"G4Polyhedra::BoundingLimits()",
"GeomMgt0001",
659 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
663 return exist = (pMin < pMax) ?
true :
false;
672 std::vector<G4int> iout;
684 if (area < 0.)
std::reverse(contourRZ.begin(),contourRZ.end());
690 message <<
"Triangulation of RZ contour has failed for solid: "
692 <<
"\nExtent has been calculated using boundary box";
710 std::vector<const G4ThreeVectorList *> polygons;
711 polygons.resize(ksteps+1);
712 for (
G4int k=0; k<ksteps+1; ++k) {
719 G4int ntria = triangles.size()/3;
720 for (
G4int i=0; i<ntria; ++i)
725 for (
G4int k=0; k<ksteps+1; ++k)
728 G4ThreeVectorList::iterator iter = ptr->begin();
729 iter->set(triangles[i3+0].
x()*cosCur,
730 triangles[i3+0].
x()*sinCur,
731 triangles[i3+0].
y());
733 iter->set(triangles[i3+1].
x()*cosCur,
734 triangles[i3+1].
x()*sinCur,
735 triangles[i3+1].
y());
737 iter->set(triangles[i3+2].
x()*cosCur,
738 triangles[i3+2].
x()*sinCur,
739 triangles[i3+2].
y());
742 sinCur = sinCur*cosStep + cosCur*sinStep;
743 cosCur = cosCur*cosStep - sinTmp*sinStep;
749 if (!benv.
CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax))
continue;
750 if (emin < pMin) pMin = emin;
751 if (emax > pMax) pMax =
emax;
752 if (eminlim > pMin && emaxlim < pMax)
break;
755 for (
G4int k=0; k<ksteps+1; ++k) {
delete polygons[k]; polygons[k]=0;}
756 return (pMin < pMax);
793 G4int oldprc = os.precision(16);
794 os <<
"-----------------------------------------------------------\n"
795 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
796 <<
" ===================================================\n"
797 <<
" Solid type: G4Polyhedra\n"
800 <<
" ending phi angle : " <<
endPhi/
degree <<
" degrees \n"
801 <<
" number of sides : " <<
numSide <<
" \n";
806 os <<
" number of Z planes: " << numPlanes <<
"\n"
808 for (i=0; i<numPlanes; i++)
810 os <<
" Z plane " << i <<
": "
813 os <<
" Tangent distances to inner surface (Rmin): \n";
814 for (i=0; i<numPlanes; i++)
816 os <<
" Z plane " << i <<
": "
819 os <<
" Tangent distances to outer surface (Rmax): \n";
820 for (i=0; i<numPlanes; i++)
822 os <<
" Z plane " << i <<
": "
826 os <<
" number of RZ points: " <<
numCorner <<
"\n"
827 <<
" RZ values (corners): \n";
833 os <<
"-----------------------------------------------------------\n";
834 os.precision(oldprc);
849 G4double lambda1, lambda2, chose,aOne,aTwo;
860 if( (chose>=0.) && (chose < aOne) )
864 return (p2+lambda1*v+lambda2*w);
869 return (p0+lambda1*t+lambda2*u);
888 return (p2 + lambda1*w + lambda2*v);
900 G4double chose, totArea=0., Achose1, Achose2,
901 rad1, rad2, sinphi1, sinphi2, cosphi1, cosphi2;
903 area, aTop=0., aBottom=0., zVal=0.;
906 std::vector<G4double> aVector1;
907 std::vector<G4double> aVector2;
908 std::vector<G4double> aVector3;
916 for(j=0; j<numPlanes-1; j++)
922 area = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+b)*cosksi;
923 aVector1.push_back(area);
926 for(j=0; j<numPlanes-1; j++)
932 area = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+b)*cosksi;
933 aVector2.push_back(area);
936 for(j=0; j<numPlanes-1; j++)
947 else { aVector3.push_back(0.); }
950 for(j=0; j<numPlanes-1; j++)
952 totArea +=
numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
962 aTop = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+b)*cosksi;
970 aBottom = std::sqrt(l2-
sqr((a-b)*cosksi))*(a+b)*cosksi;
974 Achose2 =
numSide*(aVector1[0]+aVector2[0])+2.*aVector3[0];
977 if( (chose >= 0.) && (chose < aTop + aBottom) )
980 rang = std::floor((chose-
startPhi)/ksi-0.01);
981 if(rang<0) { rang=0; }
982 rang = std::fabs(rang);
983 sinphi1 = std::sin(
startPhi+rang*ksi);
984 sinphi2 = std::sin(
startPhi+(rang+1)*ksi);
985 cosphi1 = std::cos(
startPhi+rang*ksi);
986 cosphi2 = std::cos(
startPhi+(rang+1)*ksi);
988 if(chose>=0. && chose<aTop)
1008 for (j=0; j<numPlanes-1; j++)
1010 if( ((chose >= Achose1) && (chose < Achose2)) || (j == numPlanes-2) )
1014 Achose1 +=
numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
1015 Achose2 = Achose1 +
numSide*(aVector1[j+1]+aVector2[j+1])
1025 totArea =
numSide*(aVector1[j]+aVector2[j])+2.*aVector3[j];
1028 if( (chose>=0.) && (chose<
numSide*aVector1[j]) )
1031 rang = std::floor((chose-
startPhi)/ksi-0.01);
1032 if(rang<0) { rang=0; }
1033 rang = std::fabs(rang);
1036 sinphi1 = std::sin(
startPhi+rang*ksi);
1037 sinphi2 = std::sin(
startPhi+(rang+1)*ksi);
1038 cosphi1 = std::cos(
startPhi+rang*ksi);
1039 cosphi2 = std::cos(
startPhi+(rang+1)*ksi);
1051 else if ( (chose >=
numSide*aVector1[j])
1052 && (chose <=
numSide*(aVector1[j]+aVector2[j])) )
1055 rang = std::floor((chose-
startPhi)/ksi-0.01);
1056 if(rang<0) { rang=0; }
1057 rang = std::fabs(rang);
1060 sinphi1 = std::sin(
startPhi+rang*ksi);
1061 sinphi2 = std::sin(
startPhi+(rang+1)*ksi);
1062 cosphi1 = std::cos(
startPhi+rang*ksi);
1063 cosphi2 = std::cos(
startPhi+(rang+1)*ksi);
1077 if( (chose>=0.) && (chose < 1.) )
1152 typedef G4int int4[4];
1159 std::vector<G4bool> chopped(
numCorner,
false);
1160 std::vector<G4int*> triQuads;
1163 while (remaining >= 3)
1168 G4int iStepper = iStarter;
1171 if (A < 0) { A = iStepper; }
1172 else if (
B < 0) {
B = iStepper; }
1173 else if (
C < 0) {
C = iStepper; }
1176 if (++iStepper >=
numCorner) iStepper = 0;
1178 while (chopped[iStepper]);
1180 while (
C < 0 && iStepper != iStarter);
1195 triQuads.push_back(tq);
1203 if (++iStarter >=
numCorner) { iStarter = 0; }
1205 while (chopped[iStarter]);
1213 faces_vec =
new int4[nFaces];
1217 for (
G4int iEnd = 0; iEnd < 2; ++iEnd)
1219 for (
size_t i = 0; i < triQuads.size(); ++i)
1232 a = triQuads[i][0] + addition;
1233 b = triQuads[i][2] + addition;
1234 c = triQuads[i][1] + addition;
1237 G4int bc = std::abs(c - b);
1238 G4int ca = std::abs(a - c);
1239 faces_vec[iface][0] = (ab == 1 || ab ==
d)? a: -a;
1240 faces_vec[iface][1] = (bc == 1 || bc ==
d)? b: -b;
1241 faces_vec[iface][2] = (ca == 1 || ca ==
d)? c: -c;
1242 faces_vec[iface][3] = 0;
1249 xyz =
new double3[nNodes];
1257 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1258 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1260 if (iCorner < numCorner - 1)
1262 faces_vec[iface][0] = ixyz + 1;
1263 faces_vec[iface][1] = ixyz + numCorner + 1;
1264 faces_vec[iface][2] = ixyz + numCorner + 2;
1265 faces_vec[iface][3] = ixyz + 2;
1269 faces_vec[iface][0] = ixyz + 1;
1270 faces_vec[iface][1] = ixyz + numCorner + 1;
1271 faces_vec[iface][2] = ixyz + 2;
1272 faces_vec[iface][3] = ixyz - numCorner + 2;
1284 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1285 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1294 xyz =
new double3[nNodes];
1295 faces_vec =
new int4[nFaces];
1300 G4int ixyz = 0, iface = 0;
1305 xyz[ixyz][0] =
corners[iCorner].
r * std::cos(phi);
1306 xyz[ixyz][1] =
corners[iCorner].
r * std::sin(phi);
1308 if (iSide < numSide - 1)
1310 if (iCorner < numCorner - 1)
1312 faces_vec[iface][0] = ixyz + 1;
1313 faces_vec[iface][1] = ixyz + numCorner + 1;
1314 faces_vec[iface][2] = ixyz + numCorner + 2;
1315 faces_vec[iface][3] = ixyz + 2;
1319 faces_vec[iface][0] = ixyz + 1;
1320 faces_vec[iface][1] = ixyz + numCorner + 1;
1321 faces_vec[iface][2] = ixyz + 2;
1322 faces_vec[iface][3] = ixyz - numCorner + 2;
1327 if (iCorner < numCorner - 1)
1329 faces_vec[iface][0] = ixyz + 1;
1330 faces_vec[iface][1] = ixyz + numCorner - nFaces + 1;
1331 faces_vec[iface][2] = ixyz + numCorner - nFaces + 2;
1332 faces_vec[iface][3] = ixyz + 2;
1336 faces_vec[iface][0] = ixyz + 1;
1337 faces_vec[iface][1] = ixyz - nFaces + numCorner + 1;
1338 faces_vec[iface][2] = ixyz - nFaces + 2;
1339 faces_vec[iface][3] = ixyz - numCorner + 2;
1350 delete [] faces_vec;
1355 message <<
"Problem creating G4Polyhedron for: " <<
GetName();
1356 G4Exception(
"G4Polyhedra::CreatePolyhedron()",
"GeomSolids1002",
1372 G4bool isConvertible=
true;
1378 std::vector<G4double>
Z;
1379 std::vector<G4double> Rmin;
1380 std::vector<G4double> Rmax;
1382 G4int countPlanes=1;
1393 Rmax.push_back (
corners[1].r);icurr=1;
1395 else if (Zprev ==
corners[numPlanes-1].z)
1397 Rmin.push_back(
corners[numPlanes-1].r);
1398 Rmax.push_back (
corners[0].r);
1404 Rmax.push_back (
corners[0].r);
1409 G4int inextr=0, inextl=0;
1410 for (
G4int i=0; i < numPlanes-2; i++)
1413 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1430 Rmin.push_back(
corners[inextl].r);
1431 Rmax.push_back(
corners[icurr].r);
1435 Rmin.push_back(
corners[inextl].r);
1444 Rmin.push_back(
corners[icurl].r);
1445 Rmax.push_back(
corners[icurr].r);
1449 Rmin.push_back(
corners[icurl].r);
1456 isConvertible=
false;
break;
1458 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1466 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1468 Rmin.push_back(
corners[inextl].r);
1469 Rmax.push_back (
corners[inextr].r);
1473 Z.push_back(Zright);
1482 Rmax.push_back(
corners[inextr].r);
1483 Rmin.push_back(
corners[icurr].r);
1487 Rmin.push_back(
corners[icurl].r + (Zright-
corners[icurl].z)/difZl
1489 Rmax.push_back(
corners[inextr].r);
1497 Rmax.push_back(
corners[inextr].r);
1498 Rmin.push_back (
corners[icurr].r);
1502 Rmax.push_back(
corners[inextr].r);
1510 isConvertible=
false;
break;
1520 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1522 Rmax.push_back(
corners[inextr].r);
1523 Rmin.push_back(
corners[inextl].r);
1535 for(
G4int j=0; j < countPlanes; j++)
1551 <<
"cannot be converted to Polyhedra with (Rmin,Rmaz,Z) parameters!";
1552 G4Exception(
"G4Polyhedra::SetOriginalParameters()",
1561 for(
G4int j=0; j < numPlanes; j++)
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
void set(double x, double y, double z)
G4int createPolyhedron(G4int Nnodes, G4int Nfaces, const G4double xyz[][3], const G4int faces[][4])
CLHEP::Hep3Vector G4ThreeVector
G4Polyhedron * CreatePolyhedron() const
G4double GetMaxExtent(const EAxis pAxis) const
std::vector< ExP01TrackerHit * > a
static const G4double kInfinity
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
G4GeometryType GetEntityType() const
void CopyStuff(const G4Polyhedra &source)
CLHEP::Hep2Vector G4TwoVector
void ScaleA(G4double scale)
void message(RunManager *runmanager)
G4double GetSinStartPhi() const
static double normal(HepRandomEngine *eptr)
static const G4double emax
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool MustBeOutside(const G4ThreeVector &p) const
virtual EInside Inside(const G4ThreeVector &p) const
void SetOriginalParameters(G4PolyhedraHistorical *pars)
G4Polyhedra & operator=(const G4Polyhedra &source)
EInside Inside(const G4ThreeVector &p) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
std::vector< G4ThreeVector > G4ThreeVectorList
G4double GetStartPhi() const
G4double GetMinExtent(const EAxis pAxis) const
G4Polyhedra(const G4String &name, G4double phiStart, G4double phiTotal, G4int numSide, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
G4bool RemoveDuplicateVertices(G4double tolerance)
G4PolyhedraHistorical * original_parameters
static constexpr double twopi
G4double GetEndPhi() const
double A(double temperature)
G4bool CrossesItself(G4double tolerance)
ThreeVector shoot(const G4int Ap, const G4int Af)
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4PolyhedraSideRZ GetCorner(const G4int index) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4int NumVertices() const
G4bool RemoveRedundantVertices(G4double tolerance)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
G4double GetCosStartPhi() const
G4int GetNumRZCorner() const
G4ThreeVector GetPointOnSurfaceGeneric() const
static constexpr double degree
G4ThreeVector GetPointOnTriangle(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4PolyhedraSideRZ * corners
G4ThreeVector GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3) const
void Create(G4double phiStart, G4double phiTotal, G4int numSide, G4ReduciblePolygon *rz)
double B(double temperature)
G4bool fRebuildPolyhedron
std::ostream & StreamInfo(std::ostream &os) const
G4EnclosingCylinder * enclosingCylinder
G4ThreeVector GetPointOnSurface() const
G4Polyhedron * fpPolyhedron
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const