54 const G4int handedness,
65 axis0min, axis1min, axis0max, axis1max),
66 fKappa(kappa), fTanStereo(tanstereo),
67 fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(
twopi)
71 G4Exception(
"G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
73 "Should swap axis0 and axis1!");
76 fInside.gp.set(kInfinity, kInfinity, kInfinity);
110 if (handedness < 0) {
111 fTanStereo = TanInnerStereo;
114 fTanStereo = TanOuterStereo;
117 fTan2Stereo = fTanStereo * fTanStereo;
123 fInside.gp.set(kInfinity, kInfinity, kInfinity);
126 SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ;
136 fR0(0.), fR02(0.), fDPhi(0.)
175 normal = normal.unit();
194 if (fInside.gp == gp) {
195 return fInside.inside;
204 return fInside.inside;
211 if (distanceToOut < -halftol) {
217 G4int areacode = GetAreaCode(p);
223 if (distanceToOut <= halftol) {
229 G4cout <<
"WARNING - G4TwistTubsHypeSide::Inside()" <<
G4endl
230 <<
" Invalid option !" <<
G4endl
231 <<
" name, areacode, distanceToOut = "
232 <<
GetName() <<
", " << std::hex << areacode << std::dec <<
", "
233 << distanceToOut <<
G4endl;
237 return fInside.inside;
302 for (i=0; i<2; i++) {
303 distance[i] = kInfinity;
306 gxx[i].
set(kInfinity, kInfinity, kInfinity);
335 if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz)) {
337 distance[0] = kInfinity;
339 isvalid[0], 0, validate, &gp, &gv);
344 G4double xxz = std::sqrt(fR02 / (vslope2 - fTan2Stereo))
345 * (vz / std::fabs(vz)) ;
347 xx[0].
set(t*v.
x(), t*v.
y(), xxz);
350 xx[0].
set(v.
x()*fR0, v.
y()*fR0, 0);
352 distance[0] = xx[0].
mag();
356 areacode[0] = GetAreaCode(xx[0]);
358 if (distance[0] >= 0) isvalid[0] =
true;
361 areacode[0] = GetAreaCode(xx[0],
false);
363 if (distance[0] >= 0) isvalid[0] =
true;
367 if (distance[0] >= 0) isvalid[0] =
true;
371 isvalid[0], 1, validate, &gp, &gv);
380 G4double b = 2.0 * ( p.
x() * v.
x() + p.
y() * v.
y() - p.
z() * v.
z() * fTan2Stereo );
381 G4double c = p.
x()*p.
x() + p.
y()*p.
y() - fR02 - p.
z()*p.
z()*fTan2Stereo;
389 xx[0] = p + distance[0]*
v;
393 areacode[0] = GetAreaCode(xx[0]);
395 if (distance[0] >= 0) isvalid[0] =
true;
398 areacode[0] = GetAreaCode(xx[0],
false);
400 if (distance[0] >= 0) isvalid[0] =
true;
404 if (distance[0] >= 0) isvalid[0] =
true;
408 isvalid[0], 1, validate, &gp, &gv);
417 isvalid[0], 0, validate, &gp, &gv);
426 G4double tmpdist[2] = {kInfinity, kInfinity};
429 G4bool tmpisvalid[2] = {
false,
false};
432 for (i=0; i<2; i++) {
433 tmpdist[i] = factor*(-b - D);
435 tmpxx[i] = p + tmpdist[i]*
v;
438 tmpareacode[i] = GetAreaCode(tmpxx[i]);
440 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
444 tmpareacode[i] = GetAreaCode(tmpxx[i],
false);
446 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
451 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
456 if (tmpdist[0] <= tmpdist[1]) {
457 distance[0] = tmpdist[0];
458 distance[1] = tmpdist[1];
463 areacode[0] = tmpareacode[0];
464 areacode[1] = tmpareacode[1];
465 isvalid[0] = tmpisvalid[0];
466 isvalid[1] = tmpisvalid[1];
468 distance[0] = tmpdist[1];
469 distance[1] = tmpdist[0];
474 areacode[0] = tmpareacode[1];
475 areacode[1] = tmpareacode[0];
476 isvalid[0] = tmpisvalid[1];
477 isvalid[1] = tmpisvalid[0];
481 isvalid[0], 2, validate, &gp, &gv);
483 isvalid[1], 2, validate, &gp, &gv);
491 isvalid[0], 0, validate, &gp, &gv);
528 for (
G4int i=0; i<2; i++) {
529 distance[i] = kInfinity;
531 gxx[i].
set(kInfinity, kInfinity, kInfinity);
544 for (
G4int i=0; i<2; i++) {
548 if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol) {
567 G4double r1 = std::sqrt(fR02 + pz * pz * fTan2Stereo);
571 if (prho > r1 + halftol) {
578 G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
579 G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
587 distance[0] = (pabsz - xx1).mag();
593 }
else if (prho < r1 - halftol) {
599 xx1.
set(r1, 0. , pz);
602 xx1.
set(t * pabsz.x(), t * pabsz.y() , pz);
617 xx2.
set(r2, 0. , 0.);
620 xx2.
set(t * pabsz.x(), t * pabsz.y() , 0.);
629 xx.
set(p.
x(), p.
y(), pz);
662 G4int phiareacode = GetAreaCodeInPhi(xx);
667 if ((phiareacode &
sAxisMin) == sAxisMin) {
670 if (isoutsideinphi) isoutside =
true;
675 if (isoutsideinphi) isoutside =
true;
687 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
689 }
else if (xx.
z() >
fAxisMax[zaxis] - ctol) {
692 if (areacode & sBoundary) areacode |=
sCorner;
695 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
703 areacode = tmpareacode;
704 }
else if ((areacode & sBoundary) !=
sBoundary) {
712 G4int phiareacode = GetAreaCodeInPhi(xx,
false);
728 if (phiareacode == sAxisMin) {
731 if (areacode & sBoundary) areacode |=
sCorner;
734 }
else if (phiareacode == sAxisMax) {
737 if (areacode & sBoundary) areacode |=
sCorner;
744 if ((areacode & sBoundary) !=
sBoundary) {
750 std::ostringstream message;
751 message <<
"Feature NOT implemented !" <<
G4endl
753 <<
" fAxis[1] = " <<
fAxis[1];
790 areacode = tmpareacode;
810 void G4TwistTubsHypeSide::SetCorners(
826 for (i=0; i<2; i++) {
828 : EndInnerRadius[i]);
837 x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
838 y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
843 x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
844 y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
849 x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
850 y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
855 x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
856 y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
861 std::ostringstream message;
862 message <<
"Feature NOT implemented !" <<
G4endl
864 <<
" fAxis[1] = " <<
fAxis[1];
874 void G4TwistTubsHypeSide::SetCorners()
878 "Method NOT implemented !");
884 void G4TwistTubsHypeSide::SetBoundaries()
895 direction = direction.
unit();
901 direction = direction.
unit();
907 direction = direction.
unit();
913 direction = direction.
unit();
917 std::ostringstream message;
918 message <<
"Feature NOT implemented !" <<
G4endl
920 <<
" fAxis[1] = " <<
fAxis[1];
921 G4Exception(
"G4TwistTubsHypeSide::SetBoundaries()",
945 for ( i = 0 ; i<
n ; i++ ) {
949 for ( j = 0 ; j<k ; j++ )
951 nnode =
GetNode(i,j,k,n,iside) ;
957 x = xmin + j*(xmax-xmin)/(k-1) ;
959 x = xmax - j*(xmax-xmin)/(k-1) ;
964 xyz[nnode][0] = p.
x() ;
965 xyz[nnode][1] = p.
y() ;
966 xyz[nnode][2] = p.
z() ;
968 if ( i<n-1 && j<k-1 ) {
970 nface =
GetFace(i,j,k,n,iside) ;
static const G4int sAxisZ
void set(double x, double y, double z)
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
static const G4int sC0Min1Max
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
static const G4int sAxisPhi
virtual G4double GetBoundaryMin(G4double phi)
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
G4bool IsBoundary(G4int areacode, G4bool testbitmode=false) const
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
virtual G4double GetBoundaryMax(G4double phi)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
static const G4int sC0Min1Min
G4GLOB_DLL std::ostream G4cout
G4bool IsOutside(G4int areacode) const
virtual EInside Inside(const G4ThreeVector &gp)
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
G4double GetDistance(G4int i) const
static const G4int sAxis1
static const G4int sC0Max1Max
static const G4int sBoundary
G4double GetRadialTolerance() const
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)
G4bool IsValid(G4int i) const
static const G4int sAxis0
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisMin
static const G4int sInside
virtual G4double GetRhoAtPZ(const G4ThreeVector &p, G4bool isglobal=false) const
G4TwistTubsHypeSide(const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, const G4int handedness, const G4double kappa, const G4double tanstereo, const G4double r0, const EAxis axis0=kPhi, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const G4int sCorner
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sAxisMax
G4ThreeVector GetXX(G4int i) 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)
virtual ~G4TwistTubsHypeSide()
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
static const G4int sC0Max1Min
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
CurrentStatus fCurStatWithV
static G4GeometryTolerance * GetInstance()
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)