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 9.0)


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