Geant4 Cross Reference

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


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // G4TwistTubsHypeSide implementation          << 
 27 //                                                 26 //
 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu <<  27 // $Id: G4TwistTubsHypeSide.cc 67011 2013-01-29 16:17:41Z gcosmo $
 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), <<  28 //
 30 //               from original version in Jupi <<  29 // 
                                                   >>  30 // --------------------------------------------------------------------
                                                   >>  31 // GEANT 4 class source file
                                                   >>  32 //
                                                   >>  33 //
                                                   >>  34 // G4TwistTubsHypeSide.cc
                                                   >>  35 //
                                                   >>  36 // Author: 
                                                   >>  37 //   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
                                                   >>  38 //
                                                   >>  39 // History:
                                                   >>  40 //   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
                                                   >>  41 //                 from original version in Jupiter-2.5.02 application.
 31 // -------------------------------------------     42 // --------------------------------------------------------------------
 32                                                    43 
 33 #include "G4TwistTubsHypeSide.hh"                  44 #include "G4TwistTubsHypeSide.hh"
 34 #include "G4PhysicalConstants.hh"                  45 #include "G4PhysicalConstants.hh"
 35 #include "G4GeometryTolerance.hh"                  46 #include "G4GeometryTolerance.hh"
 36                                                    47 
 37 //============================================     48 //=====================================================================
 38 //* constructors -----------------------------     49 //* constructors ------------------------------------------------------
 39                                                    50 
 40 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const <<  51 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String         &name,
 41                                          const <<  52                                          const G4RotationMatrix &rot,
 42                                          const <<  53                                          const G4ThreeVector    &tlate,
 43                                          const     54                                          const G4int             handedness,
 44                                          const     55                                          const G4double          kappa,
 45                                          const     56                                          const G4double          tanstereo,
 46                                          const     57                                          const G4double          r0,
 47                                          const     58                                          const EAxis             axis0,
 48                                          const     59                                          const EAxis             axis1,
 49                                                    60                                                G4double          axis0min,
 50                                                    61                                                G4double          axis1min,
 51                                                    62                                                G4double          axis0max,
 52                                                    63                                                G4double          axis1max )
 53   : G4VTwistSurface(name, rot, tlate, handedne     64   : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
 54                     axis0min, axis1min, axis0m <<  65                    axis0min, axis1min, axis0max, axis1max),
 55     fKappa(kappa), fTanStereo(tanstereo),          66     fKappa(kappa), fTanStereo(tanstereo),
 56     fTan2Stereo(tanstereo*tanstereo), fR0(r0),     67     fTan2Stereo(tanstereo*tanstereo), fR0(r0), fR02(r0*r0), fDPhi(twopi)
 57 {                                                  68 {
 58    if ( (axis0 == kZAxis) && (axis1 == kPhi) )     69    if ( (axis0 == kZAxis) && (axis1 == kPhi) )
 59    {                                               70    {
 60       G4Exception("G4TwistTubsHypeSide::G4Twis     71       G4Exception("G4TwistTubsHypeSide::G4TwistTubsHypeSide()",
 61                   "GeomSolids0002", FatalError     72                   "GeomSolids0002", FatalErrorInArgument,
 62                   "Should swap axis0 and axis1     73                   "Should swap axis0 and axis1!");
 63    }                                               74    }
                                                   >>  75    
 64    fInside.gp.set(kInfinity, kInfinity, kInfin     76    fInside.gp.set(kInfinity, kInfinity, kInfinity);
 65    fInside.inside = kOutside;                      77    fInside.inside = kOutside;
 66    fIsValidNorm = false;                           78    fIsValidNorm = false;
                                                   >>  79    
 67    SetCorners();                                   80    SetCorners();
 68    SetBoundaries();                                81    SetBoundaries();
                                                   >>  82 
 69 }                                                  83 }
 70                                                    84 
 71 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const <<  85 G4TwistTubsHypeSide::G4TwistTubsHypeSide(const G4String      &name,
 72                                                <<  86                                          G4double         EndInnerRadius[2],
 73                                                <<  87                                          G4double         EndOuterRadius[2],
 74                                                <<  88                                          G4double         DPhi,
 75                                                <<  89                                          G4double         EndPhi[2],
 76                                                <<  90                                          G4double         EndZ[2], 
 77                                                <<  91                                          G4double         InnerRadius,
 78                                                <<  92                                          G4double         OuterRadius,
 79                                                <<  93                                          G4double         Kappa,
 80                                                <<  94                                          G4double         TanInnerStereo,
 81                                                <<  95                                          G4double         TanOuterStereo,
 82                                                <<  96                                          G4int            handedness)
 83    : G4VTwistSurface(name)                         97    : G4VTwistSurface(name)
 84 {                                                  98 {
 85                                                    99 
 86    fHandedness = handedness;   // +z = +ve, -z    100    fHandedness = handedness;   // +z = +ve, -z = -ve
 87    fAxis[0]    = kPhi;                            101    fAxis[0]    = kPhi;
 88    fAxis[1]    = kZAxis;                          102    fAxis[1]    = kZAxis;
 89    fAxisMin[0] = kInfinity;         // we cann    103    fAxisMin[0] = kInfinity;         // we cannot fix boundary min of Phi, 
 90    fAxisMax[0] = kInfinity;         // because    104    fAxisMax[0] = kInfinity;         // because it depends on z.
 91    fAxisMin[1] = EndZ[0];                         105    fAxisMin[1] = EndZ[0];
 92    fAxisMax[1] = EndZ[1];                         106    fAxisMax[1] = EndZ[1];
 93    fKappa      = Kappa;                           107    fKappa      = Kappa;
 94    fDPhi       = DPhi ;                           108    fDPhi       = DPhi ;
 95                                                   109 
 96    if (handedness < 0)  // inner hyperbolic su << 110    if (handedness < 0) { // inner hyperbolic surface
 97    {                                           << 
 98       fTanStereo  = TanInnerStereo;               111       fTanStereo  = TanInnerStereo;
 99       fR0         = InnerRadius;                  112       fR0         = InnerRadius;
100    }                                           << 113    } else {              // outer hyperbolic surface
101    else                 // outer hyperbolic su << 
102    {                                           << 
103       fTanStereo  = TanOuterStereo;               114       fTanStereo  = TanOuterStereo;
104       fR0         = OuterRadius;                  115       fR0         = OuterRadius;
105    }                                              116    }
106    fTan2Stereo = fTanStereo * fTanStereo;         117    fTan2Stereo = fTanStereo * fTanStereo;
107    fR02        = fR0 * fR0;                       118    fR02        = fR0 * fR0;
108                                                   119    
109    fTrans.set(0, 0, 0);                           120    fTrans.set(0, 0, 0);
110    fIsValidNorm = false;                          121    fIsValidNorm = false;
111                                                   122 
112    fInside.gp.set(kInfinity, kInfinity, kInfin    123    fInside.gp.set(kInfinity, kInfinity, kInfinity);
113    fInside.inside = kOutside;                     124    fInside.inside = kOutside;
114                                                   125    
115    SetCorners(EndInnerRadius, EndOuterRadius,     126    SetCorners(EndInnerRadius, EndOuterRadius, DPhi, EndPhi, EndZ) ; 
116                                                   127 
117    SetBoundaries();                               128    SetBoundaries();
118 }                                                 129 }
119                                                   130 
120 //============================================    131 //=====================================================================
121 //* Fake default constructor -----------------    132 //* Fake default constructor ------------------------------------------
122                                                   133 
123 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __vo    134 G4TwistTubsHypeSide::G4TwistTubsHypeSide( __void__& a )
124   : G4VTwistSurface(a)                         << 135   : G4VTwistSurface(a), fKappa(0.), fTanStereo(0.), fTan2Stereo(0.),
                                                   >> 136     fR0(0.), fR02(0.), fDPhi(0.)
125 {                                                 137 {
126 }                                                 138 }
127                                                   139 
128 //============================================    140 //=====================================================================
129 //* destructor -------------------------------    141 //* destructor --------------------------------------------------------
130                                                   142 
131 G4TwistTubsHypeSide::~G4TwistTubsHypeSide() =  << 143 G4TwistTubsHypeSide::~G4TwistTubsHypeSide()
                                                   >> 144 {
                                                   >> 145 }
132                                                   146 
133 //============================================    147 //=====================================================================
134 //* GetNormal --------------------------------    148 //* GetNormal ---------------------------------------------------------
135                                                   149 
136 G4ThreeVector G4TwistTubsHypeSide::GetNormal(c << 150 G4ThreeVector G4TwistTubsHypeSide::GetNormal(const G4ThreeVector &tmpxx, 
137                                                   151                                                    G4bool isGlobal) 
138 {                                                 152 {
139    // GetNormal returns a normal vector at a s    153    // GetNormal returns a normal vector at a surface (or very close
140    // to surface) point at tmpxx.                 154    // to surface) point at tmpxx.
141    // If isGlobal=true, it returns the normal     155    // If isGlobal=true, it returns the normal in global coordinate.
                                                   >> 156    //
142                                                   157    
143    G4ThreeVector xx;                              158    G4ThreeVector xx;
144    if (isGlobal)                               << 159    if (isGlobal) {
145    {                                           << 
146       xx = ComputeLocalPoint(tmpxx);              160       xx = ComputeLocalPoint(tmpxx);
147       if ((xx - fCurrentNormal.p).mag() < 0.5  << 161       if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
148       {                                        << 
149          return ComputeGlobalDirection(fCurren    162          return ComputeGlobalDirection(fCurrentNormal.normal);
150       }                                           163       }
151    }                                           << 164    } else {
152    else                                        << 
153    {                                           << 
154       xx = tmpxx;                                 165       xx = tmpxx;
155       if (xx == fCurrentNormal.p)              << 166       if (xx == fCurrentNormal.p) {
156       {                                        << 
157          return fCurrentNormal.normal;            167          return fCurrentNormal.normal;
158       }                                           168       }
159    }                                              169    }
160                                                   170    
161    fCurrentNormal.p = xx;                         171    fCurrentNormal.p = xx;
162                                                   172 
163    G4ThreeVector normal( xx.x(), xx.y(), -xx.z    173    G4ThreeVector normal( xx.x(), xx.y(), -xx.z() * fTan2Stereo);
164    normal *= fHandedness;                         174    normal *= fHandedness;
165    normal = normal.unit();                        175    normal = normal.unit();
166                                                   176 
167    if (isGlobal)                               << 177    if (isGlobal) {
168    {                                           << 
169       fCurrentNormal.normal = ComputeLocalDire    178       fCurrentNormal.normal = ComputeLocalDirection(normal);
170    }                                           << 179    } else {
171    else                                        << 
172    {                                           << 
173       fCurrentNormal.normal = normal;             180       fCurrentNormal.normal = normal;
174    }                                              181    }
175    return fCurrentNormal.normal;                  182    return fCurrentNormal.normal;
176 }                                                 183 }
177                                                   184 
178 //============================================    185 //=====================================================================
179 //* Inside() ---------------------------------    186 //* Inside() ----------------------------------------------------------
180                                                   187 
181 EInside G4TwistTubsHypeSide::Inside(const G4Th << 188 EInside G4TwistTubsHypeSide::Inside(const G4ThreeVector &gp) 
182 {                                                 189 {
183    // Inside returns                              190    // Inside returns 
184    const G4double halftol                      << 191    static const G4double halftol
185      = 0.5 * G4GeometryTolerance::GetInstance(    192      = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
186                                                   193 
187    if (fInside.gp == gp)                       << 194    if (fInside.gp == gp) {
188    {                                           << 
189       return fInside.inside;                      195       return fInside.inside;
190    }                                              196    }
191    fInside.gp = gp;                               197    fInside.gp = gp;
192                                                   198    
193    G4ThreeVector p = ComputeLocalPoint(gp);       199    G4ThreeVector p = ComputeLocalPoint(gp);
194                                                   200   
195                                                   201 
196    if (p.mag() < DBL_MIN)                      << 202    if (p.mag() < DBL_MIN) {
197    {                                           << 
198       fInside.inside = kOutside;                  203       fInside.inside = kOutside;
199       return fInside.inside;                      204       return fInside.inside;
200    }                                              205    }
201                                                   206    
202    G4double rhohype = GetRhoAtPZ(p);              207    G4double rhohype = GetRhoAtPZ(p);
203    G4double distanceToOut = fHandedness * (rho    208    G4double distanceToOut = fHandedness * (rhohype - p.getRho());
204                             // +ve : inside       209                             // +ve : inside
205                                                   210 
206    if (distanceToOut < -halftol)               << 211    if (distanceToOut < -halftol) {
207    {                                           << 212 
208      fInside.inside = kOutside;                   213      fInside.inside = kOutside;
209    }                                           << 214 
210    else                                        << 215    } else {
211    {                                           << 216 
212       G4int areacode = GetAreaCode(p);            217       G4int areacode = GetAreaCode(p);
213       if (IsOutside(areacode))                 << 218       if (IsOutside(areacode)) {
214       {                                        << 
215          fInside.inside = kOutside;               219          fInside.inside = kOutside;
216       }                                        << 220       } else if (IsBoundary(areacode)) {
217       else if (IsBoundary(areacode))           << 
218       {                                        << 
219          fInside.inside = kSurface;               221          fInside.inside = kSurface;
220       }                                        << 222       } else if (IsInside(areacode)) {
221       else if (IsInside(areacode))             << 223          if (distanceToOut <= halftol) {
222       {                                        << 
223          if (distanceToOut <= halftol)         << 
224          {                                     << 
225             fInside.inside = kSurface;            224             fInside.inside = kSurface;
226          }                                     << 225          } else {
227          else                                  << 
228          {                                     << 
229             fInside.inside = kInside;             226             fInside.inside = kInside;
230          }                                        227          }
231       }                                        << 228       } else {
232       else                                     << 
233       {                                        << 
234          G4cout << "WARNING - G4TwistTubsHypeS    229          G4cout << "WARNING - G4TwistTubsHypeSide::Inside()" << G4endl
235                 << "          Invalid option !    230                 << "          Invalid option !" << G4endl
236                 << "          name, areacode,     231                 << "          name, areacode, distanceToOut = "
237                 << GetName() << ", " << std::h << 232                 << GetName() << ", " << std::hex << areacode << std::dec << ", "
238                 << std::dec << ", " << distanc << 233                 << distanceToOut << G4endl;
239       }                                           234       }
240    }                                              235    }
241                                                << 236    
242    return fInside.inside;                         237    return fInside.inside; 
243 }                                                 238 }
244                                                   239 
245 //============================================    240 //=====================================================================
246 //* DistanceToSurface ------------------------    241 //* DistanceToSurface -------------------------------------------------
247                                                   242 
248 G4int G4TwistTubsHypeSide::DistanceToSurface(c << 243 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector &gp,
249                                              c << 244                                              const G4ThreeVector &gv,
250                                                   245                                                    G4ThreeVector  gxx[],
251                                                   246                                                    G4double       distance[],
252                                                   247                                                    G4int          areacode[],
253                                                   248                                                    G4bool         isvalid[],
254                                                   249                                                    EValidate      validate)
255 {                                                 250 {
                                                   >> 251    //
256    // Decide if and where a line intersects wi    252    // Decide if and where a line intersects with a hyperbolic
257    // surface (of infinite extent)                253    // surface (of infinite extent)
258    //                                             254    //
259    // Arguments:                                  255    // Arguments:
260    //     p       - (in) Point on trajectory      256    //     p       - (in) Point on trajectory
261    //     v       - (in) Vector along trajecto    257    //     v       - (in) Vector along trajectory
262    //     r2      - (in) Square of radius at z    258    //     r2      - (in) Square of radius at z = 0
263    //     tan2phi - (in) std::tan(stereo)**2      259    //     tan2phi - (in) std::tan(stereo)**2
264    //     s       - (out) Up to two points of     260    //     s       - (out) Up to two points of intersection, where the
265    //                     intersection point i    261    //                     intersection point is p + s*v, and if there are
266    //                     two intersections, s    262    //                     two intersections, s[0] < s[1]. May be negative.
267    // Returns:                                    263    // Returns:
268    //     The number of intersections. If 0, t    264    //     The number of intersections. If 0, the trajectory misses.
269    //                                             265    //
270    //                                             266    //
271    // Equation of a line:                         267    // Equation of a line:
272    //                                             268    //
273    //       x = x0 + s*tx      y = y0 + s*ty      269    //       x = x0 + s*tx      y = y0 + s*ty      z = z0 + s*tz
274    //                                             270    //
275    // Equation of a hyperbolic surface:           271    // Equation of a hyperbolic surface:
276    //                                             272    //
277    //       x**2 + y**2 = r**2 + (z*tanPhi)**2    273    //       x**2 + y**2 = r**2 + (z*tanPhi)**2
278    //                                             274    //
279    // Solution is quadratic:                      275    // Solution is quadratic:
280    //                                             276    //
281    //  a*s**2 + b*s + c = 0                       277    //  a*s**2 + b*s + c = 0
282    //                                             278    //
283    // where:                                      279    // where:
284    //                                             280    //
285    //  a = tx**2 + ty**2 - (tz*tanPhi)**2         281    //  a = tx**2 + ty**2 - (tz*tanPhi)**2
286    //                                             282    //
287    //  b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2    283    //  b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
288    //                                             284    //
289    //  c = x0**2 + y0**2 - r**2 - (z0*tanPhi)*    285    //  c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
290    //                                             286    //
291                                                   287       
292    fCurStatWithV.ResetfDone(validate, &gp, &gv    288    fCurStatWithV.ResetfDone(validate, &gp, &gv);
293                                                   289 
294    if (fCurStatWithV.IsDone())                 << 290    if (fCurStatWithV.IsDone()) {
295    {                                           << 291       G4int i;
296       for (G4int i=0; i<fCurStatWithV.GetNXX() << 292       for (i=0; i<fCurStatWithV.GetNXX(); i++) {
297       {                                        << 
298          gxx[i] = fCurStatWithV.GetXX(i);         293          gxx[i] = fCurStatWithV.GetXX(i);
299          distance[i] = fCurStatWithV.GetDistan    294          distance[i] = fCurStatWithV.GetDistance(i);
300          areacode[i] = fCurStatWithV.GetAreaco    295          areacode[i] = fCurStatWithV.GetAreacode(i);
301          isvalid[i]  = fCurStatWithV.IsValid(i    296          isvalid[i]  = fCurStatWithV.IsValid(i);
302       }                                           297       }
303       return fCurStatWithV.GetNXX();              298       return fCurStatWithV.GetNXX();
304    }                                           << 299    } else {
305    else   // initialize                        << 300       // initialize
306    {                                           << 301       G4int i;
307       for (auto i=0; i<2; ++i)                 << 302       for (i=0; i<2; i++) {
308       {                                        << 
309          distance[i] = kInfinity;                 303          distance[i] = kInfinity;
310          areacode[i] = sOutside;                  304          areacode[i] = sOutside;
311          isvalid[i]  = false;                     305          isvalid[i]  = false;
312          gxx[i].set(kInfinity, kInfinity, kInf    306          gxx[i].set(kInfinity, kInfinity, kInfinity);
313       }                                           307       }
314    }                                              308    }
315                                                   309    
316    G4ThreeVector p = ComputeLocalPoint(gp);       310    G4ThreeVector p = ComputeLocalPoint(gp);
317    G4ThreeVector v = ComputeLocalDirection(gv)    311    G4ThreeVector v = ComputeLocalDirection(gv);
318    G4ThreeVector xx[2];                           312    G4ThreeVector xx[2]; 
319                                                   313    
320    //                                             314    //
321    // special case!  p is on origin.              315    // special case!  p is on origin.
322    //                                             316    // 
323                                                   317 
324    if (p.mag() == 0)                           << 318    if (p.mag() == 0) {
325    {                                           << 
326       // p is origin.                             319       // p is origin. 
327       // unique solution of 2-dimension questi    320       // unique solution of 2-dimension question in r-z plane 
328       // Equations:                               321       // Equations:
329       //    r^2 = fR02 + z^2*fTan2Stere0          322       //    r^2 = fR02 + z^2*fTan2Stere0
330       //    r = beta*z                            323       //    r = beta*z
331       //        where                             324       //        where 
332       //        beta = vrho / vz                  325       //        beta = vrho / vz
333       // Solution (z value of intersection poi    326       // Solution (z value of intersection point):
334       //    xxz = +- std::sqrt (fR02 / (beta^2    327       //    xxz = +- std::sqrt (fR02 / (beta^2 - fTan2Stereo))
335       //                                          328       //
336                                                   329 
337       G4double vz    = v.z();                     330       G4double vz    = v.z();
338       G4double absvz = std::fabs(vz);             331       G4double absvz = std::fabs(vz);
339       G4double vrho  = v.getRho();                332       G4double vrho  = v.getRho();       
340       G4double vslope = vrho/vz;                  333       G4double vslope = vrho/vz;
341       G4double vslope2 = vslope * vslope;         334       G4double vslope2 = vslope * vslope;
342       if (vrho == 0 || (vrho/absvz) <= (absvz* << 335       if (vrho == 0 || (vrho/absvz) <= (absvz*std::fabs(fTanStereo)/absvz)) {
343       {                                        << 
344          // vz/vrho is bigger than slope of as    336          // vz/vrho is bigger than slope of asymptonic line
345          distance[0] = kInfinity;                 337          distance[0] = kInfinity;
346          fCurStatWithV.SetCurrentStatus(0, gxx    338          fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
347                                         isvali    339                                         isvalid[0], 0, validate, &gp, &gv);
348          return 0;                                340          return 0;
349       }                                           341       }
350                                                   342        
351       if (vz != 0.0)                           << 343       if (vz) { 
352       {                                        << 
353          G4double xxz  = std::sqrt(fR02 / (vsl    344          G4double xxz  = std::sqrt(fR02 / (vslope2 - fTan2Stereo)) 
354                         * (vz / std::fabs(vz))    345                         * (vz / std::fabs(vz)) ;
355          G4double t = xxz / vz;                   346          G4double t = xxz / vz;
356          xx[0].set(t*v.x(), t*v.y(), xxz);        347          xx[0].set(t*v.x(), t*v.y(), xxz);
357       }                                        << 348       } else {
358       else                                     << 
359       {                                        << 
360          // p.z = 0 && v.z =0                     349          // p.z = 0 && v.z =0
361          xx[0].set(v.x()*fR0, v.y()*fR0, 0);      350          xx[0].set(v.x()*fR0, v.y()*fR0, 0);  // v is a unit vector.
362       }                                           351       }
363       distance[0] = xx[0].mag();                  352       distance[0] = xx[0].mag();
364       gxx[0]      = ComputeGlobalPoint(xx[0]);    353       gxx[0]      = ComputeGlobalPoint(xx[0]);
365                                                   354 
366       if (validate == kValidateWithTol)        << 355       if (validate == kValidateWithTol) {
367       {                                        << 
368          areacode[0] = GetAreaCode(xx[0]);        356          areacode[0] = GetAreaCode(xx[0]);
369          if (!IsOutside(areacode[0]))          << 357          if (!IsOutside(areacode[0])) {
370          {                                     << 
371             if (distance[0] >= 0) isvalid[0] =    358             if (distance[0] >= 0) isvalid[0] = true;
372          }                                        359          }
373       }                                        << 360       } else if (validate == kValidateWithoutTol) {
374       else if (validate == kValidateWithoutTol << 
375       {                                        << 
376          areacode[0] = GetAreaCode(xx[0], fals    361          areacode[0] = GetAreaCode(xx[0], false);
377          if (IsInside(areacode[0]))            << 362          if (IsInside(areacode[0])) {
378          {                                     << 
379             if (distance[0] >= 0) isvalid[0] =    363             if (distance[0] >= 0) isvalid[0] = true;
380          }                                        364          }
381       }                                        << 365       } else { // kDontValidate                       
382       else  // kDontValidate                   << 
383       {                                        << 
384          areacode[0] = sInside;                   366          areacode[0] = sInside;
385             if (distance[0] >= 0) isvalid[0] =    367             if (distance[0] >= 0) isvalid[0] = true;
386       }                                           368       }
387                                                   369                  
388       fCurStatWithV.SetCurrentStatus(0, gxx[0]    370       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
389                                      isvalid[0 << 371                                         isvalid[0], 1, validate, &gp, &gv);
390       return 1;                                   372       return 1;
391    }                                              373    }
392                                                   374 
393    //                                             375    //
394    // special case end.                           376    // special case end.
395    //                                             377    // 
396                                                   378 
397    G4double a = v.x()*v.x() + v.y()*v.y() - v.    379    G4double a = v.x()*v.x() + v.y()*v.y() - v.z()*v.z()*fTan2Stereo;
398    G4double b = 2.0                            << 380    G4double b = 2.0 * ( p.x() * v.x() + p.y() * v.y() - p.z() * v.z() * fTan2Stereo );
399               * ( p.x() * v.x() + p.y() * v.y( << 
400    G4double c = p.x()*p.x() + p.y()*p.y() - fR    381    G4double c = p.x()*p.x() + p.y()*p.y() - fR02 - p.z()*p.z()*fTan2Stereo;
401    G4double D = b*b - 4*a*c;          //discri    382    G4double D = b*b - 4*a*c;          //discriminant
402    G4int vout = 0;                                383    G4int vout = 0;
403                                                   384    
404    if (std::fabs(a) < DBL_MIN)                 << 385    if (std::fabs(a) < DBL_MIN) {
405    {                                           << 386       if (std::fabs(b) > DBL_MIN) {           // single solution
406       if (std::fabs(b) > DBL_MIN)            / << 387 
407       {                                        << 
408          distance[0] = -c/b;                      388          distance[0] = -c/b;
409          xx[0] = p + distance[0]*v;               389          xx[0] = p + distance[0]*v;
410          gxx[0] = ComputeGlobalPoint(xx[0]);      390          gxx[0] = ComputeGlobalPoint(xx[0]);
411                                                   391 
412          if (validate == kValidateWithTol)     << 392          if (validate == kValidateWithTol) {
413          {                                     << 
414             areacode[0] = GetAreaCode(xx[0]);     393             areacode[0] = GetAreaCode(xx[0]);
415             if (!IsOutside(areacode[0]))       << 394             if (!IsOutside(areacode[0])) {
416             {                                  << 
417                if (distance[0] >= 0) isvalid[0    395                if (distance[0] >= 0) isvalid[0] = true;
418             }                                     396             }
419          }                                     << 397          } else if (validate == kValidateWithoutTol) {
420          else if (validate == kValidateWithout << 
421          {                                     << 
422             areacode[0] = GetAreaCode(xx[0], f    398             areacode[0] = GetAreaCode(xx[0], false);
423             if (IsInside(areacode[0]))         << 399             if (IsInside(areacode[0])) {
424             {                                  << 
425                if (distance[0] >= 0) isvalid[0    400                if (distance[0] >= 0) isvalid[0] = true;
426             }                                     401             }
427          }                                     << 402          } else { // kDontValidate                       
428          else  // kDontValidate                << 
429          {                                     << 
430             areacode[0] = sInside;                403             areacode[0] = sInside;
431                if (distance[0] >= 0) isvalid[0    404                if (distance[0] >= 0) isvalid[0] = true;
432          }                                        405          }
433                                                   406                  
434          fCurStatWithV.SetCurrentStatus(0, gxx    407          fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
435                                         isvali    408                                         isvalid[0], 1, validate, &gp, &gv);
436          vout = 1;                                409          vout = 1;
437       }                                        << 410          
438       else                                     << 411       } else {
439       {                                        << 412          // if a=b=0 and c != 0, p is origin and v is parallel to asymptotic line.
440         // if a=b=0 and c != 0, p is origin an << 413          // if a=b=c=0, p is on surface and v is paralell to stereo wire. 
441         // if a=b=c=0, p is on surface and v i << 414          // return distance = infinity.
442         // return distance = infinity          << 
443                                                   415 
444          fCurStatWithV.SetCurrentStatus(0, gxx    416          fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
445                                         isvali    417                                         isvalid[0], 0, validate, &gp, &gv);
                                                   >> 418 
446          vout = 0;                                419          vout = 0;
447       }                                           420       }
448    }                                           << 421       
449    else if (D > DBL_MIN)          // double so << 422    } else if (D > DBL_MIN) {         // double solutions
450    {                                           << 423       
451       D = std::sqrt(D);                           424       D = std::sqrt(D);
452       G4double      factor = 0.5/a;               425       G4double      factor = 0.5/a;
453       G4double      tmpdist[2] = {kInfinity, k    426       G4double      tmpdist[2] = {kInfinity, kInfinity};
454       G4ThreeVector tmpxx[2] ;                    427       G4ThreeVector tmpxx[2] ;
455       G4int         tmpareacode[2] = {sOutside    428       G4int         tmpareacode[2] = {sOutside, sOutside};
456       G4bool        tmpisvalid[2]  = {false, f    429       G4bool        tmpisvalid[2]  = {false, false};
                                                   >> 430       G4int i;
457                                                   431 
458       for (auto i=0; i<2; ++i)                 << 432       for (i=0; i<2; i++) {
459       {                                        << 
460          tmpdist[i] = factor*(-b - D);            433          tmpdist[i] = factor*(-b - D);
461          D = -D;                                  434          D = -D;
462          tmpxx[i] = p + tmpdist[i]*v;             435          tmpxx[i] = p + tmpdist[i]*v;
463                                                   436         
464          if (validate == kValidateWithTol)     << 437          if (validate == kValidateWithTol) {
465          {                                     << 
466             tmpareacode[i] = GetAreaCode(tmpxx    438             tmpareacode[i] = GetAreaCode(tmpxx[i]);
467             if (!IsOutside(tmpareacode[i]))    << 439             if (!IsOutside(tmpareacode[i])) {
468             {                                  << 
469                if (tmpdist[i] >= 0) tmpisvalid    440                if (tmpdist[i] >= 0) tmpisvalid[i] = true;
470                continue;                          441                continue;
471             }                                     442             }
472          }                                     << 443          } else if (validate == kValidateWithoutTol) {
473          else if (validate == kValidateWithout << 
474          {                                     << 
475             tmpareacode[i] = GetAreaCode(tmpxx    444             tmpareacode[i] = GetAreaCode(tmpxx[i], false);
476             if (IsInside(tmpareacode[i]))      << 445             if (IsInside(tmpareacode[i])) {
477             {                                  << 
478                if (tmpdist[i] >= 0) tmpisvalid    446                if (tmpdist[i] >= 0) tmpisvalid[i] = true;
479                continue;                          447                continue;
480             }                                     448             }
481          }                                     << 449          } else { // kDontValidate
482          else  // kDontValidate                << 
483          {                                     << 
484             tmpareacode[i] = sInside;             450             tmpareacode[i] = sInside;
485                if (tmpdist[i] >= 0) tmpisvalid    451                if (tmpdist[i] >= 0) tmpisvalid[i] = true;
486             continue;                             452             continue;
487          }                                        453          }
488       }                                           454       }      
489                                                   455 
490       if (tmpdist[0] <= tmpdist[1])            << 456       if (tmpdist[0] <= tmpdist[1]) {
491       {                                        << 
492           distance[0] = tmpdist[0];               457           distance[0] = tmpdist[0];
493           distance[1] = tmpdist[1];               458           distance[1] = tmpdist[1];
494           xx[0]       = tmpxx[0];                 459           xx[0]       = tmpxx[0];
495           xx[1]       = tmpxx[1];                 460           xx[1]       = tmpxx[1];
496           gxx[0]      = ComputeGlobalPoint(tmp    461           gxx[0]      = ComputeGlobalPoint(tmpxx[0]);
497           gxx[1]      = ComputeGlobalPoint(tmp    462           gxx[1]      = ComputeGlobalPoint(tmpxx[1]);
498           areacode[0] = tmpareacode[0];           463           areacode[0] = tmpareacode[0];
499           areacode[1] = tmpareacode[1];           464           areacode[1] = tmpareacode[1];
500           isvalid[0]  = tmpisvalid[0];            465           isvalid[0]  = tmpisvalid[0];
501           isvalid[1]  = tmpisvalid[1];            466           isvalid[1]  = tmpisvalid[1];
502       }                                        << 467       } else {
503       else                                     << 
504       {                                        << 
505           distance[0] = tmpdist[1];               468           distance[0] = tmpdist[1];
506           distance[1] = tmpdist[0];               469           distance[1] = tmpdist[0];
507           xx[0]       = tmpxx[1];                 470           xx[0]       = tmpxx[1];
508           xx[1]       = tmpxx[0];                 471           xx[1]       = tmpxx[0];
509           gxx[0]      = ComputeGlobalPoint(tmp    472           gxx[0]      = ComputeGlobalPoint(tmpxx[1]);
510           gxx[1]      = ComputeGlobalPoint(tmp    473           gxx[1]      = ComputeGlobalPoint(tmpxx[0]);
511           areacode[0] = tmpareacode[1];           474           areacode[0] = tmpareacode[1];
512           areacode[1] = tmpareacode[0];           475           areacode[1] = tmpareacode[0];
513           isvalid[0]  = tmpisvalid[1];            476           isvalid[0]  = tmpisvalid[1];
514           isvalid[1]  = tmpisvalid[0];            477           isvalid[1]  = tmpisvalid[0];
515       }                                           478       }
516                                                   479          
517       fCurStatWithV.SetCurrentStatus(0, gxx[0]    480       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
518                                      isvalid[0    481                                      isvalid[0], 2, validate, &gp, &gv);
519       fCurStatWithV.SetCurrentStatus(1, gxx[1]    482       fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
520                                      isvalid[1    483                                      isvalid[1], 2, validate, &gp, &gv);
521       vout = 2;                                   484       vout = 2;
522    }                                           << 485       
523    else                                        << 486    } else {
524    {                                           << 
525       // if D<0, no solution                      487       // if D<0, no solution
526       // if D=0, just grazing the surfaces, re    488       // if D=0, just grazing the surfaces, return kInfinity
527                                                   489 
528       fCurStatWithV.SetCurrentStatus(0, gxx[0]    490       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
529                                      isvalid[0    491                                      isvalid[0], 0, validate, &gp, &gv);
530       vout = 0;                                   492       vout = 0;
531    }                                              493    }
532    return vout;                                   494    return vout;
533 }                                                 495 }
534                                                   496 
                                                   >> 497    
535 //============================================    498 //=====================================================================
536 //* DistanceToSurface ------------------------    499 //* DistanceToSurface -------------------------------------------------
537                                                   500 
538 G4int G4TwistTubsHypeSide::DistanceToSurface(c << 501 G4int G4TwistTubsHypeSide::DistanceToSurface(const G4ThreeVector &gp,
539                                                   502                                                    G4ThreeVector  gxx[],
540                                                   503                                                    G4double       distance[],
541                                                   504                                                    G4int          areacode[])
542 {                                                 505 {
543     // Find the approximate distance of a poin    506     // Find the approximate distance of a point of a hyperbolic surface.
544     // The distance must be an underestimate.     507     // The distance must be an underestimate. 
545     // It will also be nice (although not nece    508     // It will also be nice (although not necessary) that the estimate is
546     // always finite no matter how close the p    509     // always finite no matter how close the point is.
547     //                                            510     //
548     // We arranged G4Hype::ApproxDistOutside a    511     // We arranged G4Hype::ApproxDistOutside and G4Hype::ApproxDistInside
549     // for this function. See these discriptio    512     // for this function. See these discriptions.
550                                                   513     
551    const G4double halftol                      << 514    static const G4double halftol
552      = 0.5 * G4GeometryTolerance::GetInstance(    515      = 0.5 * G4GeometryTolerance::GetInstance()->GetRadialTolerance();
553                                                   516 
554    fCurStat.ResetfDone(kDontValidate, &gp);       517    fCurStat.ResetfDone(kDontValidate, &gp);
555                                                   518 
556    if (fCurStat.IsDone())                      << 519    if (fCurStat.IsDone()) {
557    {                                           << 520       for (G4int i=0; i<fCurStat.GetNXX(); i++) {
558       for (G4int i=0; i<fCurStat.GetNXX(); ++i << 
559       {                                        << 
560          gxx[i] = fCurStat.GetXX(i);              521          gxx[i] = fCurStat.GetXX(i);
561          distance[i] = fCurStat.GetDistance(i)    522          distance[i] = fCurStat.GetDistance(i);
562          areacode[i] = fCurStat.GetAreacode(i)    523          areacode[i] = fCurStat.GetAreacode(i);
563       }                                           524       }
564       return fCurStat.GetNXX();                   525       return fCurStat.GetNXX();
565    }                                           << 526    } else {
566    else  // initialize                         << 527       // initialize
567    {                                           << 528       for (G4int i=0; i<2; i++) {
568       for (auto i=0; i<2; ++i)                 << 
569       {                                        << 
570          distance[i] = kInfinity;                 529          distance[i] = kInfinity;
571          areacode[i] = sOutside;                  530          areacode[i] = sOutside;
572          gxx[i].set(kInfinity, kInfinity, kInf    531          gxx[i].set(kInfinity, kInfinity, kInfinity);
573       }                                           532       }
574    }                                              533    }
575                                                   534    
576                                                   535 
577    G4ThreeVector p = ComputeLocalPoint(gp);       536    G4ThreeVector p = ComputeLocalPoint(gp);
578    G4ThreeVector xx;                              537    G4ThreeVector xx;
579                                                   538 
580    //                                             539    //
581    // special case!                               540    // special case!
582    // If p is on surface, return distance = 0     541    // If p is on surface, return distance = 0 immediatery .
583    //                                             542    //
584    G4ThreeVector  lastgxx[2];                     543    G4ThreeVector  lastgxx[2];
585    for (auto i=0; i<2; ++i)                    << 544    for (G4int i=0; i<2; i++) {
586    {                                           << 
587       lastgxx[i] = fCurStatWithV.GetXX(i);        545       lastgxx[i] = fCurStatWithV.GetXX(i);
588    }                                              546    }
589                                                   547 
590    if ((gp - lastgxx[0]).mag() < halftol || (g << 548    if ((gp - lastgxx[0]).mag() < halftol || (gp - lastgxx[1]).mag() < halftol) {
591    {                                           << 
592       // last winner, or last poststep point i    549       // last winner, or last poststep point is on the surface.
593       xx = p;                                     550       xx = p;             
594       gxx[0] = gp;                                551       gxx[0] = gp;
595       distance[0] = 0;                            552       distance[0] = 0;      
596                                                   553 
597       G4bool isvalid = true;                      554       G4bool isvalid = true;
598       fCurStat.SetCurrentStatus(0, gxx[0], dis    555       fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
599                                 isvalid, 1, kD    556                                 isvalid, 1, kDontValidate, &gp);
600                                                   557 
601       return 1;                                   558       return 1;
                                                   >> 559 
602    }                                              560    }
603    //                                             561    //
604    // special case end                            562    // special case end
605    //                                             563    //
606                                                   564        
607    G4double prho       = p.getRho();              565    G4double prho       = p.getRho();
608    G4double pz         = std::fabs(p.z());        566    G4double pz         = std::fabs(p.z());           // use symmetry
609    G4double r1         = std::sqrt(fR02 + pz *    567    G4double r1         = std::sqrt(fR02 + pz * pz * fTan2Stereo);
610                                                   568    
611    G4ThreeVector pabsz(p.x(), p.y(), pz);         569    G4ThreeVector pabsz(p.x(), p.y(), pz);
612                                                   570     
613    if (prho > r1 + halftol)   // p is outside  << 571    if (prho > r1 + halftol) {  // p is outside of Hyperbolic surface
614    {                                           << 572 
615       // First point xx1                          573       // First point xx1
616       G4double t = r1 / prho;                     574       G4double t = r1 / prho;
617       G4ThreeVector xx1(t * pabsz.x(), t * pab    575       G4ThreeVector xx1(t * pabsz.x(), t * pabsz.y() , pz);
618                                                   576       
619       // Second point xx2                         577       // Second point xx2
620       G4double z2 = (prho * fTanStereo + pz) /    578       G4double z2 = (prho * fTanStereo + pz) / (1 + fTan2Stereo);
621       G4double r2 = std::sqrt(fR02 + z2 * z2 *    579       G4double r2 = std::sqrt(fR02 + z2 * z2 * fTan2Stereo);
622       t = r2 / prho;                              580       t = r2 / prho;
623       G4ThreeVector xx2(t * pabsz.x(), t * pab    581       G4ThreeVector xx2(t * pabsz.x(), t * pabsz.y() , z2);
624                                                   582             
625       G4double len = (xx2 - xx1).mag();           583       G4double len = (xx2 - xx1).mag();
626       if (len < DBL_MIN)                       << 584       if (len < DBL_MIN) {
627       {                                        << 
628          // xx2 = xx1?? I guess we                585          // xx2 = xx1?? I guess we
629          // must have really bracketed the nor    586          // must have really bracketed the normal
630          distance[0] = (pabsz - xx1).mag();       587          distance[0] = (pabsz - xx1).mag();
631          xx = xx1;                                588          xx = xx1;
632       }                                        << 589       } else {
633       else                                     << 
634       {                                        << 
635          distance[0] = DistanceToLine(pabsz, x    590          distance[0] = DistanceToLine(pabsz, xx1, (xx2 - xx1) , xx);
636       }                                           591       }
637                                                   592       
638    }                                           << 593    } else if (prho < r1 - halftol) { // p is inside of Hyperbolic surface.
639    else if (prho < r1 - halftol)  // p is insi << 594            
640    {                                           << 
641       // First point xx1                          595       // First point xx1
642       G4double t;                                 596       G4double t;
643       G4ThreeVector xx1;                          597       G4ThreeVector xx1;
644       if (prho < DBL_MIN)                      << 598       if (prho < DBL_MIN) {
645       {                                        << 
646          xx1.set(r1, 0. , pz);                    599          xx1.set(r1, 0. , pz);
647       }                                        << 600       } else {
648       else                                     << 
649       {                                        << 
650          t = r1 / prho;                           601          t = r1 / prho;
651          xx1.set(t * pabsz.x(), t * pabsz.y()     602          xx1.set(t * pabsz.x(), t * pabsz.y() , pz);
652       }                                           603       }
653                                                   604       
654       // dr, dz is tangential vector of Hyparb    605       // dr, dz is tangential vector of Hyparbolic surface at xx1
655       // dr = r, dz = z*tan2stereo                606       // dr = r, dz = z*tan2stereo
656       G4double dr  = pz * fTan2Stereo;            607       G4double dr  = pz * fTan2Stereo;
657       G4double dz  = r1;                          608       G4double dz  = r1;
658       G4double tanbeta   = dr / dz;               609       G4double tanbeta   = dr / dz;
659       G4double pztanbeta = pz * tanbeta;          610       G4double pztanbeta = pz * tanbeta;
660                                                   611       
661       // Second point xx2                         612       // Second point xx2 
662       // xx2 is intersection between x-axis an    613       // xx2 is intersection between x-axis and tangential vector
663       G4double r2 = r1 - pztanbeta;               614       G4double r2 = r1 - pztanbeta;
664       G4ThreeVector xx2;                          615       G4ThreeVector xx2;
665       if (prho < DBL_MIN)                      << 616       if (prho < DBL_MIN) {
666       {                                        << 
667          xx2.set(r2, 0. , 0.);                    617          xx2.set(r2, 0. , 0.);
668       }                                        << 618       } else {
669       else                                     << 
670       {                                        << 
671          t  = r2 / prho;                          619          t  = r2 / prho;
672          xx2.set(t * pabsz.x(), t * pabsz.y()     620          xx2.set(t * pabsz.x(), t * pabsz.y() , 0.);
673       }                                           621       }
674                                                   622       
675       G4ThreeVector d = xx2 - xx1;                623       G4ThreeVector d = xx2 - xx1;
676       distance[0] = DistanceToLine(pabsz, xx1,    624       distance[0] = DistanceToLine(pabsz, xx1, d, xx);
677                                                   625           
678    }                                           << 626    } else {  // p is on Hyperbolic surface.
679    else   // p is on Hyperbolic surface.       << 627    
680    {                                           << 
681       distance[0] = 0;                            628       distance[0] = 0;
682       xx.set(p.x(), p.y(), pz);                   629       xx.set(p.x(), p.y(), pz);
                                                   >> 630 
683    }                                              631    }
684                                                   632 
685    if (p.z() < 0)                              << 633    if (p.z() < 0) {
686    {                                           << 
687       G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.    634       G4ThreeVector tmpxx(xx.x(), xx.y(), -xx.z());
688       xx = tmpxx;                                 635       xx = tmpxx;
689    }                                              636    }
690                                                   637        
691    gxx[0] = ComputeGlobalPoint(xx);               638    gxx[0] = ComputeGlobalPoint(xx);
692    areacode[0]    = sInside;                      639    areacode[0]    = sInside;
693    G4bool isvalid = true;                         640    G4bool isvalid = true;
694    fCurStat.SetCurrentStatus(0, gxx[0], distan    641    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
695                              isvalid, 1, kDont    642                              isvalid, 1, kDontValidate, &gp);
696    return 1;                                      643    return 1;
697 }                                                 644 }
698                                                   645 
699 //============================================    646 //=====================================================================
700 //* GetAreaCode ------------------------------    647 //* GetAreaCode -------------------------------------------------------
701                                                   648 
702 G4int G4TwistTubsHypeSide::GetAreaCode(const G << 649 G4int G4TwistTubsHypeSide::GetAreaCode(const G4ThreeVector &xx, 
703                                              G    650                                              G4bool         withTol)
704 {                                                 651 {
705    const G4double ctol = 0.5 * kCarTolerance;  << 652    static const G4double ctol = 0.5 * kCarTolerance;
706    G4int areacode = sInside;                      653    G4int areacode = sInside;
707                                                   654 
708    if ((fAxis[0] == kPhi && fAxis[1] == kZAxis << 655    if ((fAxis[0] == kPhi && fAxis[1] == kZAxis))  {
709    {                                           << 656       //G4int phiaxis = 0;
710       G4int zaxis   = 1;                          657       G4int zaxis   = 1;
711                                                   658       
712       if (withTol)                             << 659       if (withTol) {
713       {                                        << 660 
714          G4bool isoutside      = false;           661          G4bool isoutside      = false;
715          G4int  phiareacode    = GetAreaCodeIn    662          G4int  phiareacode    = GetAreaCodeInPhi(xx);
716          G4bool isoutsideinphi = IsOutside(phi    663          G4bool isoutsideinphi = IsOutside(phiareacode);
717                                                   664 
718          // test boundary of phiaxis              665          // test boundary of phiaxis
719                                                   666 
720          if ((phiareacode & sAxisMin) == sAxis << 667          if ((phiareacode & sAxisMin) == sAxisMin) {
721          {                                     << 668 
722             areacode |= (sAxis0 & (sAxisPhi |     669             areacode |= (sAxis0 & (sAxisPhi | sAxisMin)) | sBoundary;
723             if (isoutsideinphi) isoutside = tr    670             if (isoutsideinphi) isoutside = true;
724          }                                     << 671 
725          else if ((phiareacode & sAxisMax)  == << 672          } else if ((phiareacode & sAxisMax)  == sAxisMax) {
726          {                                     << 673 
727             areacode |= (sAxis0 & (sAxisPhi |     674             areacode |= (sAxis0 & (sAxisPhi | sAxisMax)) | sBoundary;
728             if (isoutsideinphi) isoutside = tr    675             if (isoutsideinphi) isoutside = true;
                                                   >> 676 
729          }                                        677          }
730                                                   678 
731          // test boundary of zaxis                679          // test boundary of zaxis
732                                                   680 
733          if (xx.z() < fAxisMin[zaxis] + ctol)  << 681          if (xx.z() < fAxisMin[zaxis] + ctol) {
734          {                                     << 682 
735             areacode |= (sAxis1 & (sAxisZ | sA    683             areacode |= (sAxis1 & (sAxisZ | sAxisMin));
736             if   ((areacode & sBoundary) != 0) << 684             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
737             else                        areaco    685             else                        areacode |= sBoundary;
738                                                   686 
739             if (xx.z() <= fAxisMin[zaxis] - ct    687             if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
740                                                   688 
741          }                                     << 689          } else if (xx.z() > fAxisMax[zaxis] - ctol) {
742          else if (xx.z() > fAxisMax[zaxis] - c << 690 
743          {                                     << 
744             areacode |= (sAxis1 & (sAxisZ | sA    691             areacode |= (sAxis1 & (sAxisZ | sAxisMax));
745             if   ((areacode & sBoundary) != 0) << 692             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
746             else                        areaco    693             else                        areacode |= sBoundary;
747                                                   694 
748             if (xx.z() >= fAxisMax[zaxis] + ct    695             if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
749          }                                        696          }
750                                                   697 
751          // if isoutside = true, clear sInside    698          // if isoutside = true, clear sInside bit.
752          // if not on boundary, add boundary i    699          // if not on boundary, add boundary information. 
753                                                   700 
754          if (isoutside)                        << 701          if (isoutside) {
755          {                                     << 
756             G4int tmpareacode = areacode & (~s    702             G4int tmpareacode = areacode & (~sInside);
757             areacode = tmpareacode;               703             areacode = tmpareacode;
758          }                                     << 704          } else if ((areacode & sBoundary) != sBoundary) {
759          else if ((areacode & sBoundary) != sB << 
760          {                                     << 
761             areacode |= (sAxis0 & sAxisPhi) |     705             areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
762          }                                        706          }
763                                                   707 
764          return areacode;                         708          return areacode;
765       }                                        << 709       
766       else                                     << 710       } else {
767       {                                        << 711 
768          G4int phiareacode = GetAreaCodeInPhi(    712          G4int phiareacode = GetAreaCodeInPhi(xx, false);
769                                                   713          
770          // test boundary of z-axis               714          // test boundary of z-axis
771                                                   715 
772          if (xx.z() < fAxisMin[zaxis])         << 716          if (xx.z() < fAxisMin[zaxis]) {
773          {                                     << 717 
774             areacode |= (sAxis1 & (sAxisZ | sA    718             areacode |= (sAxis1 & (sAxisZ | sAxisMin)) | sBoundary;
775                                                   719 
776          }                                     << 720          } else if (xx.z() > fAxisMax[zaxis]) {
777          else if (xx.z() > fAxisMax[zaxis])    << 721 
778          {                                     << 
779             areacode |= (sAxis1 & (sAxisZ | sA    722             areacode |= (sAxis1 & (sAxisZ | sAxisMax)) | sBoundary;
                                                   >> 723 
780          }                                        724          }
781                                                   725 
782          // boundary of phi-axis                  726          // boundary of phi-axis
783                                                   727 
784          if (phiareacode == sAxisMin)          << 728          if (phiareacode == sAxisMin) {
785          {                                     << 729 
786             areacode |= (sAxis0 & (sAxisPhi |     730             areacode |= (sAxis0 & (sAxisPhi | sAxisMin));
787             if   ((areacode & sBoundary) != 0) << 731             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
788             else                        areaco    732             else                        areacode |= sBoundary; 
789                                                   733              
790          }                                     << 734          } else if (phiareacode == sAxisMax) {
791          else if (phiareacode == sAxisMax)     << 735 
792          {                                     << 
793             areacode |= (sAxis0 & (sAxisPhi |     736             areacode |= (sAxis0 & (sAxisPhi | sAxisMax));
794             if   ((areacode & sBoundary) != 0) << 737             if   (areacode & sBoundary) areacode |= sCorner;  // xx is on the corner.
795             else                        areaco << 738             else                        areacode |= sBoundary; 
                                                   >> 739            
796          }                                        740          }
797                                                   741 
798          // if not on boundary, add boundary i    742          // if not on boundary, add boundary information. 
799                                                   743 
800          if ((areacode & sBoundary) != sBounda << 744          if ((areacode & sBoundary) != sBoundary) {
801          {                                     << 
802             areacode |= (sAxis0 & sAxisPhi) |     745             areacode |= (sAxis0 & sAxisPhi) | (sAxis1 & sAxisZ);
803          }                                        746          }
804          return areacode;                         747          return areacode;
805       }                                           748       }
806    }                                           << 749    } else {
807    else                                        << 
808    {                                           << 
809       std::ostringstream message;                 750       std::ostringstream message;
810       message << "Feature NOT implemented !" <    751       message << "Feature NOT implemented !" << G4endl
811               << "        fAxis[0] = " << fAxi    752               << "        fAxis[0] = " << fAxis[0] << G4endl
812               << "        fAxis[1] = " << fAxi    753               << "        fAxis[1] = " << fAxis[1];
813       G4Exception("G4TwistTubsHypeSide::GetAre    754       G4Exception("G4TwistTubsHypeSide::GetAreaCode()",
814                   "GeomSolids0001", FatalExcep    755                   "GeomSolids0001", FatalException, message);
815    }                                              756    }
816    return areacode;                               757    return areacode;
817 }                                                 758 }
818                                                   759 
819 //============================================    760 //=====================================================================
820 //* GetAreaCodeInPhi -------------------------    761 //* GetAreaCodeInPhi --------------------------------------------------
821                                                   762 
822 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(co << 763 G4int G4TwistTubsHypeSide::GetAreaCodeInPhi(const G4ThreeVector &xx,
823                                                   764                                                   G4bool withTol)
824 {                                                 765 {
825                                                   766    
826    G4ThreeVector lowerlimit; // lower phi-boun    767    G4ThreeVector lowerlimit; // lower phi-boundary limit at z = xx.z()
827    G4ThreeVector upperlimit; // upper phi-boun    768    G4ThreeVector upperlimit; // upper phi-boundary limit at z = xx.z()
828    lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxis    769    lowerlimit = GetBoundaryAtPZ(sAxis0 & sAxisMin, xx);
829    upperlimit = GetBoundaryAtPZ(sAxis0 & sAxis    770    upperlimit = GetBoundaryAtPZ(sAxis0 & sAxisMax, xx);
830                                                   771 
831    G4int  areacode  = sInside;                    772    G4int  areacode  = sInside;
832    G4bool isoutside = false;                      773    G4bool isoutside = false; 
833                                                   774    
834    if (withTol)                                << 775    if (withTol) {
835    {                                           << 776          
836       if (AmIOnLeftSide(xx, lowerlimit) >= 0)  << 777       if (AmIOnLeftSide(xx, lowerlimit) >= 0) {        // xx is on lowerlimit
837       {                                        << 
838          areacode |= (sAxisMin | sBoundary);      778          areacode |= (sAxisMin | sBoundary);
839          if (AmIOnLeftSide(xx, lowerlimit) > 0    779          if (AmIOnLeftSide(xx, lowerlimit) > 0) isoutside = true; 
840                                                   780 
841       }                                        << 781       } else if (AmIOnLeftSide(xx, upperlimit) <= 0) { // xx is on upperlimit
842       else if (AmIOnLeftSide(xx, upperlimit) < << 
843       {                                        << 
844          areacode |= (sAxisMax | sBoundary);      782          areacode |= (sAxisMax | sBoundary);
845          if (AmIOnLeftSide(xx, upperlimit) < 0    783          if (AmIOnLeftSide(xx, upperlimit) < 0) isoutside = true; 
846       }                                           784       }
847                                                   785 
848       // if isoutside = true, clear inside bit    786       // if isoutside = true, clear inside bit.
849                                                   787 
850       if (isoutside)                           << 788       if (isoutside) {
851       {                                        << 
852          G4int tmpareacode = areacode & (~sIns    789          G4int tmpareacode = areacode & (~sInside);
853          areacode = tmpareacode;                  790          areacode = tmpareacode;
854       }                                           791       }
855    }                                           << 792 
856    else                                        << 793 
857    {                                           << 794    } else {
858       if (AmIOnLeftSide(xx, lowerlimit, false) << 795    
859       {                                        << 796       if (AmIOnLeftSide(xx, lowerlimit, false) >= 0) {
860          areacode |= (sAxisMin | sBoundary);      797          areacode |= (sAxisMin | sBoundary);
861       }                                        << 798       } else if (AmIOnLeftSide(xx, upperlimit, false) <= 0) {
862       else if (AmIOnLeftSide(xx, upperlimit, f << 
863       {                                        << 
864          areacode |= (sAxisMax | sBoundary);      799          areacode |= (sAxisMax | sBoundary);
865       }                                           800       }
866    }                                              801    }
867                                                   802 
868    return areacode;                               803    return areacode;
                                                   >> 804    
869 }                                                 805 }
870                                                   806 
871 //============================================    807 //=====================================================================
872 //* SetCorners(EndInnerRadius, EndOuterRadius,    808 //* SetCorners(EndInnerRadius, EndOuterRadius,DPhi,EndPhi,EndZ) -------
873                                                   809 
874 void G4TwistTubsHypeSide::SetCorners( G4double << 810 void G4TwistTubsHypeSide::SetCorners(
875                                       G4double << 811                                      G4double         EndInnerRadius[2],
876                                       G4double << 812                                      G4double         EndOuterRadius[2],
877                                       G4double << 813                                      G4double         DPhi,
878                                       G4double << 814                                      G4double         endPhi[2],
                                                   >> 815                                      G4double         endZ[2] 
                                                   >> 816                                      )
879 {                                                 817 {
880    // Set Corner points in local coodinate.       818    // Set Corner points in local coodinate.
881                                                   819 
882    if (fAxis[0] == kPhi && fAxis[1] == kZAxis)    820    if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
883                                                   821 
                                                   >> 822       G4int i;
884       G4double endRad[2];                         823       G4double endRad[2];
885       G4double halfdphi = 0.5*DPhi;               824       G4double halfdphi = 0.5*DPhi;
886                                                   825       
887       for (auto i=0; i<2; ++i) // i=0,1 : -ve  << 826       for (i=0; i<2; i++) { // i=0,1 : -ve z, +ve z
888       {                                        << 827          endRad[i] = (fHandedness == 1 ? EndOuterRadius[i]
889         endRad[i] = (fHandedness == 1 ? EndOut << 828                                       : EndInnerRadius[i]);
890       }                                           829       }
891                                                   830 
892       G4int zmin = 0 ;  // at -ve z               831       G4int zmin = 0 ;  // at -ve z
893       G4int zmax = 1 ;  // at +ve z               832       G4int zmax = 1 ;  // at +ve z
894                                                   833 
895       G4double x, y, z;                           834       G4double x, y, z;
896                                                   835       
897       // corner of Axis0min and Axis1min          836       // corner of Axis0min and Axis1min
898       x = endRad[zmin]*std::cos(endPhi[zmin] -    837       x = endRad[zmin]*std::cos(endPhi[zmin] - halfdphi);
899       y = endRad[zmin]*std::sin(endPhi[zmin] -    838       y = endRad[zmin]*std::sin(endPhi[zmin] - halfdphi);
900       z = endZ[zmin];                             839       z = endZ[zmin];
901       SetCorner(sC0Min1Min, x, y, z);             840       SetCorner(sC0Min1Min, x, y, z);
902                                                   841       
903       // corner of Axis0max and Axis1min          842       // corner of Axis0max and Axis1min
904       x = endRad[zmin]*std::cos(endPhi[zmin] +    843       x = endRad[zmin]*std::cos(endPhi[zmin] + halfdphi);
905       y = endRad[zmin]*std::sin(endPhi[zmin] +    844       y = endRad[zmin]*std::sin(endPhi[zmin] + halfdphi);
906       z = endZ[zmin];                             845       z = endZ[zmin];
907       SetCorner(sC0Max1Min, x, y, z);             846       SetCorner(sC0Max1Min, x, y, z);
908                                                   847       
909       // corner of Axis0max and Axis1max          848       // corner of Axis0max and Axis1max
910       x = endRad[zmax]*std::cos(endPhi[zmax] +    849       x = endRad[zmax]*std::cos(endPhi[zmax] + halfdphi);
911       y = endRad[zmax]*std::sin(endPhi[zmax] +    850       y = endRad[zmax]*std::sin(endPhi[zmax] + halfdphi);
912       z = endZ[zmax];                             851       z = endZ[zmax];
913       SetCorner(sC0Max1Max, x, y, z);             852       SetCorner(sC0Max1Max, x, y, z);
914                                                   853       
915       // corner of Axis0min and Axis1max          854       // corner of Axis0min and Axis1max
916       x = endRad[zmax]*std::cos(endPhi[zmax] -    855       x = endRad[zmax]*std::cos(endPhi[zmax] - halfdphi);
917       y = endRad[zmax]*std::sin(endPhi[zmax] -    856       y = endRad[zmax]*std::sin(endPhi[zmax] - halfdphi);
918       z = endZ[zmax];                             857       z = endZ[zmax];
919       SetCorner(sC0Min1Max, x, y, z);             858       SetCorner(sC0Min1Max, x, y, z);
920                                                   859 
921    }                                           << 860    } else {
922    else                                        << 
923    {                                           << 
924       std::ostringstream message;                 861       std::ostringstream message;
925       message << "Feature NOT implemented !" <    862       message << "Feature NOT implemented !" << G4endl
926               << "        fAxis[0] = " << fAxi    863               << "        fAxis[0] = " << fAxis[0] << G4endl
927               << "        fAxis[1] = " << fAxi    864               << "        fAxis[1] = " << fAxis[1];
928       G4Exception("G4TwistTubsHypeSide::SetCor    865       G4Exception("G4TwistTubsHypeSide::SetCorners()",
929                   "GeomSolids0001", FatalExcep    866                   "GeomSolids0001", FatalException, message);
930    }                                              867    }
931 }                                                 868 }
932                                                   869 
                                                   >> 870 
933 //============================================    871 //=====================================================================
934 //* SetCorners() -----------------------------    872 //* SetCorners() ------------------------------------------------------
935                                                   873 
936 void G4TwistTubsHypeSide::SetCorners()            874 void G4TwistTubsHypeSide::SetCorners()
937 {                                                 875 {
938    G4Exception("G4TwistTubsHypeSide::SetCorner    876    G4Exception("G4TwistTubsHypeSide::SetCorners()",
939                "GeomSolids0001", FatalExceptio    877                "GeomSolids0001", FatalException,
940                "Method NOT implemented !");       878                "Method NOT implemented !");
941 }                                                 879 }
942                                                   880 
943 //============================================    881 //=====================================================================
944 //* SetBoundaries() --------------------------    882 //* SetBoundaries() ---------------------------------------------------
945                                                   883 
946 void G4TwistTubsHypeSide::SetBoundaries()         884 void G4TwistTubsHypeSide::SetBoundaries()
947 {                                                 885 {
948    // Set direction-unit vector of phi-boundar    886    // Set direction-unit vector of phi-boundary-lines in local coodinate.
949    // sAxis0 must be kPhi.                        887    // sAxis0 must be kPhi.
950    // This fanction set lower phi-boundary and    888    // This fanction set lower phi-boundary and upper phi-boundary.
951                                                   889 
952    if (fAxis[0] == kPhi && fAxis[1] == kZAxis) << 890    if (fAxis[0] == kPhi && fAxis[1] == kZAxis) {
953    {                                           << 891 
954       G4ThreeVector direction;                    892       G4ThreeVector direction;
955       // sAxis0 & sAxisMin                        893       // sAxis0 & sAxisMin
956       direction = GetCorner(sC0Min1Max) - GetC    894       direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
957       direction = direction.unit();               895       direction = direction.unit();
958       SetBoundary(sAxis0 & (sAxisPhi | sAxisMi    896       SetBoundary(sAxis0 & (sAxisPhi | sAxisMin), direction, 
959                    GetCorner(sC0Min1Min), sAxi    897                    GetCorner(sC0Min1Min), sAxisZ);
960                                                   898 
961       // sAxis0 & sAxisMax                        899       // sAxis0 & sAxisMax
962       direction = GetCorner(sC0Max1Max) - GetC    900       direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
963       direction = direction.unit();               901       direction = direction.unit();
964       SetBoundary(sAxis0 & (sAxisPhi | sAxisMa    902       SetBoundary(sAxis0 & (sAxisPhi | sAxisMax), direction, 
965                   GetCorner(sC0Max1Min), sAxis    903                   GetCorner(sC0Max1Min), sAxisZ);
966                                                   904 
967       // sAxis1 & sAxisMin                        905       // sAxis1 & sAxisMin
968       direction = GetCorner(sC0Max1Min) - GetC    906       direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
969       direction = direction.unit();               907       direction = direction.unit();
970       SetBoundary(sAxis1 & (sAxisZ | sAxisMin)    908       SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction, 
971                    GetCorner(sC0Min1Min), sAxi    909                    GetCorner(sC0Min1Min), sAxisPhi);
972                                                   910 
973       // sAxis1 & sAxisMax                        911       // sAxis1 & sAxisMax
974       direction = GetCorner(sC0Max1Max) - GetC    912       direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
975       direction = direction.unit();               913       direction = direction.unit();
976       SetBoundary(sAxis1 & (sAxisZ | sAxisMax)    914       SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction, 
977                   GetCorner(sC0Min1Max), sAxis    915                   GetCorner(sC0Min1Max), sAxisPhi);
978    }                                           << 916    } else {
979    else                                        << 
980    {                                           << 
981       std::ostringstream message;                 917       std::ostringstream message;
982       message << "Feature NOT implemented !" <    918       message << "Feature NOT implemented !" << G4endl
983               << "        fAxis[0] = " << fAxi    919               << "        fAxis[0] = " << fAxis[0] << G4endl
984               << "        fAxis[1] = " << fAxi    920               << "        fAxis[1] = " << fAxis[1];
985       G4Exception("G4TwistTubsHypeSide::SetBou    921       G4Exception("G4TwistTubsHypeSide::SetBoundaries()",
986                   "GeomSolids0001", FatalExcep    922                   "GeomSolids0001", FatalException, message);
987    }                                              923    }
988 }                                                 924 }
989                                                   925 
990 //============================================    926 //=====================================================================
991 //* GetFacets() ------------------------------    927 //* GetFacets() -------------------------------------------------------
992                                                   928 
993 void G4TwistTubsHypeSide::GetFacets( G4int k,     929 void G4TwistTubsHypeSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
994                                      G4int fac    930                                      G4int faces[][4], G4int iside ) 
995 {                                                 931 {
                                                   >> 932 
996   G4double z ;     // the two parameters for t    933   G4double z ;     // the two parameters for the surface equation
997   G4double x,xmin,xmax ;                          934   G4double x,xmin,xmax ;
998                                                   935 
999   G4ThreeVector p ;  // a point on the surface    936   G4ThreeVector p ;  // a point on the surface, given by (z,u)
1000                                                  937 
1001   G4int nnode ;                                  938   G4int nnode ;
1002   G4int nface ;                                  939   G4int nface ;
1003                                                  940 
1004   // calculate the (n-1)*(k-1) vertices          941   // calculate the (n-1)*(k-1) vertices
1005                                                  942 
1006   for ( G4int i = 0 ; i<n ; ++i )             << 943   G4int i,j ;
1007   {                                           << 944 
                                                   >> 945   for ( i = 0 ; i<n ; i++ ) {
                                                   >> 946 
1008     z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin    947     z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;
1009                                                  948 
1010     for ( G4int j = 0 ; j<k ; ++j )           << 949     for ( j = 0 ; j<k ; j++ )
1011     {                                            950     {
1012       nnode = GetNode(i,j,k,n,iside) ;           951       nnode = GetNode(i,j,k,n,iside) ;
1013                                                  952 
1014       xmin = GetBoundaryMin(z) ;                 953       xmin = GetBoundaryMin(z) ; 
1015       xmax = GetBoundaryMax(z) ;                 954       xmax = GetBoundaryMax(z) ;
1016                                                  955 
1017       if (fHandedness < 0)  // inner hyperbol << 956       if (fHandedness < 0) { // inner hyperbolic surface
1018       {                                       << 
1019         x = xmin + j*(xmax-xmin)/(k-1) ;         957         x = xmin + j*(xmax-xmin)/(k-1) ;
1020       }                                       << 958       } else {               // outer hyperbolic surface
1021       else                  // outer hyperbol << 
1022       {                                       << 
1023         x = xmax - j*(xmax-xmin)/(k-1) ;         959         x = xmax - j*(xmax-xmin)/(k-1) ;
1024       }                                          960       }
1025                                                  961 
1026       p = SurfacePoint(x,z,true) ;  // surfac    962       p = SurfacePoint(x,z,true) ;  // surface point in global coord.system
1027                                                  963 
1028       xyz[nnode][0] = p.x() ;                    964       xyz[nnode][0] = p.x() ;
1029       xyz[nnode][1] = p.y() ;                    965       xyz[nnode][1] = p.y() ;
1030       xyz[nnode][2] = p.z() ;                    966       xyz[nnode][2] = p.z() ;
1031                                                  967 
1032       if ( i<n-1 && j<k-1 )    // clock wise  << 968       if ( i<n-1 && j<k-1 ) {   // clock wise filling
1033       {                                       << 969         
1034         nface = GetFace(i,j,k,n,iside) ;         970         nface = GetFace(i,j,k,n,iside) ;
1035   faces[nface][0] = GetEdgeVisibility(i,j,k,n << 971 
1036                         * ( GetNode(i  ,j  ,k << 972   faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i  ,j  ,k,n,iside)+1) ;  
1037   faces[nface][1] = GetEdgeVisibility(i,j,k,n << 973   faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j  ,k,n,iside)+1) ;
1038                         * ( GetNode(i+1,j  ,k << 974   faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
1039   faces[nface][2] = GetEdgeVisibility(i,j,k,n << 975   faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i  ,j+1,k,n,iside)+1) ;
1040                         * ( GetNode(i+1,j+1,k << 976 
1041   faces[nface][3] = GetEdgeVisibility(i,j,k,n << 
1042                         * ( GetNode(i  ,j+1,k << 
1043       }                                          977       }
1044     }                                            978     }
1045   }                                              979   }
1046 }                                                980 }
1047                                                  981