Geant4 Cross Reference

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


  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 // G4VTwistSurface implementation              << 
 27 //                                                 26 //
 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu <<  27 // $Id: G4VTwistSurface.cc,v 1.10 2010-07-12 15:25:37 gcosmo Exp $
 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), <<  28 // GEANT4 tag $Name: not supported by cvs2svn $
 30 //               from original version in Jupi <<  29 //
                                                   >>  30 // 
                                                   >>  31 // --------------------------------------------------------------------
                                                   >>  32 // GEANT 4 class source file
                                                   >>  33 //
                                                   >>  34 //
                                                   >>  35 // G4VTwistSurface.cc
                                                   >>  36 //
                                                   >>  37 // Author: 
                                                   >>  38 //   01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp)
                                                   >>  39 //
                                                   >>  40 // History:
                                                   >>  41 //   13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
                                                   >>  42 //                 from original version in Jupiter-2.5.02 application.
 31 // -------------------------------------------     43 // --------------------------------------------------------------------
 32                                                    44 
 33 #include <iomanip>                                 45 #include <iomanip>
 34                                                    46 
 35 #include "G4VTwistSurface.hh"                      47 #include "G4VTwistSurface.hh"
 36 #include "G4GeometryTolerance.hh"                  48 #include "G4GeometryTolerance.hh"
 37                                                    49 
 38 const G4int  G4VTwistSurface::sOutside             50 const G4int  G4VTwistSurface::sOutside        = 0x00000000;
 39 const G4int  G4VTwistSurface::sInside              51 const G4int  G4VTwistSurface::sInside         = 0x10000000;
 40 const G4int  G4VTwistSurface::sBoundary            52 const G4int  G4VTwistSurface::sBoundary       = 0x20000000;
 41 const G4int  G4VTwistSurface::sCorner              53 const G4int  G4VTwistSurface::sCorner         = 0x40000000;
 42 const G4int  G4VTwistSurface::sC0Min1Min           54 const G4int  G4VTwistSurface::sC0Min1Min      = 0x40000101; 
 43 const G4int  G4VTwistSurface::sC0Max1Min           55 const G4int  G4VTwistSurface::sC0Max1Min      = 0x40000201;
 44 const G4int  G4VTwistSurface::sC0Max1Max           56 const G4int  G4VTwistSurface::sC0Max1Max      = 0x40000202; 
 45 const G4int  G4VTwistSurface::sC0Min1Max           57 const G4int  G4VTwistSurface::sC0Min1Max      = 0x40000102; 
 46 const G4int  G4VTwistSurface::sAxisMin             58 const G4int  G4VTwistSurface::sAxisMin        = 0x00000101; 
 47 const G4int  G4VTwistSurface::sAxisMax             59 const G4int  G4VTwistSurface::sAxisMax        = 0x00000202; 
 48 const G4int  G4VTwistSurface::sAxisX               60 const G4int  G4VTwistSurface::sAxisX          = 0x00000404;
 49 const G4int  G4VTwistSurface::sAxisY               61 const G4int  G4VTwistSurface::sAxisY          = 0x00000808;
 50 const G4int  G4VTwistSurface::sAxisZ               62 const G4int  G4VTwistSurface::sAxisZ          = 0x00000C0C;
 51 const G4int  G4VTwistSurface::sAxisRho             63 const G4int  G4VTwistSurface::sAxisRho        = 0x00001010;
 52 const G4int  G4VTwistSurface::sAxisPhi             64 const G4int  G4VTwistSurface::sAxisPhi        = 0x00001414;
 53                                                    65 
 54 // mask                                            66 // mask
 55 const G4int  G4VTwistSurface::sAxis0               67 const G4int  G4VTwistSurface::sAxis0          = 0x0000FF00;
 56 const G4int  G4VTwistSurface::sAxis1               68 const G4int  G4VTwistSurface::sAxis1          = 0x000000FF;
 57 const G4int  G4VTwistSurface::sSizeMask            69 const G4int  G4VTwistSurface::sSizeMask       = 0x00000303;
 58 const G4int  G4VTwistSurface::sAxisMask            70 const G4int  G4VTwistSurface::sAxisMask       = 0x0000FCFC;
 59 const G4int  G4VTwistSurface::sAreaMask            71 const G4int  G4VTwistSurface::sAreaMask       = 0XF0000000;
 60                                                    72 
 61 //============================================     73 //=====================================================================
 62 //* constructors -----------------------------     74 //* constructors ------------------------------------------------------
 63                                                    75 
 64 G4VTwistSurface::G4VTwistSurface(const G4Strin     76 G4VTwistSurface::G4VTwistSurface(const G4String &name)
 65   : fIsValidNorm(false), fName(name)               77   : fIsValidNorm(false), fName(name)
 66 {                                                  78 {
 67                                                    79 
 68    fAxis[0]    = kUndefined;                       80    fAxis[0]    = kUndefined;
 69    fAxis[1]    = kUndefined;                       81    fAxis[1]    = kUndefined;
 70    fAxisMin[0] = kInfinity;                        82    fAxisMin[0] = kInfinity;
 71    fAxisMin[1] = kInfinity;                        83    fAxisMin[1] = kInfinity;
 72    fAxisMax[0] = kInfinity;                        84    fAxisMax[0] = kInfinity;
 73    fAxisMax[1] = kInfinity;                        85    fAxisMax[1] = kInfinity;
 74    fHandedness = 1;                                86    fHandedness = 1;
 75                                                    87 
 76    for (auto i=0; i<4; ++i)                    <<  88    for (G4int i=0; i<4; i++)
 77    {                                               89    {
 78       fCorners[i].set(kInfinity, kInfinity, kI     90       fCorners[i].set(kInfinity, kInfinity, kInfinity);
 79       fNeighbours[i] = nullptr;                <<  91       fNeighbours[i] = 0;
 80    }                                               92    }
 81                                                    93 
 82    fCurrentNormal.p.set(kInfinity, kInfinity,      94    fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
 83                                                    95    
 84    fAmIOnLeftSide.me.set(kInfinity, kInfinity,     96    fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
 85    fAmIOnLeftSide.vec.set(kInfinity, kInfinity     97    fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
 86    kCarTolerance = G4GeometryTolerance::GetIns     98    kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
 87 }                                                  99 }
 88                                                   100 
 89 G4VTwistSurface::G4VTwistSurface(const G4Strin << 101 G4VTwistSurface::G4VTwistSurface(const G4String         &name,
 90                        const G4RotationMatrix& << 102                        const G4RotationMatrix &rot,
 91                        const G4ThreeVector&    << 103                        const G4ThreeVector    &tlate,
 92                              G4int                104                              G4int             handedness,
 93                        const EAxis                105                        const EAxis             axis0 ,
 94                        const EAxis                106                        const EAxis             axis1 ,
 95                              G4double             107                              G4double          axis0min,
 96                              G4double             108                              G4double          axis1min,
 97                              G4double             109                              G4double          axis0max,
 98                              G4double             110                              G4double          axis1max )
 99    : fIsValidNorm(false), fName(name)             111    : fIsValidNorm(false), fName(name)
100 {                                                 112 {
101    fAxis[0]    = axis0;                           113    fAxis[0]    = axis0;
102    fAxis[1]    = axis1;                           114    fAxis[1]    = axis1;
103    fAxisMin[0] = axis0min;                        115    fAxisMin[0] = axis0min;
104    fAxisMin[1] = axis1min;                        116    fAxisMin[1] = axis1min;
105    fAxisMax[0] = axis0max;                        117    fAxisMax[0] = axis0max;
106    fAxisMax[1] = axis1max;                        118    fAxisMax[1] = axis1max;
107    fHandedness = handedness;                      119    fHandedness = handedness;
108    fRot        = rot;                             120    fRot        = rot;
109    fTrans      = tlate;                           121    fTrans      = tlate;
110                                                   122 
111    for (auto i=0; i<4; ++i)                    << 123    for (G4int i=0; i<4; i++)
112    {                                              124    {
113       fCorners[i].set(kInfinity, kInfinity, kI    125       fCorners[i].set(kInfinity, kInfinity, kInfinity);
114       fNeighbours[i] = nullptr;                << 126       fNeighbours[i] = 0;
115    }                                              127    }
116                                                   128 
117    fCurrentNormal.p.set(kInfinity, kInfinity,     129    fCurrentNormal.p.set(kInfinity, kInfinity, kInfinity);
118                                                   130    
119    fAmIOnLeftSide.me.set(kInfinity, kInfinity,    131    fAmIOnLeftSide.me.set(kInfinity, kInfinity, kInfinity);
120    fAmIOnLeftSide.vec.set(kInfinity, kInfinity    132    fAmIOnLeftSide.vec.set(kInfinity, kInfinity, kInfinity);
121    kCarTolerance = G4GeometryTolerance::GetIns    133    kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
122 }                                                 134 }
123                                                   135 
124 //============================================    136 //=====================================================================
125 //* Fake default constructor -----------------    137 //* Fake default constructor ------------------------------------------
126                                                   138 
127 G4VTwistSurface::G4VTwistSurface( __void__& )     139 G4VTwistSurface::G4VTwistSurface( __void__& )
128   : fHandedness(0), fIsValidNorm(false), kCarT    140   : fHandedness(0), fIsValidNorm(false), kCarTolerance(0.),
129     fName("")                                     141     fName("")
130 {                                                 142 {
131    fAxis[0] = fAxis[1] = kXAxis;                  143    fAxis[0] = fAxis[1] = kXAxis;
132    fAxisMin[0] = fAxisMin[1] = 0.;                144    fAxisMin[0] = fAxisMin[1] = 0.;
133    fAxisMax[0] = fAxisMax[1] = 0.;                145    fAxisMax[0] = fAxisMax[1] = 0.;
134    fNeighbours[0] = fNeighbours[1] = fNeighbou << 146    fNeighbours[0] = fNeighbours[1] = fNeighbours[2] = fNeighbours[3] = 0;
                                                   >> 147 }
                                                   >> 148 
                                                   >> 149 //=====================================================================
                                                   >> 150 //* destructor --------------------------------------------------------
                                                   >> 151 
                                                   >> 152 G4VTwistSurface::~G4VTwistSurface()
                                                   >> 153 {
135 }                                                 154 }
136                                                   155 
137 //============================================    156 //=====================================================================
138 //* AmIOnLeftSide ----------------------------    157 //* AmIOnLeftSide -----------------------------------------------------
139                                                   158 
140 G4int G4VTwistSurface::AmIOnLeftSide(const G4T << 159 G4int G4VTwistSurface::AmIOnLeftSide(const G4ThreeVector &me, 
141                                      const G4T << 160                                 const G4ThreeVector &vec,
142                                            G4b << 161                                       G4bool        withtol) 
143 {                                                 162 {
144    // AmIOnLeftSide returns phi-location of "m    163    // AmIOnLeftSide returns phi-location of "me"
145    // (phi relation between me and vec project    164    // (phi relation between me and vec projected on z=0 plane).
146    // If "me" is on -ve-phi-side of "vec", it     165    // If "me" is on -ve-phi-side of "vec", it returns 1.
147    // On the other hand, if "me" is on +ve-phi    166    // On the other hand, if "me" is on +ve-phi-side of "vec",
148    // it returns -1.                              167    // it returns -1.
149    // (The return value represents z-coordinat    168    // (The return value represents z-coordinate of normal vector
150    //  of me.cross(vec).)                         169    //  of me.cross(vec).)
151    // If me is on boundary of vec, return 0.      170    // If me is on boundary of vec, return 0.
152                                                   171 
153    const G4double kAngTolerance                << 172    static const G4double kAngTolerance
154      = G4GeometryTolerance::GetInstance()->Get    173      = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
155                                                   174 
156    G4RotationMatrix unitrot;                   << 175    static G4RotationMatrix unitrot;  // unit matrix
157    const G4RotationMatrix rottol    = unitrot. << 176    static const G4RotationMatrix rottol    = unitrot.rotateZ(0.5*kAngTolerance);
158    const G4RotationMatrix invrottol = unitrot. << 177    static const G4RotationMatrix invrottol = unitrot.rotateZ(-1.*kAngTolerance);
159                                                   178 
160    if (fAmIOnLeftSide.me == me                    179    if (fAmIOnLeftSide.me == me 
161        && fAmIOnLeftSide.vec == vec               180        && fAmIOnLeftSide.vec == vec
162        && fAmIOnLeftSide.withTol == withtol)   << 181        && fAmIOnLeftSide.withTol == withtol) {
163    {                                           << 
164       return fAmIOnLeftSide.amIOnLeftSide;        182       return fAmIOnLeftSide.amIOnLeftSide;
165    }                                              183    }
166                                                   184    
167    fAmIOnLeftSide.me      = me;                   185    fAmIOnLeftSide.me      = me;
168    fAmIOnLeftSide.vec     = vec;                  186    fAmIOnLeftSide.vec     = vec;
169    fAmIOnLeftSide.withTol = withtol;              187    fAmIOnLeftSide.withTol = withtol;
170                                                   188    
171    G4ThreeVector met   = (G4ThreeVector(me.x()    189    G4ThreeVector met   = (G4ThreeVector(me.x(), me.y(), 0.)).unit();
172    G4ThreeVector vect  = (G4ThreeVector(vec.x(    190    G4ThreeVector vect  = (G4ThreeVector(vec.x(), vec.y(), 0.)).unit();
173                                                   191    
174    G4ThreeVector ivect = invrottol * vect;        192    G4ThreeVector ivect = invrottol * vect;
175    G4ThreeVector rvect = rottol * vect;           193    G4ThreeVector rvect = rottol * vect;
176                                                   194 
177    G4double metcrossvect = met.x() * vect.y()     195    G4double metcrossvect = met.x() * vect.y() - met.y() * vect.x();
178                                                   196    
179    if (withtol)                                << 197    if (withtol) {
180    {                                           << 
181       if (met.x() * ivect.y() - met.y() * ivec    198       if (met.x() * ivect.y() - met.y() * ivect.x() > 0 && 
182           metcrossvect >= 0)  {                   199           metcrossvect >= 0)  {
183          fAmIOnLeftSide.amIOnLeftSide = 1;        200          fAmIOnLeftSide.amIOnLeftSide = 1;
184       } else if (met.x() * rvect.y() - met.y()    201       } else if (met.x() * rvect.y() - met.y() * rvect.x() < 0 &&
185                  metcrossvect <= 0)  {            202                  metcrossvect <= 0)  {
186          fAmIOnLeftSide.amIOnLeftSide = -1;       203          fAmIOnLeftSide.amIOnLeftSide = -1;
187       } else {                                    204       } else {
188          fAmIOnLeftSide.amIOnLeftSide = 0;        205          fAmIOnLeftSide.amIOnLeftSide = 0;
189       }                                           206       }
190    }                                           << 207    } else {
191    else                                        << 
192    {                                           << 
193       if (metcrossvect > 0) {                     208       if (metcrossvect > 0) {    
194          fAmIOnLeftSide.amIOnLeftSide = 1;        209          fAmIOnLeftSide.amIOnLeftSide = 1;
195       } else if (metcrossvect < 0 ) {             210       } else if (metcrossvect < 0 ) {
196          fAmIOnLeftSide.amIOnLeftSide = -1;       211          fAmIOnLeftSide.amIOnLeftSide = -1;
197       } else {                                    212       } else {       
198          fAmIOnLeftSide.amIOnLeftSide = 0;        213          fAmIOnLeftSide.amIOnLeftSide = 0;
199       }                                           214       }
200    }                                              215    }
201                                                   216 
202 #ifdef G4TWISTDEBUG                               217 #ifdef G4TWISTDEBUG
203    G4cout << "         === G4VTwistSurface::Am    218    G4cout << "         === G4VTwistSurface::AmIOnLeftSide() =============="
204           << G4endl;                              219           << G4endl;
205    G4cout << "             Name , returncode      220    G4cout << "             Name , returncode  : " << fName << " " 
206                        << fAmIOnLeftSide.amIOn    221                        << fAmIOnLeftSide.amIOnLeftSide <<  G4endl;
207    G4cout << "             me, vec    : " << s    222    G4cout << "             me, vec    : " << std::setprecision(14) << me 
208                                           << "    223                                           << " " << vec  << G4endl;
209    G4cout << "             met, vect  : " << m    224    G4cout << "             met, vect  : " << met << " " << vect  << G4endl;
210    G4cout << "             ivec, rvec : " << i    225    G4cout << "             ivec, rvec : " << ivect << " " << rvect << G4endl;
211    G4cout << "             met x vect : " << m    226    G4cout << "             met x vect : " << metcrossvect << G4endl;
212    G4cout << "             met x ivec : " << m    227    G4cout << "             met x ivec : " << met.cross(ivect) << G4endl;
213    G4cout << "             met x rvec : " << m    228    G4cout << "             met x rvec : " << met.cross(rvect) << G4endl;
214    G4cout << "         =======================    229    G4cout << "         =============================================="
215           << G4endl;                              230           << G4endl;
216 #endif                                            231 #endif
217                                                   232 
218    return fAmIOnLeftSide.amIOnLeftSide;           233    return fAmIOnLeftSide.amIOnLeftSide;
219 }                                                 234 }
220                                                   235 
221 //============================================    236 //=====================================================================
222 //* DistanceToBoundary -----------------------    237 //* DistanceToBoundary ------------------------------------------------
223                                                   238 
224 G4double G4VTwistSurface::DistanceToBoundary(G    239 G4double G4VTwistSurface::DistanceToBoundary(G4int areacode,
225                                              G << 240                                         G4ThreeVector &xx,
226                                        const G << 241                                         const G4ThreeVector &p) 
227 {                                                 242 {
228    // DistanceToBoundary                          243    // DistanceToBoundary 
229    //                                             244    //
230    // return distance to nearest boundary from    245    // return distance to nearest boundary from arbitrary point p 
231    // in local coodinate.                         246    // in local coodinate.
232    // Argument areacode must be one of them:      247    // Argument areacode must be one of them:
233    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,       248    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
234    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.       249    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
235    //                                             250    //
236                                                   251 
237    G4ThreeVector d;    // direction vector of     252    G4ThreeVector d;    // direction vector of the boundary
238    G4ThreeVector x0;   // reference point of t    253    G4ThreeVector x0;   // reference point of the boundary
239    G4double      dist = kInfinity;                254    G4double      dist = kInfinity;
240    G4int         boundarytype;                    255    G4int         boundarytype;
241                                                   256 
242    if (IsAxis0(areacode) && IsAxis1(areacode)) << 257    if (IsAxis0(areacode) && IsAxis1(areacode)) {
243    {                                           << 
244       std::ostringstream message;                 258       std::ostringstream message;
245       message << "Point is in the corner area.    259       message << "Point is in the corner area." << G4endl
246               << "        Point is in the corn    260               << "        Point is in the corner area. This function returns"
247               << G4endl                           261               << G4endl
248               << "        a direction vector o    262               << "        a direction vector of a boundary line." << G4endl
249               << "        areacode = " << area    263               << "        areacode = " << areacode;
250       G4Exception("G4VTwistSurface::DistanceTo    264       G4Exception("G4VTwistSurface::DistanceToBoundary()", "GeomSolids0003",
251                   FatalException, message);       265                   FatalException, message);
252    }                                           << 266    } else if (IsAxis0(areacode) || IsAxis1(areacode)) {
253    else if (IsAxis0(areacode) || IsAxis1(areac << 
254    {                                           << 
255       GetBoundaryParameters(areacode, d, x0, b    267       GetBoundaryParameters(areacode, d, x0, boundarytype);
256       if (boundarytype == sAxisPhi)            << 268       if (boundarytype == sAxisPhi) {
257       {                                        << 
258          G4double t = x0.getRho() / p.getRho()    269          G4double t = x0.getRho() / p.getRho();
259          xx.set(t*p.x(), t*p.y(), x0.z());        270          xx.set(t*p.x(), t*p.y(), x0.z());
260          dist = (xx - p).mag();                   271          dist = (xx - p).mag();
261       }                                        << 272       } else { 
262       else                                     << 
263       {                                        << 
264          // linear boundary                       273          // linear boundary
265          // sAxisX, sAxisY, sAxisZ, sAxisRho      274          // sAxisX, sAxisY, sAxisZ, sAxisRho
266          dist = DistanceToLine(p, x0, d, xx);     275          dist = DistanceToLine(p, x0, d, xx);
267       }                                           276       }
268    }                                           << 277    } else {
269    else                                        << 
270    {                                           << 
271       std::ostringstream message;                 278       std::ostringstream message;
272       message << "Bad areacode of boundary." <    279       message << "Bad areacode of boundary." << G4endl
273               << "        areacode = " << area    280               << "        areacode = " << areacode;
274       G4Exception("G4VTwistSurface::DistanceTo    281       G4Exception("G4VTwistSurface::DistanceToBoundary()", "GeomSolids0003",
275                   FatalException, message);       282                   FatalException, message);
276    }                                              283    }
277    return dist;                                   284    return dist;
278 }                                                 285 }
279                                                   286 
280 //============================================    287 //=====================================================================
281 //* DistanceToIn -----------------------------    288 //* DistanceToIn ------------------------------------------------------
282                                                   289 
283 G4double G4VTwistSurface::DistanceToIn(const G << 290 G4double G4VTwistSurface::DistanceToIn(const G4ThreeVector &gp,
284                                        const G << 291                                   const G4ThreeVector &gv,
285                                              G << 292                                         G4ThreeVector &gxxbest)
286 {                                                 293 {
287 #ifdef G4TWISTDEBUG                               294 #ifdef G4TWISTDEBUG
288    G4cout << " ~~~~ G4VTwistSurface::DistanceT << 295    G4cout << " ~~~~~ G4VTwistSurface::DistanceToIn(p,v) - Start ~~~~~" << G4endl;
289    G4cout << "      Name : " << fName << G4end    296    G4cout << "      Name : " << fName << G4endl;
290    G4cout << "      gp   : " << gp << G4endl;     297    G4cout << "      gp   : " << gp << G4endl;
291    G4cout << "      gv   : " <<  gv << G4endl;    298    G4cout << "      gv   : " <<  gv << G4endl;
292    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    299    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
293 #endif                                            300 #endif
294                                                   301    
295    G4ThreeVector gxx[G4VSURFACENXX];              302    G4ThreeVector gxx[G4VSURFACENXX];
296    G4double      distance[G4VSURFACENXX]  ;       303    G4double      distance[G4VSURFACENXX]  ; 
297    G4int         areacode[G4VSURFACENXX]  ;       304    G4int         areacode[G4VSURFACENXX]  ;
298    G4bool        isvalid[G4VSURFACENXX]   ;       305    G4bool        isvalid[G4VSURFACENXX]   ; 
299                                                   306    
300    for (G4int i = 0 ; i<G4VSURFACENXX ; ++i )  << 307    for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
301    {                                           << 
302      distance[i] = kInfinity ;                    308      distance[i] = kInfinity ;
303      areacode[i] = sOutside ;                     309      areacode[i] = sOutside ;
304      isvalid[i] = false ;                         310      isvalid[i] = false ;
305    }                                              311    }
306                                                   312 
307    G4double      bestdistance   = kInfinity;      313    G4double      bestdistance   = kInfinity;
308 #ifdef G4TWISTDEBUG                               314 #ifdef G4TWISTDEBUG
309    G4int         besti          = -1;             315    G4int         besti          = -1;  
310 #endif                                            316 #endif
311    G4ThreeVector bestgxx(kInfinity, kInfinity,    317    G4ThreeVector bestgxx(kInfinity, kInfinity, kInfinity);
312                                                   318 
313    G4int         nxx = DistanceToSurface(gp, g    319    G4int         nxx = DistanceToSurface(gp, gv, gxx, distance, areacode, 
314                                          isval    320                                          isvalid, kValidateWithTol);
315                                                   321 
316    for (G4int i=0; i<nxx; ++i)                 << 322    for (G4int i=0; i< nxx; i++) {
317    {                                           << 
318                                                   323 
319       // skip this intersection if:               324       // skip this intersection if:
320       //   - invalid intersection                 325       //   - invalid intersection
321       //   - particle goes outword the surface    326       //   - particle goes outword the surface
322                                                   327 
323       if (!isvalid[i])                         << 328       if (!isvalid[i]) {
324       {                                        << 
325          // xx[i] is sOutside or distance[i] <    329          // xx[i] is sOutside or distance[i] < 0
326          continue;                                330          continue;      
327       }                                           331       }
328                                                   332 
329       G4ThreeVector normal = GetNormal(gxx[i],    333       G4ThreeVector normal = GetNormal(gxx[i], true);
330                                                   334 
331       if ((normal * gv) >= 0)                  << 335       if ((normal * gv) >= 0) {
332       {                                        << 
333                                                   336 
334 #ifdef G4TWISTDEBUG                               337 #ifdef G4TWISTDEBUG
335          G4cout << "   G4VTwistSurface::Distan    338          G4cout << "   G4VTwistSurface::DistanceToIn(p,v): "
336                 << "particle goes outword the     339                 << "particle goes outword the surface." << G4endl;
337 #endif                                            340 #endif 
338          continue;                                341          continue; 
339       }                                           342       }
340                                                   343       
341       //                                          344       //
342       // accept this intersection if the inter    345       // accept this intersection if the intersection is inside.
343       //                                          346       //
344                                                   347 
345       if (IsInside(areacode[i]))               << 348       if (IsInside(areacode[i])) {
346       {                                        << 349          if (distance[i] < bestdistance) {
347          if (distance[i] < bestdistance)       << 
348          {                                     << 
349             bestdistance = distance[i];           350             bestdistance = distance[i];
350             bestgxx = gxx[i];                     351             bestgxx = gxx[i];
351 #ifdef G4TWISTDEBUG                               352 #ifdef G4TWISTDEBUG
352             besti   = i;                          353             besti   = i;
353             G4cout << "   G4VTwistSurface::Dis    354             G4cout << "   G4VTwistSurface::DistanceToIn(p,v): "
354                    << " areacode sInside name,    355                    << " areacode sInside name, distance = "
355                    << fName <<  " "<< bestdist    356                    << fName <<  " "<< bestdistance << G4endl;
356 #endif                                            357 #endif 
357          }                                        358          }
358                                                   359 
359       //                                          360       //
360       // else, the intersection is on boundary    361       // else, the intersection is on boundary or corner.
361       //                                          362       //
362                                                   363 
363       }                                        << 364       } else {
364       else                                     << 365 
365       {                                        << 366          G4VTwistSurface *neighbours[2];
366          G4VTwistSurface* neighbours[2];       << 
367          G4bool      isaccepted[2] = {false, f    367          G4bool      isaccepted[2] = {false, false};
368          G4int       nneighbours   = GetNeighb    368          G4int       nneighbours   = GetNeighbours(areacode[i], neighbours);
369                                                   369             
370          for (G4int j=0; j<nneighbours; ++j)   << 370          for (G4int j=0; j< nneighbours; j++) {
371          {                                     << 
372             // if on corner, nneighbours = 2.     371             // if on corner, nneighbours = 2.
373             // if on boundary, nneighbours = 1    372             // if on boundary, nneighbours = 1.
374                                                   373 
375             G4ThreeVector tmpgxx[G4VSURFACENXX    374             G4ThreeVector tmpgxx[G4VSURFACENXX];
376             G4double      tmpdist[G4VSURFACENX    375             G4double      tmpdist[G4VSURFACENXX] ;
377             G4int         tmpareacode[G4VSURFA    376             G4int         tmpareacode[G4VSURFACENXX] ;
378             G4bool        tmpisvalid[G4VSURFAC    377             G4bool        tmpisvalid[G4VSURFACENXX] ;
379                                                   378 
380             for (G4int l = 0 ; l<G4VSURFACENXX << 379             for (G4int l = 0 ; l<G4VSURFACENXX ; l++ ) {
381             {                                  << 
382               tmpdist[l] = kInfinity ;            380               tmpdist[l] = kInfinity ;
383               tmpareacode[l] = sOutside ;         381               tmpareacode[l] = sOutside ;
384               tmpisvalid[l] = false ;             382               tmpisvalid[l] = false ;
385             }                                     383             }
386                                                   384 
387             G4int tmpnxx = neighbours[j]->Dist    385             G4int tmpnxx = neighbours[j]->DistanceToSurface(
388                                           gp,     386                                           gp, gv, tmpgxx, tmpdist,
389                                           tmpa    387                                           tmpareacode, tmpisvalid,
390                                           kVal    388                                           kValidateWithTol);
391             G4ThreeVector neighbournormal;        389             G4ThreeVector neighbournormal;
392                                                   390 
393             for (G4int k=0; k< tmpnxx; ++k)    << 391             for (G4int k=0; k< tmpnxx; k++) {
394             {                                  << 392 
395                //                                 393                //  
396                // if tmpxx[k] is valid && sIns    394                // if tmpxx[k] is valid && sInside, the final winner must
397                // be neighbour surface. return    395                // be neighbour surface. return kInfinity. 
398                // else , choose tmpxx on same     396                // else , choose tmpxx on same boundary of xx, then check normal 
399                //                                 397                //  
400                                                   398 
401                if (IsInside(tmpareacode[k]))   << 399                if (IsInside(tmpareacode[k])) {
402                {                               << 400 
403 #ifdef G4TWISTDEBUG                               401 #ifdef G4TWISTDEBUG
404                   G4cout << "   G4VTwistSurfac    402                   G4cout << "   G4VTwistSurface:DistanceToIn(p,v): "
405                          << " intersection "<<    403                          << " intersection "<< tmpgxx[k] << G4endl
406                          << "   is inside of n    404                          << "   is inside of neighbour surface of " << fName 
407                          << " . returning kInf    405                          << " . returning kInfinity." << G4endl;
408                   G4cout << "~~ G4VTwistSurfac << 406                   G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~"
409                          << G4endl;               407                          << G4endl;
410                   G4cout << "      No intersec    408                   G4cout << "      No intersections " << G4endl; 
411                   G4cout << "      Name : " <<    409                   G4cout << "      Name : " << fName << G4endl; 
412                   G4cout << "~~~~~~~~~~~~~~~~~    410                   G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
413                          << G4endl;               411                          << G4endl;
414 #endif                                            412 #endif 
415                   if (tmpisvalid[k])  return k    413                   if (tmpisvalid[k])  return kInfinity;
416                   continue;                       414                   continue;
417                                                   415 
418                //                                 416                //  
419                // if tmpxx[k] is valid && sIns    417                // if tmpxx[k] is valid && sInside, the final winner must
420                // be neighbour surface. return    418                // be neighbour surface. return .  
421                //                                 419                //
422                                                   420 
423                }                               << 421                } else if (IsSameBoundary(this,areacode[i],
424                else if (IsSameBoundary(this,ar << 422                                          neighbours[j], tmpareacode[k])) { 
425                                        neighbo << 
426                {                               << 
427                   // tmpxx[k] is same boundary    423                   // tmpxx[k] is same boundary (or corner) of xx.
428                                                   424                  
429                   neighbournormal = neighbours    425                   neighbournormal = neighbours[j]->GetNormal(tmpgxx[k], true);
430                   if (neighbournormal * gv < 0    426                   if (neighbournormal * gv < 0) isaccepted[j] = true;
431                }                                  427                }
432             }                                     428             } 
433                                                   429 
434             // if nneighbours = 1, chabge isac    430             // if nneighbours = 1, chabge isaccepted[1] before
435             // exiting neighboursurface loop.     431             // exiting neighboursurface loop.  
436                                                   432 
437             if (nneighbours == 1) isaccepted[1    433             if (nneighbours == 1) isaccepted[1] = true;
438                                                   434 
439          } // neighboursurface loop end           435          } // neighboursurface loop end
440                                                   436 
441          // now, we can accept xx intersection    437          // now, we can accept xx intersection
442                                                   438 
443          if (isaccepted[0] && isaccepted[1])   << 439          if (isaccepted[0] == true && isaccepted[1] == true) {
444          {                                     << 440             if (distance[i] < bestdistance) {
445             if (distance[i] < bestdistance)    << 
446             {                                  << 
447                bestdistance = distance[i];        441                bestdistance = distance[i];
448                gxxbest = gxx[i];                  442                gxxbest = gxx[i];
449 #ifdef G4TWISTDEBUG                               443 #ifdef G4TWISTDEBUG
450                besti = i;                      << 444                besti   = i;
451                G4cout << "   G4VTwistSurface::    445                G4cout << "   G4VTwistSurface::DistanceToIn(p,v): "
452                       << " areacode sBoundary     446                       << " areacode sBoundary & sBoundary distance = "
453                       << fName  << " " << dist    447                       << fName  << " " << distance[i] << G4endl;
454 #endif                                            448 #endif 
455             }                                     449             }
456          }                                        450          }
                                                   >> 451 
457       } // else end                               452       } // else end
458    } // intersection loop end                     453    } // intersection loop end
459                                                   454 
460    gxxbest = bestgxx;                             455    gxxbest = bestgxx;
461                                                   456 
462 #ifdef G4TWISTDEBUG                               457 #ifdef G4TWISTDEBUG
463    if (besti < 0)                              << 458    if (besti < 0) {
464    {                                           << 459       G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) - return ~~~~" << G4endl;
465       G4cout << "~~~ G4VTwistSurface::Distance << 
466       G4cout << "      No intersections " << G    460       G4cout << "      No intersections " << G4endl; 
467       G4cout << "      Name : " << fName << G4    461       G4cout << "      Name : " << fName << G4endl; 
468       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    462       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
469    }                                           << 463    } else {
470    else                                        << 464       G4cout << "~~~~~ G4VTwistSurface::DistanceToIn(p,v) : return ~~~~" << G4endl;
471    {                                           << 
472       G4cout << "~~~ G4VTwistSurface::Distance << 
473       G4cout << "      Name, i  : " << fName <    465       G4cout << "      Name, i  : " << fName << " , " << besti << G4endl; 
474       G4cout << "      gxx[i]   : " << gxxbest    466       G4cout << "      gxx[i]   : " << gxxbest << G4endl; 
475       G4cout << "      bestdist : " << bestdis    467       G4cout << "      bestdist : " << bestdistance << G4endl;
476       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    468       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
477    }                                              469    } 
478                                                   470 
479 #endif                                            471 #endif
480                                                   472 
481    return bestdistance;                           473    return bestdistance;
482 }                                                 474 }
483                                                   475 
484 //============================================    476 //=====================================================================
485 //* DistanceToOut(p, v) ----------------------    477 //* DistanceToOut(p, v) -----------------------------------------------
486                                                   478 
487 G4double G4VTwistSurface::DistanceToOut(const  << 479 G4double G4VTwistSurface::DistanceToOut(const G4ThreeVector &gp,
488                                         const  << 480                                    const G4ThreeVector &gv,
489                                                << 481                                          G4ThreeVector &gxxbest)
490 {                                                 482 {
491 #ifdef G4TWISTDEBUG                               483 #ifdef G4TWISTDEBUG
492    G4cout << "~~~~~ G4VTwistSurface::DistanceT    484    G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - Start ~~~~" << G4endl;
493    G4cout << "      Name : " << fName << G4end    485    G4cout << "      Name : " << fName << G4endl;
494    G4cout << "      gp   : " << gp << G4endl;     486    G4cout << "      gp   : " << gp << G4endl;
495    G4cout << "      gv   : " <<  gv << G4endl;    487    G4cout << "      gv   : " <<  gv << G4endl;
496    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    488    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
497 #endif                                            489 #endif
498                                                   490 
499    G4ThreeVector gxx[G4VSURFACENXX];              491    G4ThreeVector gxx[G4VSURFACENXX];
500    G4double      distance[G4VSURFACENXX];      << 492    G4double      distance[G4VSURFACENXX]  ; 
501    G4int         areacode[G4VSURFACENXX];      << 493    G4int         areacode[G4VSURFACENXX]  ;
502    G4bool        isvalid[G4VSURFACENXX];       << 494    G4bool        isvalid[G4VSURFACENXX]   ; 
                                                   >> 495    G4int         i;
503                                                   496    
504    for ( G4int i = 0 ; i<G4VSURFACENXX ; ++i ) << 497    for ( i = 0 ; i<G4VSURFACENXX ; i++ )
505    {                                              498    {
506      distance[i] = kInfinity ;                    499      distance[i] = kInfinity ;
507      areacode[i] = sOutside ;                     500      areacode[i] = sOutside ;
508      isvalid[i] = false ;                         501      isvalid[i] = false ;
509    }                                              502    }
510                                                   503 
511    G4int         nxx;                             504    G4int         nxx;
512    G4double      bestdistance   = kInfinity;      505    G4double      bestdistance   = kInfinity;
513                                                   506 
514    nxx = DistanceToSurface(gp, gv, gxx, distan    507    nxx = DistanceToSurface(gp, gv, gxx, distance, areacode,
515                            isvalid, kValidateW    508                            isvalid, kValidateWithTol);
516                                                   509 
517    for (G4int i=0; i<nxx; ++i)                 << 510    for (i=0; i<nxx; i++) {
518    {                                           << 511       if (!(isvalid[i])) {
519       if (!(isvalid[i]))                       << 
520       {                                        << 
521          continue;                                512          continue;
522       }                                           513       }
523                                                   514 
524       G4ThreeVector normal = GetNormal(gxx[i],    515       G4ThreeVector normal = GetNormal(gxx[i], true);
525       if (normal * gv <= 0)                    << 516       if (normal * gv <= 0) {
526       {                                        << 
527          // particle goes toword inside of sol    517          // particle goes toword inside of solid, return kInfinity
528 #ifdef G4TWISTDEBUG                               518 #ifdef G4TWISTDEBUG
529           G4cout << "   G4VTwistSurface::Dista << 519           G4cout << "   G4VTwistSurface::DistanceToOut(p,v): normal*gv < 0, normal " 
530                  << fName << " " << normal        520                  << fName << " " << normal 
531                  << G4endl;                       521                  << G4endl;
532 #endif                                            522 #endif 
533       }                                        << 523       } else {
534       else                                     << 
535       {                                        << 
536          // gxx[i] is accepted.                   524          // gxx[i] is accepted.
537          if (distance[i] < bestdistance)       << 525          if (distance[i] < bestdistance) {
538          {                                     << 
539             bestdistance = distance[i];           526             bestdistance = distance[i];
540             gxxbest = gxx[i];                     527             gxxbest = gxx[i];
541          }                                        528          }
542       }                                           529       } 
543    }                                              530    }
544                                                   531 
545 #ifdef G4TWISTDEBUG                               532 #ifdef G4TWISTDEBUG
546    if (besti < 0)                              << 533    if (besti < 0) {
547    {                                           << 534       G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) - return ~~~" << G4endl;
548       G4cout << "~~ G4VTwistSurface::DistanceT << 
549       G4cout << "      No intersections   " <<    535       G4cout << "      No intersections   " << G4endl; 
550       G4cout << "      Name     : " << fName <    536       G4cout << "      Name     : " << fName << G4endl; 
551       G4cout << "      bestdist : " << bestdis    537       G4cout << "      bestdist : " << bestdistance << G4endl;
552       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    538       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
553    }                                           << 539    } else {
554    else                                        << 540       G4cout << "~~~~~ G4VTwistSurface::DistanceToOut(p,v) : return ~~~" << G4endl;
555    {                                           << 
556       G4cout << "~~ G4VTwistSurface::DistanceT << 
557       G4cout << "      Name, i  : " << fName <    541       G4cout << "      Name, i  : " << fName << " , " << i << G4endl; 
558       G4cout << "      gxx[i]   : " << gxxbest    542       G4cout << "      gxx[i]   : " << gxxbest << G4endl; 
559       G4cout << "      bestdist : " << bestdis    543       G4cout << "      bestdist : " << bestdistance << G4endl;
560       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    544       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
561    }                                              545    } 
562 #endif                                            546 #endif
563                                                   547 
564    return bestdistance;                           548    return bestdistance;
565 }                                                 549 }
566                                                   550 
567 //============================================    551 //=====================================================================
568 //* DistanceTo(p) ----------------------------    552 //* DistanceTo(p) -----------------------------------------------------
569                                                   553 
570 G4double G4VTwistSurface::DistanceTo(const G4T << 554 G4double G4VTwistSurface::DistanceTo(const G4ThreeVector &gp,
571                                            G4T << 555                                       G4ThreeVector &gxxbest)
572 {                                                 556 {
573 #ifdef G4TWISTDEBUG                               557 #ifdef G4TWISTDEBUG
574    G4cout << "~~~~~ G4VTwistSurface::DistanceT    558    G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - Start ~~~~~~~~~" << G4endl;
575    G4cout << "      Name : " << fName << G4end    559    G4cout << "      Name : " << fName << G4endl;
576    G4cout << "      gp   : " << gp << G4endl;     560    G4cout << "      gp   : " << gp << G4endl;
577    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    561    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
578 #endif                                            562 #endif
579                                                   563 
580                                                   564 
581    G4ThreeVector gxx[G4VSURFACENXX];              565    G4ThreeVector gxx[G4VSURFACENXX];
582    G4double      distance[G4VSURFACENXX]  ;       566    G4double      distance[G4VSURFACENXX]  ; 
583    G4int         areacode[G4VSURFACENXX]  ;       567    G4int         areacode[G4VSURFACENXX]  ;
584                                                   568    
585    for (G4int i = 0 ; i<G4VSURFACENXX ; ++i )  << 569    for (G4int i = 0 ; i<G4VSURFACENXX ; i++ ) {
586    {                                           << 
587      distance[i] = kInfinity ;                    570      distance[i] = kInfinity ;
588      areacode[i] = sOutside ;                     571      areacode[i] = sOutside ;
589    }                                              572    }
590                                                   573 
                                                   >> 574 
591    DistanceToSurface(gp, gxx, distance, areaco    575    DistanceToSurface(gp, gxx, distance, areacode);
592    gxxbest = gxx[0];                              576    gxxbest = gxx[0];
593                                                   577 
594 #ifdef G4TWISTDEBUG                               578 #ifdef G4TWISTDEBUG
595    G4cout << "~~~~~ G4VTwistSurface::DistanceT    579    G4cout << "~~~~~ G4VTwistSurface::DistanceTo(p) - return ~~~~~~~~" << G4endl;
596    G4cout << "      Name     : " << fName << G    580    G4cout << "      Name     : " << fName << G4endl; 
597    G4cout << "      gxx      : " << gxxbest <<    581    G4cout << "      gxx      : " << gxxbest << G4endl; 
598    G4cout << "      bestdist : " << distance[0    582    G4cout << "      bestdist : " << distance[0] << G4endl;
599    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    583    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << G4endl;
600 #endif                                            584 #endif
601                                                   585 
602    return distance[0];                            586    return distance[0];
603 }                                                 587 }
604                                                   588 
605 //============================================    589 //=====================================================================
606 //* IsSameBoundary ---------------------------    590 //* IsSameBoundary ----------------------------------------------------
607                                                   591 
608 G4bool                                         << 592 G4bool G4VTwistSurface::IsSameBoundary(G4VTwistSurface *surface1, G4int areacode1,
609 G4VTwistSurface::IsSameBoundary(G4VTwistSurfac << 593                                   G4VTwistSurface *surface2, G4int areacode2 ) const
610                                 G4VTwistSurfac << 
611 {                                                 594 {
612    //                                             595    //
613    // IsSameBoundary                              596    // IsSameBoundary
614    //                                             597    //
615    // checking tool whether two boundaries on     598    // checking tool whether two boundaries on different surfaces are same or not.
616    //                                             599    //
617                                                   600 
618    G4bool testbitmode = true;                     601    G4bool testbitmode = true;
619    G4bool iscorner[2] = {IsCorner(areacode1, t    602    G4bool iscorner[2] = {IsCorner(areacode1, testbitmode), 
620                          IsCorner(areacode2, t    603                          IsCorner(areacode2, testbitmode)};
621                                                   604 
622    if (iscorner[0] && iscorner[1])             << 605    if (iscorner[0] && iscorner[1]) {
623    {                                           << 
624       // on corner                                606       // on corner 
625       G4ThreeVector corner1 =                     607       G4ThreeVector corner1 = 
626            surf1->ComputeGlobalPoint(surf1->Ge << 608            surface1->ComputeGlobalPoint(surface1->GetCorner(areacode1));
627       G4ThreeVector corner2 =                     609       G4ThreeVector corner2 = 
628            surf2->ComputeGlobalPoint(surf2->Ge << 610            surface2->ComputeGlobalPoint(surface2->GetCorner(areacode2));
629                                                   611 
630       return (corner1 - corner2).mag() < kCarT << 612       if ((corner1 - corner2).mag() < kCarTolerance) {
631    }                                           << 613          return true;
632    else if ((IsBoundary(areacode1, testbitmode << 614       } else {
633             (IsBoundary(areacode2, testbitmode << 615          return false;
634    {                                           << 616       }
                                                   >> 617     
                                                   >> 618    } else if ((IsBoundary(areacode1, testbitmode) && (!iscorner[0])) &&
                                                   >> 619               (IsBoundary(areacode2, testbitmode) && (!iscorner[1]))) {
635       // on boundary                              620       // on boundary  
636       G4ThreeVector d1, d2, ld1, ld2;             621       G4ThreeVector d1, d2, ld1, ld2;
637       G4ThreeVector x01, x02, lx01, lx02;         622       G4ThreeVector x01, x02, lx01, lx02;
638       G4int         type1, type2;                 623       G4int         type1, type2;
639       surf1->GetBoundaryParameters(areacode1,  << 624       surface1->GetBoundaryParameters(areacode1, ld1, lx01, type1);
640       surf2->GetBoundaryParameters(areacode2,  << 625       surface2->GetBoundaryParameters(areacode2, ld2, lx02, type2);
641                                                   626 
642       x01 = surf1->ComputeGlobalPoint(lx01);   << 627       x01 = surface1->ComputeGlobalPoint(lx01);
643       x02 = surf2->ComputeGlobalPoint(lx02);   << 628       x02 = surface2->ComputeGlobalPoint(lx02);
644       d1  = surf1->ComputeGlobalDirection(ld1) << 629       d1  = surface1->ComputeGlobalDirection(ld1);
645       d2  = surf2->ComputeGlobalDirection(ld2) << 630       d2  = surface2->ComputeGlobalDirection(ld2);
                                                   >> 631 
                                                   >> 632       if ((x01 - x02).mag() < kCarTolerance &&
                                                   >> 633           (d1 - d2).mag() < kCarTolerance) {
                                                   >> 634         return true;
                                                   >> 635       } else {
                                                   >> 636         return false;
                                                   >> 637       }
646                                                   638 
647       return (x01 - x02).mag() < kCarTolerance << 639    } else {
648             && (d1 - d2).mag() < kCarTolerance << 
649    }                                           << 
650    else                                        << 
651    {                                           << 
652       return false;                               640       return false;
653    }                                              641    }
654 }                                                 642 }
655                                                   643 
656 //============================================    644 //=====================================================================
657 //* GetBoundaryParameters --------------------    645 //* GetBoundaryParameters ---------------------------------------------
658                                                   646 
659 void G4VTwistSurface::GetBoundaryParameters(co << 647 void G4VTwistSurface::GetBoundaryParameters(const G4int         &areacode,
660                                              G << 648                                              G4ThreeVector &d,
661                                              G << 649                                              G4ThreeVector &x0,
662                                              G << 650                                              G4int         &boundarytype) const
663 {                                                 651 {
664    // areacode must be one of them:               652    // areacode must be one of them:
665    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,       653    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
666    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.       654    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
667                                                   655    
668    for (const auto & boundary : fBoundaries)   << 656    G4int i;
669    {                                           << 657    for (i=0; i<4; i++) {
670       if (boundary.GetBoundaryParameters(areac << 658       if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0,
671       {                                        << 659                                                boundarytype)) {
672          return;                                  660          return;
673       }                                           661       }
674    }                                              662    }
675                                                   663 
676    std::ostringstream message;                    664    std::ostringstream message;
677    message << "Not registered boundary." << G4    665    message << "Not registered boundary." << G4endl
678            << "        Boundary at areacode "     666            << "        Boundary at areacode " << std::hex << areacode
679            << std::dec << G4endl                  667            << std::dec << G4endl
680            << "        is not registered.";       668            << "        is not registered.";
681    G4Exception("G4VTwistSurface::GetBoundaryPa    669    G4Exception("G4VTwistSurface::GetBoundaryParameters()", "GeomSolids0002",
682                FatalException, message);          670                FatalException, message);
683 }                                                 671 }
684                                                   672 
685 //============================================    673 //=====================================================================
686 //* GetBoundaryAtPZ --------------------------    674 //* GetBoundaryAtPZ ---------------------------------------------------
687                                                   675 
688 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ    676 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ(G4int areacode,
689                                                << 677                                           const G4ThreeVector &p) const
690 {                                                 678 {
691    // areacode must be one of them:               679    // areacode must be one of them:
692    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,       680    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
693    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.       681    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.
694                                                   682 
695    if (((areacode & sAxis0) != 0) && ((areacod << 683    if (areacode & sAxis0 && areacode & sAxis1) {
696    {                                           << 
697      std::ostringstream message;                  684      std::ostringstream message;
698      message << "Point is in the corner area."    685      message << "Point is in the corner area." << G4endl
699              << "        This function returns    686              << "        This function returns "
700              << "a direction vector of a bound    687              << "a direction vector of a boundary line." << G4endl
701              << "        areacode = " << areac    688              << "        areacode = " << areacode;
702      G4Exception("G4VTwistSurface::GetBoundary    689      G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0003",
703                  FatalException, message);        690                  FatalException, message);
704    }                                              691    }
705                                                   692 
706    G4ThreeVector d;                               693    G4ThreeVector d;
707    G4ThreeVector x0;                              694    G4ThreeVector x0;
708    G4int         boundarytype = 0;             << 695    G4int         boundarytype;
709    G4bool        found = false;                   696    G4bool        found = false;
710                                                   697    
711    for (const auto & boundary : fBoundaries)   << 698    for (G4int i=0; i<4; i++) {
712    {                                           << 699       if (fBoundaries[i].GetBoundaryParameters(areacode, d, x0, 
713       if (boundary.GetBoundaryParameters(areac << 700                                                 boundarytype)){
714       {                                        << 
715          found = true;                            701          found = true;
716          continue;                                702          continue;
717       }                                           703       }
718    }                                              704    }
719                                                   705 
720    if (!found)                                 << 706    if (!found) {
721    {                                           << 
722      std::ostringstream message;                  707      std::ostringstream message;
723      message << "Not registered boundary." <<     708      message << "Not registered boundary." << G4endl
724              << "        Boundary at areacode     709              << "        Boundary at areacode " << areacode << G4endl
725              << "        is not registered.";     710              << "        is not registered.";
726      G4Exception("G4VTwistSurface::GetBoundary    711      G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0002",
727                  FatalException, message);        712                  FatalException, message);
728    }                                              713    }
729                                                   714 
730    if (((boundarytype & sAxisPhi) == sAxisPhi)    715    if (((boundarytype & sAxisPhi) == sAxisPhi) ||
731        ((boundarytype & sAxisRho) == sAxisRho) << 716        ((boundarytype & sAxisRho) == sAxisRho)) {
732    {                                           << 
733      std::ostringstream message;                  717      std::ostringstream message;
734      message << "Not a z-depended line boundar    718      message << "Not a z-depended line boundary." << G4endl
735              << "        Boundary at areacode     719              << "        Boundary at areacode " << areacode << G4endl
736              << "        is not a z-depended l    720              << "        is not a z-depended line.";
737      G4Exception("G4VTwistSurface::GetBoundary    721      G4Exception("G4VTwistSurface::GetBoundaryAtPZ()", "GeomSolids0002",
738                  FatalException, message);        722                  FatalException, message);
739    }                                              723    }
740    return ((p.z() - x0.z()) / d.z()) * d + x0;    724    return ((p.z() - x0.z()) / d.z()) * d + x0;
741 }                                                 725 }
742                                                   726 
743 //============================================    727 //=====================================================================
744 //* SetCorner --------------------------------    728 //* SetCorner ---------------------------------------------------------
745                                                   729 
746 void G4VTwistSurface::SetCorner(G4int areacode << 730 void G4VTwistSurface::SetCorner(G4int areacode, G4double x, G4double y, G4double z)
747                                 G4double x, G4 << 
748 {                                                 731 {
749    if ((areacode & sCorner) != sCorner)        << 732    if ((areacode & sCorner) != sCorner){
750    {                                           << 
751      std::ostringstream message;                  733      std::ostringstream message;
752      message << "Area code must represents cor    734      message << "Area code must represents corner." << G4endl
753              << "        areacode " << areacod    735              << "        areacode " << areacode;
754      G4Exception("G4VTwistSurface::SetCorner()    736      G4Exception("G4VTwistSurface::SetCorner()", "GeomSolids0002",
755                  FatalException, message);        737                  FatalException, message);
756    }                                              738    }
757                                                   739 
758    if ((areacode & sC0Min1Min) == sC0Min1Min)     740    if ((areacode & sC0Min1Min) == sC0Min1Min) {
759       fCorners[0].set(x, y, z);                   741       fCorners[0].set(x, y, z);
760    } else if ((areacode & sC0Max1Min) == sC0Ma    742    } else if ((areacode & sC0Max1Min) == sC0Max1Min) {
761       fCorners[1].set(x, y, z);                   743       fCorners[1].set(x, y, z);
762    } else if ((areacode & sC0Max1Max) == sC0Ma    744    } else if ((areacode & sC0Max1Max) == sC0Max1Max) {
763       fCorners[2].set(x, y, z);                   745       fCorners[2].set(x, y, z);
764    } else if ((areacode & sC0Min1Max) == sC0Mi    746    } else if ((areacode & sC0Min1Max) == sC0Min1Max) {
765       fCorners[3].set(x, y, z);                   747       fCorners[3].set(x, y, z);
766    }                                              748    }
767 }                                                 749 }
768                                                   750 
769 //============================================    751 //=====================================================================
770 //* SetBoundaryAxis --------------------------    752 //* SetBoundaryAxis ---------------------------------------------------
771                                                   753 
772 void G4VTwistSurface::GetBoundaryAxis(G4int ar    754 void G4VTwistSurface::GetBoundaryAxis(G4int areacode, EAxis axis[]) const
773 {                                                 755 {
774    if ((areacode & sBoundary) != sBoundary) {     756    if ((areacode & sBoundary) != sBoundary) {
775      G4Exception("G4VTwistSurface::GetBoundary    757      G4Exception("G4VTwistSurface::GetBoundaryAxis()", "GeomSolids0003",
776                  FatalException, "Not located     758                  FatalException, "Not located on a boundary!");
777    }                                              759    }
778    for (G4int i=0; i<2; ++i)                   << 760    G4int i;
779    {                                           << 761    for (i=0; i<2; i++) {
                                                   >> 762 
780       G4int whichaxis = 0 ;                       763       G4int whichaxis = 0 ;
781       if (i == 0) {                               764       if (i == 0) {
782          whichaxis = sAxis0;                      765          whichaxis = sAxis0;
783       } else if (i == 1) {                        766       } else if (i == 1) {
784          whichaxis = sAxis1;                      767          whichaxis = sAxis1;
785       }                                           768       }
786                                                   769       
787       // extracted axiscode of whichaxis          770       // extracted axiscode of whichaxis
788       G4int axiscode = whichaxis & sAxisMask &    771       G4int axiscode = whichaxis & sAxisMask & areacode ; 
789       if (axiscode != 0) {                     << 772       if (axiscode) {
790          if (axiscode == (whichaxis & sAxisX))    773          if (axiscode == (whichaxis & sAxisX)) {
791             axis[i] = kXAxis;                     774             axis[i] = kXAxis;
792          } else if (axiscode == (whichaxis & s    775          } else if (axiscode == (whichaxis & sAxisY)) {
793             axis[i] = kYAxis;                     776             axis[i] = kYAxis;
794          } else if (axiscode == (whichaxis & s    777          } else if (axiscode == (whichaxis & sAxisZ)) {
795             axis[i] = kZAxis;                     778             axis[i] = kZAxis;
796          } else if (axiscode == (whichaxis & s    779          } else if (axiscode == (whichaxis & sAxisRho)) {
797             axis[i] = kRho;                       780             axis[i] = kRho;
798          } else if (axiscode == (whichaxis & s    781          } else if (axiscode == (whichaxis & sAxisPhi)) {
799             axis[i] = kPhi;                       782             axis[i] = kPhi;
800          } else {                                 783          } else {
801            std::ostringstream message;            784            std::ostringstream message;
802            message << "Not supported areacode.    785            message << "Not supported areacode." << G4endl
803                    << "        areacode " << a    786                    << "        areacode " << areacode;
804            G4Exception("G4VTwistSurface::GetBo    787            G4Exception("G4VTwistSurface::GetBoundaryAxis()", "GeomSolids0001",
805                        FatalException, message    788                        FatalException, message);
806          }                                        789          }
807       }                                           790       }
808    }                                              791    }
809 }                                                 792 }
810                                                   793 
811 //============================================    794 //=====================================================================
812 //* SetBoundaryLimit -------------------------    795 //* SetBoundaryLimit --------------------------------------------------
813                                                   796 
814 void G4VTwistSurface::GetBoundaryLimit(G4int a    797 void G4VTwistSurface::GetBoundaryLimit(G4int areacode, G4double limit[]) const
815 {                                                 798 {
816    if ((areacode & sCorner) != 0) {            << 799    if (areacode & sCorner) {
817       if ((areacode & sC0Min1Min) != 0) {      << 800       if (areacode & sC0Min1Max) {
818          limit[0] = fAxisMin[0];                  801          limit[0] = fAxisMin[0];
819          limit[1] = fAxisMin[1];                  802          limit[1] = fAxisMin[1];
820       } else if ((areacode & sC0Max1Min) != 0) << 803       } else if (areacode & sC0Max1Min) {
821          limit[0] = fAxisMax[0];                  804          limit[0] = fAxisMax[0];
822          limit[1] = fAxisMin[1];                  805          limit[1] = fAxisMin[1];
823       } else if ((areacode & sC0Max1Max) != 0) << 806       } else if (areacode & sC0Max1Max) {
824          limit[0] = fAxisMax[0];                  807          limit[0] = fAxisMax[0];
825          limit[1] = fAxisMax[1];                  808          limit[1] = fAxisMax[1];
826       } else if ((areacode & sC0Min1Max) != 0) << 809       } else if (areacode & sC0Min1Max) {
827          limit[0] = fAxisMin[0];                  810          limit[0] = fAxisMin[0];
828          limit[1] = fAxisMax[1];                  811          limit[1] = fAxisMax[1];
829       }                                           812       }
830    } else if ((areacode & sBoundary) != 0) {   << 813    } else if (areacode & sBoundary) {
831       if ((areacode & (sAxis0 | sAxisMin)) !=  << 814       if (areacode & (sAxis0 | sAxisMin)) {
832          limit[0] = fAxisMin[0];                  815          limit[0] = fAxisMin[0];
833       } else if ((areacode & (sAxis1 | sAxisMi << 816       } else if (areacode & (sAxis1 | sAxisMin)) {
834          limit[0] = fAxisMin[1];                  817          limit[0] = fAxisMin[1];
835       } else if ((areacode & (sAxis0 | sAxisMa << 818       } else if (areacode & (sAxis0 | sAxisMax)) {
836          limit[0] = fAxisMax[0];                  819          limit[0] = fAxisMax[0];
837       } else if ((areacode & (sAxis1 | sAxisMa << 820       } else if (areacode & (sAxis1 | sAxisMax)) {
838          limit[0] = fAxisMax[1];                  821          limit[0] = fAxisMax[1];
839       }                                           822       }
840    } else {                                       823    } else {
841      std::ostringstream message;                  824      std::ostringstream message;
842      message << "Not located on a boundary!" <    825      message << "Not located on a boundary!" << G4endl
843              << "          areacode " << areac    826              << "          areacode " << areacode;
844      G4Exception("G4VTwistSurface::GetBoundary    827      G4Exception("G4VTwistSurface::GetBoundaryLimit()", "GeomSolids1002",
845                  JustWarning, message);           828                  JustWarning, message);
846    }                                              829    }
847 }                                                 830 }
848                                                   831 
849 //============================================    832 //=====================================================================
850 //* SetBoundary ------------------------------    833 //* SetBoundary -------------------------------------------------------
851                                                   834 
852 void G4VTwistSurface::SetBoundary(const G4int& << 835 void G4VTwistSurface::SetBoundary(const G4int         &axiscode,
853                                   const G4Thre << 836                              const G4ThreeVector &direction,
854                                   const G4Thre << 837                              const G4ThreeVector &x0,
855                                   const G4int& << 838                              const G4int         &boundarytype)
856 {                                                 839 {
857    G4int code = (~sAxisMask) & axiscode;          840    G4int code = (~sAxisMask) & axiscode;
858    if ((code == (sAxis0 & sAxisMin)) ||           841    if ((code == (sAxis0 & sAxisMin)) ||
859        (code == (sAxis0 & sAxisMax)) ||           842        (code == (sAxis0 & sAxisMax)) ||
860        (code == (sAxis1 & sAxisMin)) ||           843        (code == (sAxis1 & sAxisMin)) ||
861        (code == (sAxis1 & sAxisMax)))          << 844        (code == (sAxis1 & sAxisMax))) {
862    {                                           << 845 
                                                   >> 846       G4int i;
863       G4bool done = false;                        847       G4bool done = false;
864       for (auto & boundary : fBoundaries)      << 848       for (i=0; i<4; i++) {
865       {                                        << 849          if (fBoundaries[i].IsEmpty()) {
866          if (boundary.IsEmpty())               << 850             fBoundaries[i].SetFields(axiscode, direction,
867          {                                     << 851                                      x0, boundarytype);
868             boundary.SetFields(axiscode, direc << 
869             done = true;                          852             done = true;
870             break;                                853             break;
871          }                                        854          }
872       }                                           855       }
873                                                   856 
874       if (!done)                               << 857       if (!done) {
875       {                                        << 
876          G4Exception("G4VTwistSurface::SetBoun    858          G4Exception("G4VTwistSurface::SetBoundary()", "GeomSolids0003",
877                       FatalException, "Number     859                       FatalException, "Number of boundary exceeding 4!");
878       }                                           860       }
879    }                                           << 861    } else {
880    else                                        << 
881    {                                           << 
882       std::ostringstream message;                 862       std::ostringstream message;
883       message << "Invalid axis-code." << G4end    863       message << "Invalid axis-code." << G4endl
884               << "        axiscode = "            864               << "        axiscode = "
885               << std::hex << axiscode << std::    865               << std::hex << axiscode << std::dec;
886       G4Exception("G4VTwistSurface::SetBoundar    866       G4Exception("G4VTwistSurface::SetBoundary()", "GeomSolids0003",
887                   FatalException, message);       867                   FatalException, message);
888    }                                              868    }
889 }                                                 869 }
890                                                   870 
891 //============================================    871 //=====================================================================
892 //* GetFace ----------------------------------    872 //* GetFace -----------------------------------------------------------
893                                                   873 
894 G4int G4VTwistSurface::GetFace( G4int i, G4int << 874 G4int G4VTwistSurface::GetFace( G4int i, G4int j, G4int m,
895                                 G4int n, G4int    875                                 G4int n, G4int iside ) 
896 {                                                 876 {
897   // this is the face mapping function            877   // this is the face mapping function
898   // (i,j) -> face number                         878   // (i,j) -> face number
899                                                   879 
900   if ( iside == 0 ) {                             880   if ( iside == 0 ) {
901     return i * ( k - 1 ) + j ;                 << 881     return i * ( m - 1 ) + j ;
902   }                                               882   }
903                                                   883 
904   else if ( iside == 1 ) {                        884   else if ( iside == 1 ) {
905     return (k-1)*(k-1) + i*(k-1) + j ;         << 885     return (m-1)*(m-1) + i*(m-1) + j ;
906   }                                               886   }
907                                                   887 
908   else if ( iside == 2 ) {                        888   else if ( iside == 2 ) {
909     return 2*(k-1)*(k-1) + i*(k-1) + j ;       << 889     return 2*(m-1)*(m-1) + i*(m-1) + j ;
910   }                                               890   }
911                                                   891 
912   else if ( iside == 3 ) {                        892   else if ( iside == 3 ) {
913     return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k- << 893     return 2*(m-1)*(m-1) + (n-1)*(m-1) + i*(m-1) + j ;
914   }                                               894   }
915                                                   895   
916   else if ( iside == 4 ) {                        896   else if ( iside == 4 ) {
917     return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*( << 897     return 2*(m-1)*(m-1) + 2*(n-1)*(m-1) + i*(m-1) + j ;
918   }                                               898   }
919                                                   899   
920   else if ( iside == 5 ) {                        900   else if ( iside == 5 ) {
921     return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*( << 901     return 2*(m-1)*(m-1) + 3*(n-1)*(m-1) + i*(m-1) + j ;
922   }                                               902   }
923                                                   903 
924   else {                                          904   else {
                                                   >> 905 
925     std::ostringstream message;                   906     std::ostringstream message;
926     message << "Not correct side number: "        907     message << "Not correct side number: "
927             << GetName() << G4endl                908             << GetName() << G4endl
928             << "iside is " << iside << " but s    909             << "iside is " << iside << " but should be "
929             << "0,1,2,3,4 or 5" << ".";           910             << "0,1,2,3,4 or 5" << ".";
930     G4Exception("G4TwistSurface::G4GetFace()",    911     G4Exception("G4TwistSurface::G4GetFace()", "GeomSolids0002",
931                 FatalException, message);         912                 FatalException, message);
                                                   >> 913 
                                                   >> 914 
932   }                                               915   }
933                                                   916 
934   return -1 ;  // wrong face                      917   return -1 ;  // wrong face
935 }                                                 918 }
936                                                   919 
937 //============================================    920 //=====================================================================
938 //* GetNode ----------------------------------    921 //* GetNode -----------------------------------------------------------
939                                                   922 
940 G4int G4VTwistSurface::GetNode( G4int i, G4int << 923 G4int G4VTwistSurface::GetNode( G4int i, G4int j, G4int m,
941                                 G4int n, G4int    924                                 G4int n, G4int iside ) 
942 {                                                 925 {
943   // this is the node mapping function            926   // this is the node mapping function
944   // (i,j) -> node number                         927   // (i,j) -> node number
945   // Depends on the side iside and the used me    928   // Depends on the side iside and the used meshing of the surface
946                                                   929 
947   if ( iside == 0 )                            << 930   if ( iside == 0 ) {
948   {                                            << 931     // lower endcap is mxm squared. 
949     // lower endcap is kxk squared.            << 932     // n = m 
950     // n = k                                   << 933     return i * m + j ;
951     return i * k + j ;                         << 
952   }                                               934   }
953                                                   935 
954   if ( iside == 1 )                            << 936   if ( iside == 1 ) {
955   {                                            << 937     // upper endcap is mxm squared. Shift by m*m
956     // upper endcap is kxk squared. Shift by k << 938     // n = m 
957     // n = k                                   << 939     return  m*m + i*m + j ;
958     return  k*k + i*k + j ;                    << 
959   }                                               940   }
960                                                   941 
961   else if ( iside == 2 )                       << 942   else if ( iside == 2 ) {
962   {                                            << 
963     // front side.                                943     // front side.
964     if      ( i == 0 )   { return       j ;  } << 944     if      ( i == 0 )     {   return       j ;  }
965     else if ( i == n-1 ) { return k*k + j ;  } << 945     else if ( i == n-1 )   {   return m*m + j ;  } 
966     else                 { return 2*k*k + 4*(i << 946     else                   {   return 2*m*m + 4*(i-1)*(m-1) + j ; }
967   }                                               947   }
968                                                   948 
969   else if ( iside == 3 )                       << 949   else if ( iside == 3 ) {
970   {                                            << 
971     // right side                                 950     // right side
972     if      ( i == 0 )   { return       (j+1)* << 951     if      ( i == 0 )     {   return       (j+1)*m - 1 ; } 
973     else if ( i == n-1 ) { return k*k + (j+1)* << 952     else if ( i == n-1 )   {   return m*m + (j+1)*m - 1 ; }
974     else                                       << 953     else                   {   return 2*m*m + 4*(i-1)*(m-1) + (m-1) + j ; }
975     {                                          << 
976       return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j << 
977     }                                          << 
978   }                                               954   }
979   else if ( iside == 4 )                       << 955   else if ( iside == 4 ) {
980   {                                            << 956     // back side.
981     // back side                               << 957     if      ( i == 0 )     {   return   m*m - 1 - j ; }               // reversed order
982     if      ( i == 0 )   { return   k*k - 1 -  << 958     else if ( i == n-1 )   {   return 2*m*m - 1 - j ; }               // reversed order 
983     else if ( i == n-1 ) { return 2*k*k - 1 -  << 959     else                   {   return 2*m*m + 4*(i-1)*(m-1) + 2*(m-1) + j ; // normal order
984     else                                       << 
985     {                                          << 
986       return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + << 
987     }                                             960     }
988   }                                               961   }
989   else if ( iside == 5 )                       << 962   else if ( iside == 5 ) {
990   {                                            << 
991     // left side                                  963     // left side 
992     if      ( i == 0 )   { return k*k   - (j+1 << 964     if      ( i == 0 )     {   return m*m   - (j+1)*m ; }             // reversed order
993     else if ( i == n-1)  { return 2*k*k - (j+1 << 965     else if ( i == n-1)    {   return 2*m*m - (j+1)*m ; }             // reverded order
994     else                                       << 966     else {
995     {                                          << 967       if ( j == m-1 )      {   return 2*m*m + 4*(i-1)*(m-1) ; }       // special case
996       if ( j == k-1 ) { return 2*k*k + 4*(i-1) << 968       else                 {   return 2*m*m + 4*(i-1)*(m-1) + 3*(m-1) + j ; }  // normal order
997       else                                     << 969     }  
998       {                                        << 
999         return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) << 
1000       }                                       << 
1001     }                                         << 
1002   }                                              970   }
1003   else                                        << 971 
1004   {                                           << 972   else {
                                                   >> 973 
1005     std::ostringstream message;                  974     std::ostringstream message;
1006     message << "Not correct side number: "       975     message << "Not correct side number: "
1007             << GetName() << G4endl               976             << GetName() << G4endl
1008             << "iside is " << iside << " but     977             << "iside is " << iside << " but should be "
1009             << "0,1,2,3,4 or 5" << ".";          978             << "0,1,2,3,4 or 5" << ".";
1010     G4Exception("G4TwistSurface::G4GetNode()"    979     G4Exception("G4TwistSurface::G4GetNode()", "GeomSolids0002",
1011                 FatalException, message);        980                 FatalException, message);
1012   }                                              981   } 
1013   return -1 ;  // wrong node                     982   return -1 ;  // wrong node 
1014 }                                                983 } 
1015                                                  984 
1016 //===========================================    985 //=====================================================================
1017 //* GetEdgeVisiblility ----------------------    986 //* GetEdgeVisiblility ------------------------------------------------
1018                                                  987 
1019 G4int G4VTwistSurface::GetEdgeVisibility( G4i << 988 G4int G4VTwistSurface::GetEdgeVisibility( G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation) 
1020                                           G4i << 
1021 {                                                989 {
                                                   >> 990 
1022   // clockwise filling         -> positive or    991   // clockwise filling         -> positive orientation
1023   // counter clockwise filling -> negative or    992   // counter clockwise filling -> negative orientation 
1024                                                  993 
1025   //                                             994   //
1026   //   d    C    c                               995   //   d    C    c
1027   //     +------+                                996   //     +------+ 
1028   //     |      |                                997   //     |      |
1029   //     |      |                                998   //     |      |
1030   //     |      |                                999   //     |      |
1031   //   D |      |B                               1000   //   D |      |B
1032   //     |      |                                1001   //     |      |
1033   //     |      |                                1002   //     |      |
1034   //     |      |                                1003   //     |      |
1035   //     +------+                                1004   //     +------+
1036   //    a   A    b                               1005   //    a   A    b
1037   //                                             1006   //    
1038   //  a = +--+    A = ---+                       1007   //  a = +--+    A = ---+
1039   //  b = --++    B = --+-                       1008   //  b = --++    B = --+-
1040   //  c = -++-    C = -+--                       1009   //  c = -++-    C = -+--
1041   //  d = ++--    D = +---                       1010   //  d = ++--    D = +---
1042                                                  1011 
1043                                                  1012 
1044   // check first invisible faces                 1013   // check first invisible faces
1045                                                  1014 
1046   if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) ) << 1015   if ( ( i>0 && i<n-2 ) && ( j>0 && j<m-2 ) ) {
1047   {                                           << 
1048     return -1 ;   // always invisible, signs:    1016     return -1 ;   // always invisible, signs:   ----
1049   }                                              1017   }
1050                                                  1018   
1051   // change first the vertex number (depends     1019   // change first the vertex number (depends on the orientation)
1052   // 0,1,2,3  -> 3,2,1,0                         1020   // 0,1,2,3  -> 3,2,1,0
1053   if ( orientation < 0 ) { number = ( 3 - num    1021   if ( orientation < 0 ) { number = ( 3 - number ) ;  }
1054                                                  1022   
1055   // check true edges                            1023   // check true edges
1056   if ( ( j>=1 && j<=k-3 ) )                   << 1024   if ( ( j>=1 && j<=m-3 ) ) {
1057   {                                           << 1025 
1058     if ( i == 0 )  {        // signs (A):  --    1026     if ( i == 0 )  {        // signs (A):  ---+  
1059       return ( number == 3 ) ? 1 : -1 ;          1027       return ( number == 3 ) ? 1 : -1 ;
1060     }                                            1028     }
1061                                                  1029     
1062     else if ( i == n-2 ) {  // signs (C):  -+    1030     else if ( i == n-2 ) {  // signs (C):  -+--
1063       return  ( number == 1 ) ? 1 : -1 ;         1031       return  ( number == 1 ) ? 1 : -1 ; 
1064     }                                            1032     }
1065                                                  1033     
1066     else                                      << 1034     else {
1067     {                                         << 
1068       std::ostringstream message;                1035       std::ostringstream message;
1069       message << "Not correct face number: "     1036       message << "Not correct face number: " << GetName() << " !";
1070       G4Exception("G4TwistSurface::G4GetEdgeV    1037       G4Exception("G4TwistSurface::G4GetEdgeVisibility()",
1071                   "GeomSolids0003", FatalExce    1038                   "GeomSolids0003", FatalException, message);
1072     }                                            1039     }
1073   }                                              1040   }
1074                                                  1041   
1075   if ( ( i>=1 && i<=n-3 ) )                   << 1042   if ( ( i>=1 && i<=n-3 ) ) {
1076   {                                           << 1043 
1077     if ( j == 0 )  {        // signs (D):  +-    1044     if ( j == 0 )  {        // signs (D):  +---
1078       return ( number == 0 ) ? 1 : -1 ;          1045       return ( number == 0 ) ? 1 : -1 ;
1079     }                                            1046     }
1080                                                  1047 
1081     else if ( j == k-2 ) {  // signs (B):  -- << 1048     else if ( j == m-2 ) {  // signs (B):  --+-
1082       return  ( number == 2 ) ? 1 : -1 ;         1049       return  ( number == 2 ) ? 1 : -1 ; 
1083     }                                            1050     }
1084                                                  1051 
1085     else                                      << 1052     else {
1086     {                                         << 
1087       std::ostringstream message;                1053       std::ostringstream message;
1088       message << "Not correct face number: "     1054       message << "Not correct face number: " << GetName() << " !";
1089       G4Exception("G4TwistSurface::G4GetEdgeV    1055       G4Exception("G4TwistSurface::G4GetEdgeVisibility()",
1090                   "GeomSolids0003", FatalExce    1056                   "GeomSolids0003", FatalException, message);
1091     }                                            1057     }
1092   }                                              1058   }
1093                                                  1059   
1094   // now the corners                             1060   // now the corners
1095   if ( i == 0 && j == 0 ) {          // signs    1061   if ( i == 0 && j == 0 ) {          // signs (a) : +--+
1096     return ( number == 0 || number == 3 ) ? 1    1062     return ( number == 0 || number == 3 ) ? 1 : -1 ;
1097   }                                              1063   }
1098   else if ( i == 0 && j == k-2 ) {   // signs << 1064   else if ( i == 0 && j == m-2 ) {   // signs (b) : --++
1099     return ( number == 2 || number == 3 ) ? 1    1065     return ( number == 2 || number == 3 ) ? 1 : -1 ;
1100   }                                              1066   }
1101   else if ( i == n-2 && j == k-2 ) { // signs << 1067   else if ( i == n-2 && j == m-2 ) { // signs (c) : -++-
1102     return ( number == 1 || number == 2 ) ? 1    1068     return ( number == 1 || number == 2 ) ? 1 : -1 ;
1103   }                                              1069   }
1104   else if ( i == n-2 && j == 0 ) {   // signs    1070   else if ( i == n-2 && j == 0 ) {   // signs (d) : ++--
1105     return ( number == 0 || number == 1 ) ? 1    1071     return ( number == 0 || number == 1 ) ? 1 : -1 ;
1106   }                                              1072   }
1107   else                                        << 1073   else {
1108   {                                           << 
1109     std::ostringstream message;                  1074     std::ostringstream message;
1110     message << "Not correct face number: " <<    1075     message << "Not correct face number: " << GetName() << " !";
1111     G4Exception("G4TwistSurface::G4GetEdgeVis    1076     G4Exception("G4TwistSurface::G4GetEdgeVisibility()",
1112                 "GeomSolids0003", FatalExcept    1077                 "GeomSolids0003", FatalException, message);
1113   }                                              1078   }
1114                                                  1079 
1115   std::ostringstream message;                    1080   std::ostringstream message;
1116   message << "Not correct face number: " << G    1081   message << "Not correct face number: " << GetName() << " !";
1117   G4Exception("G4TwistSurface::G4GetEdgeVisib    1082   G4Exception("G4TwistSurface::G4GetEdgeVisibility()", "GeomSolids0003",
1118               FatalException, message);          1083               FatalException, message);
1119                                                  1084 
1120   return 0 ;                                     1085   return 0 ;
1121 }                                                1086 }
1122                                                  1087 
1123                                                  1088 
1124 //===========================================    1089 //=====================================================================
1125 //* DebugPrint ------------------------------    1090 //* DebugPrint --------------------------------------------------------
1126                                                  1091 
1127 void G4VTwistSurface::DebugPrint() const         1092 void G4VTwistSurface::DebugPrint() const
1128 {                                                1093 {
1129    G4ThreeVector A = fRot * GetCorner(sC0Min1    1094    G4ThreeVector A = fRot * GetCorner(sC0Min1Min) + fTrans;
1130    G4ThreeVector B = fRot * GetCorner(sC0Max1    1095    G4ThreeVector B = fRot * GetCorner(sC0Max1Min) + fTrans;
1131    G4ThreeVector C = fRot * GetCorner(sC0Max1    1096    G4ThreeVector C = fRot * GetCorner(sC0Max1Max) + fTrans;
1132    G4ThreeVector D = fRot * GetCorner(sC0Min1    1097    G4ThreeVector D = fRot * GetCorner(sC0Min1Max) + fTrans;
1133                                                  1098   
1134    G4cout << "/* G4VTwistSurface::DebugPrint( << 1099    G4cout << "/* G4VTwistSurface::DebugPrint():-------------------------------"
1135           << G4endl;                             1100           << G4endl;
1136    G4cout << "/* Name = " << fName << G4endl;    1101    G4cout << "/* Name = " << fName << G4endl;
1137    G4cout << "/* Axis = " << std::hex << fAxi    1102    G4cout << "/* Axis = " << std::hex << fAxis[0] << " "
1138           << std::hex << fAxis[1]                1103           << std::hex << fAxis[1] 
1139           << " (0,1,2,3,5 = kXAxis,kYAxis,kZA    1104           << " (0,1,2,3,5 = kXAxis,kYAxis,kZAxis,kRho,kPhi)"
1140           << std::dec << G4endl;                 1105           << std::dec << G4endl;
1141    G4cout << "/* BoundaryLimit(in local) fAxi    1106    G4cout << "/* BoundaryLimit(in local) fAxis0(min, max) = ("<<fAxisMin[0] 
1142           << ", " << fAxisMax[0] << ")" << G4    1107           << ", " << fAxisMax[0] << ")" << G4endl;
1143    G4cout << "/* BoundaryLimit(in local) fAxi    1108    G4cout << "/* BoundaryLimit(in local) fAxis1(min, max) = ("<<fAxisMin[1] 
1144           << ", " << fAxisMax[1] << ")" << G4    1109           << ", " << fAxisMax[1] << ")" << G4endl;
1145    G4cout << "/* Cornar point sC0Min1Min = "     1110    G4cout << "/* Cornar point sC0Min1Min = " << A << G4endl;
1146    G4cout << "/* Cornar point sC0Max1Min = "     1111    G4cout << "/* Cornar point sC0Max1Min = " << B << G4endl;
1147    G4cout << "/* Cornar point sC0Max1Max = "     1112    G4cout << "/* Cornar point sC0Max1Max = " << C << G4endl;
1148    G4cout << "/* Cornar point sC0Min1Max = "     1113    G4cout << "/* Cornar point sC0Min1Max = " << D << G4endl;
1149    G4cout << "/*-----------------------------    1114    G4cout << "/*---------------------------------------------------------"
1150           << G4endl;                             1115           << G4endl;
1151 }                                                1116 }
1152                                                  1117 
1153 //===========================================    1118 //=====================================================================
1154 // G4VTwistSurface::CurrentStatus class          1119 // G4VTwistSurface::CurrentStatus class
1155 //===========================================    1120 //=====================================================================
1156                                                  1121 
1157 //===========================================    1122 //=====================================================================
1158 //* CurrentStatus::CurrentStatus ------------    1123 //* CurrentStatus::CurrentStatus --------------------------------------
1159                                                  1124 
1160 G4VTwistSurface::CurrentStatus::CurrentStatus    1125 G4VTwistSurface::CurrentStatus::CurrentStatus() 
1161 {                                                1126 {
1162   for (size_t i=0; i<G4VSURFACENXX; ++i)      << 1127   for (size_t i=0; i<G4VSURFACENXX; i++)
1163   {                                              1128   {
1164     fDistance[i] = kInfinity;                    1129     fDistance[i] = kInfinity;
1165     fAreacode[i] = sOutside;                     1130     fAreacode[i] = sOutside;
1166     fIsValid[i]  = false;                        1131     fIsValid[i]  = false;
1167     fXX[i].set(kInfinity, kInfinity, kInfinit    1132     fXX[i].set(kInfinity, kInfinity, kInfinity);
1168   }                                              1133   }
1169   fNXX   = 0;                                    1134   fNXX   = 0;
1170   fLastp.set(kInfinity, kInfinity, kInfinity)    1135   fLastp.set(kInfinity, kInfinity, kInfinity);
1171   fLastv.set(kInfinity, kInfinity, kInfinity)    1136   fLastv.set(kInfinity, kInfinity, kInfinity);
1172   fLastValidate = kUninitialized;                1137   fLastValidate = kUninitialized;
1173   fDone = false;                                 1138   fDone = false;
1174 }                                                1139 }
1175                                                  1140 
1176 //===========================================    1141 //=====================================================================
1177 //* CurrentStatus::~CurrentStatus -----------    1142 //* CurrentStatus::~CurrentStatus -------------------------------------
1178                                                  1143 
1179 G4VTwistSurface::CurrentStatus::~CurrentStatu    1144 G4VTwistSurface::CurrentStatus::~CurrentStatus() 
1180 = default;                                    << 1145 {
                                                   >> 1146 }
1181                                                  1147 
1182 //===========================================    1148 //=====================================================================
1183 //* CurrentStatus::SetCurrentStatus ---------    1149 //* CurrentStatus::SetCurrentStatus -----------------------------------
1184                                                  1150 
1185 void                                             1151 void
1186 G4VTwistSurface::CurrentStatus::SetCurrentSta << 1152 G4VTwistSurface::CurrentStatus::SetCurrentStatus(G4int                i, 
1187                                               << 1153                                             G4ThreeVector       &xx, 
1188                                               << 1154                                             G4double            &dist, 
1189                                               << 1155                                             G4int               &areacode, 
1190                                               << 1156                                             G4bool              &isvalid,
1191                                               << 1157                                             G4int                nxx,
1192                                               << 1158                                             EValidate            validate,
1193                                            co << 1159                                       const G4ThreeVector *p, 
1194                                            co << 1160                                       const G4ThreeVector *v)
1195 {                                                1161 {
1196   fDistance[i]  = dist;                          1162   fDistance[i]  = dist;
1197   fAreacode[i]  = areacode;                      1163   fAreacode[i]  = areacode;
1198   fIsValid[i]   = isvalid;                       1164   fIsValid[i]   = isvalid;
1199   fXX[i]        = xx;                            1165   fXX[i]        = xx;
1200   fNXX          = nxx;                           1166   fNXX          = nxx;
1201   fLastValidate = validate;                      1167   fLastValidate = validate;
1202   if (p != nullptr)                           << 1168   if (p)
1203   {                                              1169   {
1204     fLastp = *p;                                 1170     fLastp = *p;
1205   }                                              1171   }
1206   else                                           1172   else
1207   {                                              1173   {
1208     G4Exception("G4VTwistSurface::CurrentStat    1174     G4Exception("G4VTwistSurface::CurrentStatus::SetCurrentStatus()",
1209                 "GeomSolids0003", FatalExcept    1175                 "GeomSolids0003", FatalException, "SetCurrentStatus: p = 0!");
1210   }                                              1176   }
1211   if (v != nullptr)                           << 1177   if (v) 
1212   {                                              1178   {
1213     fLastv = *v;                                 1179     fLastv = *v;
1214   }                                              1180   }
1215   else                                           1181   else
1216   {                                              1182   {
1217     fLastv.set(kInfinity, kInfinity, kInfinit    1183     fLastv.set(kInfinity, kInfinity, kInfinity);
1218   }                                              1184   }
1219   fDone = true;                                  1185   fDone = true;
1220 }                                                1186 }
1221                                                  1187 
1222 //===========================================    1188 //=====================================================================
1223 //* CurrentStatus::ResetfDone ---------------    1189 //* CurrentStatus::ResetfDone -----------------------------------------
1224                                                  1190 
1225 void                                             1191 void
1226 G4VTwistSurface::CurrentStatus::ResetfDone(EV << 1192 G4VTwistSurface::CurrentStatus::ResetfDone(EValidate     validate,
1227                                      const G4 << 1193                                 const G4ThreeVector *p, 
1228                                      const G4 << 1194                                 const G4ThreeVector *v)
1229                                                  1195 
1230 {                                                1196 {
1231   if (validate == fLastValidate && p != nullp << 1197   if (validate == fLastValidate && p && *p == fLastp)
1232   {                                              1198   {
1233      if (v == nullptr || (*v == fLastv)) retu << 1199      if (!v || (v && *v == fLastv)) return;
1234   }                                              1200   }         
1235   G4ThreeVector xx(kInfinity, kInfinity, kInf    1201   G4ThreeVector xx(kInfinity, kInfinity, kInfinity);
1236   for (size_t i=0; i<G4VSURFACENXX; ++i)      << 1202   for (size_t i=0; i<G4VSURFACENXX; i++)
1237   {                                              1203   {
1238     fDistance[i] = kInfinity;                    1204     fDistance[i] = kInfinity;
1239     fAreacode[i] = sOutside;                     1205     fAreacode[i] = sOutside;
1240     fIsValid[i]  = false;                        1206     fIsValid[i]  = false;
1241     fXX[i] = xx;   // bug in old code ( was f    1207     fXX[i] = xx;   // bug in old code ( was fXX[i] =  xx[i]  )
1242   }                                              1208   }
1243   fNXX = 0;                                   << 1209   fNXX   = 0;
1244   fLastp.set(kInfinity, kInfinity, kInfinity)    1210   fLastp.set(kInfinity, kInfinity, kInfinity);
1245   fLastv.set(kInfinity, kInfinity, kInfinity)    1211   fLastv.set(kInfinity, kInfinity, kInfinity);
1246   fLastValidate = kUninitialized;                1212   fLastValidate = kUninitialized;
1247   fDone = false;                                 1213   fDone = false;
1248 }                                                1214 }
1249                                                  1215 
1250 //===========================================    1216 //=====================================================================
1251 //* CurrentStatus::DebugPrint ---------------    1217 //* CurrentStatus::DebugPrint -----------------------------------------
1252                                                  1218 
1253 void                                             1219 void
1254 G4VTwistSurface::CurrentStatus::DebugPrint()     1220 G4VTwistSurface::CurrentStatus::DebugPrint() const
1255 {                                                1221 {
1256   G4cout << "CurrentStatus::Dist0,1= " << fDi    1222   G4cout << "CurrentStatus::Dist0,1= " << fDistance[0] 
1257          << " " << fDistance[1] << " areacode    1223          << " " << fDistance[1] << " areacode = " << fAreacode[0] 
1258          << " " << fAreacode[1] << G4endl;       1224          << " " << fAreacode[1] << G4endl;
1259 }                                                1225 }
1260                                                  1226 
1261 //===========================================    1227 //=====================================================================
1262 // G4VTwistSurface::Boundary class               1228 // G4VTwistSurface::Boundary class
1263 //===========================================    1229 //=====================================================================
1264                                                  1230 
1265 //===========================================    1231 //=====================================================================
                                                   >> 1232 //* Boundary::Boundary ------------------------------------------------
                                                   >> 1233 
                                                   >> 1234 G4VTwistSurface::Boundary::Boundary()
                                                   >> 1235  : fBoundaryAcode(-1), fBoundaryType(0)
                                                   >> 1236 {
                                                   >> 1237 }
                                                   >> 1238 
                                                   >> 1239 //=====================================================================
                                                   >> 1240 //* Boundary::~Boundary -----------------------------------------------
                                                   >> 1241 
                                                   >> 1242 G4VTwistSurface::Boundary::~Boundary()
                                                   >> 1243 {
                                                   >> 1244 }
                                                   >> 1245 
                                                   >> 1246 //=====================================================================
1266 //* Boundary::SetFields ---------------------    1247 //* Boundary::SetFields -----------------------------------------------
1267                                                  1248 
1268 void                                             1249 void
1269 G4VTwistSurface::Boundary::SetFields(const G4 << 1250 G4VTwistSurface::Boundary::SetFields(const G4int         &areacode, 
1270                                 const G4Three << 1251                                 const G4ThreeVector &d, 
1271                                 const G4Three << 1252                                 const G4ThreeVector &x0, 
1272                                 const G4int&  << 1253                                 const G4int         &boundarytype)
1273 {                                                1254 {
1274   fBoundaryAcode     = areacode;                 1255   fBoundaryAcode     = areacode;
1275   fBoundaryDirection = d;                        1256   fBoundaryDirection = d;
1276   fBoundaryX0        = x0;                       1257   fBoundaryX0        = x0;
1277   fBoundaryType      = boundarytype;             1258   fBoundaryType      = boundarytype;
1278 }                                                1259 }
1279                                                  1260 
1280 //===========================================    1261 //=====================================================================
1281 //* Boundary::IsEmpty -----------------------    1262 //* Boundary::IsEmpty -------------------------------------------------
1282                                                  1263 
1283 G4bool G4VTwistSurface::Boundary::IsEmpty() c    1264 G4bool G4VTwistSurface::Boundary::IsEmpty() const 
1284 {                                                1265 {
1285   return fBoundaryAcode == -1;                << 1266   if (fBoundaryAcode == -1) return true;
                                                   >> 1267   return false;
1286 }                                                1268 }
1287                                                  1269 
1288 //===========================================    1270 //=====================================================================
1289 //* Boundary::GetBoundaryParameters ---------    1271 //* Boundary::GetBoundaryParameters -----------------------------------
1290                                                  1272 
1291 G4bool                                           1273 G4bool
1292 G4VTwistSurface::Boundary::GetBoundaryParamet << 1274 G4VTwistSurface::Boundary::GetBoundaryParameters(const G4int         &areacode, 
1293                                               << 1275                                                   G4ThreeVector &d,
1294                                               << 1276                                                   G4ThreeVector &x0, 
1295                                               << 1277                                                   G4int  &boundarytype) const 
1296 {                                                1278 {  
1297   // areacode must be one of them:               1279   // areacode must be one of them:
1298   // sAxis0 & sAxisMin, sAxis0 & sAxisMax,       1280   // sAxis0 & sAxisMin, sAxis0 & sAxisMax,
1299   // sAxis1 & sAxisMin, sAxis1 & sAxisMax        1281   // sAxis1 & sAxisMin, sAxis1 & sAxisMax
1300   //                                          << 1282   if ((areacode & sAxis0) && (areacode & sAxis1))
1301   if (((areacode & sAxis0) != 0) && ((areacod << 
1302   {                                              1283   {
1303     std::ostringstream message;                  1284     std::ostringstream message;
1304     message << "Located in the corner area."     1285     message << "Located in the corner area." << G4endl
1305             << "        This function returns    1286             << "        This function returns a direction vector of "
1306             << "a boundary line." << G4endl      1287             << "a boundary line." << G4endl
1307             << "        areacode = " << areac    1288             << "        areacode = " << areacode;
1308     G4Exception("G4VTwistSurface::Boundary::G    1289     G4Exception("G4VTwistSurface::Boundary::GetBoundaryParameters()",
1309                 "GeomSolids0003", FatalExcept    1290                 "GeomSolids0003", FatalException, message);
1310   }                                              1291   } 
1311   if ((areacode & sSizeMask) != (fBoundaryAco    1292   if ((areacode & sSizeMask) != (fBoundaryAcode & sSizeMask))
1312   {                                              1293   {
1313     return false;                                1294     return false;
1314   }                                              1295   }
1315   d  = fBoundaryDirection;                       1296   d  = fBoundaryDirection;
1316   x0 = fBoundaryX0;                              1297   x0 = fBoundaryX0;
1317   boundarytype = fBoundaryType;                  1298   boundarytype = fBoundaryType;
1318   return true;                                   1299   return true;
1319 }                                                1300 }
1320                                                  1301