62 axis0min, axis1min, axis0max, axis1max),
66 G4Exception(
"G4TwistTubsSide::G4TwistTubsSide()",
"GeomSolids0002",
101 SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
226 for (i=0; i<2; i++) {
246 if ((absvz <
DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x()) <
DBL_MIN)) {
251 isvalid[0], 0, validate, &gp, &gv);
271 distance[0] = - c / b;
272 xx[0] = p + distance[0]*v;
278 if (distance[0] >= 0) isvalid[0] =
true;
283 if (distance[0] >= 0) isvalid[0] =
true;
289 if (distance[0] >= 0) isvalid[0] =
true;
293 areacode[0], isvalid[0],
294 0, validate, &gp, &gv);
300 isvalid[0], 1, validate, &gp, &gv);
312 isvalid[0], 0, validate, &gp, &gv);
324 G4bool tmpisvalid[2] = {
false,
false};
327 for (i=0; i<2; i++) {
333 if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
335 tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
337 tmpdist[i] = factor * bminusD;
341 tmpxx[i] = p + tmpdist[i]*v;
346 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
352 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
357 if (tmpxx[i].
x() > 0) {
359 if (tmpdist[i] >= 0) tmpisvalid[i] =
true;
368 if (tmpdist[0] <= tmpdist[1]) {
369 distance[0] = tmpdist[0];
370 distance[1] = tmpdist[1];
375 areacode[0] = tmpareacode[0];
376 areacode[1] = tmpareacode[1];
377 isvalid[0] = tmpisvalid[0];
378 isvalid[1] = tmpisvalid[1];
380 distance[0] = tmpdist[1];
381 distance[1] = tmpdist[0];
386 areacode[0] = tmpareacode[1];
387 areacode[1] = tmpareacode[0];
388 isvalid[0] = tmpisvalid[1];
389 isvalid[1] = tmpisvalid[0];
393 isvalid[0], 2, validate, &gp, &gv);
395 isvalid[1], 2, validate, &gp, &gv);
399 for (
G4int k=0; k<2; k++) {
400 if (!isvalid[k])
continue;
403 * xx[k].
z() , xx[k].
z());
404 G4double deltaY = (xx[k] - xxonsurface).mag();
411 for (l=0; l<maxcount; l++) {
414 surfacenormal, xx[k]);
415 deltaY = (xx[k] - xxonsurface).mag();
416 if (deltaY > lastdeltaY) {
425 xxonsurface.set(xx[k].
x(),
426 fKappa * std::fabs(xx[k].
x()) * xx[k].
z(),
430 std::ostringstream message;
431 message <<
"Exceeded maxloop count!" <<
G4endl
432 <<
" maxloop count " << maxcount;
433 G4Exception(
"G4TwistTubsFlatSide::DistanceToSurface(p,v)",
445 isvalid[0], 0, validate, &gp, &gv);
470 for (i=0; i<2; i++) {
491 for (i=0; i<2; i++) {
495 if ((gp - lastgxx[0]).mag() < halftol
496 || (gp - lastgxx[1]).mag() < halftol) {
508 if (p.getRho() == 0) {
548 }
else if (p.z() < C.z()) {
554 }
else if (p.z() < A.z()) {
563 for (i=0; i<2; i++) {
566 B = x0[i] + ((A.z() - x0[i].z()) / d[i].
z()) * d[i];
570 D = x0[i] + ((C.z() - x0[i].z()) / d[i].
z()) * d[i];
579 G4double test = (A.z() - C.z()) * parity * pside;
604 }
else if (test < 0) {
636 if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
637 xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad);
647 if (distToACB * distToCAD > 0 && distToACB < 0) {
653 if (distToACB * distToCAD > 0) {
656 if (distToACB <= distToCAD) {
657 distance[0] = distToACB;
660 distance[0] = distToCAD;
667 distance[0] = distToACB;
670 distance[0] = distToCAD;
709 if (distToanm * distTocmn > 0 && distToanm < 0) {
712 "Point p is behind the surfaces.");
716 if (std::fabs(distToanm) <= halftol) {
720 }
else if (std::fabs(distTocmn) <= halftol) {
726 if (distToanm <= distTocmn) {
771 if (xx.x() <
fAxisMin[xaxis] + ctol) {
773 if (xx.x() <=
fAxisMin[xaxis] - ctol) isoutside =
true;
775 }
else if (xx.x() >
fAxisMax[xaxis] - ctol) {
777 if (xx.x() >=
fAxisMax[xaxis] + ctol) isoutside =
true;
782 if (xx.z() <
fAxisMin[zaxis] + ctol) {
787 if (xx.z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
789 }
else if (xx.z() >
fAxisMax[zaxis] - ctol) {
794 if (xx.z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
802 areacode = tmpareacode;
813 }
else if (xx.x() >
fAxisMax[xaxis]) {
824 }
else if (xx.z() >
fAxisMax[zaxis]) {
838 "Feature NOT implemented !");
862 x = endInnerRad[zmin]*std::cos(endPhi[zmin]);
863 y = endInnerRad[zmin]*std::sin(endPhi[zmin]);
868 x = endOuterRad[zmin]*std::cos(endPhi[zmin]);
869 y = endOuterRad[zmin]*std::sin(endPhi[zmin]);
874 x = endOuterRad[zmax]*std::cos(endPhi[zmax]);
875 y = endOuterRad[zmax]*std::sin(endPhi[zmax]);
880 x = endInnerRad[zmax]*std::cos(endPhi[zmax]);
881 y = endInnerRad[zmax]*std::sin(endPhi[zmax]);
886 std::ostringstream message;
887 message <<
"Feature NOT implemented !" <<
G4endl
889 <<
" fAxis[1] = " <<
fAxis[1];
902 "Method NOT implemented !");
918 direction = direction.unit();
924 direction = direction.unit();
930 direction = direction.unit();
936 direction = direction.unit();
941 std::ostringstream message;
942 message <<
"Feature NOT implemented !" <<
G4endl
944 <<
" fAxis[1] = " <<
fAxis[1];
969 for ( i = 0 ; i<
n ; i++ )
974 for ( j = 0 ; j<k ; j++ ) {
976 nnode =
GetNode(i,j,k,n,iside) ;
982 x = xmin + j*(xmax-xmin)/(k-1) ;
984 x = xmax - j*(xmax-xmin)/(k-1) ;
989 xyz[nnode][0] = p.x() ;
990 xyz[nnode][1] = p.y() ;
991 xyz[nnode][2] = p.z() ;
993 if ( i<n-1 && j<k-1 ) {
995 nface =
GetFace(i,j,k,n,iside) ;
static const G4int sAxisZ
G4int GetAreacode(G4int i) const
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
static const G4int sC0Min1Max
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
static const G4double kInfinity
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
CLHEP::HepRotation G4RotationMatrix
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4SurfCurNormal fCurrentNormal
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector GetCorner(G4int areacode) const
static const G4int sOutside
virtual G4ThreeVector GetBoundaryAtPZ(G4int areacode, const G4ThreeVector &p) const
virtual ~G4TwistTubsSide()
double B(double temperature)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
static const G4int sAxisX
static double normal(HepRandomEngine *eptr)
static const G4int sC0Min1Min
virtual void SetCorners()
virtual void SetBoundaries()
virtual void GetBoundaryParameters(const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
double A(double temperature)
G4bool IsOutside(G4int areacode) const
G4double GetDistance(G4int i) const
static const G4int sAxis1
static const G4int sC0Max1Max
static const G4int sBoundary
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)
G4bool IsValid(G4int i) const
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
static const G4int sAxis0
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
static const G4int sAxisMin
virtual G4double GetBoundaryMax(G4double phi)
virtual G4ThreeVector SurfacePoint(G4double, G4double, G4bool isGlobal=false)
static const G4int sInside
virtual G4double GetBoundaryMin(G4double phi)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static const G4int sCorner
static const G4double factor
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sAxisMax
const G4double x[NPOINTSGL]
G4ThreeVector GetXX(G4int i) const
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
G4TwistTubsSide(const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0=kXAxis, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
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)
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 G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
CurrentStatus fCurStatWithV
virtual G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &A, const G4ThreeVector &B, const G4ThreeVector &C, const G4ThreeVector &D, const G4int parity, G4ThreeVector &xx, G4ThreeVector &n)
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)