G4TwistBoxSide.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4TwistBoxSide.cc 67011 2013-01-29 16:17:41Z gcosmo $
00028 //
00029 // 
00030 // --------------------------------------------------------------------
00031 // GEANT 4 class source file
00032 //
00033 //
00034 // G4TwistBoxSide.cc
00035 //
00036 // Author:
00037 //
00038 //   18/03/2005 - O.Link (Oliver.Link@cern.ch)
00039 //
00040 // --------------------------------------------------------------------
00041 
00042 #include <cmath>
00043 
00044 #include "G4TwistBoxSide.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4JTPolynomialSolver.hh"
00047 
00048 //=====================================================================
00049 //* constructors ------------------------------------------------------
00050 
00051 G4TwistBoxSide::G4TwistBoxSide(const G4String     &name,
00052                            G4double      PhiTwist,    // twist angle
00053                            G4double      pDz,         // half z lenght
00054                            G4double      pTheta,      // direction between end planes
00055                            G4double      pPhi,        // defined by polar and azimutal angles.
00056                            G4double      pDy1,        // half y length at -pDz
00057                            G4double      pDx1,        // half x length at -pDz,-pDy
00058                            G4double      pDx2,        // half x length at -pDz,+pDy
00059                            G4double      pDy2,        // half y length at +pDz
00060                            G4double      pDx3,        // half x length at +pDz,-pDy
00061                            G4double      pDx4,        // half x length at +pDz,+pDy
00062                            G4double      pAlph,       // tilt angle at +pDz
00063                            G4double      AngleSide    // parity
00064                                                ) : G4VTwistSurface(name)
00065 {  
00066   
00067                  
00068   fAxis[0]    = kYAxis; // in local coordinate system
00069   fAxis[1]    = kZAxis;
00070   fAxisMin[0] = -kInfinity ;  // Y Axis boundary
00071   fAxisMax[0] = kInfinity ;   //   depends on z !!
00072   fAxisMin[1] = -pDz ;      // Z Axis boundary
00073   fAxisMax[1] = pDz ;
00074   
00075   fDx1  = pDx1 ;
00076   fDx2  = pDx2 ;  // box
00077   fDx3  = pDx3 ;
00078   fDx4  = pDx4 ;  // box
00079 
00080   // this is an overhead. But the parameter naming scheme fits to the other surfaces.
00081  
00082   if ( ! (fDx1 == fDx2 && fDx3 == fDx4 ) ) {
00083     std::ostringstream message;
00084     message << "TwistedTrapBoxSide is not used as a the side of a box: "
00085             << GetName() << G4endl
00086             << "        Not a box !";
00087     G4Exception("G4TwistBoxSide::G4TwistBoxSide()", "GeomSolids0002",
00088                 FatalException, message);
00089   }
00090  
00091   fDy1   = pDy1 ;
00092   fDy2   = pDy2 ;
00093 
00094   fDz   = pDz ;
00095 
00096   fAlph = pAlph  ;
00097   fTAlph = std::tan(fAlph) ;
00098 
00099   fTheta = pTheta ;
00100   fPhi   = pPhi ;
00101 
00102  // precalculate frequently used parameters
00103   fDx4plus2  = fDx4 + fDx2 ;
00104   fDx4minus2 = fDx4 - fDx2 ;
00105   fDx3plus1  = fDx3 + fDx1 ; 
00106   fDx3minus1 = fDx3 - fDx1 ;
00107   fDy2plus1  = fDy2 + fDy1 ;
00108   fDy2minus1 = fDy2 - fDy1 ;
00109 
00110   fa1md1 = 2*fDx2 - 2*fDx1  ; 
00111   fa2md2 = 2*fDx4 - 2*fDx3 ;
00112 
00113 
00114   fPhiTwist = PhiTwist ;     // dphi
00115   fAngleSide = AngleSide ;  // 0,90,180,270 deg
00116 
00117   fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi)  ;  // dx in surface equation
00118   fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi)  ;  // dy in surface equation
00119   
00120   fRot.rotateZ( AngleSide ) ; 
00121   
00122   fTrans.set(0, 0, 0);  // No Translation
00123   fIsValidNorm = false;
00124   
00125   SetCorners() ;
00126   SetBoundaries() ;
00127 
00128 }
00129 
00130 
00131 //=====================================================================
00132 //* Fake default constructor ------------------------------------------
00133 
00134 G4TwistBoxSide::G4TwistBoxSide( __void__& a )
00135   : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.), 
00136     fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.), 
00137     fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.), 
00138     fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), 
00139     fa2md2(0.)
00140 {
00141 }
00142 
00143 
00144 //=====================================================================
00145 //* destructor --------------------------------------------------------
00146 
00147 G4TwistBoxSide::~G4TwistBoxSide()
00148 {
00149 }
00150 
00151 //=====================================================================
00152 //* GetNormal ---------------------------------------------------------
00153 
00154 G4ThreeVector G4TwistBoxSide::GetNormal(const G4ThreeVector &tmpxx, 
00155                                                 G4bool isGlobal) 
00156 {
00157    // GetNormal returns a normal vector at a surface (or very close
00158    // to surface) point at tmpxx.
00159    // If isGlobal=true, it returns the normal in global coordinate.
00160    //
00161 
00162    G4ThreeVector xx;
00163    if (isGlobal) {
00164       xx = ComputeLocalPoint(tmpxx);
00165       if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
00166          return ComputeGlobalDirection(fCurrentNormal.normal);
00167       }
00168    } else {
00169       xx = tmpxx;
00170       if (xx == fCurrentNormal.p) {
00171          return fCurrentNormal.normal;
00172       }
00173    }
00174 
00175    G4double phi ;
00176    G4double u ;
00177 
00178    GetPhiUAtX(xx,phi,u) ;   // phi,u for point xx close to surface 
00179 
00180    G4ThreeVector normal =  NormAng(phi,u) ;  // the normal vector at phi,u
00181 
00182 #ifdef G4TWISTDEBUG
00183    G4cout  << "normal vector = " << normal << G4endl ;
00184    G4cout << "phi = " << phi << " , u = " << u << G4endl ;
00185 #endif
00186 
00187    //    normal = normal/normal.mag() ;
00188 
00189    if (isGlobal) {
00190       fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
00191    } else {
00192       fCurrentNormal.normal = normal.unit();
00193    }
00194    return fCurrentNormal.normal;
00195 }
00196 
00197 //=====================================================================
00198 //* DistanceToSurface -------------------------------------------------
00199 
00200 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector &gp,
00201                                         const G4ThreeVector &gv,
00202                                               G4ThreeVector  gxx[],
00203                                               G4double       distance[],
00204                                               G4int          areacode[],
00205                                               G4bool         isvalid[],
00206                                               EValidate      validate)
00207 {
00208 
00209   static const G4double ctol = 0.5 * kCarTolerance;
00210   static const G4double pihalf = pi/2 ;
00211 
00212   G4bool IsParallel = false ;
00213   G4bool IsConverged =  false ;
00214 
00215   G4int nxx = 0 ;  // number of physical solutions
00216 
00217   fCurStatWithV.ResetfDone(validate, &gp, &gv);
00218 
00219   if (fCurStatWithV.IsDone()) {
00220     G4int i;
00221     for (i=0; i<fCurStatWithV.GetNXX(); i++) {
00222       gxx[i] = fCurStatWithV.GetXX(i);
00223       distance[i] = fCurStatWithV.GetDistance(i);
00224       areacode[i] = fCurStatWithV.GetAreacode(i);
00225       isvalid[i]  = fCurStatWithV.IsValid(i);
00226     }
00227     return fCurStatWithV.GetNXX();
00228   } else {
00229    
00230    // initialize
00231     G4int i;
00232     for (i=0; i<G4VSURFACENXX ; i++) {
00233       distance[i] = kInfinity;
00234       areacode[i] = sOutside;
00235       isvalid[i]  = false;
00236       gxx[i].set(kInfinity, kInfinity, kInfinity);
00237     }
00238   }
00239 
00240   G4ThreeVector p = ComputeLocalPoint(gp);
00241   G4ThreeVector v = ComputeLocalDirection(gv);
00242   
00243 #ifdef G4TWISTDEBUG
00244   G4cout << "Local point p = " << p << G4endl ;
00245   G4cout << "Local direction v = " << v << G4endl ; 
00246 #endif
00247 
00248   G4double phi=0.,u=0.,q=0.;  // parameters
00249 
00250   // temporary variables
00251 
00252   G4double      tmpdist = kInfinity ;
00253   G4ThreeVector tmpxx;
00254   G4int         tmpareacode = sOutside ;
00255   G4bool        tmpisvalid  = false ;
00256 
00257   std::vector<Intersection> xbuf ;
00258   Intersection xbuftmp ;
00259   
00260   // prepare some variables for the intersection finder
00261 
00262   G4double L = 2*fDz ;
00263 
00264   G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
00265   G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
00266 
00267   // special case vz = 0
00268 
00269   if ( v.z() == 0. ) {         
00270 
00271     if ( (std::fabs(p.z()) <= L) && fPhiTwist ) {  // intersection possible in z
00272       
00273       phi = p.z() * fPhiTwist / L ;  // phi is determined by the z-position 
00274 
00275       q = (2.* fPhiTwist*((v.x() - fTAlph*v.y())*std::cos(phi)
00276                              + (fTAlph*v.x() + v.y())*std::sin(phi)));
00277       
00278       if (q) {
00279         u = (2*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x()
00280                 + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y())
00281              + (fDx4plus2*fPhiTwist + 2*fDx4minus2*phi)
00282              * (v.y()*std::cos(phi) - v.x()*std::sin(phi))) / q;
00283       }
00284 
00285       xbuftmp.phi = phi ;
00286       xbuftmp.u = u ;
00287       xbuftmp.areacode = sOutside ;
00288       xbuftmp.distance = kInfinity ;
00289       xbuftmp.isvalid = false ;
00290 
00291       xbuf.push_back(xbuftmp) ;  // store it to xbuf
00292 
00293     }
00294 
00295     else {                        // no intersection possible
00296 
00297       distance[0] = kInfinity;
00298       gxx[0].set(kInfinity,kInfinity,kInfinity);
00299       isvalid[0] = false ;
00300       areacode[0] = sOutside ;
00301       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
00302                                      areacode[0], isvalid[0],
00303                                      0, validate, &gp, &gv);
00304       
00305       return 0;
00306 
00307 
00308     }  // end std::fabs(p.z() <= L 
00309 
00310   } // end v.z() == 0
00311   
00312 
00313   // general solution for non-zero vz
00314 
00315   else {
00316 
00317     G4double c[8],srd[7],si[7] ;  
00318 
00319     c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz + fDx4plus2*fPhiTwist*v.z()) ;
00320     c[6] = 28800*(phiyz + 2*fDz*v.x() - (fdeltaX + fDx4minus2)*v.z() + fTAlph*(phixz - 2*fDz*v.y() + fdeltaY*v.z())) ;
00321     c[5] = -1200*(10*phixz - 48*fDz*v.y() + 24*fdeltaY*v.z() + fDx4plus2*fPhiTwist*v.z() - 2*fTAlph*(5*phiyz + 24*fDz*v.x() - 12*fdeltaX*v.z())) ;
00322     c[4] = -2400*(phiyz + 10*fDz*v.x() - 5*fdeltaX*v.z() + fDx4minus2*v.z() + fTAlph*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())) ;
00323     c[3] = 24*(2*phixz - 200*fDz*v.y() + 100*fdeltaY*v.z() - fDx4plus2*fPhiTwist*v.z() - 2*fTAlph*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z())) ;
00324     c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y() + 3*(fdeltaX + fDx4minus2 - fdeltaY*fTAlph)*v.z()) ;
00325     c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x() - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()) ;
00326     c[0] = 36*(2* fDz*(v.x() - fTAlph*v.y()) - fdeltaX*v.z() + fdeltaY*fTAlph*v.z()) ;
00327 
00328 
00329 #ifdef G4TWISTDEBUG
00330     G4cout << "coef = " << c[0] << " " 
00331            <<  c[1] << " "  
00332            <<  c[2] << " "  
00333            <<  c[3] << " "  
00334            <<  c[4] << " "  
00335            <<  c[5] << " "  
00336            <<  c[6] << " "  
00337            <<  c[7] << G4endl ;
00338 #endif    
00339 
00340     G4JTPolynomialSolver trapEq ;
00341     G4int num = trapEq.FindRoots(c,7,srd,si);
00342   
00343 
00344     for (G4int i = 0 ; i<num ; i++ ) {  // loop over all mathematical solutions
00345       if ( (si[i]==0.0) && fPhiTwist ) {  // only real solutions
00346 #ifdef G4TWISTDEBUG
00347         G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
00348 #endif
00349         phi = std::fmod(srd[i] , pihalf)  ;
00350 
00351         u   = (2*phiyz + 4*fDz*phi*v.y() - 2*fdeltaY*phi*v.z()
00352              - fDx4plus2*fPhiTwist*v.z()*std::sin(phi)
00353              - 2*fDx4minus2*phi*v.z()*std::sin(phi))
00354              /(2*fPhiTwist*v.z()*std::cos(phi)
00355                + 2*fPhiTwist*fTAlph*v.z()*std::sin(phi)) ;
00356 
00357         xbuftmp.phi = phi ;
00358         xbuftmp.u = u ;
00359         xbuftmp.areacode = sOutside ;
00360         xbuftmp.distance = kInfinity ;
00361         xbuftmp.isvalid = false ;
00362         
00363         xbuf.push_back(xbuftmp) ;  // store it to xbuf
00364       
00365 #ifdef G4TWISTDEBUG
00366         G4cout << "solution " << i << " = " << phi << " , " << u  << G4endl ;
00367 #endif
00368 
00369       }  // end if real solution
00370     }  // end loop i
00371     
00372   }    // end general case
00373 
00374 
00375   nxx = xbuf.size() ;  // save the number of  solutions
00376 
00377   G4ThreeVector xxonsurface  ;       // point on surface
00378   G4ThreeVector surfacenormal  ;     // normal vector  
00379   G4double deltaX  ;                 // distance between intersection point and point on surface
00380   G4double theta  ;                  // angle between track and surfacenormal
00381   G4double factor ;                  // a scaling factor
00382   G4int maxint = 30 ;                // number of iterations
00383 
00384 
00385   for ( size_t k = 0 ; k<xbuf.size() ; k++ ) {
00386 
00387 #ifdef G4TWISTDEBUG
00388     G4cout << "Solution " << k << " : " 
00389            << "reconstructed phiR = " << xbuf[k].phi
00390            << ", uR = " << xbuf[k].u << G4endl ; 
00391 #endif
00392     
00393     phi = xbuf[k].phi ;  // get the stored values for phi and u
00394     u = xbuf[k].u ;
00395 
00396     IsConverged = false ;   // no convergence at the beginning
00397     
00398     for ( G4int i = 1 ; i<maxint ; i++ ) {
00399       
00400       xxonsurface = SurfacePoint(phi,u) ;
00401       surfacenormal = NormAng(phi,u) ;
00402       tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 
00403       deltaX = ( tmpxx - xxonsurface ).mag() ; 
00404       theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
00405       if ( theta < 0.001 ) { 
00406         factor = 50 ;
00407         IsParallel = true ;
00408       }
00409       else {
00410         factor = 1 ;
00411       }
00412 
00413 #ifdef G4TWISTDEBUG
00414       G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ;
00415       G4cout << "X = " << tmpxx << G4endl ;
00416 #endif
00417       
00418       GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced
00419       
00420 #ifdef G4TWISTDEBUG
00421       G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
00422 #endif
00423       
00424       if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
00425       
00426     }  // end iterative loop (i)
00427     
00428 
00429     // new code  21.09.05 O.Link
00430     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 
00431 
00432 #ifdef G4TWISTDEBUG
00433     G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
00434     G4cout << "distance = " << tmpdist << G4endl ;
00435     G4cout << "local X = " << tmpxx << G4endl ;
00436 #endif
00437     
00438     tmpisvalid = false ;  // init 
00439 
00440     if ( IsConverged ) {
00441       
00442       if (validate == kValidateWithTol) {
00443         tmpareacode = GetAreaCode(tmpxx);
00444         if (!IsOutside(tmpareacode)) {
00445           if (tmpdist >= 0) tmpisvalid = true;
00446         }
00447       } else if (validate == kValidateWithoutTol) {
00448         tmpareacode = GetAreaCode(tmpxx, false);
00449         if (IsInside(tmpareacode)) {
00450           if (tmpdist >= 0) tmpisvalid = true;
00451         }
00452       } else { // kDontValidate
00453         G4Exception("G4TwistBoxSide::DistanceToSurface()",
00454                     "GeomSolids0001", FatalException,
00455                     "Feature NOT implemented !");
00456       }
00457 
00458     } 
00459     else {
00460       tmpdist = kInfinity;     // no convergence after 10 steps 
00461       tmpisvalid = false ;     // solution is not vaild
00462     }  
00463 
00464 
00465     // store the found values 
00466     xbuf[k].xx = tmpxx ;
00467     xbuf[k].distance = tmpdist ;
00468     xbuf[k].areacode = tmpareacode ;
00469     xbuf[k].isvalid = tmpisvalid ;
00470 
00471 
00472   }  // end loop over physical solutions (variable k)
00473 
00474 
00475   std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
00476 
00477 #ifdef G4TWISTDEBUG
00478   G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
00479   G4cout << G4endl << G4endl ;
00480 #endif
00481 
00482 
00483   // erase identical intersection (within kCarTolerance) 
00484   xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ;
00485 
00486 
00487   // add guesses
00488 
00489   G4int nxxtmp = xbuf.size() ;
00490 
00491   if ( nxxtmp<2 || IsParallel  ) {
00492 
00493     // positive end
00494 #ifdef G4TWISTDEBUG
00495     G4cout << "add guess at +z/2 .. " << G4endl ;
00496 #endif
00497 
00498     phi = fPhiTwist/2 ;
00499     u   = 0 ;
00500 
00501     
00502      
00503     xbuftmp.phi = phi ;
00504     xbuftmp.u = u ;
00505     xbuftmp.areacode = sOutside ;
00506     xbuftmp.distance = kInfinity ;
00507     xbuftmp.isvalid = false ;
00508     
00509     xbuf.push_back(xbuftmp) ;  // store it to xbuf
00510 
00511 
00512 #ifdef G4TWISTDEBUG
00513     G4cout << "add guess at -z/2 .. " << G4endl ;
00514 #endif
00515 
00516     phi = -fPhiTwist/2 ;
00517     u   = 0 ;
00518 
00519     xbuftmp.phi = phi ;
00520     xbuftmp.u = u ;
00521     xbuftmp.areacode = sOutside ;
00522     xbuftmp.distance = kInfinity ;
00523     xbuftmp.isvalid = false ;
00524     
00525     xbuf.push_back(xbuftmp) ;  // store it to xbuf
00526 
00527     for ( size_t k = nxxtmp ; k<xbuf.size() ; k++ ) {
00528 
00529 #ifdef G4TWISTDEBUG
00530       G4cout << "Solution " << k << " : " 
00531              << "reconstructed phiR = " << xbuf[k].phi
00532              << ", uR = " << xbuf[k].u << G4endl ; 
00533 #endif
00534       
00535       phi = xbuf[k].phi ;  // get the stored values for phi and u
00536       u   = xbuf[k].u ;
00537 
00538       IsConverged = false ;   // no convergence at the beginning
00539       
00540       for ( G4int i = 1 ; i<maxint ; i++ ) {
00541         
00542         xxonsurface = SurfacePoint(phi,u) ;
00543         surfacenormal = NormAng(phi,u) ;
00544         tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 
00545         deltaX = ( tmpxx - xxonsurface ).mag() ; 
00546         theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
00547         if ( theta < 0.001 ) { 
00548           factor = 50 ;    
00549         }
00550         else {
00551           factor = 1 ;
00552         }
00553         
00554 #ifdef G4TWISTDEBUG
00555         G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ;
00556         G4cout << "X = " << tmpxx << G4endl ;
00557 #endif
00558 
00559         GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced
00560       
00561 #ifdef G4TWISTDEBUG
00562         G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
00563 #endif
00564       
00565         if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
00566       
00567       }  // end iterative loop (i)
00568     
00569 
00570     // new code  21.09.05 O.Link
00571     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 
00572 
00573 #ifdef G4TWISTDEBUG
00574       G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
00575       G4cout << "distance = " << tmpdist << G4endl ;
00576       G4cout << "local X = " << tmpxx << G4endl ;
00577 #endif
00578 
00579       tmpisvalid = false ;  // init 
00580 
00581       if ( IsConverged ) {
00582 
00583         if (validate == kValidateWithTol) {
00584           tmpareacode = GetAreaCode(tmpxx);
00585           if (!IsOutside(tmpareacode)) {
00586             if (tmpdist >= 0) tmpisvalid = true;
00587           }
00588         } else if (validate == kValidateWithoutTol) {
00589           tmpareacode = GetAreaCode(tmpxx, false);
00590           if (IsInside(tmpareacode)) {
00591             if (tmpdist >= 0) tmpisvalid = true;
00592           }
00593         } else { // kDontValidate
00594           G4Exception("G4TwistedBoxSide::DistanceToSurface()",
00595                       "GeomSolids0001", FatalException,
00596                       "Feature NOT implemented !");
00597         }
00598         
00599       } 
00600       else {
00601         tmpdist = kInfinity;     // no convergence after 10 steps 
00602         tmpisvalid = false ;     // solution is not vaild
00603       }  
00604         
00605         
00606       // store the found values 
00607       xbuf[k].xx = tmpxx ;
00608       xbuf[k].distance = tmpdist ;
00609       xbuf[k].areacode = tmpareacode ;
00610       xbuf[k].isvalid = tmpisvalid ;
00611 
00612 
00613     }  // end loop over physical solutions 
00614 
00615 
00616   }  // end less than 2 solutions
00617 
00618 
00619   // sort again
00620   std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
00621 
00622   // erase identical intersection (within kCarTolerance) 
00623   xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ;
00624 
00625 #ifdef G4TWISTDEBUG
00626   G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
00627   G4cout << G4endl << G4endl ;
00628 #endif
00629 
00630   nxx = xbuf.size() ;   // determine number of solutions again.
00631 
00632   for ( size_t i = 0 ; i<xbuf.size() ; i++ ) {
00633     
00634     distance[i] = xbuf[i].distance;
00635     gxx[i]      = ComputeGlobalPoint(xbuf[i].xx);
00636     areacode[i] = xbuf[i].areacode ;
00637     isvalid[i]  = xbuf[i].isvalid ;
00638     
00639     fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
00640                                      isvalid[i], nxx, validate, &gp, &gv);
00641 
00642 #ifdef G4TWISTDEBUG
00643     G4cout << "element Nr. " << i 
00644            << ", local Intersection = " << xbuf[i].xx 
00645            << ", distance = " << xbuf[i].distance 
00646            << ", u = " << xbuf[i].u 
00647            << ", phi = " << xbuf[i].phi 
00648            << ", isvalid = " << xbuf[i].isvalid 
00649            << G4endl ;
00650 #endif
00651 
00652   }  // end for( i ) loop
00653 
00654     
00655 #ifdef G4TWISTDEBUG
00656   G4cout << "G4TwistBoxSide finished " << G4endl ;
00657   G4cout << nxx << " possible physical solutions found" << G4endl ;
00658   for ( G4int k= 0 ; k< nxx ; k++ ) {
00659     G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
00660     G4cout << "distance = " << distance[k] << G4endl ;
00661     G4cout << "isvalid = " << isvalid[k] << G4endl ;
00662   }
00663 #endif
00664 
00665   return nxx ;
00666     
00667 }
00668 
00669 
00670 //=====================================================================
00671 //* DistanceToSurface -------------------------------------------------
00672 
00673 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector &gp,
00674                                                 G4ThreeVector  gxx[],
00675                                                 G4double       distance[],
00676                                                 G4int          areacode[])
00677 {  
00678   // to do
00679 
00680   static const G4double ctol = 0.5 * kCarTolerance;
00681 
00682   fCurStat.ResetfDone(kDontValidate, &gp);
00683 
00684    if (fCurStat.IsDone()) {
00685       G4int i;
00686       for (i=0; i<fCurStat.GetNXX(); i++) {
00687          gxx[i] = fCurStat.GetXX(i);
00688          distance[i] = fCurStat.GetDistance(i);
00689          areacode[i] = fCurStat.GetAreacode(i);
00690       }
00691       return fCurStat.GetNXX();
00692    } else {
00693       // initialize
00694       G4int i;
00695       for (i=0; i<G4VSURFACENXX; i++) {
00696          distance[i] = kInfinity;
00697          areacode[i] = sOutside;
00698          gxx[i].set(kInfinity, kInfinity, kInfinity);
00699       }
00700    }
00701    
00702    G4ThreeVector p = ComputeLocalPoint(gp);
00703    G4ThreeVector xx;  // intersection point
00704    G4ThreeVector xxonsurface ; // interpolated intersection point 
00705 
00706    // the surfacenormal at that surface point
00707    G4double phiR = 0  ; // 
00708    G4double uR = 0 ;
00709 
00710    G4ThreeVector surfacenormal ; 
00711    G4double deltaX ;
00712    
00713    G4int maxint = 20 ;
00714 
00715    for ( G4int i = 1 ; i<maxint ; i++ ) {
00716 
00717      xxonsurface = SurfacePoint(phiR,uR) ;
00718      surfacenormal = NormAng(phiR,uR) ;
00719      distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX
00720      deltaX = ( xx - xxonsurface ).mag() ; 
00721 
00722 #ifdef G4TWISTDEBUG
00723      G4cout << "i = " << i << ", distance = " << distance[0] << ", " << deltaX << G4endl ;
00724      G4cout << "X = " << xx << G4endl ;
00725 #endif
00726 
00727      // the new point xx is accepted and phi/psi replaced
00728      GetPhiUAtX(xx, phiR, uR) ;
00729      
00730      if ( deltaX <= ctol ) { break ; }
00731 
00732    }
00733 
00734    // check validity of solution ( valid phi,psi ) 
00735 
00736    G4double halfphi = 0.5*fPhiTwist ;
00737    G4double uMax = GetBoundaryMax(phiR) ;
00738 
00739    if (  phiR > halfphi ) phiR =  halfphi ;
00740    if ( phiR < -halfphi ) phiR = -halfphi ;
00741    if ( uR > uMax ) uR = uMax ;
00742    if ( uR < -uMax ) uR = -uMax ;
00743 
00744    xxonsurface = SurfacePoint(phiR,uR) ;
00745    distance[0] = (  p - xx ).mag() ;
00746    if ( distance[0] <= ctol ) { distance[0] = 0 ; } 
00747 
00748    // end of validity 
00749 
00750 #ifdef G4TWISTDEBUG
00751    G4cout << "refined solution "  << phiR << " , " << uR << " , " <<  G4endl ;
00752    G4cout << "distance = " << distance[0] << G4endl ;
00753    G4cout << "X = " << xx << G4endl ;
00754 #endif
00755 
00756    G4bool isvalid = true;
00757    gxx[0]      = ComputeGlobalPoint(xx);
00758    
00759 #ifdef G4TWISTDEBUG
00760    G4cout << "intersection Point found: " << gxx[0] << G4endl ;
00761    G4cout << "distance = " << distance[0] << G4endl ;
00762 #endif
00763 
00764    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00765                             isvalid, 1, kDontValidate, &gp);
00766    return 1;
00767    
00768 
00769 }
00770 
00771 
00772 //=====================================================================
00773 //* GetAreaCode -------------------------------------------------------
00774 
00775 G4int G4TwistBoxSide::GetAreaCode(const G4ThreeVector &xx, 
00776                                           G4bool withTol)
00777 {
00778    // We must use the function in local coordinate system.
00779    // See the description of DistanceToSurface(p,v).
00780    
00781    static const G4double ctol = 0.5 * kCarTolerance;
00782 
00783    G4double phi ;
00784    G4double yprime ;
00785    GetPhiUAtX(xx, phi,yprime ) ;
00786 
00787    G4double fYAxisMax =  GetBoundaryMax(phi) ;   // Boundaries are symmetric
00788    G4double fYAxisMin =  - fYAxisMax ;
00789 
00790 #ifdef G4TWISTDEBUG
00791    G4cout << "GetAreaCode: phi = " << phi << G4endl ;
00792    G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
00793    G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
00794 #endif
00795 
00796    G4int areacode = sInside;
00797    
00798    if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
00799 
00800       G4int zaxis = 1;
00801       
00802       if (withTol) {
00803 
00804         G4bool isoutside   = false;
00805         
00806         // test boundary of yaxis
00807 
00808          if (yprime < fYAxisMin + ctol) {
00809             areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary; 
00810             if (yprime <= fYAxisMin - ctol) isoutside = true;
00811 
00812          } else if (yprime > fYAxisMax - ctol) {
00813             areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
00814             if (yprime >= fYAxisMax + ctol)  isoutside = true;
00815          }
00816 
00817          // test boundary of z-axis
00818 
00819          if (xx.z() < fAxisMin[zaxis] + ctol) {
00820             areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 
00821 
00822             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00823             else                        areacode |= sBoundary;
00824             if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
00825 
00826          } else if (xx.z() > fAxisMax[zaxis] - ctol) {
00827             areacode |= (sAxis1 & (sAxisZ | sAxisMax));
00828 
00829             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00830             else                        areacode |= sBoundary; 
00831             if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
00832          }
00833 
00834          // if isoutside = true, clear inside bit.             
00835          // if not on boundary, add axis information.             
00836          
00837          if (isoutside) {
00838             G4int tmpareacode = areacode & (~sInside);
00839             areacode = tmpareacode;
00840          } else if ((areacode & sBoundary) != sBoundary) {
00841             areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
00842          }           
00843          
00844       } else {
00845 
00846          // boundary of y-axis
00847 
00848          if (yprime < fYAxisMin ) {
00849             areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
00850          } else if (yprime > fYAxisMax) {
00851             areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
00852          }
00853          
00854          // boundary of z-axis
00855 
00856          if (xx.z() < fAxisMin[zaxis]) {
00857             areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00858             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00859             else                        areacode |= sBoundary; 
00860            
00861          } else if (xx.z() > fAxisMax[zaxis]) {
00862             areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
00863             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
00864             else                        areacode |= sBoundary; 
00865          }
00866 
00867          if ((areacode & sBoundary) != sBoundary) {
00868             areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
00869          }           
00870       }
00871       return areacode;
00872    } else {
00873       G4Exception("G4TwistBoxSide::GetAreaCode()",
00874                   "GeomSolids0001", FatalException,
00875                   "Feature NOT implemented !");
00876    }
00877    return areacode;
00878 }
00879 
00880 //=====================================================================
00881 //* SetCorners() ------------------------------------------------------
00882 
00883 void G4TwistBoxSide::SetCorners()
00884 {
00885 
00886   // Set Corner points in local coodinate.   
00887 
00888   if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
00889     
00890     G4double x, y, z;
00891 
00892     // corner of Axis0min and Axis1min
00893 
00894     x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std::cos(fPhiTwist/2.) - fDy1*std::sin(fPhiTwist/2.) ;
00895     y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.) + (-fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
00896     z = -fDz ;
00897 
00898     SetCorner(sC0Min1Min, x, y, z);
00899 
00900     // corner of Axis0max and Axis1min
00901 
00902     x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
00903     y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
00904     z = -fDz ;
00905 
00906     SetCorner(sC0Max1Min, x, y, z);
00907       
00908     // corner of Axis0max and Axis1max
00909 
00910     x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
00911     y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
00912     z = fDz ;
00913 
00914     SetCorner(sC0Max1Max, x, y, z);
00915       
00916     // corner of Axis0min and Axis1max
00917 
00918     x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std::cos(fPhiTwist/2.) + fDy2*std::sin(fPhiTwist/2.) ;
00919     y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.) + (fDx4 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
00920     z = fDz ;
00921 
00922     SetCorner(sC0Min1Max, x, y, z);
00923 
00924   } else {
00925 
00926     G4Exception("G4TwistBoxSide::SetCorners()",
00927                 "GeomSolids0001", FatalException,
00928                 "Method NOT implemented !");
00929   }
00930 }
00931 
00932 //=====================================================================
00933 //* SetBoundaries() ---------------------------------------------------
00934 
00935 void G4TwistBoxSide::SetBoundaries()
00936 {
00937    // Set direction-unit vector of boundary-lines in local coodinate. 
00938    //   
00939 
00940   G4ThreeVector direction;
00941    
00942   if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
00943       
00944     // sAxis0 & sAxisMin
00945     direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
00946     direction = direction.unit();
00947     SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction, 
00948                 GetCorner(sC0Min1Min), sAxisZ) ;
00949       
00950       // sAxis0 & sAxisMax
00951     direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
00952     direction = direction.unit();
00953     SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction, 
00954                 GetCorner(sC0Max1Min), sAxisZ);
00955     
00956     // sAxis1 & sAxisMin
00957     direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
00958     direction = direction.unit();
00959     SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 
00960                 GetCorner(sC0Min1Min), sAxisY);
00961     
00962     // sAxis1 & sAxisMax
00963     direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
00964     direction = direction.unit();
00965     SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 
00966                 GetCorner(sC0Min1Max), sAxisY);
00967     
00968   } else {
00969     
00970   G4Exception("G4TwistBoxSide::SetCorners()",
00971               "GeomSolids0001", FatalException,
00972               "Feature NOT implemented !");
00973   }
00974 }
00975 
00976 
00977 void G4TwistBoxSide::GetPhiUAtX( G4ThreeVector p, G4double &phi, G4double &u) 
00978 {
00979   // find closest point XX on surface for a given point p
00980   // X0 is a point on the surface,  d is the direction ( both for a fixed z = pz)
00981   
00982   // phi is given by the z coordinate of p
00983 
00984   phi = p.z()/(2*fDz)*fPhiTwist ;
00985 
00986   u =  -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4minus2*phi) + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi - fPhiTwist*(fTAlph*p.x() + p.y()))* std::cos(phi) + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*phi + fPhiTwist*(p.x() - fTAlph*p.y()))*std::sin(phi))/(2.*(fPhiTwist + fPhiTwist*fTAlph*fTAlph)) ;
00987 
00988 
00989 }
00990 
00991 
00992 G4ThreeVector G4TwistBoxSide::ProjectPoint(const G4ThreeVector &p, 
00993                                                     G4bool isglobal) 
00994 {
00995   // Get Rho at p.z() on Hyperbolic Surface.
00996   G4ThreeVector tmpp;
00997   if (isglobal) {
00998      tmpp = fRot.inverse()*p - fTrans;
00999   } else {
01000      tmpp = p;
01001   }
01002 
01003   G4double phi ;
01004   G4double u ;
01005 
01006   GetPhiUAtX( tmpp, phi, u ) ;  // calculate (phi, u) for a point p close the surface
01007   
01008   G4ThreeVector xx = SurfacePoint(phi,u) ;  // transform back to cartesian coordinates
01009 
01010   if (isglobal) {
01011      return (fRot * xx + fTrans);
01012   } else {
01013      return xx;
01014   }
01015 }
01016 
01017 void G4TwistBoxSide::GetFacets( G4int k, G4int n,  G4double xyz[][3], G4int faces[][4], G4int iside ) 
01018 {
01019 
01020   G4double phi ;
01021   G4double b ;   
01022 
01023   G4double z, u ;     // the two parameters for the surface equation
01024   G4ThreeVector p ;  // a point on the surface, given by (z,u)
01025 
01026   G4int nnode ;
01027   G4int nface ;
01028 
01029   // calculate the (n-1)*(k-1) vertices
01030 
01031   G4int i,j ;
01032 
01033   for ( i = 0 ; i<n ; i++ ) {
01034 
01035     z = -fDz+i*(2.*fDz)/(n-1) ;
01036     phi = z*fPhiTwist/(2*fDz) ;
01037     b = GetValueB(phi) ;
01038 
01039     for ( j = 0 ; j<k ; j++ ) {
01040 
01041       nnode = GetNode(i,j,k,n,iside) ;
01042       u = -b/2 +j*b/(k-1) ;
01043       p = SurfacePoint(phi,u,true) ;  // surface point in global coordinate system
01044 
01045       xyz[nnode][0] = p.x() ;
01046       xyz[nnode][1] = p.y() ;
01047       xyz[nnode][2] = p.z() ;
01048 
01049       if ( i<n-1 && j<k-1 ) {   // conterclock wise filling
01050         
01051         nface = GetFace(i,j,k,n,iside) ;
01052         faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i  ,j  ,k,n,iside)+1) ;  // fortran numbering
01053         faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i  ,j+1,k,n,iside)+1) ;
01054         faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ;
01055         faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j  ,k,n,iside)+1) ;
01056 
01057       }
01058     }
01059   }
01060 }

Generated on Mon May 27 17:50:04 2013 for Geant4 by  doxygen 1.4.7