36 const double zPlane[],
37 const double rInner[],
38 const double rOuter[])
43 Init(phiStart, phiTotal, numZPlanes, zPlane, rInner, rOuter);
66 std::ostringstream message;
67 message <<
"Polycone " <<
GetName() <<
"cannot be converted" << std::endl
68 <<
"to Polycone with (Rmin,Rmaz,Z) parameters! Use GenericPolycone" ;
74 std::cout <<
"INFO: Converting polycone " <<
GetName() << std::endl
75 <<
"to optimized polycone with (Rmin,Rmaz,Z) parameters !"
82 for (
int i = 0; i < num; i++)
90 Init(phiStart, phiTotal, num, Z, R1, R2);
106 const double zPlane[],
107 const double rInner[],
108 const double rOuter[])
127 endPhi = phiStart+phiTotal;
144 double RMaxextent=rOuter[0];
145 for (
int j=1; j < numZPlanes; j++)
147 if (rOuter[j] > RMaxextent) RMaxextent=rOuter[j];
148 if (rInner[j]>rOuter[j])
150 std::ostringstream message;
151 message <<
"Cannot create Polycone with rInner > rOuter for the same Z"
153 <<
" rInner > rOuter for the same Z !" << std::endl
154 <<
" rMin[" << j <<
"] = " << rInner[j]
155 <<
" -- rMax[" << j <<
"] = " << rOuter[j];
161 double prevZ = zPlane[0], prevRmax = 0, prevRmin = 0;
163 if (zPlane[1] < zPlane[0])dirZ = -1;
167 for (i = 0; i < numZPlanes; i++)
169 if ((i < numZPlanes - 1) && (zPlane[i] == zPlane[i + 1]))
171 if ((rInner[i] > rOuter[i + 1])
172 || (rInner[i + 1] > rOuter[i]))
175 std::ostringstream message;
176 message <<
"Cannot create a Polycone with no contiguous segments."
178 <<
"Segments are not contiguous !" << std::endl
179 <<
" rMin[" << i <<
"] = " << rInner[i]
180 <<
" -- rMax[" << i + 1 <<
"] = " << rOuter[i + 1] << std::endl
181 <<
" rMin[" << i + 1 <<
"] = " << rInner[i + 1]
182 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
190 double rMin = rInner[i];
191 double rMax = rOuter[i];
192 double z = zPlane[i];
196 if (((z > prevZ)&&(dirZ>0))||((z < prevZ)&&(dirZ<0)))
198 if (dirZ*(z-prevZ)< 0)
200 std::ostringstream message;
201 message <<
"Cannot create a Polycone with different Z directions.Use GenericPolycone."
203 <<
"ZPlane is changing direction !" << std::endl
204 <<
" zPlane[0] = " << zPlane[0]
205 <<
" -- zPlane[1] = " << zPlane[1] << std::endl
206 <<
" zPlane[" << i - 1 <<
"] = " << zPlane[i - 1]
207 <<
" -- rPlane[" << i <<
"] = " << zPlane[i];
215 double dz = (z - prevZ)*dirZ / 2;
217 bool tubular = (rMin == prevRmin && prevRmax == rMax);
223 solid =
new UTubs(
"", rMin, rMax, dz, phiStart, phiTotal);
227 solid =
new UCons(
"", prevRmin, prevRmax, rMin, rMax, dz, phiStart, phiTotal);
237 int zi =
fZs.size() - 1;
238 double shift =
fZs[zi - 1] + 0.5 * (
fZs[zi] -
fZs[zi - 1]);
241 section.
shift = shift;
243 section.
solid = solid;
246 if (rMax < RMaxextent) { section.
convex =
false;}
247 else { section.
convex =
true;}
251 if ((rMax<prevRmax)||(rMax < RMaxextent)||(prevRmax < RMaxextent))
252 { section.
convex =
false;}
263 else fZs.push_back(z);
281 double mxy = rz->
Amax();
284 double r = rz->
Amax();
288 if (rz->
Amin() < 0.0)
290 std::ostringstream message;
291 message <<
"Illegal input parameters - " <<
GetName() << std::endl
292 <<
" All R values must be >= 0 !";
353 int oldprc = os.precision(16);
354 os <<
"-----------------------------------------------------------\n"
355 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
356 <<
" ===================================================\n"
357 <<
" Solid type: UPolycone\n"
363 os <<
" number of Z planes: " << numPlanes <<
"\n"
365 for (i = 0; i < numPlanes; i++)
367 os <<
" Z plane " << i <<
": "
370 os <<
" Tangent distances to inner surface (Rmin): \n";
371 for (i = 0; i < numPlanes; i++)
373 os <<
" Z plane " << i <<
": "
376 os <<
" Tangent distances to outer surface (Rmax): \n";
377 for (i = 0; i < numPlanes; i++)
379 os <<
" Z plane " << i <<
": "
382 os <<
"-----------------------------------------------------------\n";
383 os.precision(oldprc);
395 double rMinPlus, rMaxPlus , rMinMinus, rMaxMinus;
397 static const double halfTolerance =
fgTolerance * 0.5;
399 double r2 = p.
x() * p.
x() + p.
y() * p.
y();
405 rMaxPlus = tubs->
GetRMax()+ halfTolerance;
406 rMinPlus = tubs->
GetRMin() + halfTolerance;
407 rMinMinus = tubs->
GetRMin() - halfTolerance;
408 rMaxMinus = tubs->
GetRMax() - halfTolerance;
420 double ratio = (ps.z() + dz) / (2 * dz);
421 rMinPlus = rMin1 + (rMin2 - rMin1) * ratio + halfTolerance;
422 rMaxPlus = rMax1 + (rMax2 - rMax1) * ratio + halfTolerance;
423 rMinMinus = rMinPlus - 2*halfTolerance;
424 rMaxMinus = rMaxPlus - 2*halfTolerance;
434 if (rMinMinus <= 0) { pos =
eSurface; }
436 && ((ps.z() < -dz+halfTolerance) || (ps.z() > dz-halfTolerance)))
443 if ( (r2 < rMinMinus*rMinMinus) || (r2 > rMaxPlus*rMaxPlus) )
447 if ( (r2 < rMinPlus*rMinPlus) || (r2 > rMaxMinus*rMaxMinus) )
454 if ( (ps.z() < -dz+halfTolerance) || (ps.z() > dz-halfTolerance) )
460 double phi = std::atan2(p.
y(), p.
x());
469 if ( (ps.z() < -dz+halfTolerance) || (ps.z() > dz-halfTolerance) )
473 if ( (r2 < rMinPlus*rMinPlus) || (r2 > rMaxMinus*rMaxMinus) )
501 static const double htolerance = 0.5 *
fgTolerance;
510 if (index > 0 && p.
z() -
fZs[index] < htolerance)
512 nextSection = index - 1;
517 nextSection = index + 1;
561 pb = p + idistance * v;
564 int increment = (v.z() > 0) ? 1 : -1;
571 pb.z() -= section.
shift;
576 pb.z() += section.
shift;
585 UVector3&
n,
bool& convex,
double )
const
600 if ( indexLow != indexHigh )
618 double totalDistance = 0;
619 int increment = (v.
z() > 0) ? 1 : -1;
620 bool convexloc =
true;
629 if (totalDistance != 0||(istep < 2))
631 pn = p + (totalDistance ) * v;
655 if( (index == 0) && ( increment > 0 ) ) { index1 += increment; }
656 if( (index ==
fMaxSection) && (increment<0) ) { index1 += increment;}
658 UVector3 pte = p+(totalDistance+distance)*v;
660 pte.
z() -= section1.shift;
661 if (section1.solid->Inside(pte) ==
eOutside)
669 if((convexloc) && (section.
convex))
679 totalDistance += distance;
687 pn = p + (totalDistance) * v;
691 double dz1 = std::fabs(pn.
z()-
fZs[index]);
692 double dz2 = std::fabs(pn.
z()-
fZs[index+1]);
693 if(dz1 < halfTolerance) { convex=
false; }
694 if(dz2 < halfTolerance) { convex=
false; }
700 if (std::fabs(pn.
z()-
fZs[
fMaxSection]) < halfTolerance) { convex=
false; }
706 if (std::fabs(pn.
z()-
fZs[1]) < halfTolerance) { convex=
false;}
713 return totalDistance;
722 double rho=std::sqrt(p.
x()*p.
x()+p.
y()*p.
y());
724 double safeDown = p.
z()-
fZs[index];
725 double safeUp =
fZs[index+1]-p.
z();
730 if (minSafety < 1e-6)
return 0;
734 double dz1 =
fZs[i] - p.
z();
735 double dz2 =
fZs[i+1] - p.
z();
737 if (safeR < 0.) { safeUp=dz1;
break; }
738 if (dz1 < dz2) { safeR = std::sqrt(safeR*safeR+dz1*dz1); }
739 else {safeR = std::sqrt(safeR*safeR+dz1*dz1); }
740 if (safeR < dz1) { safeUp = safeR;
break; }
741 if (safeR < dz2) { safeUp = safeR;
break; }
747 for (
int i = index - 1; i >= 0; --i)
749 double dz1 = p.
z()-
fZs[i+1];
750 double dz2 = p.
z()-
fZs[i];
752 if (safeR < 0.) { safeDown=dz1;
break; }
753 if(dz1 < dz2) { safeR = std::sqrt(safeR*safeR+dz1*dz1); }
754 else { safeR = std::sqrt(safeR*safeR+dz1*dz1); }
755 if (safeR < dz1) { safeDown = safeR;
break; }
756 if (safeR < dz2) { safeDown = safeR;
break; }
760 if (safeUp < minSafety) minSafety=safeUp;
761 if (safeDown < minSafety) minSafety=safeDown;
775 double zbase =
fZs[index + 1];
778 double dz =
fZs[i] - zbase;
779 if (dz >= minSafety)
break;
781 if (safety < minSafety) minSafety = safety;
784 zbase =
fZs[index - 1];
785 for (
int i = index - 1; i >= 0; --i)
787 double dz = zbase -
fZs[i];
788 if (dz >= minSafety)
break;
790 if (safety < minSafety) minSafety = safety;
803 if (index > 0 && p.
z() -
fZs[index] < htolerance)
805 nextSection = index - 1;
810 nextSection = index + 1;
850 if (pos !=
eSurface) index = nextSection;
863 aMin.
Set(-r, -r,
fZs.front());
864 aMax.
Set(r, r,
fZs.back());
892 double Area = 0, totArea = 0;
897 std::vector<double> areas;
898 std::vector<UVector3> points;
903 for (i = 0; i < numPlanes - 1; i++)
928 areas.push_back(Area);
935 totArea += (areas[0] + areas[numPlanes]);
952 double fRmin2,
double fRmax2,
953 double zOne,
double zTwo,
954 double& totArea)
const
958 double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
959 double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
960 double fDz = (zTwo - zOne) / 2., afDz = std::fabs(fDz);
963 rone = (fRmax1 - fRmax2) / (2.*fDz);
964 rtwo = (fRmin1 - fRmin2) / (2.*fDz);
965 if (fRmax1 == fRmax2)
971 qone = fDz * (fRmax1 + fRmax2) / (fRmax1 - fRmax2);
973 if (fRmin1 == fRmin2)
979 qtwo = fDz * (fRmin1 + fRmin2) / (fRmin1 - fRmin2);
983 Afive = fDz * (fRmax1 - fRmin1 + fRmax2 - fRmin2);
984 totArea = Aone + Atwo + 2.*Afive;
987 cosu = std::cos(phi);
988 sinu = std::sin(phi);
996 if ((chose >= 0) && (chose < Aone))
998 if (fRmax1 != fRmax2)
1001 point =
UVector3(rone * cosu * (qone - zRand),
1002 rone * sinu * (qone - zRand), zRand);
1006 point =
UVector3(fRmax1 * cosu, fRmax1 * sinu,
1011 else if (chose >= Aone && chose < Aone + Atwo)
1013 if (fRmin1 != fRmin2)
1016 point =
UVector3(rtwo * cosu * (qtwo - zRand),
1017 rtwo * sinu * (qtwo - zRand), zRand);
1022 point =
UVector3(fRmin1 * cosu, fRmin1 * sinu,
1026 else if ((chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive))
1029 rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
1030 rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
1032 point =
UVector3(rRand1 * std::cos(startPhi),
1033 rRand1 * std::sin(startPhi), zRand);
1038 rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
1039 rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
1042 rRand1 * std::sin(
endPhi), zRand);
1046 return point + offset;
1056 double zOne,
double zTwo,
1057 double& totArea)
const
1059 double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1060 aOne, aTwo, aFou, rRand, fDz, fSPhi, fDPhi;
1061 fDz = std::fabs(0.5 * (zTwo - zOne));
1065 aOne = 2.*fDz * fDPhi * fRMax;
1066 aTwo = 2.*fDz * fDPhi * fRMin;
1067 aFou = 2.*fDz * (fRMax - fRMin);
1068 totArea = aOne + aTwo + 2.*aFou;
1070 cosphi = std::cos(phi);
1071 sinphi = std::sin(phi);
1078 if ((chose >= 0) && (chose < aOne))
1080 xRand = fRMax * cosphi;
1081 yRand = fRMax * sinphi;
1083 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1085 else if ((chose >= aOne) && (chose < aOne + aTwo))
1087 xRand = fRMin * cosphi;
1088 yRand = fRMin * sinphi;
1090 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1092 else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aFou))
1094 xRand = rRand * std::cos(fSPhi + fDPhi);
1095 yRand = rRand * std::sin(fSPhi + fDPhi);
1097 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1102 xRand = rRand * std::cos(fSPhi + fDPhi);
1103 yRand = rRand * std::sin(fSPhi + fDPhi);
1105 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1115 double fRMin2,
double fRMax2,
1118 double xRand, yRand, phi, cosphi, sinphi, rRand1, rRand2, A1, Atot, rCh;
1120 cosphi = std::cos(phi);
1121 sinphi = std::sin(phi);
1123 if (fRMin1 == fRMin2)
1131 A1 = std::fabs(fRMin2 * fRMin2 - fRMin1 * fRMin1);
1133 if (fRMax1 == fRMax2)
1141 Atot = A1 + std::fabs(fRMax2 * fRMax2 - fRMax1 * fRMax1);
1150 xRand = rRand1 * cosphi;
1151 yRand = rRand1 * sinphi;
1153 return UVector3(xRand, yRand, zOne);
1163 double fRMin2,
double fRMax2,
1164 double zOne,
double zTwo,
1165 double& totArea)
const
1171 if ((fRMin1 == fRMin2) && (fRMax1 == fRMax2))
1175 return GetPointOnCone(fRMin1, fRMax1, fRMin2, fRMax2, zOne, zTwo, totArea);
1184 double Area = 0, totArea = 0, Achose1 = 0, Achose2 = 0, phi, cosphi, sinphi, rRand;
1189 cosphi = std::cos(phi);
1190 sinphi = std::sin(phi);
1196 std::vector<double> areas;
1197 std::vector<UVector3> points;
1202 for (i = 0; i < numPlanes - 1; i++)
1227 areas.push_back(Area);
1234 totArea += (areas[0] + areas[numPlanes]);
1237 if ((chose >= 0.) && (chose < areas[0]))
1239 return UVector3(rRand * cosphi, rRand * sinphi,
1243 for (i = 0; i < numPlanes - 1; i++)
1245 Achose1 += areas[i];
1246 Achose2 = (Achose1 + areas[i + 1]);
1247 if (chose >= Achose1 && chose < Achose2)
1262 return UVector3(rRand * cosphi, rRand * sinphi,
1272 : fStartAngle(0.), fOpeningAngle(0.), fNumZPlanes(0),
1273 fZValues(0), Rmin(0), Rmax(0)
1303 if (&right ==
this)
return *
this;
1341 if (
this == &source)
return *
this;
1390 bool isConvertible =
true;
1391 double Zmax = rz->
Bmax();
1396 std::vector<double> Z;
1397 std::vector<double> Rmin;
1398 std::vector<double> Rmax;
1400 int countPlanes = 1;
1407 double Zprev = Z[0];
1414 else if (Zprev ==
corners[numPlanes - 1].z)
1416 Rmin.push_back(
corners[numPlanes - 1].r);
1418 icurl = numPlanes - 1;
1428 int inextr = 0, inextl = 0;
1429 for (
int i = 0; i < numPlanes - 2; i++)
1432 inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1452 Rmin.push_back(
corners[inextl].r);
1453 Rmax.push_back(
corners[icurr].r);
1457 Rmin.push_back(
corners[inextl].r);
1466 Rmin.push_back(
corners[icurl].r);
1467 Rmax.push_back(
corners[icurr].r);
1471 Rmin.push_back(
corners[icurl].r);
1478 isConvertible =
false;
1481 icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1489 icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1491 Rmin.push_back(
corners[inextl].r);
1492 Rmax.push_back(
corners[inextr].r);
1496 Z.push_back(Zright);
1505 Rmax.push_back(
corners[inextr].r);
1506 Rmin.push_back (
corners[icurr].r);
1510 Rmin.push_back (
corners[icurl].r + (Zright-
corners[icurl].z)/difZl
1512 Rmax.push_back(
corners[inextr].r);
1520 Rmax.push_back(
corners[inextr].r);
1521 Rmin.push_back(
corners[icurr].r);
1525 Rmax.push_back(
corners[inextr].r);
1526 Rmin.push_back(
corners[icurl].r + (Zright -
corners[icurl].z) / difZl
1533 isConvertible =
false;
1544 inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1548 Rmax.push_back(
corners[inextr].r);
1549 Rmin.push_back(
corners[inextl].r);
1553 Rmax.push_back(
corners[inextr].r);
1554 Rmin.push_back(
corners[inextl].r);
1566 for (
int j = 0; j < countPlanes; j++)
1579 std::ostringstream message;
1580 message <<
"Polycone " <<
GetName() << std::endl
1581 <<
"cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1591 for (
int j = 0; j < numPlanes; j++)
1601 return isConvertible;
1616 double* Z, *R1, *R2;
1618 Z =
new double[num];
1619 R1 =
new double[num];
1620 R2 =
new double[num];
1621 for (
int i = 0; i < num; i++)
double SafetyFromOutside(const UVector3 &p) const
double DistanceToIn(const UVector3 &p, const UVector3 &v, double aPstep=UUtils::kInfinity) const
std::vector< double > Rmin
static double frTolerance
std::vector< double > fZValues
UVector3 GetPointOnTubs(double fRMin, double fRMax, double zOne, double zTwo, double &totArea) const
bool NormalSection(int index, const UVector3 &p, UVector3 &n) const
int GetSection(double z) const
const std::string & GetName() const
static double fgTolerance
void CopyStuff(const UPolycone &source)
std::vector< UPolyconeSection > fSections
UVector3 GetPointOnCut(double fRMin1, double fRMax1, double fRMin2, double fRMax2, double zOne, double zTwo, double &totArea) const
double SafetyFromOutsideSection(int index, const double rho, const UVector3 &p) const
std::ostream & StreamInfo(std::ostream &os) const
UPolyconeHistorical & operator=(const UPolyconeHistorical &right)
virtual bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const =0
virtual EnumInside Inside(const UVector3 &aPoint) const =0
VUSolid::EnumInside Inside(const UVector3 &p) const
std::vector< double > Rmax
bool Normal(const UVector3 &aPoint, UVector3 &aNormal) const
static int BinarySearch(const std::vector< T > &vec, T value)
static const double kInfinity
virtual double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const =0
static const double kTwoPi
UVector3 GetPointOnSurface() const
void Init(double phiStart, double phiTotal, int numZPlanes, const double zPlane[], const double rInner[], const double rOuter[])
UPolyconeHistorical * fOriginalParameters
UPolyconeSideRZ * corners
UGeometryType GetEntityType() const
double DistanceToIn(const UVector3 &aPoint, const UVector3 &aDirection, double aPstep=UUtils::kInfinity) const
static const G4double minSafety
double GetZHalfLength() const
void Exception(const char *originOfException, const char *exceptionCode, UExceptionSeverity severity, int level, const char *description)
void SetOriginalParameters()
VUSolid::EnumInside InsideSection(int index, const UVector3 &p) const
void Set(double xx, double yy, double zz)
double SafetyFromOutside(const UVector3 &aPoint, bool aAccurate=false) const
double SafetyFromInside(const UVector3 &aPoint, bool aAccurate=false) const
void Set(double dx, double dy, double dz)
std::string UGeometryType
EnumInside Inside(const UVector3 &aPoint) const
UVector3 GetPointOnRing(double fRMin, double fRMax, double fRMin2, double fRMax2, double zOne) const
UVector3 GetPointOnCone(double fRmin1, double fRmax1, double fRmin2, double fRmax2, double zOne, double zTwo, double &totArea) const
virtual double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const =0
UPolycone & operator=(const UPolycone &source)
double SafetyFromInsideSection(int index, const double rho, const UVector3 &p) const
std::vector< double > fZs
UEnclosingCylinder * enclosingCylinder
void Extent(UVector3 &aMin, UVector3 &aMax) const
double Random(double min=0.0, double max=1.0)
double DistanceToOut(const UVector3 &aPoint, const UVector3 &aDirection, UVector3 &aNormalVector, bool &aConvex, double aPstep=UUtils::kInfinity) const
UPolycone(const std::string &name)
static const G4double pos
virtual double Capacity()=0