84 message <<
"TwistedTrapBoxSide is not used as a the side of a box: "
87 G4Exception(
"G4TwistBoxSide::G4TwistBoxSide()",
"GeomSolids0002",
135 :
G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.),
136 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.),
137 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.),
138 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.),
212 G4bool IsParallel = false ;
213 G4bool IsConverged = false ;
255 G4bool tmpisvalid = false ;
257 std::vector<Intersection> xbuf ;
276 + (
fTAlph*v.
x() + v.
y())*std::sin(phi)));
282 * (v.
y()*std::cos(phi) - v.
x()*std::sin(phi))) / q;
291 xbuf.push_back(xbuftmp) ;
302 areacode[0], isvalid[0],
303 0, validate, &gp, &gv);
330 G4cout <<
"coef = " << c[0] <<
" "
344 for (
G4int i = 0 ; i<num ; i++ ) {
347 G4cout <<
"Solution " << i <<
" : " << srd[i] <<
G4endl ;
349 phi = std::fmod(srd[i] , pihalf) ;
363 xbuf.push_back(xbuftmp) ;
366 G4cout <<
"solution " << i <<
" = " << phi <<
" , " << u <<
G4endl ;
385 for (
size_t k = 0 ; k<xbuf.size() ; k++ ) {
388 G4cout <<
"Solution " << k <<
" : "
389 <<
"reconstructed phiR = " << xbuf[k].phi
390 <<
", uR = " << xbuf[k].u <<
G4endl ;
396 IsConverged = false ;
398 for (
G4int i = 1 ; i<maxint ; i++ ) {
401 surfacenormal =
NormAng(phi,u) ;
403 deltaX = ( tmpxx - xxonsurface ).mag() ;
404 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
405 if ( theta < 0.001 ) {
414 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
421 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
424 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
430 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
433 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
445 if (tmpdist >= 0) tmpisvalid =
true;
450 if (tmpdist >= 0) tmpisvalid =
true;
455 "Feature NOT implemented !");
467 xbuf[k].distance = tmpdist ;
468 xbuf[k].areacode = tmpareacode ;
469 xbuf[k].isvalid = tmpisvalid ;
478 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
484 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
489 G4int nxxtmp = xbuf.size() ;
491 if ( nxxtmp<2 || IsParallel ) {
509 xbuf.push_back(xbuftmp) ;
525 xbuf.push_back(xbuftmp) ;
527 for (
size_t k = nxxtmp ; k<xbuf.size() ; k++ ) {
530 G4cout <<
"Solution " << k <<
" : "
531 <<
"reconstructed phiR = " << xbuf[k].phi
532 <<
", uR = " << xbuf[k].u <<
G4endl ;
538 IsConverged = false ;
540 for (
G4int i = 1 ; i<maxint ; i++ ) {
543 surfacenormal =
NormAng(phi,u) ;
545 deltaX = ( tmpxx - xxonsurface ).mag() ;
546 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
547 if ( theta < 0.001 ) {
555 G4cout <<
"Step i = " << i <<
", distance = " << tmpdist <<
", " << deltaX <<
G4endl ;
562 G4cout <<
"approximated phi = " << phi <<
", u = " << u <<
G4endl ;
565 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
571 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ;
574 G4cout <<
"refined solution " << phi <<
" , " << u <<
G4endl ;
586 if (tmpdist >= 0) tmpisvalid =
true;
591 if (tmpdist >= 0) tmpisvalid =
true;
594 G4Exception(
"G4TwistedBoxSide::DistanceToSurface()",
596 "Feature NOT implemented !");
608 xbuf[k].distance = tmpdist ;
609 xbuf[k].areacode = tmpareacode ;
610 xbuf[k].isvalid = tmpisvalid ;
623 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() ,
EqualIntersection ) , xbuf.end() ) ;
626 G4cout << G4endl <<
"list xbuf after sorting : " <<
G4endl ;
632 for (
size_t i = 0 ; i<xbuf.size() ; i++ ) {
634 distance[i] = xbuf[i].distance;
636 areacode[i] = xbuf[i].areacode ;
637 isvalid[i] = xbuf[i].isvalid ;
640 isvalid[i], nxx, validate, &gp, &gv);
643 G4cout <<
"element Nr. " << i
644 <<
", local Intersection = " << xbuf[i].xx
645 <<
", distance = " << xbuf[i].distance
646 <<
", u = " << xbuf[i].u
647 <<
", phi = " << xbuf[i].phi
648 <<
", isvalid = " << xbuf[i].isvalid
657 G4cout << nxx <<
" possible physical solutions found" <<
G4endl ;
658 for (
G4int k= 0 ; k< nxx ; k++ ) {
659 G4cout <<
"global intersection Point found: " << gxx[k] <<
G4endl ;
715 for (
G4int i = 1 ; i<maxint ; i++ ) {
718 surfacenormal =
NormAng(phiR,uR) ;
720 deltaX = ( xx - xxonsurface ).mag() ;
723 G4cout <<
"i = " << i <<
", distance = " << distance[0] <<
", " << deltaX <<
G4endl ;
730 if ( deltaX <= ctol ) { break ; }
739 if ( phiR > halfphi ) phiR = halfphi ;
740 if ( phiR < -halfphi ) phiR = -halfphi ;
741 if ( uR > uMax ) uR = uMax ;
742 if ( uR < -uMax ) uR = -uMax ;
745 distance[0] = ( p -
xx ).mag() ;
746 if ( distance[0] <= ctol ) { distance[0] = 0 ; }
751 G4cout <<
"refined solution " << phiR <<
" , " << uR <<
" , " <<
G4endl ;
760 G4cout <<
"intersection Point found: " << gxx[0] <<
G4endl ;
792 G4cout <<
"GetAreaCode: yprime = " << yprime <<
G4endl ;
793 G4cout <<
"Intervall is " << fYAxisMin <<
" to " << fYAxisMax <<
G4endl ;
808 if (yprime < fYAxisMin + ctol) {
810 if (yprime <= fYAxisMin - ctol) isoutside =
true;
812 }
else if (yprime > fYAxisMax - ctol) {
814 if (yprime >= fYAxisMax + ctol) isoutside =
true;
824 if (xx.
z() <=
fAxisMin[zaxis] - ctol) isoutside =
true;
826 }
else if (xx.
z() >
fAxisMax[zaxis] - ctol) {
831 if (xx.
z() >=
fAxisMax[zaxis] + ctol) isoutside =
true;
839 areacode = tmpareacode;
848 if (yprime < fYAxisMin ) {
850 }
else if (yprime > fYAxisMax) {
875 "Feature NOT implemented !");
928 "Method NOT implemented !");
946 direction = direction.
unit();
952 direction = direction.
unit();
958 direction = direction.
unit();
964 direction = direction.
unit();
972 "Feature NOT implemented !");
1033 for ( i = 0 ; i<
n ; i++ ) {
1035 z = -
fDz+i*(2.*
fDz)/(n-1) ;
1039 for ( j = 0 ; j<k ; j++ ) {
1041 nnode =
GetNode(i,j,k,n,iside) ;
1042 u = -b/2 +j*b/(k-1) ;
1045 xyz[nnode][0] = p.
x() ;
1046 xyz[nnode][1] = p.
y() ;
1047 xyz[nnode][2] = p.
z() ;
1049 if ( i<n-1 && j<k-1 ) {
1051 nface =
GetFace(i,j,k,n,iside) ;
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
G4ThreeVector NormAng(G4double phi, G4double u)
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
std::vector< ExP01TrackerHit * > a
G4ThreeVector ProjectPoint(const G4ThreeVector &p, G4bool isglobal=false)
static const G4double kInfinity
virtual void SetCorners()
virtual void SetBoundaries()
virtual void SetBoundary(const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
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)
static double normal(HepRandomEngine *eptr)
virtual ~G4TwistBoxSide()
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
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
virtual G4String GetName() const
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
G4double GetValueB(G4double phi)
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 G4double GetBoundaryMax(G4double phi)
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 GetCorner(G4int areacode) const
G4bool DistanceSort(const Intersection &a, const Intersection &b)
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
virtual G4ThreeVector GetNormal(const G4ThreeVector &xx, G4bool isGlobal=false)
static constexpr double pi
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
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
G4double DistanceToPlane(const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4ThreeVector GetXX(G4int i) const
HepRotation inverse() const
virtual G4int DistanceToSurface(const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
static const G4int sAxisY
G4TwistBoxSide(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)
G4ThreeVector ComputeGlobalDirection(const G4ThreeVector &lp) const
virtual void GetFacets(G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
G4bool IsValid(G4int i) const
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const