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 ;
void set(double x, double y, double z)
virtual G4ThreeVector SurfacePoint(G4double phi, G4double u, G4bool isGlobal=false)
static const G4double kInfinity
G4double GetDistance(G4int i) const
G4ThreeVector ComputeGlobalPoint(const G4ThreeVector &lp) const
G4int GetAreacode(G4int i) const
G4double DistanceToPlaneWithV(const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
static const G4int sOutside
virtual G4int GetAreaCode(const G4ThreeVector &xx, G4bool withTol=true)
G4GLOB_DLL std::ostream G4cout
G4bool IsValid(G4int i) const
void SetCurrentStatus(G4int i, G4ThreeVector &xx, G4double &dist, G4int &areacode, G4bool &isvalid, G4int nxx, EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u)
G4bool IsOutside(G4int areacode) const
static const G4double factor
G4ThreeVector GetXX(G4int i) const
G4ThreeVector ComputeLocalPoint(const G4ThreeVector &gp) const
G4bool EqualIntersection(const Intersection &a, const Intersection &b)
G4ThreeVector ComputeLocalDirection(const G4ThreeVector &gp) const
G4bool IsInside(G4int areacode, G4bool testbitmode=false) const
G4bool DistanceSort(const Intersection &a, const Intersection &b)
G4int FindRoots(G4double *op, G4int degree, G4double *zeror, G4double *zeroi)
CurrentStatus fCurStatWithV
G4ThreeVector NormAng(G4double phi, G4double u)
void ResetfDone(EValidate validate, const G4ThreeVector *p, const G4ThreeVector *v=0)