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" ;
76 std::cout <<
"INFO: Converting polycone " <<
GetName() << std::endl
77 <<
"to optimized polycone with (Rmin,Rmaz,Z) parameters !"
84 for (
int i = 0; i < num; i++)
92 Init(phiStart, phiTotal, num, Z, R1, R2);
108 const double zPlane[],
109 const double rInner[],
110 const double rOuter[])
129 endPhi = phiStart+phiTotal;
143 double RMaxextent=rOuter[0];
144 for (
int j=1; j < numZPlanes; j++)
146 if (rOuter[j] > RMaxextent) RMaxextent=rOuter[j];
147 if (rInner[j]>rOuter[j])
149 std::ostringstream message;
150 message <<
"Cannot create Polycone with rInner > rOuter for the same Z"
152 <<
" rInner > rOuter for the same Z !" << std::endl
153 <<
" rMin[" << j <<
"] = " << rInner[j]
154 <<
" -- rMax[" << j <<
"] = " << rOuter[j];
160 double prevZ = 0, prevRmax = 0, prevRmin = 0;
162 if (zPlane[1] < zPlane[0])dirZ = -1;
166 for (i = 0; i < numZPlanes; i++)
168 if ((i < numZPlanes - 1) && (zPlane[i] == zPlane[i + 1]))
170 if ((rInner[i] > rOuter[i + 1])
171 || (rInner[i + 1] > rOuter[i]))
174 std::ostringstream message;
175 message <<
"Cannot create a Polycone with no contiguous segments."
177 <<
"Segments are not contiguous !" << std::endl
178 <<
" rMin[" << i <<
"] = " << rInner[i]
179 <<
" -- rMax[" << i + 1 <<
"] = " << rOuter[i + 1] << std::endl
180 <<
" rMin[" << i + 1 <<
"] = " << rInner[i + 1]
181 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
189 double rMin = rInner[i];
190 double rMax = rOuter[i];
191 double z = zPlane[i];
199 std::ostringstream message;
200 message <<
"Cannot create a Polycone with different Z directions.Use GenericPolycone."
202 <<
"ZPlane is changing direction !" << std::endl
203 <<
" zPlane[0] = " << zPlane[0]
204 <<
" -- zPlane[1] = " << zPlane[1] << std::endl
205 <<
" zPlane[" << i - 1 <<
"] = " << zPlane[i - 1]
206 <<
" -- rPlane[" << i <<
"] = " << zPlane[i];
214 double dz = (z - prevZ) / 2;
216 bool tubular = (rMin == prevRmin && prevRmax == rMax);
222 solid =
new UTubs(
"", rMin, rMax, dz, phiStart, phiTotal);
226 solid =
new UCons(
"", prevRmin, prevRmax, rMin, rMax, dz, phiStart, phiTotal);
236 int zi =
fZs.size() - 1;
237 double shift =
fZs[zi - 1] + 0.5 * (
fZs[zi] -
fZs[zi - 1]);
240 section.
shift = shift;
242 section.
solid = solid;
245 if (rMax < RMaxextent) { section.
convex =
false;}
246 else { section.
convex =
true;}
250 if ((rMax<prevRmax)||(rMax < RMaxextent)||(prevRmax < RMaxextent))
251 { section.
convex =
false;}
262 else fZs.push_back(z);
280 double mxy = rz->
Amax();
283 double r = rz->
Amax();
287 if (rz->
Amin() < 0.0)
289 std::ostringstream message;
290 message <<
"Illegal input parameters - " <<
GetName() << std::endl
291 <<
" All R values must be >= 0 !";
363 int oldprc = os.precision(16);
364 os <<
"-----------------------------------------------------------\n"
365 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
366 <<
" ===================================================\n"
367 <<
" Solid type: UPolycone3\n"
373 os <<
" number of Z planes: " << numPlanes <<
"\n"
375 for (i = 0; i < numPlanes; i++)
377 os <<
" Z plane " << i <<
": "
380 os <<
" Tangent distances to inner surface (Rmin): \n";
381 for (i = 0; i < numPlanes; i++)
383 os <<
" Z plane " << i <<
": "
386 os <<
" Tangent distances to outer surface (Rmax): \n";
387 for (i = 0; i < numPlanes; i++)
389 os <<
" Z plane " << i <<
": "
392 os <<
"-----------------------------------------------------------\n";
393 os.precision(oldprc);
406 double rMinPlus, rMaxPlus;
413 rMinPlus = tubs->
GetRMin() + halfTolerance;
414 rMaxPlus = tubs->
GetRMax() + halfTolerance;
427 double ratio = (ps.
z + dz) / (2 * dz);
428 rMinPlus = rMin1 + (rMin2 - rMin1) * ratio + halfTolerance;
429 rMaxPlus = rMax1 + (rMax2 - rMax1) * ratio + halfTolerance;
435 double r2 = p.
x * p.
x + p.
y * p.
y;
437 if (r2 < rMinMinus * rMinMinus || r2 > rMaxPlus * rMaxPlus)
return eOutside;
438 if (r2 < rMinPlus * rMinPlus || r2 > rMaxMinus * rMaxMinus)
return eSurface;
442 if (ps.
z < -dz + halfTolerance || ps.
z > dz - halfTolerance)
447 if (r2 < 1e-10)
return eInside;
449 double phi = std::atan2(p.
y, p.
x);
456 if (ps.
z < -dz + halfTolerance || ps.
z > dz - halfTolerance)
477 static const double htolerance = 0.5 *
fgTolerance;
486 if (index > 0 && p.
z -
fZs[index] < htolerance)
488 nextSection = index - 1;
493 nextSection = index + 1;
541 pb = p + idistance * v;
544 int increment = (v.z > 0) ? 1 : -1;
564 UVector3&
n,
bool& convex,
double )
const
579 if ( indexLow != indexHigh )
584 else{index=indexLow;}
587 else{index=indexLow;}
588 double totalDistance = 0;
589 int increment = (v.
z > 0) ? 1 : -1;
598 if (totalDistance != 0)
600 pn = p + (totalDistance ) * v;
607 else pn.
z -= section.
shift;
612 {
int index1 = index;
613 if(( index > 0) && ( index <
fMaxSection )){index1 += increment;}
615 if((index == 0) && ( increment > 0 ))index1 += increment;
616 if((index ==
fMaxSection) && (increment<0 ))index1 += increment;
618 UVector3 pte = p+(totalDistance+distance)*v;
620 pte.
z -= section1.shift;
621 if (section1.solid->Inside(pte) ==
eOutside)
628 if((convexloc) && (section.
convex)) {convexloc=
true;}
629 else{convexloc=
false;}
632 totalDistance += distance;
642 pn = p + (totalDistance) * v;
646 double dz1 = std::fabs(pn.
z-
fZs[index]);
647 double dz2 = std::fabs(pn.
z-
fZs[index+1]);
648 if(dz1 < halfTolerance)convex=
false;
649 if(dz2 < halfTolerance)convex=
false;
655 if(index<=0){
if(std::fabs(pn.
z-
fZs[1]) < halfTolerance)convex=
false;}
659 return totalDistance;
667 double rho=std::sqrt(p.
x*p.
x+p.
y*p.
y);
669 double safeDown = p.
z-
fZs[index];
670 double safeUp =
fZs[index+1]-p.
z;
672 double minSafety =safeR;
675 if (minSafety < 1e-6)
return 0;
679 double dz1 =
fZs[i] - p.
z;
680 double dz2 =
fZs[i+1] - p.
z;
682 if (safeR < 0.) { safeUp=dz1;
break; }
683 if (dz1 < dz2) { safeR = std::sqrt(safeR*safeR+dz1*dz1); }
684 else {safeR = std::sqrt(safeR*safeR+dz1*dz1); }
685 if (safeR < dz1) { safeUp = safeR;
break; }
686 if (safeR < dz2) { safeUp = safeR;
break; }
692 for (
int i = index - 1; i >= 0; --i)
694 double dz1 = p.
z-
fZs[i+1];
695 double dz2 = p.
z-
fZs[i];
697 if (safeR < 0.) { safeDown=dz1;
break; }
698 if(dz1 < dz2) { safeR = std::sqrt(safeR*safeR+dz1*dz1); }
699 else { safeR = std::sqrt(safeR*safeR+dz1*dz1); }
700 if (safeR < dz1) { safeDown = safeR;
break; }
701 if (safeR < dz2) { safeDown = safeR;
break; }
705 if (safeUp < minSafety) minSafety=safeUp;
706 if (safeDown < minSafety) minSafety=safeDown;
718 if (minSafety < 1e-6)
return minSafety;
720 double zbase =
fZs[index + 1];
723 double dz =
fZs[i] - zbase;
724 if (dz >= minSafety)
break;
726 if (safety < minSafety) minSafety = safety;
729 zbase =
fZs[index - 1];
730 for (
int i = index - 1; i >= 0; --i)
732 double dz = zbase -
fZs[i];
733 if (dz >= minSafety)
break;
735 if (safety < minSafety) minSafety = safety;
748 if (index > 0 && p.
z -
fZs[index] < htolerance)
750 nextSection = index - 1;
755 nextSection = index + 1;
814 if (pos !=
eSurface) index = nextSection;
828 aMin.
Set(-r, -r,
fZs.front());
829 aMax.
Set(r, r,
fZs.back());
857 double Area = 0, totArea = 0;
862 std::vector<double> areas;
863 std::vector<UVector3> points;
868 for (i = 0; i < numPlanes - 1; i++)
893 areas.push_back(Area);
900 totArea += (areas[0] + areas[numPlanes]);
917 double fRmin2,
double fRmax2,
918 double zOne,
double zTwo,
919 double& totArea)
const
923 double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
924 double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
925 double fDz = (zTwo - zOne) / 2., afDz = std::fabs(fDz);
928 rone = (fRmax1 - fRmax2) / (2.*fDz);
929 rtwo = (fRmin1 - fRmin2) / (2.*fDz);
930 if (fRmax1 == fRmax2)
936 qone = fDz * (fRmax1 + fRmax2) / (fRmax1 - fRmax2);
938 if (fRmin1 == fRmin2)
944 qtwo = fDz * (fRmin1 + fRmin2) / (fRmin1 - fRmin2);
948 Afive = fDz * (fRmax1 - fRmin1 + fRmax2 - fRmin2);
949 totArea = Aone + Atwo + 2.*Afive;
952 cosu = std::cos(phi);
953 sinu = std::sin(phi);
961 if ((chose >= 0) && (chose < Aone))
963 if (fRmax1 != fRmax2)
966 point =
UVector3(rone * cosu * (qone - zRand),
967 rone * sinu * (qone - zRand), zRand);
971 point =
UVector3(fRmax1 * cosu, fRmax1 * sinu,
976 else if (chose >= Aone && chose < Aone + Atwo)
978 if (fRmin1 != fRmin2)
981 point =
UVector3(rtwo * cosu * (qtwo - zRand),
982 rtwo * sinu * (qtwo - zRand), zRand);
987 point =
UVector3(fRmin1 * cosu, fRmin1 * sinu,
991 else if ((chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive))
994 rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
995 rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
997 point =
UVector3(rRand1 * std::cos(startPhi),
998 rRand1 * std::sin(startPhi), zRand);
1003 rmin = fRmin2 - ((zRand - fDz) / (2.*fDz)) * (fRmin1 - fRmin2);
1004 rmax = fRmax2 - ((zRand - fDz) / (2.*fDz)) * (fRmax1 - fRmax2);
1007 rRand1 * std::sin(
endPhi), zRand);
1011 return point + offset;
1021 double zOne,
double zTwo,
1022 double& totArea)
const
1024 double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1025 aOne, aTwo, aFou, rRand, fDz, fSPhi, fDPhi;
1026 fDz = std::fabs(0.5 * (zTwo - zOne));
1030 aOne = 2.*fDz * fDPhi * fRMax;
1031 aTwo = 2.*fDz * fDPhi * fRMin;
1032 aFou = 2.*fDz * (fRMax - fRMin);
1033 totArea = aOne + aTwo + 2.*aFou;
1035 cosphi = std::cos(phi);
1036 sinphi = std::sin(phi);
1043 if ((chose >= 0) && (chose < aOne))
1045 xRand = fRMax * cosphi;
1046 yRand = fRMax * sinphi;
1048 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1050 else if ((chose >= aOne) && (chose < aOne + aTwo))
1052 xRand = fRMin * cosphi;
1053 yRand = fRMin * sinphi;
1055 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1057 else if ((chose >= aOne + aTwo) && (chose < aOne + aTwo + aFou))
1059 xRand = rRand * std::cos(fSPhi + fDPhi);
1060 yRand = rRand * std::sin(fSPhi + fDPhi);
1062 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1067 xRand = rRand * std::cos(fSPhi + fDPhi);
1068 yRand = rRand * std::sin(fSPhi + fDPhi);
1070 return UVector3(xRand, yRand, zRand + 0.5 * (zTwo + zOne));
1080 double fRMin2,
double fRMax2,
1083 double xRand, yRand, phi, cosphi, sinphi, rRand1, rRand2, A1, Atot, rCh;
1085 cosphi = std::cos(phi);
1086 sinphi = std::sin(phi);
1088 if (fRMin1 == fRMin2)
1096 A1 = std::fabs(fRMin2 * fRMin2 - fRMin1 * fRMin1);
1098 if (fRMax1 == fRMax2)
1106 Atot = A1 + std::fabs(fRMax2 * fRMax2 - fRMax1 * fRMax1);
1115 xRand = rRand1 * cosphi;
1116 yRand = rRand1 * sinphi;
1118 return UVector3(xRand, yRand, zOne);
1128 double fRMin2,
double fRMax2,
1129 double zOne,
double zTwo,
1130 double& totArea)
const
1136 if ((fRMin1 == fRMin2) && (fRMax1 == fRMax2))
1140 return GetPointOnCone(fRMin1, fRMax1, fRMin2, fRMax2, zOne, zTwo, totArea);
1149 double Area = 0, totArea = 0, Achose1 = 0, Achose2 = 0, phi, cosphi, sinphi, rRand;
1154 cosphi = std::cos(phi);
1155 sinphi = std::sin(phi);
1161 std::vector<double> areas;
1162 std::vector<UVector3> points;
1167 for (i = 0; i < numPlanes - 1; i++)
1192 areas.push_back(Area);
1199 totArea += (areas[0] + areas[numPlanes]);
1202 if ((chose >= 0.) && (chose < areas[0]))
1204 return UVector3(rRand * cosphi, rRand * sinphi,
1208 for (i = 0; i < numPlanes - 1; i++)
1210 Achose1 += areas[i];
1211 Achose2 = (Achose1 + areas[i + 1]);
1212 if (chose >= Achose1 && chose < Achose2)
1227 return UVector3(rRand * cosphi, rRand * sinphi,
1237 : fStartAngle(0.), fOpeningAngle(0.), fNumZPlanes(0),
1238 fZValues(0), Rmin(0), Rmax(0)
1268 if (&right ==
this)
return *
this;
1308 if (
this == &source)
return *
this;
1357 bool isConvertible =
true;
1358 double Zmax = rz->
Bmax();
1363 std::vector<double> Z;
1364 std::vector<double> Rmin;
1365 std::vector<double> Rmax;
1367 int countPlanes = 1;
1374 double Zprev = Z[0];
1381 else if (Zprev ==
corners[numPlanes - 1].z)
1383 Rmin.push_back(
corners[numPlanes - 1].r);
1385 icurl = numPlanes - 1;
1395 int inextr = 0, inextl = 0;
1396 for (
int i = 0; i < numPlanes - 2; i++)
1399 inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1419 Rmin.push_back(
corners[icurl].r);
1420 Rmax.push_back(Rmax[countPlanes - 2]);
1421 Rmax[countPlanes - 2] =
corners[icurl].
r;
1425 Rmin.push_back(
corners[inextl].r);
1426 Rmax.push_back(
corners[icurl].r);
1431 Rmin.push_back(
corners[inextl].r);
1432 Rmax.push_back(
corners[icurr].r + (Zleft -
corners[icurr].z) / difZr
1437 isConvertible =
false;
1440 icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1448 icurl = (icurl == 0) ? numPlanes - 1 : icurl - 1;
1450 Rmin.push_back(
corners[inextl].r);
1451 Rmax.push_back(
corners[inextr].r);
1455 Z.push_back(Zright);
1464 Rmin.push_back(
corners[icurr].r);
1465 Rmax.push_back(
corners[inextr].r);
1469 Rmin.push_back(
corners[inextr].r);
1470 Rmax.push_back(
corners[icurr].r);
1471 Rmax[countPlanes - 2] =
corners[inextr].
r;
1479 Rmax.push_back(
corners[inextr].r);
1480 Rmin.push_back(
corners[icurr].r);
1484 Rmax.push_back(
corners[inextr].r);
1485 Rmin.push_back(
corners[icurl].r + (Zright -
corners[icurl].z) / difZl
1492 isConvertible =
false;
1503 inextl = (icurl <= 0) ? numPlanes - 1 : icurl - 1;
1507 Rmax.push_back(
corners[inextr].r);
1508 Rmin.push_back(
corners[inextl].r);
1512 Rmax.push_back(
corners[inextr].r);
1513 Rmin.push_back(
corners[inextl].r);
1525 for (
int j = 0; j < countPlanes; j++)
1538 std::ostringstream message;
1539 message <<
"Polycone " <<
GetName() << std::endl
1540 <<
"cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1542 Warning, 1,
"can not convert");
1550 for (
int j = 0; j < numPlanes; j++)
1560 return isConvertible;
1575 double* Z, *R1, *R2;
1577 Z =
new double[num];
1578 R1 =
new double[num];
1579 R2 =
new double[num];
1580 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
double GetZHalfLength() const
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)
void Exception(const char *originOfException, const char *exceptionCode, ExceptionSeverity severity, int level, const char *description)
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