125 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
126 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
127 fAngleSide(0.), fDx4plus2(0.), fDx4minus2(0.), fDx3plus1(0.), fDx3minus1(0.),
128 fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), fa2md2(0.), fdeltaX(0.),
211 G4bool IsParallel = false ;
212 G4bool IsConverged = false ;
255 G4bool tmpisvalid = false ;
257 std::vector<Intersection> xbuf ;
272 if ( std::fabs(p.
z()) <= L )
279 *(v.
y()*std::cos(phi) - v.
x()*std::sin(phi))))
289 xbuf.push_back(xbuftmp) ;
298 areacode[0], isvalid[0],
299 0, validate, &gp, &gv);
338 + 4*
fDy1*(9*phixz - 9*fTAlph*phiyz - 56*
fDz*fTAlph*v.
x()
346 G4cout <<
"coef = " << c[0] <<
" "
359 for (
G4int i = 0 ; i<num ; i++ )
364 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
366 phi = std::fmod(srd[i] , pihalf) ;
378 xbuf.push_back(xbuftmp) ;
381 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
396 for (
size_t k = 0 ; k<xbuf.size() ; k++ )
399 G4cout <<
"Solution " << k <<
" : "
400 <<
"reconstructed phiR = " << xbuf[k].phi
401 <<
", uR = " << xbuf[k].u <<
G4endl ;
407 IsConverged = false ;
409 for (
G4int i = 1 ; i<maxint ; i++ )
412 surfacenormal =
NormAng(phi,u) ;
415 deltaX = ( tmpxx - xxonsurface ).mag() ;
416 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
428 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist
429 <<
", " << deltaX <<
G4endl ;
437 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
440 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
444 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0 ; }
447 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
461 if (tmpdist >= 0) tmpisvalid =
true;
469 if (tmpdist >= 0) { tmpisvalid =
true; }
474 G4Exception(
"G4TwistTrapAlphaSide::DistanceToSurface()",
476 "Feature NOT implemented !");
488 xbuf[k].distance = tmpdist ;
489 xbuf[k].areacode = tmpareacode ;
490 xbuf[k].isvalid = tmpisvalid ;
497 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
509 G4int nxxtmp = xbuf.size() ;
511 if ( nxxtmp<2 || IsParallel )
527 xbuf.push_back(xbuftmp) ;
542 xbuf.push_back(xbuftmp) ;
544 for (
size_t k = nxxtmp ; k<xbuf.size() ; k++ )
548 G4cout <<
"Solution " << k <<
" : "
549 <<
"reconstructed phiR = " << xbuf[k].phi
550 <<
", uR = " << xbuf[k].u <<
G4endl ;
556 IsConverged = false ;
558 for (
G4int i = 1 ; i<maxint ; i++ )
561 surfacenormal =
NormAng(phi,u) ;
563 deltaX = ( tmpxx - xxonsurface ).mag();
564 theta = std::fabs(std::acos(v*surfacenormal) - pihalf);
575 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist
576 <<
", " << deltaX << G4endl
577 <<
"X = " << tmpxx <<
G4endl ;
584 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
587 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
591 if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0; }
594 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
608 if (tmpdist >= 0) { tmpisvalid =
true; }
616 if (tmpdist >= 0) { tmpisvalid =
true; }
621 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
623 "Feature NOT implemented !");
635 xbuf[k].distance = tmpdist ;
636 xbuf[k].areacode = tmpareacode ;
637 xbuf[k].isvalid = tmpisvalid ;
650 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
656 for (
size_t i = 0 ; i<xbuf.size() ; i++ )
658 distance[i] = xbuf[i].distance;
660 areacode[i] = xbuf[i].areacode ;
661 isvalid[i] = xbuf[i].isvalid ;
664 isvalid[i], nxx, validate, &gp, &gv);
666 G4cout <<
"element Nr. " << i
667 <<
", local Intersection = " << xbuf[i].xx
668 <<
", distance = " << xbuf[i].distance
669 <<
", u = " << xbuf[i].u
670 <<
", phi = " << xbuf[i].phi
671 <<
", isvalid = " << xbuf[i].isvalid
679 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
680 for (
G4int k= 0 ; k< nxx ; k++ )
682 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
738 for (
G4int i = 1 ; i<20 ; i++ )
741 surfacenormal =
NormAng(phiR,uR) ;
743 deltaX = ( xx - xxonsurface ).mag() ;
746 G4cout <<
"i = " << i <<
", distance = " << distance[0]
747 <<
", " << deltaX <<
G4endl
748 <<
"X = " << xx <<
G4endl ;
755 if ( deltaX <= ctol ) { break ; }
762 if ( phiR > halfphi ) { phiR = halfphi ; }
763 if ( phiR < -halfphi ) { phiR = -halfphi ; }
764 if ( uR > uMax ) { uR = uMax ; }
765 if ( uR < -uMax ) { uR = -uMax ; }
768 distance[0] = ( p -
xx ).mag() ;
769 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
774 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
783 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
813 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
814 G4cout <<
"Intervall is " << fYAxisMin <<
" to " << fYAxisMax <<
G4endl ;
829 if (yprime < fYAxisMin + ctol)
832 if (yprime <= fYAxisMin - ctol) { isoutside =
true; }
835 else if (yprime > fYAxisMax - ctol)
838 if (yprime >= fYAxisMax + ctol) { isoutside =
true; }
852 if (xx.
z() <=
fAxisMin[zaxis] - ctol) { isoutside =
true; }
854 else if (xx.
z() >
fAxisMax[zaxis] - ctol)
862 if (xx.
z() >=
fAxisMax[zaxis] + ctol) { isoutside =
true; }
871 areacode = tmpareacode;
883 if (yprime < fYAxisMin )
887 else if (yprime > fYAxisMax)
922 "Feature NOT implemented !");
992 "Method NOT implemented !");
1010 direction = direction.
unit();
1016 direction = direction.
unit();
1022 direction = direction.
unit();
1028 direction = direction.
unit();
1037 "Feature NOT implemented !");
1060 -
fPhiTwist*(fTAlph*p.
x() + p.
y())))*std::cos(phi)
1065 /
fDy1 - 4*std::sin(phi)))
1066 *(std::fabs(((
fa1md1 + 4*
fDy1*fTAlph)*std::cos(phi))
1067 /
fDy1 - 4*std::sin(phi)))
1068 + (std::fabs(4*std::cos(phi)
1070 * (std::fabs(4*std::cos(phi)
1130 for (
G4int i = 0 ; i<
n ; i++ )
1132 z = -
fDz+i*(2.*
fDz)/(n-1) ;
1136 for (
G4int j = 0 ; j<k ; j++ )
1138 nnode =
GetNode(i,j,k,n,iside) ;
1139 u = -b/2 +j*b/(k-1) ;
1142 xyz[nnode][0] = p.
x() ;
1143 xyz[nnode][1] = p.
y() ;
1144 xyz[nnode][2] = p.
z() ;
1146 if ( i<n-1 && j<k-1 )
1148 nface =
GetFace(i,j,k,n,iside) ;
1150 * (
GetNode(i ,j ,k,n,iside)+1) ;
1152 * (
GetNode(i ,j+1,k,n,iside)+1) ;
1154 * (
GetNode(i+1,j+1,k,n,iside)+1) ;
1156 * (
GetNode(i+1,j ,k,n,iside)+1) ;
virtual ~G4TwistTrapAlphaSide()
G4int GetAreacode(G4int i) const
static const G4int sC0Min1Max
static constexpr double L
void set(double x, double y, double z)
CurrentStatus fCurStatWithV
static const G4int sInside
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
std::vector< ExP01TrackerHit * > a
static const G4double kInfinity
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
virtual void SetCorners()
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
virtual G4double GetBoundaryMin(G4double phi)
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)
static double normal(HepRandomEngine *eptr)
G4double GetValueB(G4double phi)
virtual G4double GetBoundaryMax(G4double phi)
G4bool IsOutside(G4int areacode) const
G4int GetEdgeVisibility(G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
static const G4int sAxisZ
static const G4int sC0Min1Min
G4ThreeVector ProjectPoint(const G4ThreeVector &p, G4bool isglobal=false)
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
static const G4int sAxis1
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
G4int GetFace(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sC0Max1Max
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
void SetCorner(G4int areacode, G4double x, G4double y, G4double z)
G4ThreeVector NormAng(G4double phi, G4double u)
G4ThreeVector GetCorner(G4int areacode) const
G4bool DistanceSort(const Intersection &a, const Intersection &b)
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4double GetDistance(G4int i) const
static const G4int sAxis0
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
HepRotation & rotateZ(double delta)
G4GLOB_DLL std::ostream G4cout
static const G4int sOutside
static constexpr double pi
G4SurfCurNormal fCurrentNormal
static const G4int sAxisMin
static const G4int sBoundary
G4TwistTrapAlphaSide(const G4String &name, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph, G4double AngleSide)
G4int GetNode(G4int i, G4int j, G4int m, G4int n, G4int iside)
static const G4int sC0Max1Min
virtual void SetBoundaries()
static const G4int sAxisMax
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector GetXX(G4int i) const
HepRotation inverse() const
static const G4int sAxisY
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
G4bool IsValid(G4int i) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const