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 7.1.p1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // G4VTwistSurface implementation                 
 27 //                                                
 28 // 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepbu    
 29 // 13-Nov-2003 - O.Link (Oliver.Link@cern.ch),    
 30 //               from original version in Jupi    
 31 // -------------------------------------------    
 32                                                   
 33 #include <iomanip>                                
 34                                                   
 35 #include "G4VTwistSurface.hh"                     
 36 #include "G4GeometryTolerance.hh"                 
 37                                                   
 38 const G4int  G4VTwistSurface::sOutside            
 39 const G4int  G4VTwistSurface::sInside             
 40 const G4int  G4VTwistSurface::sBoundary           
 41 const G4int  G4VTwistSurface::sCorner             
 42 const G4int  G4VTwistSurface::sC0Min1Min          
 43 const G4int  G4VTwistSurface::sC0Max1Min          
 44 const G4int  G4VTwistSurface::sC0Max1Max          
 45 const G4int  G4VTwistSurface::sC0Min1Max          
 46 const G4int  G4VTwistSurface::sAxisMin            
 47 const G4int  G4VTwistSurface::sAxisMax            
 48 const G4int  G4VTwistSurface::sAxisX              
 49 const G4int  G4VTwistSurface::sAxisY              
 50 const G4int  G4VTwistSurface::sAxisZ              
 51 const G4int  G4VTwistSurface::sAxisRho            
 52 const G4int  G4VTwistSurface::sAxisPhi            
 53                                                   
 54 // mask                                           
 55 const G4int  G4VTwistSurface::sAxis0              
 56 const G4int  G4VTwistSurface::sAxis1              
 57 const G4int  G4VTwistSurface::sSizeMask           
 58 const G4int  G4VTwistSurface::sAxisMask           
 59 const G4int  G4VTwistSurface::sAreaMask           
 60                                                   
 61 //============================================    
 62 //* constructors -----------------------------    
 63                                                   
 64 G4VTwistSurface::G4VTwistSurface(const G4Strin    
 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, kI    
 79       fNeighbours[i] = nullptr;                   
 80    }                                              
 81                                                   
 82    fCurrentNormal.p.set(kInfinity, kInfinity,     
 83                                                   
 84    fAmIOnLeftSide.me.set(kInfinity, kInfinity,    
 85    fAmIOnLeftSide.vec.set(kInfinity, kInfinity    
 86    kCarTolerance = G4GeometryTolerance::GetIns    
 87 }                                                 
 88                                                   
 89 G4VTwistSurface::G4VTwistSurface(const G4Strin    
 90                        const G4RotationMatrix&    
 91                        const G4ThreeVector&       
 92                              G4int                
 93                        const EAxis                
 94                        const EAxis                
 95                              G4double             
 96                              G4double             
 97                              G4double             
 98                              G4double             
 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, kI    
114       fNeighbours[i] = nullptr;                   
115    }                                              
116                                                   
117    fCurrentNormal.p.set(kInfinity, kInfinity,     
118                                                   
119    fAmIOnLeftSide.me.set(kInfinity, kInfinity,    
120    fAmIOnLeftSide.vec.set(kInfinity, kInfinity    
121    kCarTolerance = G4GeometryTolerance::GetIns    
122 }                                                 
123                                                   
124 //============================================    
125 //* Fake default constructor -----------------    
126                                                   
127 G4VTwistSurface::G4VTwistSurface( __void__& )     
128   : fHandedness(0), fIsValidNorm(false), kCarT    
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] = fNeighbou    
135 }                                                 
136                                                   
137 //============================================    
138 //* AmIOnLeftSide ----------------------------    
139                                                   
140 G4int G4VTwistSurface::AmIOnLeftSide(const G4T    
141                                      const G4T    
142                                            G4b    
143 {                                                 
144    // AmIOnLeftSide returns phi-location of "m    
145    // (phi relation between me and vec project    
146    // If "me" is on -ve-phi-side of "vec", it     
147    // On the other hand, if "me" is on +ve-phi    
148    // it returns -1.                              
149    // (The return value represents z-coordinat    
150    //  of me.cross(vec).)                         
151    // If me is on boundary of vec, return 0.      
152                                                   
153    const G4double kAngTolerance                   
154      = G4GeometryTolerance::GetInstance()->Get    
155                                                   
156    G4RotationMatrix unitrot;                      
157    const G4RotationMatrix rottol    = unitrot.    
158    const G4RotationMatrix invrottol = unitrot.    
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()    
172    G4ThreeVector vect  = (G4ThreeVector(vec.x(    
173                                                   
174    G4ThreeVector ivect = invrottol * vect;        
175    G4ThreeVector rvect = rottol * vect;           
176                                                   
177    G4double metcrossvect = met.x() * vect.y()     
178                                                   
179    if (withtol)                                   
180    {                                              
181       if (met.x() * ivect.y() - met.y() * ivec    
182           metcrossvect >= 0)  {                   
183          fAmIOnLeftSide.amIOnLeftSide = 1;        
184       } else if (met.x() * rvect.y() - met.y()    
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::Am    
204           << G4endl;                              
205    G4cout << "             Name , returncode      
206                        << fAmIOnLeftSide.amIOn    
207    G4cout << "             me, vec    : " << s    
208                                           << "    
209    G4cout << "             met, vect  : " << m    
210    G4cout << "             ivec, rvec : " << i    
211    G4cout << "             met x vect : " << m    
212    G4cout << "             met x ivec : " << m    
213    G4cout << "             met x rvec : " << m    
214    G4cout << "         =======================    
215           << G4endl;                              
216 #endif                                            
217                                                   
218    return fAmIOnLeftSide.amIOnLeftSide;           
219 }                                                 
220                                                   
221 //============================================    
222 //* DistanceToBoundary -----------------------    
223                                                   
224 G4double G4VTwistSurface::DistanceToBoundary(G    
225                                              G    
226                                        const G    
227 {                                                 
228    // DistanceToBoundary                          
229    //                                             
230    // return distance to nearest boundary from    
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     
238    G4ThreeVector x0;   // reference point of t    
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.    
246               << "        Point is in the corn    
247               << G4endl                           
248               << "        a direction vector o    
249               << "        areacode = " << area    
250       G4Exception("G4VTwistSurface::DistanceTo    
251                   FatalException, message);       
252    }                                              
253    else if (IsAxis0(areacode) || IsAxis1(areac    
254    {                                              
255       GetBoundaryParameters(areacode, d, x0, b    
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." <    
273               << "        areacode = " << area    
274       G4Exception("G4VTwistSurface::DistanceTo    
275                   FatalException, message);       
276    }                                              
277    return dist;                                   
278 }                                                 
279                                                   
280 //============================================    
281 //* DistanceToIn -----------------------------    
282                                                   
283 G4double G4VTwistSurface::DistanceToIn(const G    
284                                        const G    
285                                              G    
286 {                                                 
287 #ifdef G4TWISTDEBUG                               
288    G4cout << " ~~~~ G4VTwistSurface::DistanceT    
289    G4cout << "      Name : " << fName << G4end    
290    G4cout << "      gp   : " << gp << G4endl;     
291    G4cout << "      gv   : " <<  gv << G4endl;    
292    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
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,    
312                                                   
313    G4int         nxx = DistanceToSurface(gp, g    
314                                          isval    
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] <    
326          continue;                                
327       }                                           
328                                                   
329       G4ThreeVector normal = GetNormal(gxx[i],    
330                                                   
331       if ((normal * gv) >= 0)                     
332       {                                           
333                                                   
334 #ifdef G4TWISTDEBUG                               
335          G4cout << "   G4VTwistSurface::Distan    
336                 << "particle goes outword the     
337 #endif                                            
338          continue;                                
339       }                                           
340                                                   
341       //                                          
342       // accept this intersection if the inter    
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::Dis    
354                    << " areacode sInside name,    
355                    << fName <<  " "<< bestdist    
356 #endif                                            
357          }                                        
358                                                   
359       //                                          
360       // else, the intersection is on boundary    
361       //                                          
362                                                   
363       }                                           
364       else                                        
365       {                                           
366          G4VTwistSurface* neighbours[2];          
367          G4bool      isaccepted[2] = {false, f    
368          G4int       nneighbours   = GetNeighb    
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[G4VSURFACENX    
377             G4int         tmpareacode[G4VSURFA    
378             G4bool        tmpisvalid[G4VSURFAC    
379                                                   
380             for (G4int l = 0 ; l<G4VSURFACENXX    
381             {                                     
382               tmpdist[l] = kInfinity ;            
383               tmpareacode[l] = sOutside ;         
384               tmpisvalid[l] = false ;             
385             }                                     
386                                                   
387             G4int tmpnxx = neighbours[j]->Dist    
388                                           gp,     
389                                           tmpa    
390                                           kVal    
391             G4ThreeVector neighbournormal;        
392                                                   
393             for (G4int k=0; k< tmpnxx; ++k)       
394             {                                     
395                //                                 
396                // if tmpxx[k] is valid && sIns    
397                // be neighbour surface. return    
398                // else , choose tmpxx on same     
399                //                                 
400                                                   
401                if (IsInside(tmpareacode[k]))      
402                {                                  
403 #ifdef G4TWISTDEBUG                               
404                   G4cout << "   G4VTwistSurfac    
405                          << " intersection "<<    
406                          << "   is inside of n    
407                          << " . returning kInf    
408                   G4cout << "~~ G4VTwistSurfac    
409                          << G4endl;               
410                   G4cout << "      No intersec    
411                   G4cout << "      Name : " <<    
412                   G4cout << "~~~~~~~~~~~~~~~~~    
413                          << G4endl;               
414 #endif                                            
415                   if (tmpisvalid[k])  return k    
416                   continue;                       
417                                                   
418                //                                 
419                // if tmpxx[k] is valid && sIns    
420                // be neighbour surface. return    
421                //                                 
422                                                   
423                }                                  
424                else if (IsSameBoundary(this,ar    
425                                        neighbo    
426                {                                  
427                   // tmpxx[k] is same boundary    
428                                                   
429                   neighbournormal = neighbours    
430                   if (neighbournormal * gv < 0    
431                }                                  
432             }                                     
433                                                   
434             // if nneighbours = 1, chabge isac    
435             // exiting neighboursurface loop.     
436                                                   
437             if (nneighbours == 1) isaccepted[1    
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::    
452                       << " areacode sBoundary     
453                       << fName  << " " << dist    
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::Distance    
466       G4cout << "      No intersections " << G    
467       G4cout << "      Name : " << fName << G4    
468       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
469    }                                              
470    else                                           
471    {                                              
472       G4cout << "~~~ G4VTwistSurface::Distance    
473       G4cout << "      Name, i  : " << fName <    
474       G4cout << "      gxx[i]   : " << gxxbest    
475       G4cout << "      bestdist : " << bestdis    
476       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
477    }                                              
478                                                   
479 #endif                                            
480                                                   
481    return bestdistance;                           
482 }                                                 
483                                                   
484 //============================================    
485 //* DistanceToOut(p, v) ----------------------    
486                                                   
487 G4double G4VTwistSurface::DistanceToOut(const     
488                                         const     
489                                                   
490 {                                                 
491 #ifdef G4TWISTDEBUG                               
492    G4cout << "~~~~~ G4VTwistSurface::DistanceT    
493    G4cout << "      Name : " << fName << G4end    
494    G4cout << "      gp   : " << gp << G4endl;     
495    G4cout << "      gv   : " <<  gv << G4endl;    
496    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
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, distan    
515                            isvalid, kValidateW    
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],    
525       if (normal * gv <= 0)                       
526       {                                           
527          // particle goes toword inside of sol    
528 #ifdef G4TWISTDEBUG                               
529           G4cout << "   G4VTwistSurface::Dista    
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::DistanceT    
549       G4cout << "      No intersections   " <<    
550       G4cout << "      Name     : " << fName <    
551       G4cout << "      bestdist : " << bestdis    
552       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
553    }                                              
554    else                                           
555    {                                              
556       G4cout << "~~ G4VTwistSurface::DistanceT    
557       G4cout << "      Name, i  : " << fName <    
558       G4cout << "      gxx[i]   : " << gxxbest    
559       G4cout << "      bestdist : " << bestdis    
560       G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
561    }                                              
562 #endif                                            
563                                                   
564    return bestdistance;                           
565 }                                                 
566                                                   
567 //============================================    
568 //* DistanceTo(p) ----------------------------    
569                                                   
570 G4double G4VTwistSurface::DistanceTo(const G4T    
571                                            G4T    
572 {                                                 
573 #ifdef G4TWISTDEBUG                               
574    G4cout << "~~~~~ G4VTwistSurface::DistanceT    
575    G4cout << "      Name : " << fName << G4end    
576    G4cout << "      gp   : " << gp << G4endl;     
577    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
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, areaco    
592    gxxbest = gxx[0];                              
593                                                   
594 #ifdef G4TWISTDEBUG                               
595    G4cout << "~~~~~ G4VTwistSurface::DistanceT    
596    G4cout << "      Name     : " << fName << G    
597    G4cout << "      gxx      : " << gxxbest <<    
598    G4cout << "      bestdist : " << distance[0    
599    G4cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
600 #endif                                            
601                                                   
602    return distance[0];                            
603 }                                                 
604                                                   
605 //============================================    
606 //* IsSameBoundary ---------------------------    
607                                                   
608 G4bool                                            
609 G4VTwistSurface::IsSameBoundary(G4VTwistSurfac    
610                                 G4VTwistSurfac    
611 {                                                 
612    //                                             
613    // IsSameBoundary                              
614    //                                             
615    // checking tool whether two boundaries on     
616    //                                             
617                                                   
618    G4bool testbitmode = true;                     
619    G4bool iscorner[2] = {IsCorner(areacode1, t    
620                          IsCorner(areacode2, t    
621                                                   
622    if (iscorner[0] && iscorner[1])                
623    {                                              
624       // on corner                                
625       G4ThreeVector corner1 =                     
626            surf1->ComputeGlobalPoint(surf1->Ge    
627       G4ThreeVector corner2 =                     
628            surf2->ComputeGlobalPoint(surf2->Ge    
629                                                   
630       return (corner1 - corner2).mag() < kCarT    
631    }                                              
632    else if ((IsBoundary(areacode1, testbitmode    
633             (IsBoundary(areacode2, testbitmode    
634    {                                              
635       // on boundary                              
636       G4ThreeVector d1, d2, ld1, ld2;             
637       G4ThreeVector x01, x02, lx01, lx02;         
638       G4int         type1, type2;                 
639       surf1->GetBoundaryParameters(areacode1,     
640       surf2->GetBoundaryParameters(areacode2,     
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(co    
660                                              G    
661                                              G    
662                                              G    
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(areac    
671       {                                           
672          return;                                  
673       }                                           
674    }                                              
675                                                   
676    std::ostringstream message;                    
677    message << "Not registered boundary." << G4    
678            << "        Boundary at areacode "     
679            << std::dec << G4endl                  
680            << "        is not registered.";       
681    G4Exception("G4VTwistSurface::GetBoundaryPa    
682                FatalException, message);          
683 }                                                 
684                                                   
685 //============================================    
686 //* GetBoundaryAtPZ --------------------------    
687                                                   
688 G4ThreeVector G4VTwistSurface::GetBoundaryAtPZ    
689                                                   
690 {                                                 
691    // areacode must be one of them:               
692    // sAxis0 & sAxisMin, sAxis0 & sAxisMax,       
693    // sAxis1 & sAxisMin, sAxis1 & sAxisMax.       
694                                                   
695    if (((areacode & sAxis0) != 0) && ((areacod    
696    {                                              
697      std::ostringstream message;                  
698      message << "Point is in the corner area."    
699              << "        This function returns    
700              << "a direction vector of a bound    
701              << "        areacode = " << areac    
702      G4Exception("G4VTwistSurface::GetBoundary    
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(areac    
714       {                                           
715          found = true;                            
716          continue;                                
717       }                                           
718    }                                              
719                                                   
720    if (!found)                                    
721    {                                              
722      std::ostringstream message;                  
723      message << "Not registered boundary." <<     
724              << "        Boundary at areacode     
725              << "        is not registered.";     
726      G4Exception("G4VTwistSurface::GetBoundary    
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 boundar    
735              << "        Boundary at areacode     
736              << "        is not a z-depended l    
737      G4Exception("G4VTwistSurface::GetBoundary    
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, G4    
748 {                                                 
749    if ((areacode & sCorner) != sCorner)           
750    {                                              
751      std::ostringstream message;                  
752      message << "Area code must represents cor    
753              << "        areacode " << areacod    
754      G4Exception("G4VTwistSurface::SetCorner()    
755                  FatalException, message);        
756    }                                              
757                                                   
758    if ((areacode & sC0Min1Min) == sC0Min1Min)     
759       fCorners[0].set(x, y, z);                   
760    } else if ((areacode & sC0Max1Min) == sC0Ma    
761       fCorners[1].set(x, y, z);                   
762    } else if ((areacode & sC0Max1Max) == sC0Ma    
763       fCorners[2].set(x, y, z);                   
764    } else if ((areacode & sC0Min1Max) == sC0Mi    
765       fCorners[3].set(x, y, z);                   
766    }                                              
767 }                                                 
768                                                   
769 //============================================    
770 //* SetBoundaryAxis --------------------------    
771                                                   
772 void G4VTwistSurface::GetBoundaryAxis(G4int ar    
773 {                                                 
774    if ((areacode & sBoundary) != sBoundary) {     
775      G4Exception("G4VTwistSurface::GetBoundary    
776                  FatalException, "Not located     
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 &    
789       if (axiscode != 0) {                        
790          if (axiscode == (whichaxis & sAxisX))    
791             axis[i] = kXAxis;                     
792          } else if (axiscode == (whichaxis & s    
793             axis[i] = kYAxis;                     
794          } else if (axiscode == (whichaxis & s    
795             axis[i] = kZAxis;                     
796          } else if (axiscode == (whichaxis & s    
797             axis[i] = kRho;                       
798          } else if (axiscode == (whichaxis & s    
799             axis[i] = kPhi;                       
800          } else {                                 
801            std::ostringstream message;            
802            message << "Not supported areacode.    
803                    << "        areacode " << a    
804            G4Exception("G4VTwistSurface::GetBo    
805                        FatalException, message    
806          }                                        
807       }                                           
808    }                                              
809 }                                                 
810                                                   
811 //============================================    
812 //* SetBoundaryLimit -------------------------    
813                                                   
814 void G4VTwistSurface::GetBoundaryLimit(G4int a    
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)) !=     
832          limit[0] = fAxisMin[0];                  
833       } else if ((areacode & (sAxis1 | sAxisMi    
834          limit[0] = fAxisMin[1];                  
835       } else if ((areacode & (sAxis0 | sAxisMa    
836          limit[0] = fAxisMax[0];                  
837       } else if ((areacode & (sAxis1 | sAxisMa    
838          limit[0] = fAxisMax[1];                  
839       }                                           
840    } else {                                       
841      std::ostringstream message;                  
842      message << "Not located on a boundary!" <    
843              << "          areacode " << areac    
844      G4Exception("G4VTwistSurface::GetBoundary    
845                  JustWarning, message);           
846    }                                              
847 }                                                 
848                                                   
849 //============================================    
850 //* SetBoundary ------------------------------    
851                                                   
852 void G4VTwistSurface::SetBoundary(const G4int&    
853                                   const G4Thre    
854                                   const G4Thre    
855                                   const G4int&    
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, direc    
869             done = true;                          
870             break;                                
871          }                                        
872       }                                           
873                                                   
874       if (!done)                                  
875       {                                           
876          G4Exception("G4VTwistSurface::SetBoun    
877                       FatalException, "Number     
878       }                                           
879    }                                              
880    else                                           
881    {                                              
882       std::ostringstream message;                 
883       message << "Invalid axis-code." << G4end    
884               << "        axiscode = "            
885               << std::hex << axiscode << std::    
886       G4Exception("G4VTwistSurface::SetBoundar    
887                   FatalException, message);       
888    }                                              
889 }                                                 
890                                                   
891 //============================================    
892 //* GetFace ----------------------------------    
893                                                   
894 G4int G4VTwistSurface::GetFace( G4int i, G4int    
895                                 G4int n, G4int    
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-    
914   }                                               
915                                                   
916   else if ( iside == 4 ) {                        
917     return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*(    
918   }                                               
919                                                   
920   else if ( iside == 5 ) {                        
921     return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*(    
922   }                                               
923                                                   
924   else {                                          
925     std::ostringstream message;                   
926     message << "Not correct side number: "        
927             << GetName() << G4endl                
928             << "iside is " << iside << " but s    
929             << "0,1,2,3,4 or 5" << ".";           
930     G4Exception("G4TwistSurface::G4GetFace()",    
931                 FatalException, message);         
932   }                                               
933                                                   
934   return -1 ;  // wrong face                      
935 }                                                 
936                                                   
937 //============================================    
938 //* GetNode ----------------------------------    
939                                                   
940 G4int G4VTwistSurface::GetNode( G4int i, G4int    
941                                 G4int n, G4int    
942 {                                                 
943   // this is the node mapping function            
944   // (i,j) -> node number                         
945   // Depends on the side iside and the used me    
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    
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    
967   }                                               
968                                                   
969   else if ( iside == 3 )                          
970   {                                               
971     // right side                                 
972     if      ( i == 0 )   { return       (j+1)*    
973     else if ( i == n-1 ) { return k*k + (j+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 -     
983     else if ( i == n-1 ) { return 2*k*k - 1 -     
984     else                                          
985     {                                             
986       return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) +    
987     }                                             
988   }                                               
989   else if ( iside == 5 )                          
990   {                                               
991     // left side                                  
992     if      ( i == 0 )   { return k*k   - (j+1    
993     else if ( i == n-1)  { return 2*k*k - (j+1    
994     else                                          
995     {                                             
996       if ( j == k-1 ) { return 2*k*k + 4*(i-1)    
997       else                                        
998       {                                           
999         return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1)    
1000       }                                          
1001     }                                            
1002   }                                              
1003   else                                           
1004   {                                              
1005     std::ostringstream message;                  
1006     message << "Not correct side number: "       
1007             << GetName() << G4endl               
1008             << "iside is " << iside << " but     
1009             << "0,1,2,3,4 or 5" << ".";          
1010     G4Exception("G4TwistSurface::G4GetNode()"    
1011                 FatalException, message);        
1012   }                                              
1013   return -1 ;  // wrong node                     
1014 }                                                
1015                                                  
1016 //===========================================    
1017 //* GetEdgeVisiblility ----------------------    
1018                                                  
1019 G4int G4VTwistSurface::GetEdgeVisibility( G4i    
1020                                           G4i    
1021 {                                                
1022   // clockwise filling         -> positive or    
1023   // counter clockwise filling -> negative or    
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     
1052   // 0,1,2,3  -> 3,2,1,0                         
1053   if ( orientation < 0 ) { number = ( 3 - num    
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: "     
1070       G4Exception("G4TwistSurface::G4GetEdgeV    
1071                   "GeomSolids0003", FatalExce    
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: "     
1089       G4Exception("G4TwistSurface::G4GetEdgeV    
1090                   "GeomSolids0003", FatalExce    
1091     }                                            
1092   }                                              
1093                                                  
1094   // now the corners                             
1095   if ( i == 0 && j == 0 ) {          // signs    
1096     return ( number == 0 || number == 3 ) ? 1    
1097   }                                              
1098   else if ( i == 0 && j == k-2 ) {   // signs    
1099     return ( number == 2 || number == 3 ) ? 1    
1100   }                                              
1101   else if ( i == n-2 && j == k-2 ) { // signs    
1102     return ( number == 1 || number == 2 ) ? 1    
1103   }                                              
1104   else if ( i == n-2 && j == 0 ) {   // signs    
1105     return ( number == 0 || number == 1 ) ? 1    
1106   }                                              
1107   else                                           
1108   {                                              
1109     std::ostringstream message;                  
1110     message << "Not correct face number: " <<    
1111     G4Exception("G4TwistSurface::G4GetEdgeVis    
1112                 "GeomSolids0003", FatalExcept    
1113   }                                              
1114                                                  
1115   std::ostringstream message;                    
1116   message << "Not correct face number: " << G    
1117   G4Exception("G4TwistSurface::G4GetEdgeVisib    
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(sC0Min1    
1130    G4ThreeVector B = fRot * GetCorner(sC0Max1    
1131    G4ThreeVector C = fRot * GetCorner(sC0Max1    
1132    G4ThreeVector D = fRot * GetCorner(sC0Min1    
1133                                                  
1134    G4cout << "/* G4VTwistSurface::DebugPrint(    
1135           << G4endl;                             
1136    G4cout << "/* Name = " << fName << G4endl;    
1137    G4cout << "/* Axis = " << std::hex << fAxi    
1138           << std::hex << fAxis[1]                
1139           << " (0,1,2,3,5 = kXAxis,kYAxis,kZA    
1140           << std::dec << G4endl;                 
1141    G4cout << "/* BoundaryLimit(in local) fAxi    
1142           << ", " << fAxisMax[0] << ")" << G4    
1143    G4cout << "/* BoundaryLimit(in local) fAxi    
1144           << ", " << fAxisMax[1] << ")" << G4    
1145    G4cout << "/* Cornar point sC0Min1Min = "     
1146    G4cout << "/* Cornar point sC0Max1Min = "     
1147    G4cout << "/* Cornar point sC0Max1Max = "     
1148    G4cout << "/* Cornar point sC0Min1Max = "     
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, kInfinit    
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::~CurrentStatu    
1180 = default;                                       
1181                                                  
1182 //===========================================    
1183 //* CurrentStatus::SetCurrentStatus ---------    
1184                                                  
1185 void                                             
1186 G4VTwistSurface::CurrentStatus::SetCurrentSta    
1187                                                  
1188                                                  
1189                                                  
1190                                                  
1191                                                  
1192                                                  
1193                                            co    
1194                                            co    
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::CurrentStat    
1209                 "GeomSolids0003", FatalExcept    
1210   }                                              
1211   if (v != nullptr)                              
1212   {                                              
1213     fLastv = *v;                                 
1214   }                                              
1215   else                                           
1216   {                                              
1217     fLastv.set(kInfinity, kInfinity, kInfinit    
1218   }                                              
1219   fDone = true;                                  
1220 }                                                
1221                                                  
1222 //===========================================    
1223 //* CurrentStatus::ResetfDone ---------------    
1224                                                  
1225 void                                             
1226 G4VTwistSurface::CurrentStatus::ResetfDone(EV    
1227                                      const G4    
1228                                      const G4    
1229                                                  
1230 {                                                
1231   if (validate == fLastValidate && p != nullp    
1232   {                                              
1233      if (v == nullptr || (*v == fLastv)) retu    
1234   }                                              
1235   G4ThreeVector xx(kInfinity, kInfinity, kInf    
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 f    
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()     
1255 {                                                
1256   G4cout << "CurrentStatus::Dist0,1= " << fDi    
1257          << " " << fDistance[1] << " areacode    
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 G4    
1270                                 const G4Three    
1271                                 const G4Three    
1272                                 const G4int&     
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() c    
1284 {                                                
1285   return fBoundaryAcode == -1;                   
1286 }                                                
1287                                                  
1288 //===========================================    
1289 //* Boundary::GetBoundaryParameters ---------    
1290                                                  
1291 G4bool                                           
1292 G4VTwistSurface::Boundary::GetBoundaryParamet    
1293                                                  
1294                                                  
1295                                                  
1296 {                                                
1297   // areacode must be one of them:               
1298   // sAxis0 & sAxisMin, sAxis0 & sAxisMax,       
1299   // sAxis1 & sAxisMin, sAxis1 & sAxisMax        
1300   //                                             
1301   if (((areacode & sAxis0) != 0) && ((areacod    
1302   {                                              
1303     std::ostringstream message;                  
1304     message << "Located in the corner area."     
1305             << "        This function returns    
1306             << "a boundary line." << G4endl      
1307             << "        areacode = " << areac    
1308     G4Exception("G4VTwistSurface::Boundary::G    
1309                 "GeomSolids0003", FatalExcept    
1310   }                                              
1311   if ((areacode & sSizeMask) != (fBoundaryAco    
1312   {                                              
1313     return false;                                
1314   }                                              
1315   d  = fBoundaryDirection;                       
1316   x0 = fBoundaryX0;                              
1317   boundarytype = fBoundaryType;                  
1318   return true;                                   
1319 }                                                
1320