Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/specific/src/G4TwistTrapParallelSide.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/G4TwistTrapParallelSide.cc (Version 11.3.0) and /geometry/solids/specific/src/G4TwistTrapParallelSide.cc (Version 10.5.p1)


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