Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TwistTrapAlphaSide.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/G4TwistTrapAlphaSide.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TwistTrapAlphaSide.cc (Version 9.0.p2)


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