76 : fIsValidNorm(false),
fName(name)
87 for (
G4int i=0; i<4; i++)
110 : fIsValidNorm(false),
fName(name)
122 for (
G4int i=0; i<4; i++)
194 G4double metcrossvect = met.
x() * vect.y() - met.
y() * vect.x();
197 if (met.
x() * ivect.
y() - met.
y() * ivect.
x() > 0 &&
200 }
else if (met.
x() * rvect.
y() - met.
y() * rvect.
x() < 0 &&
207 if (metcrossvect > 0) {
209 }
else if (metcrossvect < 0 ) {
217 G4cout <<
" === G4VTwistSurface::AmIOnLeftSide() =============="
221 G4cout <<
" me, vec : " << std::setprecision(14) << me
223 G4cout <<
" met, vect : " << met <<
" " << vect <<
G4endl;
224 G4cout <<
" ivec, rvec : " << ivect <<
" " << rvect <<
G4endl;
228 G4cout <<
" =============================================="
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();
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;
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 = "
366 G4bool isaccepted[2] = {
false,
false};
369 for (
G4int j=0; j< nneighbours; j++) {
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 <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
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 = "
458 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" <<
G4endl;
461 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
463 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" <<
G4endl;
467 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
483 G4cout <<
"~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" <<
G4endl;
487 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
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;
543 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
557 G4cout <<
"~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" <<
G4endl;
560 G4cout <<
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" <<
G4endl;
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++) {
664 message <<
"Not registered boundary." <<
G4endl
665 <<
" Boundary at areacode " << std::hex << areacode
667 <<
" is not registered.";
668 G4Exception(
"G4VTwistSurface::GetBoundaryParameters()",
"GeomSolids0002",
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++) {
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)) {
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){
733 message <<
"Area code must represents corner." <<
G4endl
734 <<
" areacode " << areacode;
735 G4Exception(
"G4VTwistSurface::SetCorner()",
"GeomSolids0002",
741 }
else if ((areacode &
sC0Max1Min) == sC0Max1Min) {
743 }
else if ((areacode &
sC0Max1Max) == sC0Max1Max) {
745 }
else if ((areacode &
sC0Min1Max) == sC0Min1Max) {
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)) {
784 message <<
"Not supported areacode." <<
G4endl
785 <<
" areacode " << areacode;
786 G4Exception(
"G4VTwistSurface::GetBoundaryAxis()",
"GeomSolids0001",
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++) {
857 G4Exception(
"G4VTwistSurface::SetBoundary()",
"GeomSolids0003",
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 ;
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 ; }
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 ;
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 ;
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 ;
1074 message <<
"Not correct face number: " <<
GetName() <<
" !";
1075 G4Exception(
"G4TwistSurface::G4GetEdgeVisibility()",
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 <<
"/*---------------------------------------------------------"
1161 fDistance[i] = dist;
1162 fAreacode[i] = areacode;
1163 fIsValid[i] = isvalid;
1166 fLastValidate = validate;
1173 G4Exception(
"G4VTwistSurface::CurrentStatus::SetCurrentStatus()",
1196 if (validate == fLastValidate && p && *p == fLastp)
1198 if (!v || (*v == fLastv))
return;
1205 fIsValid[i] =
false;
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
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 sC0Min1Max
void set(double x, double y, double z)
static const G4int sInside
static const G4int sAreaMask
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisRho
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector fCorners[4]
static const G4double kInfinity
G4bool IsCorner(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisX
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
G4bool IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1, G4VTwistSurface *surface2, G4int areacode2) const
void message(RunManager *runmanager)
static const G4int sCorner
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
G4double GetAngularTolerance() const
G4bool GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
static double normal(HepRandomEngine *eptr)
const G4double kCarTolerance
G4SurfSideQuery fAmIOnLeftSide
static const G4int sAxisMask
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)=0
virtual ~G4VTwistSurface()
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal)=0
G4int fAreacode[G4VSURFACENXX]
static const G4int sAxisZ
static const G4int sC0Min1Min
G4double fDistance[G4VSURFACENXX]
virtual G4String GetName() const
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
static const G4int sAxis1
void GetBoundaryLimit(G4int areacode, G4double limit[]) const
void SetFields(const G4int &areacode, const G4ThreeVector &d, const G4ThreeVector &x0, const G4int &boundarytype)
static const G4int sSizeMask
virtual G4double DistanceTo(const G4ThreeVector &gp, G4ThreeVector &gxx)
void GetBoundaryAxis(G4int areacode, EAxis axis[]) const
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sC0Max1Max
G4double DistanceToLine(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
G4bool IsAxis1(G4int areacode) const
static const G4int sAxisPhi
double A(double temperature)
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
virtual G4double DistanceToIn(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
G4ThreeVector GetCorner(G4int areacode) const
Hep3Vector cross(const Hep3Vector &) const
G4VTwistSurface(const G4String &name)
G4bool fIsValid[G4VSURFACENXX]
G4VTwistSurface ** GetNeighbours()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
virtual G4double DistanceToBoundary(G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
G4bool IsBoundary(G4int areacode, G4bool testbitmode=false) const
G4bool IsAxis0(G4int areacode) const
G4VTwistSurface * fNeighbours[4]
G4ThreeVector fXX[G4VSURFACENXX]
static const G4int sAxis0
HepRotation & rotateZ(double delta)
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
G4GLOB_DLL std::ostream G4cout
static G4GeometryTolerance * GetInstance()
static const G4int sOutside
G4SurfCurNormal fCurrentNormal
static const G4int sAxisMin
static const G4int sBoundary
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sC0Max1Min
static const G4int sAxisMax
double B(double temperature)
virtual G4double DistanceToOut(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
G4double GetSurfaceTolerance() const
static const G4int sAxisY
virtual G4int AmIOnLeftSide(const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const