53 #if !defined(G4GEOM_USE_UEXTRUDEDSOLID)
75 const std::vector<G4TwoVector>& polygon,
76 const std::vector<ZSection>& zsections)
79 fNz(zsections.size()),
81 fGeometryType(
"G4ExtrudedSolid"),
91 message <<
"Number of vertices in polygon < 3 - " << pName;
92 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
99 message <<
"Number of z-sides < 2 - " << pName;
100 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
106 if ( zsections[i].fZ > zsections[i+1].fZ )
109 message <<
"Z-sections have to be ordered by z value (z0 < z1 < z2...) - "
111 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
117 message <<
"Z-sections with the same z position are not supported - "
119 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0001",
130 std::vector<G4int> removedVertices;
133 if (removedVertices.size() != 0)
135 G4int nremoved = removedVertices.size();
137 message <<
"The following "<< nremoved
138 <<
" vertices have been removed from polygon in " << pName
139 <<
"\nas collinear or coincident with other vertices: "
140 << removedVertices[0];
141 for (
G4int i=1; i<nremoved; ++i) message <<
", " << removedVertices[i];
142 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids1001",
150 message <<
"Number of vertices in polygon after removal < 3 - " << pName;
151 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
175 message <<
"Making facets failed - " << pName;
176 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0003",
198 const std::vector<G4TwoVector>& polygon,
206 fGeometryType(
"G4ExtrudedSolid"),
216 message <<
"Number of vertices in polygon < 3 - " << pName;
217 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
227 std::vector<G4int> removedVertices;
230 if (removedVertices.size() != 0)
232 G4int nremoved = removedVertices.size();
234 message <<
"The following "<< nremoved
235 <<
" vertices have been removed from polygon in " << pName
236 <<
"\nas collinear or coincident with other vertices: "
237 << removedVertices[0];
238 for (
G4int i=1; i<nremoved; ++i) message <<
", " << removedVertices[i];
239 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids1001",
247 message <<
"Number of vertices in polygon after removal < 3 - " << pName;
248 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0002",
273 message <<
"Making facets failed - " << pName;
274 G4Exception(
"G4ExtrudedSolid::G4ExtrudedSolid()",
"GeomSolids0003",
283 if ((scale1 == 1) && (scale2 == 1)
295 fGeometryType(
"G4ExtrudedSolid"), fSolidType(0)
305 fPolygon(rhs.fPolygon), fZSections(rhs.fZSections),
306 fTriangles(rhs.fTriangles), fIsConvex(rhs.fIsConvex),
307 fGeometryType(rhs.fGeometryType),
308 fSolidType(rhs.fSolidType), fPlanes(rhs.fPlanes),
309 fLines(rhs.fLines), fLengths(rhs.fLengths),
310 fKScales(rhs.fKScales), fScale0s(rhs.fScale0s),
311 fKOffsets(rhs.fKOffsets), fOffset0s(rhs.fOffset0s)
321 if (
this == &rhs) {
return *
this; }
369 G4double kscale = (scale2 - scale1)/(z2 - z1);
370 G4double scale0 = scale2 - kscale*(z2 - z1)/2.0;
389 for (
G4int i=0, k=Nv-1; i<Nv; k=i++)
403 for (
G4int i=0, k=Nv-1; i<Nv; k=i++)
460 return (p2 - poffset)/pscale;
471 if ( l1.
x() == l2.
x() )
486 G4bool squareComp = (dy*dy < (1+slope*slope)
522 return ( (p1.
x() - l1.
x()) * (l2.
y() - l1.
y())
523 - (l2.
x() - l1.
x()) * (p1.
y() - l1.
y()) )
524 * ( (p2.
x() - l1.
x()) * (l2.
y() - l1.
y())
525 - (l2.
x() - l1.
x()) * (p2.
y() - l1.
y()) ) > 0;
540 if ( ( p.
x() < a.
x() && p.
x() < b.
x() && p.
x() < c.
x() ) ||
541 ( p.
x() > a.
x() && p.
x() > b.
x() && p.
x() > c.
x() ) ||
542 ( p.
y() < a.
y() && p.
y() < b.
y() && p.
y() < c.
y() ) ||
543 ( p.
y() > a.
y() && p.
y() > b.
y() && p.
y() > c.
y() ) )
return false;
555 return inside || onEdge;
572 if ( result < 0 ) result += 2*
pi;
585 std::vector<G4ThreeVector> vertices;
593 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
595 if ( cross.
z() > 0.0 )
603 vertices[1] = vertices[2];
619 std::vector<G4ThreeVector> vertices;
627 = (vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]);
629 if ( cross.
z() < 0.0 )
637 vertices[1] = vertices[2];
651 typedef std::pair < G4TwoVector, G4int > Vertex;
653 static const G4double kAngTolerance =
658 std::vector< Vertex > verticesToBeDone;
661 verticesToBeDone.push_back(Vertex(
fPolygon[i], i));
663 std::vector< Vertex > ears;
665 std::vector< Vertex >::iterator
c1 = verticesToBeDone.begin();
666 std::vector< Vertex >::iterator
c2 = c1+1;
667 std::vector< Vertex >::iterator c3 = c1+2;
668 while ( verticesToBeDone.size()>2 )
685 while ( angle >= (
pi-kAngTolerance) )
694 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
700 angle =
GetAngle(c2->first, c3->first, c1->first);
705 if ( counter > fNv) {
706 G4Exception(
"G4ExtrudedSolid::AddGeneralPolygonFacets",
708 "Triangularisation has failed.");
714 std::vector< Vertex >::iterator it;
715 for ( it=verticesToBeDone.begin(); it != verticesToBeDone.end(); ++it )
719 if ( it == c1 || it == c2 || it == c3 ) {
continue; }
721 if (
IsPointInside(c1->first, c2->first, c3->first, it->first) )
731 if ( c3 == verticesToBeDone.end() ) { c3 = verticesToBeDone.begin(); }
747 if ( ! result ) {
return false; }
750 if ( ! result ) {
return false; }
752 std::vector<G4int> triangle(3);
753 triangle[0] = c1->second;
754 triangle[1] = c2->second;
755 triangle[2] = c3->second;
760 verticesToBeDone.erase(c2);
761 c1 = verticesToBeDone.begin();
783 if ( ! good ) {
return false; }
789 if ( ! good ) {
return false; }
791 std::vector<G4int> triangle(3);
803 if ( ! good ) {
return false; }
810 if ( ! good ) {
return false; }
812 std::vector<G4int> triangle1(3);
818 std::vector<G4int> triangle2(3);
827 if ( ! good ) {
return false; }
832 for (
G4int iz = 0; iz <
fNz-1; ++iz )
836 G4int j = (i+1) % fNv;
840 if ( ! good ) {
return false; }
877 for (
G4int i=0; i<np; ++i)
880 if (dd > dist) { dist = dd; }
930 G4int j = (i+1) % fNv;
942 std::vector< std::vector<G4int> >::const_iterator it =
fTriangles.begin();
947 fPolygon[(*it)[2]], pscaled) ) { inside =
true; }
949 }
while ( (inside ==
false) && (it !=
fTriangles.end()) );
1009 if (ix*ix + iy*iy > sqrCarToleranceHalf)
continue;
1015 if (kx*kx + ky*ky > sqrCarToleranceHalf)
continue;
1020 if (dd*dd > sqrCarToleranceHalf)
continue;
1040 else if (nsurf != 0)
1050 G4int oldprc = message.precision(16);
1051 message <<
"Point p is not on surface (!?) of solid: "
1053 message <<
"Position:\n";
1054 message <<
" p.x() = " << p.
x()/
mm <<
" mm\n";
1055 message <<
" p.y() = " << p.
y()/
mm <<
" mm\n";
1056 message <<
" p.z() = " << p.
z()/
mm <<
" mm";
1057 G4cout.precision(oldprc) ;
1058 G4Exception(
"G4TesselatedSolid::SurfaceNormal(p)",
"GeomSolids1002",
1094 if (tmp < dd) { dd =
tmp; iside = i; }
1101 if (tmp < dd) { dd =
tmp; iside = i; }
1107 if (tmp < dd) { dd =
tmp; iside = i; }
1120 if (
std::max(dz0,dz1) > 0) iregion = 1;
1123 if (!in) iregion += 2;
1131 if (ddz0 <= ddz1 && ddz0 <= dd)
return G4ThreeVector(0, 0,-1);
1132 if (ddz1 <= ddz0 && ddz1 <= dd)
return G4ThreeVector(0, 0, 1);
1146 if (dzmax*dzmax > dd)
return G4ThreeVector(0,0,(dz0 > dz1) ? -1 : 1);
1174 G4double ddz = (invz < 0) ? dz : -dz;
1181 G4double txmin = tzmin, txmax = tzmax;
1182 for (
G4int i=0; i<np; ++i)
1190 if (txmin < tmp) { txmin =
tmp; }
1195 if (txmax > tmp) { txmax =
tmp; }
1201 G4double tmin = txmin, tmax = txmax;
1225 for (
G4int i=0; i<np; ++i)
1228 if (dd > dist) dist = dd;
1230 return (dist > 0) ? dist : 0.;
1238 return (distz > 0) ? distz : 0;
1244 if (distz > 0) dd += distz*distz;
1245 return std::sqrt(dd);
1262 G4bool getnorm = calcNorm;
1263 if (getnorm) *validNorm =
true;
1269 if (getnorm) n->
set(0,0,-1);
1274 if (getnorm) n->
set(0,0,1);
1289 G4int iside = (vz < 0) ? -4 : -2;
1294 for (
G4int i=0; i<np; ++i)
1306 if (tmax > tmp) { tmax =
tmp; iside = i; }
1315 { n->
set(0, 0, iside + 3); }
1331 if (validNorm) { *validNorm =
fIsConvex; }
1346 for (
G4int i=0; i<np; ++i)
1349 if (dd > dist) dist = dd;
1351 return (dist < 0) ? -dist : 0.;
1357 if (distz >= 0 || (!in))
return 0;
1378 if (x < xmin0) xmin0 =
x;
1379 if (x > xmax0) xmax0 =
x;
1381 if (y < ymin0) ymin0 =
y;
1382 if (y > ymax0) ymax0 =
y;
1389 for (
G4int i=0; i<nsect; ++i)
1395 xmin =
std::min(xmin,xmin0*scale+dx);
1396 xmax =
std::max(xmax,xmax0*scale+dx);
1397 ymin =
std::min(ymin,ymin0*scale+dy);
1398 ymax =
std::max(ymax,ymax0*scale+dy);
1404 pMin.
set(xmin,ymin,zmin);
1405 pMax.
set(xmax,ymax,zmax);
1409 if (pMin.
x() >= pMax.
x() || pMin.
y() >= pMax.
y() || pMin.
z() >= pMax.
z())
1412 message <<
"Bad bounding box (min >= max) for solid: "
1414 <<
"\npMin = " << pMin
1415 <<
"\npMax = " << pMax;
1438 #ifdef G4BBOX_EXTENT
1439 if (
true)
return bbox.
CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1443 return exist = (pMin < pMax) ?
true :
false;
1458 message <<
"Triangulation of the base polygon has failed for solid: "
1460 <<
"\nExtent has been calculated using boundary box";
1468 std::vector<const G4ThreeVectorList *> polygons;
1469 polygons.resize(nsect);
1475 G4int ntria = triangles.size()/3;
1476 for (
G4int i=0; i<ntria; ++i)
1479 for (
G4int k=0; k<nsect; ++k)
1488 G4ThreeVectorList::iterator iter = ptr->begin();
1489 G4double x0 = triangles[i3+0].x()*scale+dx;
1490 G4double y0 = triangles[i3+0].y()*scale+dy;
1505 if (!benv.
CalculateExtent(pAxis,pVoxelLimit,pTransform,emin,emax))
continue;
1506 if (emin < pMin) pMin = emin;
1507 if (emax > pMax) pMax =
emax;
1508 if (eminlim > pMin && emaxlim < pMax)
break;
1511 for (
G4int k=0; k<nsect; ++k) {
delete polygons[k]; polygons[k]=0;}
1512 return (pMin < pMax);
1519 G4int oldprc = os.precision(16);
1520 os <<
"-----------------------------------------------------------\n"
1521 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
1522 <<
" ===================================================\n"
1526 { os <<
" Convex polygon; list of vertices:" <<
G4endl; }
1528 { os <<
" Concave polygon; list of vertices:" <<
G4endl; }
1532 os << std::setw(5) <<
"#" << i
1537 os <<
" Sections:" <<
G4endl;
1562 os.precision(oldprc);
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4double GetMinZExtent() const
void set(double x, double y, double z)
ZSection GetZSection(G4int index) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4TwoVector GetVertex(G4int index) const
std::vector< G4double > fKScales
G4bool IsSameLineSegment(const G4TwoVector &p, const G4TwoVector &l1, const G4TwoVector &l2) const
CLHEP::Hep3Vector G4ThreeVector
G4double GetMaxExtent(const EAxis pAxis) const
std::vector< ExP01TrackerHit * > a
G4ThreeVector ApproxSurfaceNormal(const G4ThreeVector &p) const
G4double GetMaxXExtent() const
void ComputeProjectionParameters()
static const G4double kInfinity
G4TwoVector ProjectPoint(const G4ThreeVector &point) const
Float_t y1[n_points_granero]
static constexpr double mm
G4VFacet * MakeUpFacet(G4int ind1, G4int ind2, G4int ind3) const
Float_t x1[n_points_granero]
G4TessellatedSolid & operator=(const G4TessellatedSolid &right)
G4bool IsPointInside(const G4TwoVector &a, const G4TwoVector &b, const G4TwoVector &c, const G4TwoVector &p) const
std::vector< ZSection > fZSections
CLHEP::Hep2Vector G4TwoVector
void message(RunManager *runmanager)
G4double GetAngularTolerance() const
void BoundingLimits(G4ThreeVector &pMin, G4ThreeVector &pMax) const
G4ExtrudedSolid(const G4String &pName, const std::vector< G4TwoVector > &polygon, const std::vector< ZSection > &zsections)
G4double GetAngle(const G4TwoVector &p0, const G4TwoVector &pa, const G4TwoVector &pb) const
static const G4double emax
G4double kCarToleranceHalf
std::vector< std::vector< G4int > > fTriangles
std::ostream & StreamInfo(std::ostream &os) const
G4bool IsSameSide(const G4TwoVector &p1, const G4TwoVector &p2, const G4TwoVector &l1, const G4TwoVector &l2) const
G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
Float_t y2[n_points_geant4]
void SetSolidClosed(const G4bool t)
G4bool AddFacet(G4VFacet *aFacet)
G4GeometryType fGeometryType
std::vector< G4double > fLengths
G4double GetMinXExtent() const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool BoundingBoxVsVoxelLimits(const EAxis pAxis, const G4VoxelLimits &pVoxelLimits, const G4Transform3D &pTransform3D, G4double &pMin, G4double &pMax) const
G4ExtrudedSolid & operator=(const G4ExtrudedSolid &rhs)
void ComputeLateralPlanes()
std::vector< G4ThreeVector > G4ThreeVectorList
G4double GetMaxZExtent() const
std::vector< G4TwoVector > fOffset0s
G4double GetMinExtent(const EAxis pAxis) const
virtual G4double DistanceToOut(const G4ThreeVector &p) const
G4bool AddGeneralPolygonFacets()
std::vector< line > fLines
G4int GetNofZSections() const
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< G4TwoVector > fKOffsets
G4double GetMaxYExtent() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4GeometryType GetEntityType() const
virtual ~G4ExtrudedSolid()
std::vector< plane > fPlanes
G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
EInside Inside(const G4ThreeVector &p) const
G4GLOB_DLL std::ostream G4cout
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const
G4bool PointInPolygon(const G4ThreeVector &p) const
G4bool IsSameLine(const G4TwoVector &p, const G4TwoVector &l1, const G4TwoVector &l2) const
static G4GeometryTolerance * GetInstance()
G4double DistanceToPolygonSqr(const G4ThreeVector &p) const
static constexpr double pi
std::vector< G4double > fScale0s
G4double GetMinYExtent() const
G4int GetNofVertices() const
std::vector< G4TwoVector > fPolygon
Float_t x2[n_points_geant4]
G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VFacet * MakeDownFacet(G4int ind1, G4int ind2, G4int ind3) const