Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TwistBoxSide.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/solids/specific/src/G4TwistBoxSide.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TwistBoxSide.cc (Version 11.2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4TwistBoxSide implementation                   26 // G4TwistBoxSide implementation
 27 //                                                 27 //
 28 // Author: 18/03/2005 - O.Link (Oliver.Link@ce     28 // Author: 18/03/2005 - O.Link (Oliver.Link@cern.ch)
 29 // -------------------------------------------     29 // --------------------------------------------------------------------
 30                                                    30 
 31 #include <cmath>                                   31 #include <cmath>
 32                                                    32 
 33 #include "G4TwistBoxSide.hh"                       33 #include "G4TwistBoxSide.hh"
 34 #include "G4PhysicalConstants.hh"                  34 #include "G4PhysicalConstants.hh"
 35 #include "G4JTPolynomialSolver.hh"                 35 #include "G4JTPolynomialSolver.hh"
 36                                                    36 
 37 //============================================     37 //=====================================================================
 38 //* constructors -----------------------------     38 //* constructors ------------------------------------------------------
 39                                                    39 
 40 G4TwistBoxSide::G4TwistBoxSide(const G4String&     40 G4TwistBoxSide::G4TwistBoxSide(const G4String&     name,
 41                            G4double      PhiTw     41                            G4double      PhiTwist,    // twist angle
 42                            G4double      pDz,      42                            G4double      pDz,         // half z lenght
 43                            G4double      pThet     43                            G4double      pTheta,      // direction between end planes
 44                            G4double      pPhi,     44                            G4double      pPhi,        // defined by polar and azimutal angles.
 45                            G4double      pDy1,     45                            G4double      pDy1,        // half y length at -pDz
 46                            G4double      pDx1,     46                            G4double      pDx1,        // half x length at -pDz,-pDy
 47                            G4double      pDx2,     47                            G4double      pDx2,        // half x length at -pDz,+pDy
 48                            G4double      pDy2,     48                            G4double      pDy2,        // half y length at +pDz
 49                            G4double      pDx3,     49                            G4double      pDx3,        // half x length at +pDz,-pDy
 50                            G4double      pDx4,     50                            G4double      pDx4,        // half x length at +pDz,+pDy
 51                            G4double      pAlph     51                            G4double      pAlph,       // tilt angle at +pDz
 52                            G4double      Angle     52                            G4double      AngleSide    // parity
 53                                                    53                                                ) : G4VTwistSurface(name)
 54 {                                                  54 {  
 55                                                    55   
 56                                                    56                  
 57   fAxis[0]    = kYAxis; // in local coordinate     57   fAxis[0]    = kYAxis; // in local coordinate system
 58   fAxis[1]    = kZAxis;                            58   fAxis[1]    = kZAxis;
 59   fAxisMin[0] = -kInfinity ;  // Y Axis bounda     59   fAxisMin[0] = -kInfinity ;  // Y Axis boundary
 60   fAxisMax[0] = kInfinity ;   //   depends on      60   fAxisMax[0] = kInfinity ;   //   depends on z !!
 61   fAxisMin[1] = -pDz ;      // Z Axis boundary     61   fAxisMin[1] = -pDz ;      // Z Axis boundary
 62   fAxisMax[1] = pDz ;                              62   fAxisMax[1] = pDz ;
 63                                                    63   
 64   fDx1  = pDx1 ;                                   64   fDx1  = pDx1 ;
 65   fDx2  = pDx2 ;  // box                           65   fDx2  = pDx2 ;  // box
 66   fDx3  = pDx3 ;                                   66   fDx3  = pDx3 ;
 67   fDx4  = pDx4 ;  // box                           67   fDx4  = pDx4 ;  // box
 68                                                    68 
 69   // this is an overhead. But the parameter na     69   // this is an overhead. But the parameter naming scheme fits to the other surfaces.
 70                                                    70  
 71   if ( fDx1 != fDx2 || fDx3 != fDx4 )              71   if ( fDx1 != fDx2 || fDx3 != fDx4 )
 72   {                                                72   {
 73     std::ostringstream message;                    73     std::ostringstream message;
 74     message << "TwistedTrapBoxSide is not used     74     message << "TwistedTrapBoxSide is not used as a the side of a box: "
 75             << GetName() << G4endl                 75             << GetName() << G4endl
 76             << "        Not a box !";              76             << "        Not a box !";
 77     G4Exception("G4TwistBoxSide::G4TwistBoxSid     77     G4Exception("G4TwistBoxSide::G4TwistBoxSide()", "GeomSolids0002",
 78                 FatalException, message);          78                 FatalException, message);
 79   }                                                79   }
 80                                                    80  
 81   fDy1   = pDy1 ;                                  81   fDy1   = pDy1 ;
 82   fDy2   = pDy2 ;                                  82   fDy2   = pDy2 ;
 83                                                    83 
 84   fDz   = pDz ;                                    84   fDz   = pDz ;
 85                                                    85 
 86   fAlph = pAlph  ;                                 86   fAlph = pAlph  ;
 87   fTAlph = std::tan(fAlph) ;                       87   fTAlph = std::tan(fAlph) ;
 88                                                    88 
 89   fTheta = pTheta ;                                89   fTheta = pTheta ;
 90   fPhi   = pPhi ;                                  90   fPhi   = pPhi ;
 91                                                    91 
 92   // precalculate frequently used parameters       92   // precalculate frequently used parameters
 93                                                    93 
 94   fDx4plus2  = fDx4 + fDx2 ;                       94   fDx4plus2  = fDx4 + fDx2 ;
 95   fDx4minus2 = fDx4 - fDx2 ;                       95   fDx4minus2 = fDx4 - fDx2 ;
 96   fDx3plus1  = fDx3 + fDx1 ;                       96   fDx3plus1  = fDx3 + fDx1 ; 
 97   fDx3minus1 = fDx3 - fDx1 ;                       97   fDx3minus1 = fDx3 - fDx1 ;
 98   fDy2plus1  = fDy2 + fDy1 ;                       98   fDy2plus1  = fDy2 + fDy1 ;
 99   fDy2minus1 = fDy2 - fDy1 ;                       99   fDy2minus1 = fDy2 - fDy1 ;
100                                                   100 
101   fa1md1 = 2*fDx2 - 2*fDx1 ;                      101   fa1md1 = 2*fDx2 - 2*fDx1 ; 
102   fa2md2 = 2*fDx4 - 2*fDx3 ;                      102   fa2md2 = 2*fDx4 - 2*fDx3 ;
103                                                   103 
104                                                   104 
105   fPhiTwist = PhiTwist ;     // dphi              105   fPhiTwist = PhiTwist ;     // dphi
106   fAngleSide = AngleSide ;  // 0,90,180,270 de    106   fAngleSide = AngleSide ;  // 0,90,180,270 deg
107                                                   107 
108   fdeltaX = 2*fDz * std::tan(fTheta) * std::co    108   fdeltaX = 2*fDz * std::tan(fTheta) * std::cos(fPhi); // dx in surface equation
109   fdeltaY = 2*fDz * std::tan(fTheta) * std::si    109   fdeltaY = 2*fDz * std::tan(fTheta) * std::sin(fPhi); // dy in surface equation
110                                                   110   
111   fRot.rotateZ( AngleSide ) ;                     111   fRot.rotateZ( AngleSide ) ; 
112                                                   112   
113   fTrans.set(0, 0, 0);  // No Translation         113   fTrans.set(0, 0, 0);  // No Translation
114   fIsValidNorm = false;                           114   fIsValidNorm = false;
115                                                   115   
116   SetCorners();                                   116   SetCorners();
117   SetBoundaries();                                117   SetBoundaries();
118 }                                                 118 }
119                                                   119 
120                                                   120 
121 //============================================    121 //=====================================================================
122 //* Fake default constructor -----------------    122 //* Fake default constructor ------------------------------------------
123                                                   123 
124 G4TwistBoxSide::G4TwistBoxSide( __void__& a )     124 G4TwistBoxSide::G4TwistBoxSide( __void__& a )
125   : G4VTwistSurface(a)                            125   : G4VTwistSurface(a)
126 {                                                 126 {
127 }                                                 127 }
128                                                   128 
129                                                   129 
130 //============================================    130 //=====================================================================
131 //* destructor -------------------------------    131 //* destructor --------------------------------------------------------
132                                                   132 
133 G4TwistBoxSide::~G4TwistBoxSide() = default;      133 G4TwistBoxSide::~G4TwistBoxSide() = default;
134                                                   134 
135 //============================================    135 //=====================================================================
136 //* GetNormal --------------------------------    136 //* GetNormal ---------------------------------------------------------
137                                                   137 
138 G4ThreeVector G4TwistBoxSide::GetNormal(const     138 G4ThreeVector G4TwistBoxSide::GetNormal(const G4ThreeVector& tmpxx, 
139                                                   139                                               G4bool isGlobal) 
140 {                                                 140 {
141    // GetNormal returns a normal vector at a s    141    // GetNormal returns a normal vector at a surface (or very close
142    // to surface) point at tmpxx.                 142    // to surface) point at tmpxx.
143    // If isGlobal=true, it returns the normal     143    // If isGlobal=true, it returns the normal in global coordinate.
144    //                                             144    //
145                                                   145 
146    G4ThreeVector xx;                              146    G4ThreeVector xx;
147    if (isGlobal)                                  147    if (isGlobal)
148    {                                              148    {
149       xx = ComputeLocalPoint(tmpxx);              149       xx = ComputeLocalPoint(tmpxx);
150       if ((xx - fCurrentNormal.p).mag() < 0.5     150       if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
151       {                                           151       {
152          return ComputeGlobalDirection(fCurren    152          return ComputeGlobalDirection(fCurrentNormal.normal);
153       }                                           153       }
154    }                                              154    }
155    else                                           155    else
156    {                                              156    {
157       xx = tmpxx;                                 157       xx = tmpxx;
158       if (xx == fCurrentNormal.p)                 158       if (xx == fCurrentNormal.p)
159       {                                           159       {
160          return fCurrentNormal.normal;            160          return fCurrentNormal.normal;
161       }                                           161       }
162    }                                              162    }
163                                                   163 
164    G4double phi ;                                 164    G4double phi ;
165    G4double u ;                                   165    G4double u ;
166                                                   166 
167    GetPhiUAtX(xx,phi,u) ;   // phi,u for point    167    GetPhiUAtX(xx,phi,u) ;   // phi,u for point xx close to surface 
168                                                   168 
169    G4ThreeVector normal = NormAng(phi,u) ;  //    169    G4ThreeVector normal = NormAng(phi,u) ;  // the normal vector at phi,u
170                                                   170 
171 #ifdef G4TWISTDEBUG                               171 #ifdef G4TWISTDEBUG
172    G4cout  << "normal vector = " << normal <<     172    G4cout  << "normal vector = " << normal << G4endl ;
173    G4cout << "phi = " << phi << " , u = " << u    173    G4cout << "phi = " << phi << " , u = " << u << G4endl ;
174 #endif                                            174 #endif
175                                                   175 
176    //    normal = normal/normal.mag() ;           176    //    normal = normal/normal.mag() ;
177                                                   177 
178    if (isGlobal)                                  178    if (isGlobal)
179    {                                              179    {
180       fCurrentNormal.normal = ComputeGlobalDir    180       fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
181    }                                              181    }
182    else                                           182    else
183    {                                              183    {
184       fCurrentNormal.normal = normal.unit();      184       fCurrentNormal.normal = normal.unit();
185    }                                              185    }
186    return fCurrentNormal.normal;                  186    return fCurrentNormal.normal;
187 }                                                 187 }
188                                                   188 
189 //============================================    189 //=====================================================================
190 //* DistanceToSurface ------------------------    190 //* DistanceToSurface -------------------------------------------------
191                                                   191 
192 G4int G4TwistBoxSide::DistanceToSurface(const     192 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector& gp,
193                                         const     193                                         const G4ThreeVector& gv,
194                                                   194                                               G4ThreeVector  gxx[],
195                                                   195                                               G4double       distance[],
196                                                   196                                               G4int          areacode[],
197                                                   197                                               G4bool         isvalid[],
198                                                   198                                               EValidate      validate)
199 {                                                 199 {
200                                                   200 
201   static const G4double pihalf = pi/2 ;           201   static const G4double pihalf = pi/2 ;
202   const G4double ctol = 0.5 * kCarTolerance;      202   const G4double ctol = 0.5 * kCarTolerance;
203                                                   203 
204   G4bool IsParallel = false ;                     204   G4bool IsParallel = false ;
205   G4bool IsConverged = false ;                    205   G4bool IsConverged = false ;
206                                                   206 
207   G4int nxx = 0 ;  // number of physical solut    207   G4int nxx = 0 ;  // number of physical solutions
208                                                   208 
209   fCurStatWithV.ResetfDone(validate, &gp, &gv)    209   fCurStatWithV.ResetfDone(validate, &gp, &gv);
210                                                   210 
211   if (fCurStatWithV.IsDone())                     211   if (fCurStatWithV.IsDone())
212   {                                               212   {
213     for (G4int i=0; i<fCurStatWithV.GetNXX();     213     for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
214     {                                             214     {
215       gxx[i] = fCurStatWithV.GetXX(i);            215       gxx[i] = fCurStatWithV.GetXX(i);
216       distance[i] = fCurStatWithV.GetDistance(    216       distance[i] = fCurStatWithV.GetDistance(i);
217       areacode[i] = fCurStatWithV.GetAreacode(    217       areacode[i] = fCurStatWithV.GetAreacode(i);
218       isvalid[i]  = fCurStatWithV.IsValid(i);     218       isvalid[i]  = fCurStatWithV.IsValid(i);
219     }                                             219     }
220     return fCurStatWithV.GetNXX();                220     return fCurStatWithV.GetNXX();
221   }                                               221   }
222   else  // initialize                             222   else  // initialize
223   {                                               223   {
224     for (G4int i=0; i<G4VSURFACENXX ; ++i)        224     for (G4int i=0; i<G4VSURFACENXX ; ++i)
225     {                                             225     {
226       distance[i] = kInfinity;                    226       distance[i] = kInfinity;
227       areacode[i] = sOutside;                     227       areacode[i] = sOutside;
228       isvalid[i]  = false;                        228       isvalid[i]  = false;
229       gxx[i].set(kInfinity, kInfinity, kInfini    229       gxx[i].set(kInfinity, kInfinity, kInfinity);
230     }                                             230     }
231   }                                               231   }
232                                                   232 
233   G4ThreeVector p = ComputeLocalPoint(gp);        233   G4ThreeVector p = ComputeLocalPoint(gp);
234   G4ThreeVector v = ComputeLocalDirection(gv);    234   G4ThreeVector v = ComputeLocalDirection(gv);
235                                                   235   
236 #ifdef G4TWISTDEBUG                               236 #ifdef G4TWISTDEBUG
237   G4cout << "Local point p = " << p << G4endl     237   G4cout << "Local point p = " << p << G4endl ;
238   G4cout << "Local direction v = " << v << G4e    238   G4cout << "Local direction v = " << v << G4endl ; 
239 #endif                                            239 #endif
240                                                   240 
241   G4double phi=0.,u=0.,q=0.;  // parameters       241   G4double phi=0.,u=0.,q=0.;  // parameters
242                                                   242 
243   // temporary variables                          243   // temporary variables
244                                                   244 
245   G4double      tmpdist = kInfinity ;             245   G4double      tmpdist = kInfinity ;
246   G4ThreeVector tmpxx;                            246   G4ThreeVector tmpxx;
247   G4int         tmpareacode = sOutside ;          247   G4int         tmpareacode = sOutside ;
248   G4bool        tmpisvalid  = false ;             248   G4bool        tmpisvalid  = false ;
249                                                   249 
250   std::vector<Intersection> xbuf ;                250   std::vector<Intersection> xbuf ;
251   Intersection xbuftmp ;                          251   Intersection xbuftmp ;
252                                                   252   
253   // prepare some variables for the intersecti    253   // prepare some variables for the intersection finder
254                                                   254 
255   G4double L = 2*fDz ;                            255   G4double L = 2*fDz ;
256                                                   256 
257   G4double phixz = fPhiTwist * ( p.x() * v.z()    257   G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
258   G4double phiyz = fPhiTwist * ( p.y() * v.z()    258   G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
259                                                   259 
260   // special case vz = 0                          260   // special case vz = 0
261                                                   261 
262   if ( v.z() == 0. )                              262   if ( v.z() == 0. )
263   {                                               263   {
264     if ( (std::fabs(p.z()) <= L) && (fPhiTwist    264     if ( (std::fabs(p.z()) <= L) && (fPhiTwist != 0.0) )  // intersection possible in z
265     {                                             265     {
266       phi = p.z() * fPhiTwist / L ;  // phi is    266       phi = p.z() * fPhiTwist / L ;  // phi is determined by the z-position 
267                                                   267 
268       q = (2.* fPhiTwist*((v.x() - fTAlph*v.y(    268       q = (2.* fPhiTwist*((v.x() - fTAlph*v.y())*std::cos(phi)
269                              + (fTAlph*v.x() +    269                              + (fTAlph*v.x() + v.y())*std::sin(phi)));
270                                                   270       
271       if (q != 0.0)                               271       if (q != 0.0)
272       {                                           272       {
273         u = (2*(-(fdeltaY*phi*v.x()) + fPhiTwi    273         u = (2*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x()
274                 + fdeltaX*phi*v.y() - fPhiTwis    274                 + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y())
275              + (fDx4plus2*fPhiTwist + 2*fDx4mi    275              + (fDx4plus2*fPhiTwist + 2*fDx4minus2*phi)
276              * (v.y()*std::cos(phi) - v.x()*st    276              * (v.y()*std::cos(phi) - v.x()*std::sin(phi))) / q;
277       }                                           277       }
278       xbuftmp.phi = phi ;                         278       xbuftmp.phi = phi ;
279       xbuftmp.u = u ;                             279       xbuftmp.u = u ;
280       xbuftmp.areacode = sOutside ;               280       xbuftmp.areacode = sOutside ;
281       xbuftmp.distance = kInfinity ;              281       xbuftmp.distance = kInfinity ;
282       xbuftmp.isvalid = false ;                   282       xbuftmp.isvalid = false ;
283                                                   283 
284       xbuf.push_back(xbuftmp) ;  // store it t    284       xbuf.push_back(xbuftmp) ;  // store it to xbuf
285     }                                             285     }
286     else                         // no interse    286     else                         // no intersection possible
287     {                                             287     {
288       distance[0] = kInfinity;                    288       distance[0] = kInfinity;
289       gxx[0].set(kInfinity,kInfinity,kInfinity    289       gxx[0].set(kInfinity,kInfinity,kInfinity);
290       isvalid[0] = false ;                        290       isvalid[0] = false ;
291       areacode[0] = sOutside ;                    291       areacode[0] = sOutside ;
292       fCurStatWithV.SetCurrentStatus(0, gxx[0]    292       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
293                                      areacode[    293                                      areacode[0], isvalid[0],
294                                      0, valida    294                                      0, validate, &gp, &gv);
295                                                   295       
296       return 0;                                   296       return 0;
297     }  // end std::fabs(p.z() <= L                297     }  // end std::fabs(p.z() <= L
298   } // end v.z() == 0                             298   } // end v.z() == 0
299   else  // general solution for non-zero vz       299   else  // general solution for non-zero vz
300   {                                               300   {
301     G4double c[8],srd[7],si[7] ;                  301     G4double c[8],srd[7],si[7] ;  
302                                                   302 
303     c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz +    303     c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz + fDx4plus2*fPhiTwist*v.z()) ;
304     c[6] = 28800*(phiyz + 2*fDz*v.x() - (fdelt    304     c[6] = 28800*(phiyz + 2*fDz*v.x() - (fdeltaX + fDx4minus2)*v.z() + fTAlph*(phixz - 2*fDz*v.y() + fdeltaY*v.z())) ;
305     c[5] = -1200*(10*phixz - 48*fDz*v.y() + 24    305     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())) ;
306     c[4] = -2400*(phiyz + 10*fDz*v.x() - 5*fde    306     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())) ;
307     c[3] = 24*(2*phixz - 200*fDz*v.y() + 100*f    307     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())) ;
308     c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6*    308     c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y() + 3*(fdeltaX + fDx4minus2 - fdeltaY*fTAlph)*v.z()) ;
309     c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fD    309     c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x() - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()) ;
310     c[0] = 36*(2* fDz*(v.x() - fTAlph*v.y()) -    310     c[0] = 36*(2* fDz*(v.x() - fTAlph*v.y()) - fdeltaX*v.z() + fdeltaY*fTAlph*v.z()) ;
311                                                   311 
312 #ifdef G4TWISTDEBUG                               312 #ifdef G4TWISTDEBUG
313     G4cout << "coef = " << c[0] << " "            313     G4cout << "coef = " << c[0] << " " 
314            <<  c[1] << " "                        314            <<  c[1] << " "  
315            <<  c[2] << " "                        315            <<  c[2] << " "  
316            <<  c[3] << " "                        316            <<  c[3] << " "  
317            <<  c[4] << " "                        317            <<  c[4] << " "  
318            <<  c[5] << " "                        318            <<  c[5] << " "  
319            <<  c[6] << " "                        319            <<  c[6] << " "  
320            <<  c[7] << G4endl ;                   320            <<  c[7] << G4endl ;
321 #endif                                            321 #endif    
322                                                   322 
323     G4JTPolynomialSolver trapEq ;                 323     G4JTPolynomialSolver trapEq ;
324     G4int num = trapEq.FindRoots(c,7,srd,si);     324     G4int num = trapEq.FindRoots(c,7,srd,si);
325                                                   325 
326     for (G4int i = 0 ; i<num ; ++i )    // loo    326     for (G4int i = 0 ; i<num ; ++i )    // loop over all mathematical solutions
327     {                                             327     {
328       if ( (si[i]==0.0) && (fPhiTwist != 0.0)     328       if ( (si[i]==0.0) && (fPhiTwist != 0.0) )    // only real solutions
329       {                                           329       {
330 #ifdef G4TWISTDEBUG                               330 #ifdef G4TWISTDEBUG
331         G4cout << "Solution " << i << " : " <<    331         G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
332 #endif                                            332 #endif
333         phi = std::fmod(srd[i], pihalf) ;         333         phi = std::fmod(srd[i], pihalf) ;
334                                                   334 
335         u   = (2*phiyz + 4*fDz*phi*v.y() - 2*f    335         u   = (2*phiyz + 4*fDz*phi*v.y() - 2*fdeltaY*phi*v.z()
336              - fDx4plus2*fPhiTwist*v.z()*std::    336              - fDx4plus2*fPhiTwist*v.z()*std::sin(phi)
337              - 2*fDx4minus2*phi*v.z()*std::sin    337              - 2*fDx4minus2*phi*v.z()*std::sin(phi))
338              /(2*fPhiTwist*v.z()*std::cos(phi)    338              /(2*fPhiTwist*v.z()*std::cos(phi)
339                + 2*fPhiTwist*fTAlph*v.z()*std:    339                + 2*fPhiTwist*fTAlph*v.z()*std::sin(phi)) ;
340                                                   340 
341         xbuftmp.phi = phi ;                       341         xbuftmp.phi = phi ;
342         xbuftmp.u = u ;                           342         xbuftmp.u = u ;
343         xbuftmp.areacode = sOutside ;             343         xbuftmp.areacode = sOutside ;
344         xbuftmp.distance = kInfinity ;            344         xbuftmp.distance = kInfinity ;
345         xbuftmp.isvalid = false ;                 345         xbuftmp.isvalid = false ;
346                                                   346         
347         xbuf.push_back(xbuftmp) ;  // store it    347         xbuf.push_back(xbuftmp) ;  // store it to xbuf
348                                                   348       
349 #ifdef G4TWISTDEBUG                               349 #ifdef G4TWISTDEBUG
350         G4cout << "solution " << i << " = " <<    350         G4cout << "solution " << i << " = " << phi << " , " << u  << G4endl ;
351 #endif                                            351 #endif
352       }  // end if real solution                  352       }  // end if real solution
353     }  // end loop i                              353     }  // end loop i
354   }  // end general case                          354   }  // end general case
355                                                   355 
356                                                   356 
357   nxx = (G4int)xbuf.size() ;  // save the numb    357   nxx = (G4int)xbuf.size() ;  // save the number of  solutions
358                                                   358 
359   G4ThreeVector xxonsurface  ;       // point     359   G4ThreeVector xxonsurface  ;       // point on surface
360   G4ThreeVector surfacenormal  ;     // normal    360   G4ThreeVector surfacenormal  ;     // normal vector  
361   G4double deltaX  ;                 // distan    361   G4double deltaX  ;                 // distance between intersection point and
362                                      // point     362                                      // point on surface
363   G4double theta  ;                  // angle     363   G4double theta  ;                  // angle between track and surfacenormal
364   G4double factor ;                  // a scal    364   G4double factor ;                  // a scaling factor
365   G4int maxint = 30 ;                // number    365   G4int maxint = 30 ;                // number of iterations
366                                                   366 
367                                                   367 
368   for (auto & k : xbuf)                           368   for (auto & k : xbuf)
369   {                                               369   {
370 #ifdef G4TWISTDEBUG                               370 #ifdef G4TWISTDEBUG
371     G4cout << "Solution " << k << " : "           371     G4cout << "Solution " << k << " : " 
372            << "reconstructed phiR = " << xbuf[    372            << "reconstructed phiR = " << xbuf[k].phi
373            << ", uR = " << xbuf[k].u << G4endl    373            << ", uR = " << xbuf[k].u << G4endl ; 
374 #endif                                            374 #endif
375                                                   375     
376     phi = k.phi ;  // get the stored values fo    376     phi = k.phi ;  // get the stored values for phi and u
377     u = k.u ;                                     377     u = k.u ;
378                                                   378 
379     IsConverged = false ;   // no convergence     379     IsConverged = false ;   // no convergence at the beginning
380                                                   380     
381     for ( G4int i = 1 ; i<maxint ; ++i )          381     for ( G4int i = 1 ; i<maxint ; ++i )
382     {                                             382     {
383       xxonsurface = SurfacePoint(phi,u) ;         383       xxonsurface = SurfacePoint(phi,u) ;
384       surfacenormal = NormAng(phi,u) ;            384       surfacenormal = NormAng(phi,u) ;
385       tmpdist = DistanceToPlaneWithV(p, v, xxo    385       tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 
386       deltaX = ( tmpxx - xxonsurface ).mag() ;    386       deltaX = ( tmpxx - xxonsurface ).mag() ; 
387       theta = std::fabs(std::acos(v*surfacenor    387       theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
388       if ( theta < 0.001 )                        388       if ( theta < 0.001 )
389       {                                           389       { 
390         factor = 50 ;                             390         factor = 50 ;
391         IsParallel = true ;                       391         IsParallel = true ;
392       }                                           392       }
393       else                                        393       else
394       {                                           394       {
395         factor = 1 ;                              395         factor = 1 ;
396       }                                           396       }
397                                                   397 
398 #ifdef G4TWISTDEBUG                               398 #ifdef G4TWISTDEBUG
399       G4cout << "Step i = " << i << ", distanc    399       G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ;
400       G4cout << "X = " << tmpxx << G4endl ;       400       G4cout << "X = " << tmpxx << G4endl ;
401 #endif                                            401 #endif
402                                                   402       
403       GetPhiUAtX(tmpxx, phi, u) ; // the new p    403       GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced
404                                                   404       
405 #ifdef G4TWISTDEBUG                               405 #ifdef G4TWISTDEBUG
406       G4cout << "approximated phi = " << phi <    406       G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
407 #endif                                            407 #endif
408                                                   408       
409       if ( deltaX <= factor*ctol ) { IsConverg    409       if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
410                                                   410       
411     }  // end iterative loop (i)                  411     }  // end iterative loop (i)
412                                                   412 
413     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0    413     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0. ; 
414                                                   414 
415 #ifdef G4TWISTDEBUG                               415 #ifdef G4TWISTDEBUG
416     G4cout << "refined solution "  << phi << "    416     G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
417     G4cout << "distance = " << tmpdist << G4en    417     G4cout << "distance = " << tmpdist << G4endl ;
418     G4cout << "local X = " << tmpxx << G4endl     418     G4cout << "local X = " << tmpxx << G4endl ;
419 #endif                                            419 #endif
420                                                   420     
421     tmpisvalid = false ;  // init                 421     tmpisvalid = false ;  // init 
422                                                   422 
423     if ( IsConverged )                            423     if ( IsConverged )
424     {                                             424     {
425       if (validate == kValidateWithTol)           425       if (validate == kValidateWithTol)
426       {                                           426       {
427         tmpareacode = GetAreaCode(tmpxx);         427         tmpareacode = GetAreaCode(tmpxx);
428         if (!IsOutside(tmpareacode))              428         if (!IsOutside(tmpareacode))
429         {                                         429         {
430           if (tmpdist >= 0) tmpisvalid = true;    430           if (tmpdist >= 0) tmpisvalid = true;
431         }                                         431         }
432       }                                           432       }
433       else if (validate == kValidateWithoutTol    433       else if (validate == kValidateWithoutTol)
434       {                                           434       {
435         tmpareacode = GetAreaCode(tmpxx, false    435         tmpareacode = GetAreaCode(tmpxx, false);
436         if (IsInside(tmpareacode))                436         if (IsInside(tmpareacode))
437         {                                         437         {
438           if (tmpdist >= 0) tmpisvalid = true;    438           if (tmpdist >= 0) tmpisvalid = true;
439         }                                         439         }
440       }                                           440       }
441       else   // kDontValidate                     441       else   // kDontValidate
442       {                                           442       {
443         G4Exception("G4TwistBoxSide::DistanceT    443         G4Exception("G4TwistBoxSide::DistanceToSurface()",
444                     "GeomSolids0001", FatalExc    444                     "GeomSolids0001", FatalException,
445                     "Feature NOT implemented !    445                     "Feature NOT implemented !");
446       }                                           446       }
447     }                                             447     }
448     else                                          448     else
449     {                                             449     {
450       tmpdist = kInfinity;     // no convergen    450       tmpdist = kInfinity;     // no convergence after 10 steps 
451       tmpisvalid = false ;     // solution is     451       tmpisvalid = false ;     // solution is not vaild
452     }                                             452     }  
453                                                   453 
454     // store the found values                     454     // store the found values 
455     k.xx = tmpxx ;                                455     k.xx = tmpxx ;
456     k.distance = tmpdist ;                        456     k.distance = tmpdist ;
457     k.areacode = tmpareacode ;                    457     k.areacode = tmpareacode ;
458     k.isvalid = tmpisvalid ;                      458     k.isvalid = tmpisvalid ;
459                                                   459 
460   }  // end loop over physical solutions (vari    460   }  // end loop over physical solutions (variable k)
461                                                   461 
462   std::sort(xbuf.begin() , xbuf.end(), Distanc    462   std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
463                                                   463 
464 #ifdef G4TWISTDEBUG                               464 #ifdef G4TWISTDEBUG
465   G4cout << G4endl << "list xbuf after sorting    465   G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
466   G4cout << G4endl << G4endl ;                    466   G4cout << G4endl << G4endl ;
467 #endif                                            467 #endif
468                                                   468 
469   // erase identical intersection (within kCar    469   // erase identical intersection (within kCarTolerance) 
470   xbuf.erase(std::unique(xbuf.begin(),xbuf.end    470   xbuf.erase(std::unique(xbuf.begin(),xbuf.end(),EqualIntersection),xbuf.end());
471                                                   471 
472   // add guesses                                  472   // add guesses
473                                                   473 
474   auto nxxtmp = (G4int)xbuf.size() ;              474   auto nxxtmp = (G4int)xbuf.size() ;
475                                                   475 
476   if ( nxxtmp<2 || IsParallel  )                  476   if ( nxxtmp<2 || IsParallel  )
477   {                                               477   {
478     // positive end                               478     // positive end
479 #ifdef G4TWISTDEBUG                               479 #ifdef G4TWISTDEBUG
480     G4cout << "add guess at +z/2 .. " << G4end    480     G4cout << "add guess at +z/2 .. " << G4endl ;
481 #endif                                            481 #endif
482                                                   482 
483     phi = fPhiTwist/2 ;                           483     phi = fPhiTwist/2 ;
484     u   = 0. ;                                    484     u   = 0. ;
485     xbuftmp.phi = phi ;                           485     xbuftmp.phi = phi ;
486     xbuftmp.u = u ;                               486     xbuftmp.u = u ;
487     xbuftmp.areacode = sOutside ;                 487     xbuftmp.areacode = sOutside ;
488     xbuftmp.distance = kInfinity ;                488     xbuftmp.distance = kInfinity ;
489     xbuftmp.isvalid = false ;                     489     xbuftmp.isvalid = false ;
490                                                   490     
491     xbuf.push_back(xbuftmp) ;  // store it to     491     xbuf.push_back(xbuftmp) ;  // store it to xbuf
492                                                   492 
493 #ifdef G4TWISTDEBUG                               493 #ifdef G4TWISTDEBUG
494     G4cout << "add guess at -z/2 .. " << G4end    494     G4cout << "add guess at -z/2 .. " << G4endl ;
495 #endif                                            495 #endif
496                                                   496 
497     phi = -fPhiTwist/2 ;                          497     phi = -fPhiTwist/2 ;
498     u   = 0. ;                                    498     u   = 0. ;
499                                                   499 
500     xbuftmp.phi = phi ;                           500     xbuftmp.phi = phi ;
501     xbuftmp.u = u ;                               501     xbuftmp.u = u ;
502     xbuftmp.areacode = sOutside ;                 502     xbuftmp.areacode = sOutside ;
503     xbuftmp.distance = kInfinity ;                503     xbuftmp.distance = kInfinity ;
504     xbuftmp.isvalid = false ;                     504     xbuftmp.isvalid = false ;
505                                                   505     
506     xbuf.push_back(xbuftmp) ;  // store it to     506     xbuf.push_back(xbuftmp) ;  // store it to xbuf
507                                                   507 
508     for ( std::size_t k = nxxtmp ; k<xbuf.size    508     for ( std::size_t k = nxxtmp ; k<xbuf.size() ; ++k )
509     {                                             509     {
510 #ifdef G4TWISTDEBUG                               510 #ifdef G4TWISTDEBUG
511       G4cout << "Solution " << k << " : "         511       G4cout << "Solution " << k << " : " 
512              << "reconstructed phiR = " << xbu    512              << "reconstructed phiR = " << xbuf[k].phi
513              << ", uR = " << xbuf[k].u << G4en    513              << ", uR = " << xbuf[k].u << G4endl ; 
514 #endif                                            514 #endif
515                                                   515       
516       phi = xbuf[k].phi ;  // get the stored v    516       phi = xbuf[k].phi ;  // get the stored values for phi and u
517       u   = xbuf[k].u ;                           517       u   = xbuf[k].u ;
518                                                   518 
519       IsConverged = false ;   // no convergenc    519       IsConverged = false ;   // no convergence at the beginning
520                                                   520       
521       for ( G4int i = 1 ; i<maxint ; ++i )        521       for ( G4int i = 1 ; i<maxint ; ++i )
522       {                                           522       {
523         xxonsurface = SurfacePoint(phi,u) ;       523         xxonsurface = SurfacePoint(phi,u) ;
524         surfacenormal = NormAng(phi,u) ;          524         surfacenormal = NormAng(phi,u) ;
525         tmpdist = DistanceToPlaneWithV(p, v, x    525         tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
526         deltaX = ( tmpxx - xxonsurface ).mag()    526         deltaX = ( tmpxx - xxonsurface ).mag() ; 
527         theta = std::fabs(std::acos(v*surfacen    527         theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
528         if ( theta < 0.001 )                      528         if ( theta < 0.001 )
529         {                                         529         { 
530           factor = 50 ;                           530           factor = 50 ;    
531         }                                         531         }
532         else                                      532         else
533         {                                         533         {
534           factor = 1 ;                            534           factor = 1 ;
535         }                                         535         }
536                                                   536         
537 #ifdef G4TWISTDEBUG                               537 #ifdef G4TWISTDEBUG
538         G4cout << "Step i = " << i << ", dista    538         G4cout << "Step i = " << i << ", distance = "
539                << tmpdist << ", " << deltaX <<    539                << tmpdist << ", " << deltaX << G4endl ;
540         G4cout << "X = " << tmpxx << G4endl ;     540         G4cout << "X = " << tmpxx << G4endl ;
541 #endif                                            541 #endif
542                                                   542 
543         GetPhiUAtX(tmpxx, phi, u) ;               543         GetPhiUAtX(tmpxx, phi, u) ;
544           // the new point xx is accepted and     544           // the new point xx is accepted and phi/u replaced
545                                                   545       
546 #ifdef G4TWISTDEBUG                               546 #ifdef G4TWISTDEBUG
547         G4cout << "approximated phi = " << phi    547         G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
548 #endif                                            548 #endif
549                                                   549       
550         if ( deltaX <= factor*ctol ) { IsConve    550         if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
551                                                   551       
552       }  // end iterative loop (i)                552       }  // end iterative loop (i)
553                                                   553 
554     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0    554     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0. ; 
555                                                   555 
556 #ifdef G4TWISTDEBUG                               556 #ifdef G4TWISTDEBUG
557       G4cout << "refined solution "  << phi <<    557       G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
558       G4cout << "distance = " << tmpdist << G4    558       G4cout << "distance = " << tmpdist << G4endl ;
559       G4cout << "local X = " << tmpxx << G4end    559       G4cout << "local X = " << tmpxx << G4endl ;
560 #endif                                            560 #endif
561                                                   561 
562       tmpisvalid = false ;  // init               562       tmpisvalid = false ;  // init 
563                                                   563 
564       if ( IsConverged )                          564       if ( IsConverged )
565       {                                           565       {
566         if (validate == kValidateWithTol)         566         if (validate == kValidateWithTol)
567         {                                         567         {
568           tmpareacode = GetAreaCode(tmpxx);       568           tmpareacode = GetAreaCode(tmpxx);
569           if (!IsOutside(tmpareacode))            569           if (!IsOutside(tmpareacode))
570           {                                       570           {
571             if (tmpdist >= 0) tmpisvalid = tru    571             if (tmpdist >= 0) tmpisvalid = true;
572           }                                       572           }
573         }                                         573         }
574         else if (validate == kValidateWithoutT    574         else if (validate == kValidateWithoutTol)
575         {                                         575         {
576           tmpareacode = GetAreaCode(tmpxx, fal    576           tmpareacode = GetAreaCode(tmpxx, false);
577           if (IsInside(tmpareacode))              577           if (IsInside(tmpareacode))
578           {                                       578           {
579             if (tmpdist >= 0) tmpisvalid = tru    579             if (tmpdist >= 0) tmpisvalid = true;
580           }                                       580           }
581         }                                         581         }
582         else   // kDontValidate                   582         else   // kDontValidate
583         {                                         583         {
584           G4Exception("G4TwistedBoxSide::Dista    584           G4Exception("G4TwistedBoxSide::DistanceToSurface()",
585                       "GeomSolids0001", FatalE    585                       "GeomSolids0001", FatalException,
586                       "Feature NOT implemented    586                       "Feature NOT implemented !");
587         }                                         587         }
588       }                                           588       }
589       else                                        589       else
590       {                                           590       {
591         tmpdist = kInfinity;     // no converg    591         tmpdist = kInfinity;     // no convergence after 10 steps 
592         tmpisvalid = false ;     // solution i    592         tmpisvalid = false ;     // solution is not vaild
593       }                                           593       }  
594                                                   594 
595       // store the found values                   595       // store the found values 
596       xbuf[k].xx = tmpxx ;                        596       xbuf[k].xx = tmpxx ;
597       xbuf[k].distance = tmpdist ;                597       xbuf[k].distance = tmpdist ;
598       xbuf[k].areacode = tmpareacode ;            598       xbuf[k].areacode = tmpareacode ;
599       xbuf[k].isvalid = tmpisvalid ;              599       xbuf[k].isvalid = tmpisvalid ;
600     }  // end loop over physical solutions        600     }  // end loop over physical solutions 
601   }  // end less than 2 solutions                 601   }  // end less than 2 solutions
602                                                   602 
603   // sort again                                   603   // sort again
604   std::sort(xbuf.begin() , xbuf.end(), Distanc    604   std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
605                                                   605 
606   // erase identical intersection (within kCar    606   // erase identical intersection (within kCarTolerance) 
607   xbuf.erase(std::unique(xbuf.begin(),xbuf.end    607   xbuf.erase(std::unique(xbuf.begin(),xbuf.end(),EqualIntersection),xbuf.end());
608                                                   608 
609 #ifdef G4TWISTDEBUG                               609 #ifdef G4TWISTDEBUG
610   G4cout << G4endl << "list xbuf after sorting    610   G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
611   G4cout << G4endl << G4endl ;                    611   G4cout << G4endl << G4endl ;
612 #endif                                            612 #endif
613                                                   613 
614   nxx = (G4int)xbuf.size() ;   // determine nu    614   nxx = (G4int)xbuf.size() ;   // determine number of solutions again.
615                                                   615 
616   for ( G4int i = 0 ; i<(G4int)xbuf.size() ; +    616   for ( G4int i = 0 ; i<(G4int)xbuf.size() ; ++i )
617   {                                               617   {
618     distance[i] = xbuf[i].distance;               618     distance[i] = xbuf[i].distance;
619     gxx[i]      = ComputeGlobalPoint(xbuf[i].x    619     gxx[i]      = ComputeGlobalPoint(xbuf[i].xx);
620     areacode[i] = xbuf[i].areacode ;              620     areacode[i] = xbuf[i].areacode ;
621     isvalid[i]  = xbuf[i].isvalid ;               621     isvalid[i]  = xbuf[i].isvalid ;
622                                                   622     
623     fCurStatWithV.SetCurrentStatus(i, gxx[i],     623     fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
624                                      isvalid[i    624                                      isvalid[i], nxx, validate, &gp, &gv);
625                                                   625 
626 #ifdef G4TWISTDEBUG                               626 #ifdef G4TWISTDEBUG
627     G4cout << "element Nr. " << i                 627     G4cout << "element Nr. " << i 
628            << ", local Intersection = " << xbu    628            << ", local Intersection = " << xbuf[i].xx 
629            << ", distance = " << xbuf[i].dista    629            << ", distance = " << xbuf[i].distance 
630            << ", u = " << xbuf[i].u               630            << ", u = " << xbuf[i].u 
631            << ", phi = " << xbuf[i].phi           631            << ", phi = " << xbuf[i].phi 
632            << ", isvalid = " << xbuf[i].isvali    632            << ", isvalid = " << xbuf[i].isvalid 
633            << G4endl ;                            633            << G4endl ;
634 #endif                                            634 #endif
635                                                   635 
636   }  // end for( i ) loop                         636   }  // end for( i ) loop
637                                                   637 
638 #ifdef G4TWISTDEBUG                               638 #ifdef G4TWISTDEBUG
639   G4cout << "G4TwistBoxSide finished " << G4en    639   G4cout << "G4TwistBoxSide finished " << G4endl ;
640   G4cout << nxx << " possible physical solutio    640   G4cout << nxx << " possible physical solutions found" << G4endl ;
641   for ( G4int k= 0 ; k< nxx ; ++k )               641   for ( G4int k= 0 ; k< nxx ; ++k )
642   {                                               642   {
643     G4cout << "global intersection Point found    643     G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
644     G4cout << "distance = " << distance[k] <<     644     G4cout << "distance = " << distance[k] << G4endl ;
645     G4cout << "isvalid = " << isvalid[k] << G4    645     G4cout << "isvalid = " << isvalid[k] << G4endl ;
646   }                                               646   }
647 #endif                                            647 #endif
648                                                   648 
649   return nxx ;                                    649   return nxx ;
650 }                                                 650 }
651                                                   651 
652 //============================================    652 //=====================================================================
653 //* DistanceToSurface ------------------------    653 //* DistanceToSurface -------------------------------------------------
654                                                   654 
655 G4int G4TwistBoxSide::DistanceToSurface(const     655 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector& gp,
656                                                   656                                               G4ThreeVector  gxx[],
657                                                   657                                               G4double       distance[],
658                                                   658                                               G4int          areacode[])
659 {                                                 659 {  
660   const G4double ctol = 0.5 * kCarTolerance;      660   const G4double ctol = 0.5 * kCarTolerance;
661                                                   661 
662   fCurStat.ResetfDone(kDontValidate, &gp);        662   fCurStat.ResetfDone(kDontValidate, &gp);
663                                                   663 
664    if (fCurStat.IsDone())                         664    if (fCurStat.IsDone())
665    {                                              665    {
666       for (G4int i=0; i<fCurStat.GetNXX(); ++i    666       for (G4int i=0; i<fCurStat.GetNXX(); ++i)
667       {                                           667       {
668          gxx[i] = fCurStat.GetXX(i);              668          gxx[i] = fCurStat.GetXX(i);
669          distance[i] = fCurStat.GetDistance(i)    669          distance[i] = fCurStat.GetDistance(i);
670          areacode[i] = fCurStat.GetAreacode(i)    670          areacode[i] = fCurStat.GetAreacode(i);
671       }                                           671       }
672       return fCurStat.GetNXX();                   672       return fCurStat.GetNXX();
673    }                                              673    }
674    else  // initialize                            674    else  // initialize
675    {                                              675    {
676       for (G4int i=0; i<G4VSURFACENXX; ++i)       676       for (G4int i=0; i<G4VSURFACENXX; ++i)
677       {                                           677       {
678          distance[i] = kInfinity;                 678          distance[i] = kInfinity;
679          areacode[i] = sOutside;                  679          areacode[i] = sOutside;
680          gxx[i].set(kInfinity, kInfinity, kInf    680          gxx[i].set(kInfinity, kInfinity, kInfinity);
681       }                                           681       }
682    }                                              682    }
683                                                   683    
684    G4ThreeVector p = ComputeLocalPoint(gp);       684    G4ThreeVector p = ComputeLocalPoint(gp);
685    G4ThreeVector xx;  // intersection point       685    G4ThreeVector xx;  // intersection point
686    G4ThreeVector xxonsurface ; // interpolated    686    G4ThreeVector xxonsurface ; // interpolated intersection point 
687                                                   687 
688    // the surfacenormal at that surface point     688    // the surfacenormal at that surface point
689    G4double phiR = 0. ;                           689    G4double phiR = 0. ;
690    G4double uR = 0. ;                             690    G4double uR = 0. ;
691                                                   691 
692    G4ThreeVector surfacenormal ;                  692    G4ThreeVector surfacenormal ; 
693    G4double deltaX ;                              693    G4double deltaX ;
694                                                   694    
695    G4int maxint = 20 ;                            695    G4int maxint = 20 ;
696                                                   696 
697    for ( G4int i = 1 ; i<maxint ; ++i )           697    for ( G4int i = 1 ; i<maxint ; ++i )
698    {                                              698    {
699      xxonsurface = SurfacePoint(phiR,uR) ;        699      xxonsurface = SurfacePoint(phiR,uR) ;
700      surfacenormal = NormAng(phiR,uR) ;           700      surfacenormal = NormAng(phiR,uR) ;
701      distance[0] = DistanceToPlane(p, xxonsurf    701      distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX
702      deltaX = ( xx - xxonsurface ).mag() ;        702      deltaX = ( xx - xxonsurface ).mag() ; 
703                                                   703 
704 #ifdef G4TWISTDEBUG                               704 #ifdef G4TWISTDEBUG
705      G4cout << "i = " << i << ", distance = "     705      G4cout << "i = " << i << ", distance = " << distance[0]
706             << ", " << deltaX << G4endl ;         706             << ", " << deltaX << G4endl ;
707      G4cout << "X = " << xx << G4endl ;           707      G4cout << "X = " << xx << G4endl ;
708 #endif                                            708 #endif
709                                                   709 
710      // the new point xx is accepted and phi/p    710      // the new point xx is accepted and phi/psi replaced
711      GetPhiUAtX(xx, phiR, uR) ;                   711      GetPhiUAtX(xx, phiR, uR) ;
712                                                   712      
713      if ( deltaX <= ctol ) { break ; }            713      if ( deltaX <= ctol ) { break ; }
714    }                                              714    }
715                                                   715 
716    // check validity of solution ( valid phi,p    716    // check validity of solution ( valid phi,psi ) 
717                                                   717 
718    G4double halfphi = 0.5*fPhiTwist ;             718    G4double halfphi = 0.5*fPhiTwist ;
719    G4double uMax = GetBoundaryMax(phiR) ;         719    G4double uMax = GetBoundaryMax(phiR) ;
720                                                   720 
721    if (  phiR > halfphi ) phiR =  halfphi ;       721    if (  phiR > halfphi ) phiR =  halfphi ;
722    if ( phiR < -halfphi ) phiR = -halfphi ;       722    if ( phiR < -halfphi ) phiR = -halfphi ;
723    if ( uR > uMax ) uR = uMax ;                   723    if ( uR > uMax ) uR = uMax ;
724    if ( uR < -uMax ) uR = -uMax ;                 724    if ( uR < -uMax ) uR = -uMax ;
725                                                   725 
726    xxonsurface = SurfacePoint(phiR,uR) ;          726    xxonsurface = SurfacePoint(phiR,uR) ;
727    distance[0] = (  p - xx ).mag() ;              727    distance[0] = (  p - xx ).mag() ;
728    if ( distance[0] <= ctol ) { distance[0] =     728    if ( distance[0] <= ctol ) { distance[0] = 0 ; } 
729                                                   729 
730    // end of validity                             730    // end of validity 
731                                                   731 
732 #ifdef G4TWISTDEBUG                               732 #ifdef G4TWISTDEBUG
733    G4cout << "refined solution "  << phiR << "    733    G4cout << "refined solution "  << phiR << " , " << uR << " , " <<  G4endl ;
734    G4cout << "distance = " << distance[0] << G    734    G4cout << "distance = " << distance[0] << G4endl ;
735    G4cout << "X = " << xx << G4endl ;             735    G4cout << "X = " << xx << G4endl ;
736 #endif                                            736 #endif
737                                                   737 
738    G4bool isvalid = true;                         738    G4bool isvalid = true;
739    gxx[0] = ComputeGlobalPoint(xx);               739    gxx[0] = ComputeGlobalPoint(xx);
740                                                   740    
741 #ifdef G4TWISTDEBUG                               741 #ifdef G4TWISTDEBUG
742    G4cout << "intersection Point found: " << g    742    G4cout << "intersection Point found: " << gxx[0] << G4endl ;
743    G4cout << "distance = " << distance[0] << G    743    G4cout << "distance = " << distance[0] << G4endl ;
744 #endif                                            744 #endif
745                                                   745 
746    fCurStat.SetCurrentStatus(0, gxx[0], distan    746    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
747                             isvalid, 1, kDontV    747                             isvalid, 1, kDontValidate, &gp);
748    return 1;                                      748    return 1;
749 }                                                 749 }
750                                                   750 
751 //============================================    751 //=====================================================================
752 //* GetAreaCode ------------------------------    752 //* GetAreaCode -------------------------------------------------------
753                                                   753 
754 G4int G4TwistBoxSide::GetAreaCode(const G4Thre    754 G4int G4TwistBoxSide::GetAreaCode(const G4ThreeVector& xx, 
755                                         G4bool    755                                         G4bool withTol)
756 {                                                 756 {
757    // We must use the function in local coordi    757    // We must use the function in local coordinate system.
758    // See the description of DistanceToSurface    758    // See the description of DistanceToSurface(p,v).
759                                                   759    
760    const G4double ctol = 0.5 * kCarTolerance;     760    const G4double ctol = 0.5 * kCarTolerance;
761                                                   761 
762    G4double phi ;                                 762    G4double phi ;
763    G4double yprime ;                              763    G4double yprime ;
764    GetPhiUAtX(xx, phi,yprime ) ;                  764    GetPhiUAtX(xx, phi,yprime ) ;
765                                                   765 
766    G4double fYAxisMax =  GetBoundaryMax(phi) ;    766    G4double fYAxisMax =  GetBoundaryMax(phi) ;   // Boundaries are symmetric
767    G4double fYAxisMin =  - fYAxisMax ;            767    G4double fYAxisMin =  - fYAxisMax ;
768                                                   768 
769 #ifdef G4TWISTDEBUG                               769 #ifdef G4TWISTDEBUG
770    G4cout << "GetAreaCode: phi = " << phi << G    770    G4cout << "GetAreaCode: phi = " << phi << G4endl ;
771    G4cout << "GetAreaCode: yprime = " << yprim    771    G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
772    G4cout << "Intervall is " << fYAxisMin << "    772    G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
773 #endif                                            773 #endif
774                                                   774 
775    G4int areacode = sInside;                      775    G4int areacode = sInside;
776                                                   776    
777    if (fAxis[0] == kYAxis && fAxis[1] == kZAxi    777    if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
778    {                                              778    {
779       G4int zaxis = 1;                            779       G4int zaxis = 1;
780                                                   780       
781       if (withTol)                                781       if (withTol)
782       {                                           782       {
783         G4bool isoutside = false;                 783         G4bool isoutside = false;
784                                                   784         
785         // test boundary of yaxis                 785         // test boundary of yaxis
786                                                   786 
787          if (yprime < fYAxisMin + ctol)           787          if (yprime < fYAxisMin + ctol)
788          {                                        788          {
789             areacode |= (sAxis0 & (sAxisY | sA    789             areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary; 
790             if (yprime <= fYAxisMin - ctol) is    790             if (yprime <= fYAxisMin - ctol) isoutside = true;
791                                                   791 
792          }                                        792          }
793          else if (yprime > fYAxisMax - ctol)      793          else if (yprime > fYAxisMax - ctol)
794          {                                        794          {
795             areacode |= (sAxis0 & (sAxisY | sA    795             areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
796             if (yprime >= fYAxisMax + ctol)  i    796             if (yprime >= fYAxisMax + ctol)  isoutside = true;
797          }                                        797          }
798                                                   798 
799          // test boundary of z-axis               799          // test boundary of z-axis
800                                                   800 
801          if (xx.z() < fAxisMin[zaxis] + ctol)     801          if (xx.z() < fAxisMin[zaxis] + ctol)
802          {                                        802          {
803             areacode |= (sAxis1 & (sAxisZ | sA    803             areacode |= (sAxis1 & (sAxisZ | sAxisMin)); 
804                                                   804 
805             if   ((areacode & sBoundary) != 0)    805             if   ((areacode & sBoundary) != 0) areacode |= sCorner;  // xx on corner.
806             else                        areaco    806             else                        areacode |= sBoundary;
807             if (xx.z() <= fAxisMin[zaxis] - ct    807             if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
808                                                   808 
809          }                                        809          }
810          else if (xx.z() > fAxisMax[zaxis] - c    810          else if (xx.z() > fAxisMax[zaxis] - ctol)
811          {                                        811          {
812             areacode |= (sAxis1 & (sAxisZ | sA    812             areacode |= (sAxis1 & (sAxisZ | sAxisMax));
813                                                   813 
814             if   ((areacode & sBoundary) != 0)    814             if   ((areacode & sBoundary) != 0) areacode |= sCorner;  // xx on corner.
815             else                        areaco    815             else                        areacode |= sBoundary; 
816             if (xx.z() >= fAxisMax[zaxis] + ct    816             if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
817          }                                        817          }
818                                                   818 
819          // if isoutside = true, clear inside     819          // if isoutside = true, clear inside bit.             
820          // if not on boundary, add axis infor    820          // if not on boundary, add axis information.             
821                                                   821          
822          if (isoutside)                           822          if (isoutside)
823          {                                        823          {
824             G4int tmpareacode = areacode & (~s    824             G4int tmpareacode = areacode & (~sInside);
825             areacode = tmpareacode;               825             areacode = tmpareacode;
826          }                                        826          }
827          else if ((areacode & sBoundary) != sB    827          else if ((areacode & sBoundary) != sBoundary)
828          {                                        828          {
829             areacode |= (sAxis0 & sAxisY) | (s    829             areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
830          }                                        830          }           
831                                                   831          
832       }                                           832       }
833       else                                        833       else
834       {                                           834       {
835          // boundary of y-axis                    835          // boundary of y-axis
836                                                   836 
837          if (yprime < fYAxisMin )                 837          if (yprime < fYAxisMin )
838          {                                        838          {
839             areacode |= (sAxis0 & (sAxisY | sA    839             areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
840          }                                        840          }
841          else if (yprime > fYAxisMax)             841          else if (yprime > fYAxisMax)
842          {                                        842          {
843             areacode |= (sAxis0 & (sAxisY | sA    843             areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
844          }                                        844          }
845                                                   845          
846          // boundary of z-axis                    846          // boundary of z-axis
847                                                   847 
848          if (xx.z() < fAxisMin[zaxis])            848          if (xx.z() < fAxisMin[zaxis])
849          {                                        849          {
850             areacode |= (sAxis1 & (sAxisZ | sA    850             areacode |= (sAxis1 & (sAxisZ | sAxisMin));
851             if   ((areacode & sBoundary) != 0)    851             if   ((areacode & sBoundary) != 0) areacode |= sCorner;  // xx on corner.
852             else                        areaco    852             else                        areacode |= sBoundary; 
853                                                   853            
854          }                                        854          }
855          else if (xx.z() > fAxisMax[zaxis])       855          else if (xx.z() > fAxisMax[zaxis])
856          {                                        856          {
857             areacode |= (sAxis1 & (sAxisZ | sA    857             areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
858             if   ((areacode & sBoundary) != 0)    858             if   ((areacode & sBoundary) != 0) areacode |= sCorner;  // xx on corner.
859             else                        areaco    859             else                        areacode |= sBoundary; 
860          }                                        860          }
861                                                   861 
862          if ((areacode & sBoundary) != sBounda    862          if ((areacode & sBoundary) != sBoundary)
863          {                                        863          {
864             areacode |= (sAxis0 & sAxisY) | (s    864             areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
865          }                                        865          }           
866       }                                           866       }
867       return areacode;                            867       return areacode;
868    }                                              868    }
869    else                                           869    else
870    {                                              870    {
871       G4Exception("G4TwistBoxSide::GetAreaCode    871       G4Exception("G4TwistBoxSide::GetAreaCode()",
872                   "GeomSolids0001", FatalExcep    872                   "GeomSolids0001", FatalException,
873                   "Feature NOT implemented !")    873                   "Feature NOT implemented !");
874    }                                              874    }
875    return areacode;                               875    return areacode;
876 }                                                 876 }
877                                                   877 
878 //============================================    878 //=====================================================================
879 //* SetCorners() -----------------------------    879 //* SetCorners() ------------------------------------------------------
880                                                   880 
881 void G4TwistBoxSide::SetCorners()                 881 void G4TwistBoxSide::SetCorners()
882 {                                                 882 {
883                                                   883 
884   // Set Corner points in local coodinate.        884   // Set Corner points in local coodinate.   
885                                                   885 
886   if (fAxis[0] == kYAxis && fAxis[1] == kZAxis    886   if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
887   {                                               887   {
888     G4double x, y, z;                             888     G4double x, y, z;
889                                                   889 
890     // corner of Axis0min and Axis1min            890     // corner of Axis0min and Axis1min
891                                                   891 
892     x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std    892     x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std::cos(fPhiTwist/2.) - fDy1*std::sin(fPhiTwist/2.) ;
893     y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/    893     y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.) + (-fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
894     z = -fDz ;                                    894     z = -fDz ;
895                                                   895 
896     SetCorner(sC0Min1Min, x, y, z);               896     SetCorner(sC0Min1Min, x, y, z);
897                                                   897 
898     // corner of Axis0max and Axis1min            898     // corner of Axis0max and Axis1min
899                                                   899 
900     x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std    900     x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
901     y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/    901     y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
902     z = -fDz ;                                    902     z = -fDz ;
903                                                   903 
904     SetCorner(sC0Max1Min, x, y, z);               904     SetCorner(sC0Max1Min, x, y, z);
905                                                   905       
906     // corner of Axis0max and Axis1max            906     // corner of Axis0max and Axis1max
907                                                   907 
908     x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std:    908     x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
909     y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2    909     y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
910     z = fDz ;                                     910     z = fDz ;
911                                                   911 
912     SetCorner(sC0Max1Max, x, y, z);               912     SetCorner(sC0Max1Max, x, y, z);
913                                                   913       
914     // corner of Axis0min and Axis1max            914     // corner of Axis0min and Axis1max
915                                                   915 
916     x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std:    916     x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std::cos(fPhiTwist/2.) + fDy2*std::sin(fPhiTwist/2.) ;
917     y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2    917     y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.) + (fDx4 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
918     z = fDz ;                                     918     z = fDz ;
919                                                   919 
920     SetCorner(sC0Min1Max, x, y, z);               920     SetCorner(sC0Min1Max, x, y, z);
921                                                   921 
922   }                                               922   }
923   else                                            923   else
924   {                                               924   {
925     G4Exception("G4TwistBoxSide::SetCorners()"    925     G4Exception("G4TwistBoxSide::SetCorners()",
926                 "GeomSolids0001", FatalExcepti    926                 "GeomSolids0001", FatalException,
927                 "Method NOT implemented !");      927                 "Method NOT implemented !");
928   }                                               928   }
929 }                                                 929 }
930                                                   930 
931 //============================================    931 //=====================================================================
932 //* SetBoundaries() --------------------------    932 //* SetBoundaries() ---------------------------------------------------
933                                                   933 
934 void G4TwistBoxSide::SetBoundaries()              934 void G4TwistBoxSide::SetBoundaries()
935 {                                                 935 {
936    // Set direction-unit vector of boundary-li    936    // Set direction-unit vector of boundary-lines in local coodinates. 
937                                                   937 
938   G4ThreeVector direction;                        938   G4ThreeVector direction;
939                                                   939    
940   if (fAxis[0] == kYAxis && fAxis[1] == kZAxis    940   if (fAxis[0] == kYAxis && fAxis[1] == kZAxis)
941   {                                               941   {
942     // sAxis0 & sAxisMin                          942     // sAxis0 & sAxisMin
943     direction = GetCorner(sC0Min1Max) - GetCor    943     direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
944     direction = direction.unit();                 944     direction = direction.unit();
945     SetBoundary(sAxis0 & (sAxisY | sAxisMin),     945     SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction, 
946                 GetCorner(sC0Min1Min), sAxisZ)    946                 GetCorner(sC0Min1Min), sAxisZ) ;
947                                                   947       
948       // sAxis0 & sAxisMax                        948       // sAxis0 & sAxisMax
949     direction = GetCorner(sC0Max1Max) - GetCor    949     direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
950     direction = direction.unit();                 950     direction = direction.unit();
951     SetBoundary(sAxis0 & (sAxisY | sAxisMax),     951     SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction, 
952                 GetCorner(sC0Max1Min), sAxisZ)    952                 GetCorner(sC0Max1Min), sAxisZ);
953                                                   953     
954     // sAxis1 & sAxisMin                          954     // sAxis1 & sAxisMin
955     direction = GetCorner(sC0Max1Min) - GetCor    955     direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
956     direction = direction.unit();                 956     direction = direction.unit();
957     SetBoundary(sAxis1 & (sAxisZ | sAxisMin),     957     SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 
958                 GetCorner(sC0Min1Min), sAxisY)    958                 GetCorner(sC0Min1Min), sAxisY);
959                                                   959     
960     // sAxis1 & sAxisMax                          960     // sAxis1 & sAxisMax
961     direction = GetCorner(sC0Max1Max) - GetCor    961     direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
962     direction = direction.unit();                 962     direction = direction.unit();
963     SetBoundary(sAxis1 & (sAxisZ | sAxisMax),     963     SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 
964                 GetCorner(sC0Min1Max), sAxisY)    964                 GetCorner(sC0Min1Max), sAxisY);
965                                                   965     
966   }                                               966   }
967   else                                            967   else
968   {                                               968   {
969     G4Exception("G4TwistBoxSide::SetCorners()"    969     G4Exception("G4TwistBoxSide::SetCorners()",
970                 "GeomSolids0001", FatalExcepti    970                 "GeomSolids0001", FatalException,
971                 "Feature NOT implemented !");     971                 "Feature NOT implemented !");
972   }                                               972   }
973 }                                                 973 }
974                                                   974 
975                                                   975 
976 void G4TwistBoxSide::GetPhiUAtX( const G4Three    976 void G4TwistBoxSide::GetPhiUAtX( const G4ThreeVector& p, G4double& phi, G4double& u) 
977 {                                                 977 {
978   // find closest point XX on surface for a gi    978   // find closest point XX on surface for a given point p
979   // X0 is a point on the surface,  d is the d    979   // X0 is a point on the surface,  d is the direction
980   // ( both for a fixed z = pz)                   980   // ( both for a fixed z = pz)
981                                                   981   
982   // phi is given by the z coordinate of p        982   // phi is given by the z coordinate of p
983                                                   983 
984   phi = p.z()/(2*fDz)*fPhiTwist ;                 984   phi = p.z()/(2*fDz)*fPhiTwist ;
985                                                   985 
986   u =  -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4m    986   u =  -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4minus2*phi)
987          + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi    987          + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi
988          - fPhiTwist*(fTAlph*p.x() + p.y()))*     988          - fPhiTwist*(fTAlph*p.x() + p.y()))* std::cos(phi)
989          + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*    989          + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*phi + fPhiTwist*(p.x()
990          - fTAlph*p.y()))*std::sin(phi))/(2.*(    990          - fTAlph*p.y()))*std::sin(phi))/(2.*(fPhiTwist
991          + fPhiTwist*fTAlph*fTAlph)) ;            991          + fPhiTwist*fTAlph*fTAlph)) ;
992 }                                                 992 }
993                                                   993 
994                                                   994 
995 G4ThreeVector G4TwistBoxSide::ProjectPoint(con    995 G4ThreeVector G4TwistBoxSide::ProjectPoint(const G4ThreeVector& p, 
996                                                   996                                                  G4bool isglobal) 
997 {                                                 997 {
998   // Get Rho at p.z() on Hyperbolic Surface.      998   // Get Rho at p.z() on Hyperbolic Surface.
999   G4ThreeVector tmpp;                             999   G4ThreeVector tmpp;
1000   if (isglobal)                                  1000   if (isglobal)
1001   {                                              1001   {
1002      tmpp = fRot.inverse()*p - fTrans;           1002      tmpp = fRot.inverse()*p - fTrans;
1003   }                                              1003   }
1004   else                                           1004   else
1005   {                                              1005   {
1006      tmpp = p;                                   1006      tmpp = p;
1007   }                                              1007   }
1008                                                  1008 
1009   G4double phi ;                                 1009   G4double phi ;
1010   G4double u ;                                   1010   G4double u ;
1011                                                  1011 
1012   GetPhiUAtX( tmpp, phi, u ) ;                   1012   GetPhiUAtX( tmpp, phi, u ) ;
1013     // calculate (phi, u) for a point p close    1013     // calculate (phi, u) for a point p close the surface
1014                                                  1014   
1015   G4ThreeVector xx = SurfacePoint(phi,u) ;       1015   G4ThreeVector xx = SurfacePoint(phi,u) ;
1016     // transform back to cartesian coordinate    1016     // transform back to cartesian coordinates
1017                                                  1017 
1018   if (isglobal)                                  1018   if (isglobal)
1019   {                                              1019   {
1020      return (fRot * xx + fTrans);                1020      return (fRot * xx + fTrans);
1021   }                                              1021   }
1022   else                                           1022   else
1023   {                                              1023   {
1024      return xx;                                  1024      return xx;
1025   }                                              1025   }
1026 }                                                1026 }
1027                                                  1027 
1028 void G4TwistBoxSide::GetFacets( G4int k, G4in    1028 void G4TwistBoxSide::GetFacets( G4int k, G4int n,  G4double xyz[][3],
1029                                 G4int faces[]    1029                                 G4int faces[][4], G4int iside ) 
1030 {                                                1030 {
1031   G4double phi ;                                 1031   G4double phi ;
1032   G4double b ;                                   1032   G4double b ;   
1033                                                  1033 
1034   G4double z, u ;     // the two parameters f    1034   G4double z, u ;     // the two parameters for the surface equation
1035   G4ThreeVector p ;  // a point on the surfac    1035   G4ThreeVector p ;  // a point on the surface, given by (z,u)
1036                                                  1036 
1037   G4int nnode ;                                  1037   G4int nnode ;
1038   G4int nface ;                                  1038   G4int nface ;
1039                                                  1039 
1040   // calculate the (n-1)*(k-1) vertices          1040   // calculate the (n-1)*(k-1) vertices
1041                                                  1041 
1042   G4int i,j ;                                    1042   G4int i,j ;
1043                                                  1043 
1044   for ( i = 0 ; i<n ; ++i )                      1044   for ( i = 0 ; i<n ; ++i )
1045   {                                              1045   {
1046     z = -fDz+i*(2.*fDz)/(n-1) ;                  1046     z = -fDz+i*(2.*fDz)/(n-1) ;
1047     phi = z*fPhiTwist/(2*fDz) ;                  1047     phi = z*fPhiTwist/(2*fDz) ;
1048     b = GetValueB(phi) ;                         1048     b = GetValueB(phi) ;
1049                                                  1049 
1050     for ( j = 0 ; j<k ; ++j )                    1050     for ( j = 0 ; j<k ; ++j )
1051     {                                            1051     {
1052       nnode = GetNode(i,j,k,n,iside) ;           1052       nnode = GetNode(i,j,k,n,iside) ;
1053       u = -b/2 +j*b/(k-1) ;                      1053       u = -b/2 +j*b/(k-1) ;
1054       p = SurfacePoint(phi,u,true) ;             1054       p = SurfacePoint(phi,u,true) ;
1055        // surface point in global coordinate     1055        // surface point in global coordinate system
1056                                                  1056 
1057       xyz[nnode][0] = p.x() ;                    1057       xyz[nnode][0] = p.x() ;
1058       xyz[nnode][1] = p.y() ;                    1058       xyz[nnode][1] = p.y() ;
1059       xyz[nnode][2] = p.z() ;                    1059       xyz[nnode][2] = p.z() ;
1060                                                  1060 
1061       if ( i<n-1 && j<k-1 )     // contercloc    1061       if ( i<n-1 && j<k-1 )     // conterclock wise filling
1062       {                                          1062       {
1063         nface = GetFace(i,j,k,n,iside) ;         1063         nface = GetFace(i,j,k,n,iside) ;
1064         faces[nface][0] = GetEdgeVisibility(i    1064         faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i  ,j  ,k,n,iside)+1) ;  // fortran numbering
1065         faces[nface][1] = GetEdgeVisibility(i    1065         faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i  ,j+1,k,n,iside)+1) ;
1066         faces[nface][2] = GetEdgeVisibility(i    1066         faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ;
1067         faces[nface][3] = GetEdgeVisibility(i    1067         faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j  ,k,n,iside)+1) ;
1068       }                                          1068       }
1069     }                                            1069     }
1070   }                                              1070   }
1071 }                                                1071 }
1072                                                  1072