76 : fIsValidNorm(false),
fName(name)
87 for (
G4int i=0; i<4; i++)
89 fCorners[i].
set(kInfinity, kInfinity, kInfinity);
95 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
96 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
110 : fIsValidNorm(false),
fName(name)
122 for (
G4int i=0; i<4; i++)
124 fCorners[i].
set(kInfinity, kInfinity, kInfinity);
130 fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
131 fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
139 : fHandedness(0), fIsValidNorm(false), kCarTolerance(0.),
145 fNeighbours[0] = fNeighbours[1] = fNeighbours[2] = fNeighbours[3] = 0;
178 if (fAmIOnLeftSide.me == me
179 && fAmIOnLeftSide.vec == vec
180 && fAmIOnLeftSide.withTol == withtol) {
181 return fAmIOnLeftSide.amIOnLeftSide;
184 fAmIOnLeftSide.me = me;
185 fAmIOnLeftSide.vec = vec;
186 fAmIOnLeftSide.withTol = withtol;
194 G4double metcrossvect = met.
x() * vect.y() - met.
y() * vect.x();
197 if (met.
x() * ivect.
y() - met.
y() * ivect.
x() > 0 &&
199 fAmIOnLeftSide.amIOnLeftSide = 1;
200 }
else if (met.
x() * rvect.
y() - met.
y() * rvect.
x() < 0 &&
202 fAmIOnLeftSide.amIOnLeftSide = -1;
204 fAmIOnLeftSide.amIOnLeftSide = 0;
207 if (metcrossvect > 0) {
208 fAmIOnLeftSide.amIOnLeftSide = 1;
209 }
else if (metcrossvect < 0 ) {
210 fAmIOnLeftSide.amIOnLeftSide = -1;
212 fAmIOnLeftSide.amIOnLeftSide = 0;
217 G4cout <<
" === G4VTwistSurface::AmIOnLeftSide() =============="
219 G4cout <<
" Name , returncode : " << fName <<
" "
220 << fAmIOnLeftSide.amIOnLeftSide <<
G4endl;
221 G4cout <<
" me, vec : " << std::setprecision(14) << me
223 G4cout <<
" met, vect : " << met <<
" " << vect <<
G4endl;
224 G4cout <<
" ivec, rvec : " << ivect <<
" " << rvect <<
G4endl;
228 G4cout <<
" =============================================="
232 return fAmIOnLeftSide.amIOnLeftSide;
257 std::ostringstream message;
258 message <<
"Point is in the corner area." <<
G4endl
259 <<
" Point is in the corner area. This function returns"
261 <<
" a direction vector of a boundary line." <<
G4endl
262 <<
" areacode = " << areacode;
263 G4Exception(
"G4VTwistSurface::DistanceToBoundary()",
"GeomSolids0003",
269 xx.
set(t*p.
x(), t*p.
y(), x0.
z());
270 dist = (xx -
p).mag();
277 std::ostringstream message;
278 message <<
"Bad areacode of boundary." <<
G4endl
279 <<
" areacode = " << areacode;
280 G4Exception(
"G4VTwistSurface::DistanceToBoundary()",
"GeomSolids0003",
294 G4cout <<
" ~~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~~" <<
G4endl;
298 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
307 distance[i] = kInfinity ;
321 for (
G4int i=0; i< nxx; i++) {
334 if ((normal * gv) >= 0) {
337 G4cout <<
" G4VTwistSurface::DistanceToIn(p,v): "
338 <<
"particle goes outword the surface." <<
G4endl;
348 if (distance[i] < bestdistance) {
349 bestdistance = distance[i];
353 G4cout <<
" G4VTwistSurface::DistanceToIn(p,v): "
354 <<
" areacode sInside name, distance = "
355 << fName <<
" "<< bestdistance <<
G4endl;
366 G4bool isaccepted[2] = {
false,
false};
369 for (
G4int j=0; j< nneighbours; j++) {
379 tmpdist[l] = kInfinity ;
381 tmpisvalid[l] = false ;
385 gp, gv, tmpgxx, tmpdist,
386 tmpareacode, tmpisvalid,
390 for (
G4int k=0; k< tmpnxx; k++) {
401 G4cout <<
" G4VTwistSurface:DistanceToIn(p,v): "
402 <<
" intersection "<< tmpgxx[k] << G4endl
403 <<
" is inside of neighbour surface of " << fName
404 <<
" . returning kInfinity." <<
G4endl;
405 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~"
409 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
412 if (tmpisvalid[k])
return kInfinity;
421 neighbours[j], tmpareacode[k])) {
424 neighbournormal = neighbours[j]->
GetNormal(tmpgxx[k],
true);
425 if (neighbournormal * gv < 0) isaccepted[j] =
true;
432 if (nneighbours == 1) isaccepted[1] =
true;
438 if (isaccepted[0] ==
true && isaccepted[1] ==
true) {
439 if (distance[i] < bestdistance) {
440 bestdistance = distance[i];
444 G4cout <<
" G4VTwistSurface::DistanceToIn(p,v): "
445 <<
" areacode sBoundary & sBoundary distance = "
446 << fName <<
" " << distance[i] <<
G4endl;
458 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" <<
G4endl;
461 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
463 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" <<
G4endl;
464 G4cout <<
" Name, i : " << fName <<
" , " << besti <<
G4endl;
467 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
483 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" <<
G4endl;
487 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
498 distance[i] = kInfinity ;
509 for (i=0; i<nxx; i++) {
515 if (normal * gv <= 0) {
518 G4cout <<
" G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0, normal "
519 << fName <<
" " << normal
524 if (distance[i] < bestdistance) {
525 bestdistance = distance[i];
533 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~~" <<
G4endl;
537 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
539 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~~" <<
G4endl;
540 G4cout <<
" Name, i : " << fName <<
" , " << i <<
G4endl;
543 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
557 G4cout <<
"~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" <<
G4endl;
560 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
569 distance[i] = kInfinity ;
578 G4cout <<
"~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" <<
G4endl;
582 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
600 G4bool testbitmode =
true;
604 if (iscorner[0] && iscorner[1]) {
617 }
else if ((
IsBoundary(areacode1, testbitmode) && (!iscorner[0])) &&
618 (
IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) {
649 G4int &boundarytype)
const
656 for (i=0; i<4; i++) {
663 std::ostringstream message;
664 message <<
"Not registered boundary." <<
G4endl
665 <<
" Boundary at areacode " << std::hex << areacode
667 <<
" is not registered.";
668 G4Exception(
"G4VTwistSurface::GetBoundaryParameters()",
"GeomSolids0002",
683 std::ostringstream message;
684 message <<
"Point is in the corner area." <<
G4endl
685 <<
" This function returns "
686 <<
"a direction vector of a boundary line." <<
G4endl
687 <<
" areacode = " << areacode;
688 G4Exception(
"G4VTwistSurface::GetBoundaryAtPZ()",
"GeomSolids0003",
697 for (
G4int i=0; i<4; i++) {
706 std::ostringstream message;
707 message <<
"Not registered boundary." <<
G4endl
708 <<
" Boundary at areacode " << areacode <<
G4endl
709 <<
" is not registered.";
710 G4Exception(
"G4VTwistSurface::GetBoundaryAtPZ()",
"GeomSolids0002",
714 if (((boundarytype &
sAxisPhi) == sAxisPhi) ||
715 ((boundarytype &
sAxisRho) == sAxisRho)) {
716 std::ostringstream message;
717 message <<
"Not a z-depended line boundary." <<
G4endl
718 <<
" Boundary at areacode " << areacode <<
G4endl
719 <<
" is not a z-depended line.";
720 G4Exception(
"G4VTwistSurface::GetBoundaryAtPZ()",
"GeomSolids0002",
723 return ((p.
z() - x0.
z()) / d.
z()) * d + x0;
731 if ((areacode &
sCorner) != sCorner){
732 std::ostringstream message;
733 message <<
"Area code must represents corner." <<
G4endl
734 <<
" areacode " << areacode;
735 G4Exception(
"G4VTwistSurface::SetCorner()",
"GeomSolids0002",
740 fCorners[0].
set(x, y, z);
741 }
else if ((areacode &
sC0Max1Min) == sC0Max1Min) {
742 fCorners[1].
set(x, y, z);
743 }
else if ((areacode &
sC0Max1Max) == sC0Max1Max) {
744 fCorners[2].
set(x, y, z);
745 }
else if ((areacode &
sC0Min1Max) == sC0Min1Max) {
746 fCorners[3].
set(x, y, z);
755 if ((areacode &
sBoundary) != sBoundary) {
756 G4Exception(
"G4VTwistSurface::GetBoundaryAxis()",
"GeomSolids0003",
760 for (i=0; i<2; i++) {
762 G4int whichaxis = 0 ;
772 if (axiscode == (whichaxis &
sAxisX)) {
774 }
else if (axiscode == (whichaxis &
sAxisY)) {
776 }
else if (axiscode == (whichaxis &
sAxisZ)) {
778 }
else if (axiscode == (whichaxis &
sAxisRho)) {
780 }
else if (axiscode == (whichaxis &
sAxisPhi)) {
783 std::ostringstream message;
784 message <<
"Not supported areacode." <<
G4endl
785 <<
" areacode " << areacode;
786 G4Exception(
"G4VTwistSurface::GetBoundaryAxis()",
"GeomSolids0001",
823 std::ostringstream message;
824 message <<
"Not located on a boundary!" <<
G4endl
825 <<
" areacode " << areacode;
826 G4Exception(
"G4VTwistSurface::GetBoundaryLimit()",
"GeomSolids1002",
837 const G4int &boundarytype)
847 for (i=0; i<4; i++) {
848 if (fBoundaries[i].IsEmpty()) {
849 fBoundaries[i].
SetFields(axiscode, direction,
857 G4Exception(
"G4VTwistSurface::SetBoundary()",
"GeomSolids0003",
861 std::ostringstream message;
862 message <<
"Invalid axis-code." <<
G4endl
864 << std::hex << axiscode << std::dec;
865 G4Exception(
"G4VTwistSurface::SetBoundary()",
"GeomSolids0003",
880 return i * ( k - 1 ) + j ;
883 else if ( iside == 1 ) {
884 return (k-1)*(k-1) + i*(k-1) + j ;
887 else if ( iside == 2 ) {
888 return 2*(k-1)*(k-1) + i*(k-1) + j ;
891 else if ( iside == 3 ) {
892 return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k-1) + j ;
895 else if ( iside == 4 ) {
896 return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*(k-1) + j ;
899 else if ( iside == 5 ) {
900 return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*(k-1) + j ;
905 std::ostringstream message;
906 message <<
"Not correct side number: "
908 <<
"iside is " << iside <<
" but should be "
909 <<
"0,1,2,3,4 or 5" <<
".";
910 G4Exception(
"G4TwistSurface::G4GetFace()",
"GeomSolids0002",
938 return k*k + i*k + j ;
941 else if ( iside == 2 ) {
943 if ( i == 0 ) {
return j ; }
944 else if ( i == n-1 ) {
return k*k + j ; }
945 else {
return 2*k*k + 4*(i-1)*(k-1) + j ; }
948 else if ( iside == 3 ) {
950 if ( i == 0 ) {
return (j+1)*k - 1 ; }
951 else if ( i == n-1 ) {
return k*k + (j+1)*k - 1 ; }
952 else {
return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j ; }
954 else if ( iside == 4 ) {
956 if ( i == 0 ) {
return k*k - 1 - j ; }
957 else if ( i == n-1 ) {
return 2*k*k - 1 - j ; }
958 else {
return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + j ;
961 else if ( iside == 5 ) {
963 if ( i == 0 ) {
return k*k - (j+1)*k ; }
964 else if ( i == n-1) {
return 2*k*k - (j+1)*k ; }
966 if ( j == k-1 ) {
return 2*k*k + 4*(i-1)*(k-1) ; }
967 else {
return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) + j ; }
973 std::ostringstream message;
974 message <<
"Not correct side number: "
976 <<
"iside is " << iside <<
" but should be "
977 <<
"0,1,2,3,4 or 5" <<
".";
978 G4Exception(
"G4TwistSurface::G4GetNode()",
"GeomSolids0002",
1014 if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) ) {
1020 if ( orientation < 0 ) { number = ( 3 - number ) ; }
1023 if ( ( j>=1 && j<=k-3 ) ) {
1026 return ( number == 3 ) ? 1 : -1 ;
1029 else if ( i == n-2 ) {
1030 return ( number == 1 ) ? 1 : -1 ;
1034 std::ostringstream message;
1035 message <<
"Not correct face number: " <<
GetName() <<
" !";
1036 G4Exception(
"G4TwistSurface::G4GetEdgeVisibility()",
1041 if ( ( i>=1 && i<=n-3 ) ) {
1044 return ( number == 0 ) ? 1 : -1 ;
1047 else if ( j == k-2 ) {
1048 return ( number == 2 ) ? 1 : -1 ;
1052 std::ostringstream message;
1053 message <<
"Not correct face number: " <<
GetName() <<
" !";
1054 G4Exception(
"G4TwistSurface::G4GetEdgeVisibility()",
1060 if ( i == 0 && j == 0 ) {
1061 return ( number == 0 || number == 3 ) ? 1 : -1 ;
1063 else if ( i == 0 && j == k-2 ) {
1064 return ( number == 2 || number == 3 ) ? 1 : -1 ;
1066 else if ( i == n-2 && j == k-2 ) {
1067 return ( number == 1 || number == 2 ) ? 1 : -1 ;
1069 else if ( i == n-2 && j == 0 ) {
1070 return ( number == 0 || number == 1 ) ? 1 : -1 ;
1073 std::ostringstream message;
1074 message <<
"Not correct face number: " <<
GetName() <<
" !";
1075 G4Exception(
"G4TwistSurface::G4GetEdgeVisibility()",
1079 std::ostringstream message;
1080 message <<
"Not correct face number: " <<
GetName() <<
" !";
1081 G4Exception(
"G4TwistSurface::G4GetEdgeVisibility()",
"GeomSolids0003",
1098 G4cout <<
"/* G4VTwistSurface::DebugPrint():-------------------------------"
1101 G4cout <<
"/* Axis = " << std::hex <<
fAxis[0] <<
" "
1102 << std::hex <<
fAxis[1]
1103 <<
" (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)"
1105 G4cout <<
"/* BoundaryLimit(in local) fAxis0(min, max) = ("<<
fAxisMin[0]
1107 G4cout <<
"/* BoundaryLimit(in local) fAxis1(min, max) = ("<<
fAxisMin[1]
1109 G4cout <<
"/* Cornar point sC0Min1Min = " << A <<
G4endl;
1110 G4cout <<
"/* Cornar point sC0Max1Min = " << B <<
G4endl;
1111 G4cout <<
"/* Cornar point sC0Max1Max = " << C <<
G4endl;
1112 G4cout <<
"/* Cornar point sC0Min1Max = " << D <<
G4endl;
1113 G4cout <<
"/*---------------------------------------------------------"
1128 fDistance[i] = kInfinity;
1130 fIsValid[i] =
false;
1131 fXX[i].
set(kInfinity, kInfinity, kInfinity);
1134 fLastp.
set(kInfinity, kInfinity, kInfinity);
1135 fLastv.
set(kInfinity, kInfinity, kInfinity);
1161 fDistance[i] = dist;
1162 fAreacode[i] = areacode;
1163 fIsValid[i] = isvalid;
1166 fLastValidate = validate;
1173 G4Exception(
"G4VTwistSurface::CurrentStatus::SetCurrentStatus()",
1182 fLastv.set(kInfinity, kInfinity, kInfinity);
1196 if (validate == fLastValidate && p && *p == fLastp)
1198 if (!v || (v && *v == fLastv))
return;
1203 fDistance[i] = kInfinity;
1205 fIsValid[i] =
false;
1209 fLastp.set(kInfinity, kInfinity, kInfinity);
1210 fLastv.set(kInfinity, kInfinity, kInfinity);
1221 G4cout <<
"CurrentStatus::Dist0,1= " << fDistance[0]
1222 <<
" " << fDistance[1] <<
" areacode = " << fAreacode[0]
1223 <<
" " << fAreacode[1] <<
G4endl;
1234 : fBoundaryAcode(-1), fBoundaryType(0)
1252 const G4int &boundarytype)
1254 fBoundaryAcode = areacode;
1255 fBoundaryDirection =
d;
1257 fBoundaryType = boundarytype;
1265 if (fBoundaryAcode == -1)
return true;
1276 G4int &boundarytype)
const
1283 std::ostringstream message;
1284 message <<
"Located in the corner area." <<
G4endl
1285 <<
" This function returns a direction vector of "
1286 <<
"a boundary line." <<
G4endl
1287 <<
" areacode = " << areacode;
1288 G4Exception(
"G4VTwistSurface::Boundary::GetBoundaryParameters()",
1291 if ((areacode &
sSizeMask) != (fBoundaryAcode & sSizeMask))
1295 d = fBoundaryDirection;
1297 boundarytype = fBoundaryType;
static const G4int sAxisZ
void set(double x, double y, double z)
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
static const G4int sAxisPhi
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
G4SurfCurNormal fCurrentNormal
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)=0
G4double GetSurfaceTolerance() const
G4ThreeVector GetCorner(G4int areacode) const
G4bool GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
G4bool IsAxis0(G4int areacode) const
static const G4int sOutside
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
static const G4int sAreaMask
G4bool IsBoundary(G4int areacode, G4bool testbitmode=false) const
G4bool IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1, G4VTwistSurface *surface2, G4int areacode2) const
static const G4int sAxisMask
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
static const G4int sAxisX
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
static const G4int sC0Min1Min
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
G4VTwistSurface ** GetNeighbours()
G4GLOB_DLL std::ostream G4cout
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
static const G4int sAxis1
static const G4int sC0Max1Max
static const G4int sBoundary
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
static const G4int sAxis0
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisMin
static const G4int sAxisY
G4bool IsCorner(G4int areacode, G4bool testbitmode=false) const
static const G4int sInside
G4VTwistSurface(const G4String &name)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const G4int sCorner
void GetBoundaryAxis(G4int areacode, EAxis axis[]) const
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
void SetFields(const G4int &areacode, const G4ThreeVector &d, const G4ThreeVector &x0, const G4int &boundarytype)
static const G4int sAxisMax
void GetBoundaryLimit(G4int areacode, G4double limit[]) const
virtual G4String GetName() const
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
HepRotation & rotateZ(double delta)
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
Hep3Vector cross(const Hep3Vector &) const
static const G4int sSizeMask
G4bool IsAxis1(G4int areacode) const
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
virtual ~G4VTwistSurface()
static const G4int sC0Max1Min
G4double GetAngularTolerance() const
static G4GeometryTolerance * GetInstance()
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
static const G4int sAxisRho