42 #if !defined(G4GEOM_USE_UPOLYCONE)
53 using namespace CLHEP;
80 for (i=0; i<numZPlanes; i++)
82 if(rInner[i]>rOuter[i])
85 std::ostringstream message;
86 message <<
"Cannot create a Polycone with rInner > rOuter for the same Z"
88 <<
" rInner > rOuter for the same Z !" <<
G4endl
89 <<
" rMin[" << i <<
"] = " << rInner[i]
90 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
91 G4Exception(
"G4Polycone::G4Polycone()",
"GeomSolids0002",
94 if (( i < numZPlanes-1) && ( zPlane[i] == zPlane[i+1] ))
96 if( (rInner[i] > rOuter[i+1])
97 ||(rInner[i+1] > rOuter[i]) )
100 std::ostringstream message;
101 message <<
"Cannot create a Polycone with no contiguous segments."
103 <<
" Segments are not contiguous !" <<
G4endl
104 <<
" rMin[" << i <<
"] = " << rInner[i]
105 <<
" -- rMax[" << i+1 <<
"] = " << rOuter[i+1] <<
G4endl
106 <<
" rMin[" << i+1 <<
"] = " << rInner[i+1]
107 <<
" -- rMax[" << i <<
"] = " << rOuter[i];
108 G4Exception(
"G4Polycone::G4Polycone()",
"GeomSolids0002",
126 Create( phiStart, phiTotal, rz );
146 Create( phiStart, phiTotal, rz );
155 std::ostringstream message;
156 message <<
"Polycone " <<
GetName() <<
"cannot be converted" <<
G4endl
157 <<
"to Polycone with (Rmin,Rmaz,Z) parameters!";
158 G4Exception(
"G4Polycone::G4Polycone()",
"GeomSolids0002",
164 <<
"to optimized polycone with (Rmin,Rmaz,Z) parameters !"
183 if (rz->
Amin() < 0.0)
185 std::ostringstream message;
186 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
187 <<
" All R values must be >= 0 !";
188 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
198 std::ostringstream message;
199 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
200 <<
" R/Z cross section is zero or near zero: " << rzArea;
201 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
208 std::ostringstream message;
209 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
210 <<
" Too few unique R/Z values !";
211 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
217 std::ostringstream message;
218 message <<
"Illegal input parameters - " <<
GetName() <<
G4endl
219 <<
" R/Z segments cross !";
220 G4Exception(
"G4Polycone::Create()",
"GeomSolids0002",
230 if (phiTotal <= 0 || phiTotal > twopi-1E-10)
246 endPhi = phiStart+phiTotal;
264 next->
r = iterRZ.
GetA();
265 next->
z = iterRZ.
GetB();
266 }
while( ++next, iterRZ.
Next() );
298 if (corner->
z > next->
z)
315 }
while( prev=corner, corner=next, corner >
corners );
344 :
G4VCSGfaceted(a), startPhi(0.), endPhi(0.), phiIsOpen(false),
345 numCorner(0), corners(0),original_parameters(0),
377 if (
this == &source)
return *
this;
546 G4int oldprc = os.precision(16);
547 os <<
"-----------------------------------------------------------\n"
548 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
549 <<
" ===================================================\n"
550 <<
" Solid type: G4Polycone\n"
553 <<
" ending phi angle : " <<
endPhi/
degree <<
" degrees \n";
557 os <<
" number of Z planes: " << numPlanes <<
"\n"
559 for (i=0; i<numPlanes; i++)
561 os <<
" Z plane " << i <<
": "
564 os <<
" Tangent distances to inner surface (Rmin): \n";
565 for (i=0; i<numPlanes; i++)
567 os <<
" Z plane " << i <<
": "
570 os <<
" Tangent distances to outer surface (Rmax): \n";
571 for (i=0; i<numPlanes; i++)
573 os <<
" Z plane " << i <<
": "
577 os <<
" number of RZ points: " <<
numCorner <<
"\n"
578 <<
" RZ values (corners): \n";
584 os <<
"-----------------------------------------------------------\n";
585 os.precision(oldprc);
603 G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
604 G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
605 G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
608 rone = (fRmax1-fRmax2)/(2.*fDz);
609 rtwo = (fRmin1-fRmin2)/(2.*fDz);
610 if(fRmax1==fRmax2){qone=0.;}
612 qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
614 if(fRmin1==fRmin2){qtwo=0.;}
616 qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
618 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(
sqr(fRmin1-fRmin2)+
sqr(zTwo-zOne));
619 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(
sqr(fRmax1-fRmax2)+
sqr(zTwo-zOne));
620 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
621 totArea = Aone+Atwo+2.*Afive;
624 cosu = std::cos(phi);
625 sinu = std::sin(phi);
628 if( (startPhi == 0) && (
endPhi == twopi) ) { Afive = 0; }
630 if( (chose >= 0) && (chose < Aone) )
636 rone*sinu*(qone-zRand), zRand);
645 else if(chose >= Aone && chose < Aone + Atwo)
651 rtwo*sinu*(qtwo-zRand), zRand);
660 else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
663 rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
664 rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
667 rRand1*std::sin(startPhi), zRand);
672 rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
673 rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
676 rRand1*std::sin(
endPhi), zRand);
693 G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
694 aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
695 fDz = std::fabs(0.5*(zTwo-zOne));
699 aOne = 2.*fDz*fDPhi*fRMax;
700 aTwo = 2.*fDz*fDPhi*fRMin;
701 aFou = 2.*fDz*(fRMax-fRMin);
702 totArea = aOne+aTwo+2.*aFou;
704 cosphi = std::cos(phi);
705 sinphi = std::sin(phi);
708 if(startPhi == 0 &&
endPhi == twopi)
712 if( (chose >= 0) && (chose < aOne) )
714 xRand = fRMax*cosphi;
715 yRand = fRMax*sinphi;
719 else if( (chose >= aOne) && (chose < aOne + aTwo) )
721 xRand = fRMin*cosphi;
722 yRand = fRMin*sinphi;
726 else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
728 xRand = rRand*std::cos(fSPhi+fDPhi);
729 yRand = rRand*std::sin(fSPhi+fDPhi);
736 xRand = rRand*std::cos(fSPhi+fDPhi);
737 yRand = rRand*std::sin(fSPhi+fDPhi);
752 G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
754 cosphi = std::cos(phi);
755 sinphi = std::sin(phi);
759 rRand1 = fRMin1; A1=0.;
764 A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
768 rRand2=fRMax1; Atot=A1;
773 Atot = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
777 if(rCh>A1) { rRand1=rRand2; }
779 xRand = rRand1*cosphi;
780 yRand = rRand1*sinphi;
799 if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
803 return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
812 G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
817 cosphi = std::cos(phi);
818 sinphi = std::sin(phi);
824 std::vector<G4double> areas;
825 std::vector<G4ThreeVector> points;
830 for(i=0; i<numPlanes-1; i++)
855 areas.push_back(Area);
862 totArea += (areas[0]+areas[numPlanes]);
865 if( (chose>=0.) && (chose<areas[0]) )
871 for (i=0; i<numPlanes-1; i++)
874 Achose2 = (Achose1+areas[i+1]);
875 if(chose>=Achose1 && chose<Achose2)
915 G4bool isConvertible=
true;
921 std::vector<G4double> Z;
922 std::vector<G4double> Rmin;
923 std::vector<G4double> Rmax;
936 Rmax.push_back (
corners[1].r);icurr=1;
938 else if (Zprev ==
corners[numPlanes-1].z)
940 Rmin.push_back(
corners[numPlanes-1].r);
952 G4int inextr=0, inextl=0;
953 for (
G4int i=0; i < numPlanes-2; i++)
956 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
958 if((
corners[inextr].z >= Zmax) & (
corners[inextl].
z >= Zmax)) {
break; }
973 Rmin.push_back(
corners[inextl].r);
974 Rmax.push_back(
corners[icurr].r);
978 Rmin.push_back(
corners[inextl].r);
987 Rmin.push_back(
corners[icurl].r);
988 Rmax.push_back(
corners[icurr].r);
992 Rmin.push_back(
corners[icurl].r);
999 isConvertible=
false;
break;
1001 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1009 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1011 Rmin.push_back(
corners[inextl].r);
1012 Rmax.push_back(
corners[inextr].r);
1016 Z.push_back(Zright);
1025 Rmax.push_back(
corners[inextr].r);
1026 Rmin.push_back(
corners[icurr].r);
1030 Rmin.push_back(
corners[icurl].r + (Zright-
corners[icurl].z)/difZl
1032 Rmax.push_back(
corners[inextr].r);
1040 Rmax.push_back(
corners[inextr].r);
1041 Rmin.push_back (
corners[icurr].r);
1045 Rmax.push_back(
corners[inextr].r);
1053 isConvertible=
false;
break;
1063 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1067 Rmax.push_back(
corners[inextr].r);
1068 Rmin.push_back(
corners[inextl].r);
1072 Rmax.push_back(
corners[inextr].r);
1073 Rmin.push_back(
corners[inextl].r);
1085 for(
G4int j=0; j < countPlanes; j++)
1099 std::ostringstream message;
1101 <<
"cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1102 G4Exception(
"G4Polycone::SetOriginalParameters()",
"GeomSolids0002",
1110 for(
G4int j=0; j < numPlanes; j++)
1120 return isConvertible;
G4bool CrossesItself(G4double tolerance)
G4Polyhedron * CreatePolyhedron() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4VCSGfaceted & operator=(const G4VCSGfaceted &source)
void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
G4bool fRebuildPolyhedron
G4Polyhedron * fpPolyhedron
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetPointOnSurface() const
G4ThreeVector GetPointOnRing(G4double fRMin, G4double fRMax, G4double fRMin2, G4double fRMax2, G4double zOne) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4bool MustBeOutside(const G4ThreeVector &p) const
G4PolyconeSideRZ * corners
std::ostream & StreamInfo(std::ostream &os) const
G4ThreeVector GetPointOnCone(G4double fRmin1, G4double fRmax1, G4double fRmin2, G4double fRmax2, G4double zOne, G4double zTwo, G4double &totArea) const
G4GeometryType GetEntityType() const
G4ThreeVector GetPointOnCut(G4double fRMin1, G4double fRMax1, G4double fRMin2, G4double fRMax2, G4double zOne, G4double zTwo, G4double &totArea) const
G4bool RemoveDuplicateVertices(G4double tolerance)
G4bool RemoveRedundantVertices(G4double tolerance)
G4Polycone & operator=(const G4Polycone &source)
G4GLOB_DLL std::ostream G4cout
G4int NumVertices() const
void SetOriginalParameters(G4PolyconeHistorical *pars)
EInside Inside(const G4ThreeVector &p) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
virtual void ComputeDimensions(G4Box &, const G4int, const G4VPhysicalVolume *) const
void CopyStuff(const G4Polycone &source)
G4EnclosingCylinder * enclosingCylinder
G4ThreeVector GetPointOnTubs(G4double fRMin, G4double fRMax, G4double zOne, G4double zTwo, G4double &totArea) const
static const double degree
G4bool BisectedBy(G4double a1, G4double b1, G4double a2, G4double b2, G4double tolerance)
G4Polycone(const G4String &name, G4double phiStart, G4double phiTotal, G4int numZPlanes, const G4double zPlane[], const G4double rInner[], const G4double rOuter[])
void Create(G4double phiStart, G4double phiTotal, G4ReduciblePolygon *rz)
G4PolyconeHistorical * original_parameters
virtual EInside Inside(const G4ThreeVector &p) const
G4bool ShouldMiss(const G4ThreeVector &p, const G4ThreeVector &v) const