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;
545 G4int oldprc = os.precision(16);
546 os <<
"-----------------------------------------------------------\n"
547 <<
" *** Dump for solid - " <<
GetName() <<
" ***\n"
548 <<
" ===================================================\n"
549 <<
" Solid type: G4Polycone\n"
552 <<
" ending phi angle : " <<
endPhi/
degree <<
" degrees \n";
556 os <<
" number of Z planes: " << numPlanes <<
"\n"
558 for (i=0; i<numPlanes; i++)
560 os <<
" Z plane " << i <<
": "
563 os <<
" Tangent distances to inner surface (Rmin): \n";
564 for (i=0; i<numPlanes; i++)
566 os <<
" Z plane " << i <<
": "
569 os <<
" Tangent distances to outer surface (Rmax): \n";
570 for (i=0; i<numPlanes; i++)
572 os <<
" Z plane " << i <<
": "
576 os <<
" number of RZ points: " <<
numCorner <<
"\n"
577 <<
" RZ values (corners): \n";
583 os <<
"-----------------------------------------------------------\n";
584 os.precision(oldprc);
602 G4double Aone, Atwo, Afive, phi, zRand, fDPhi, cosu, sinu;
603 G4double rRand1, rmin, rmax, chose, rone, rtwo, qone, qtwo;
604 G4double fDz=(zTwo-zOne)/2., afDz=std::fabs(fDz);
607 rone = (fRmax1-fRmax2)/(2.*fDz);
608 rtwo = (fRmin1-fRmin2)/(2.*fDz);
609 if(fRmax1==fRmax2){qone=0.;}
611 qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2);
613 if(fRmin1==fRmin2){qtwo=0.;}
615 qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2);
617 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*(
sqr(fRmin1-fRmin2)+
sqr(zTwo-zOne));
618 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*(
sqr(fRmax1-fRmax2)+
sqr(zTwo-zOne));
619 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
620 totArea = Aone+Atwo+2.*Afive;
623 cosu = std::cos(phi);
624 sinu = std::sin(phi);
627 if( (startPhi == 0) && (
endPhi == twopi) ) { Afive = 0; }
629 if( (chose >= 0) && (chose < Aone) )
635 rone*sinu*(qone-zRand), zRand);
644 else if(chose >= Aone && chose < Aone + Atwo)
650 rtwo*sinu*(qtwo-zRand), zRand);
659 else if( (chose >= Aone + Atwo + Afive) && (chose < Aone + Atwo + 2.*Afive) )
662 rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
663 rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
666 rRand1*std::sin(startPhi), zRand);
671 rmin = fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2);
672 rmax = fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2);
675 rRand1*std::sin(
endPhi), zRand);
692 G4double xRand,yRand,zRand,phi,cosphi,sinphi,chose,
693 aOne,aTwo,aFou,rRand,fDz,fSPhi,fDPhi;
694 fDz = std::fabs(0.5*(zTwo-zOne));
698 aOne = 2.*fDz*fDPhi*fRMax;
699 aTwo = 2.*fDz*fDPhi*fRMin;
700 aFou = 2.*fDz*(fRMax-fRMin);
701 totArea = aOne+aTwo+2.*aFou;
703 cosphi = std::cos(phi);
704 sinphi = std::sin(phi);
707 if(startPhi == 0 &&
endPhi == twopi)
711 if( (chose >= 0) && (chose < aOne) )
713 xRand = fRMax*cosphi;
714 yRand = fRMax*sinphi;
718 else if( (chose >= aOne) && (chose < aOne + aTwo) )
720 xRand = fRMin*cosphi;
721 yRand = fRMin*sinphi;
725 else if( (chose >= aOne+aTwo) && (chose <aOne+aTwo+aFou) )
727 xRand = rRand*std::cos(fSPhi+fDPhi);
728 yRand = rRand*std::sin(fSPhi+fDPhi);
735 xRand = rRand*std::cos(fSPhi+fDPhi);
736 yRand = rRand*std::sin(fSPhi+fDPhi);
751 G4double xRand,yRand,phi,cosphi,sinphi,rRand1,rRand2,A1,Atot,rCh;
753 cosphi = std::cos(phi);
754 sinphi = std::sin(phi);
758 rRand1 = fRMin1; A1=0.;
763 A1=std::fabs(fRMin2*fRMin2-fRMin1*fRMin1);
767 rRand2=fRMax1; Atot=A1;
772 Atot = A1+std::fabs(fRMax2*fRMax2-fRMax1*fRMax1);
776 if(rCh>A1) { rRand1=rRand2; }
778 xRand = rRand1*cosphi;
779 yRand = rRand1*sinphi;
798 if( (fRMin1 == fRMin2) && (fRMax1 == fRMax2) )
802 return GetPointOnCone(fRMin1,fRMax1,fRMin2,fRMax2,zOne,zTwo,totArea);
811 G4double Area=0,totArea=0,Achose1=0,Achose2=0,phi,cosphi,sinphi,rRand;
816 cosphi = std::cos(phi);
817 sinphi = std::sin(phi);
823 std::vector<G4double> areas;
824 std::vector<G4ThreeVector> points;
829 for(i=0; i<numPlanes-1; i++)
854 areas.push_back(Area);
861 totArea += (areas[0]+areas[numPlanes]);
864 if( (chose>=0.) && (chose<areas[0]) )
870 for (i=0; i<numPlanes-1; i++)
873 Achose2 = (Achose1+areas[i+1]);
874 if(chose>=Achose1 && chose<Achose2)
914 G4bool isConvertible=
true;
920 std::vector<G4double> Z;
921 std::vector<G4double> Rmin;
922 std::vector<G4double> Rmax;
935 Rmax.push_back (
corners[1].r);icurr=1;
937 else if (Zprev ==
corners[numPlanes-1].z)
939 Rmin.push_back(
corners[numPlanes-1].r);
951 G4int inextr=0, inextl=0;
952 for (
G4int i=0; i < numPlanes-2; i++)
955 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
957 if((
corners[inextr].z >= Zmax) & (
corners[inextl].
z >= Zmax)) {
break; }
972 Rmin.push_back(
corners[icurl].r);
973 Rmax.push_back(Rmax[countPlanes-2]);
974 Rmax[countPlanes-2]=
corners[icurl].
r;
978 Rmin.push_back(
corners[inextl].r);
979 Rmax.push_back(
corners[icurl].r);
984 Rmin.push_back(
corners[inextl].r);
990 isConvertible=
false;
break;
992 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1000 icurl=(icurl == 0)? numPlanes-1 : icurl-1;
1002 Rmin.push_back(
corners[inextl].r);
1003 Rmax.push_back (
corners[inextr].r);
1007 Z.push_back(Zright);
1016 Rmin.push_back(
corners[icurr].r);
1017 Rmax.push_back(
corners[inextr].r);
1021 Rmin.push_back(
corners[inextr].r);
1022 Rmax.push_back(
corners[icurr].r);
1023 Rmax[countPlanes-2]=
corners[inextr].
r;
1031 Rmax.push_back(
corners[inextr].r);
1032 Rmin.push_back (
corners[icurr].r);
1036 Rmax.push_back(
corners[inextr].r);
1044 isConvertible=
false;
break;
1054 inextl=(icurl <= 0)? numPlanes-1 : icurl-1;
1058 Rmax.push_back(
corners[inextr].r);
1059 Rmin.push_back(
corners[inextl].r);
1063 Rmax.push_back(
corners[inextr].r);
1064 Rmin.push_back(
corners[inextl].r);
1076 for(
G4int j=0; j < countPlanes; j++)
1090 std::ostringstream message;
1092 <<
"cannot be converted to Polycone with (Rmin,Rmaz,Z) parameters!";
1093 G4Exception(
"G4Polycone::SetOriginalParameters()",
"GeomSolids0002",
1101 for(
G4int j=0; j < numPlanes; j++)
1111 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)
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
virtual G4Polyhedron * GetPolyhedron() const