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 ]

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