source: trunk/source/geometry/solids/CSG/src/G4Tubs.cc @ 850

Last change on this file since 850 was 850, checked in by garnier, 16 years ago

geant4.8.2 beta
File size: 57.0 KB
Line 
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//
27// $Id: G4Tubs.cc,v 1.68 2008/06/23 13:37:39 gcosmo Exp $
28// GEANT4 tag $Name: HEAD $
29//
30//
31// class G4Tubs
32//
33// History:
34//
35// 02.08.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for negative value under sqrt
36//                      for the case: p on the surface and v is tangent to the surface
37// 11.05.07 T.Nikitina: bug fixed in DistanceToOut(p,v,..) for phi < 2pi
38// 03.05.05 V.Grichine: SurfaceNormal(p) according to J. Apostolakis proposal
39// 16.03.05 V.Grichine: SurfaceNormal(p) with edges/corners for boolean
40// 20.07.01 V.Grichine: bug fixed in Inside(p)
41// 20.02.01 V.Grichine: bug fixed in Inside(p) and CalculateExtent was
42//                      simplified base on G4Box::CalculateExtent
43// 07.12.00 V.Grichine: phi-section algorithm was changed in Inside(p)
44// 28.11.00 V.Grichine: bug fixed in Inside(p)
45// 31.10.00 V.Grichine: assign sr, sphi in Distance ToOut(p,v,...)
46// 08.08.00 V.Grichine: more stable roots of 2-equation in DistanceToOut(p,v,..)
47// 02.08.00 V.Grichine: point is outside check in Distance ToOut(p)
48// 17.05.00 V.Grichine: bugs (#76,#91) fixed in Distance ToOut(p,v,...)
49// 31.03.00 V.Grichine: bug fixed in Inside(p)
50// 19.11.99 V.Grichine: side = kNull in DistanceToOut(p,v,...)
51// 13.10.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
52// 28.05.99 V.Grichine: bugs fixed in DistanceToOut(p,v,...)
53// 25.05.99 V.Grichine: bugs fixed in DistanceToIn(p,v)
54// 23.03.99 V.Grichine: bug fixed in DistanceToIn(p,v)
55// 09.10.98 V.Grichine: modifications in DistanceToOut(p,v,...)
56// 18.06.98 V.Grichine: n-normalisation in DistanceToOut(p,v)
57//
58// 1994-95  P.Kent:     implementation
59//
60/////////////////////////////////////////////////////////////////////////
61
62#include "G4Tubs.hh"
63
64#include "G4VoxelLimits.hh"
65#include "G4AffineTransform.hh"
66#include "G4GeometryTolerance.hh"
67
68#include "G4VPVParameterisation.hh"
69
70#include "Randomize.hh"
71
72#include "meshdefs.hh"
73
74#include "G4VGraphicsScene.hh"
75#include "G4Polyhedron.hh"
76#include "G4NURBS.hh"
77#include "G4NURBStube.hh"
78#include "G4NURBScylinder.hh"
79#include "G4NURBStubesector.hh"
80
81using namespace CLHEP;
82
83/////////////////////////////////////////////////////////////////////////
84//
85// Constructor - check parameters, convert angles so 0<sphi+dpshi<=2_PI
86//             - note if pdphi>2PI then reset to 2PI
87
88G4Tubs::G4Tubs( const G4String &pName,
89                      G4double pRMin, G4double pRMax,
90                      G4double pDz,
91                      G4double pSPhi, G4double pDPhi )
92  : G4CSGSolid(pName)
93{
94
95  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
96  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
97
98  if (pDz>0) // Check z-len
99  {
100    fDz = pDz ;
101  }
102  else
103  {
104    G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
105           << "        Negative Z half-length ! - "
106           << pDz << G4endl;
107    G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
108                "Invalid Z half-length");
109  }
110  if ( pRMin < pRMax && pRMin >= 0 ) // Check radii
111  {
112    fRMin = pRMin ; 
113    fRMax = pRMax ;
114  }
115  else
116  {
117    G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
118           << "        Invalid values for radii !" << G4endl
119           << "        pRMin = " << pRMin << ", pRMax = " << pRMax << G4endl;
120    G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
121                "Invalid radii.");
122  }
123  if ( pDPhi >= twopi ) // Check angles
124  {
125    fDPhi=twopi;
126  }
127  else
128  {
129    if ( pDPhi > 0 )
130    {
131      fDPhi = pDPhi;
132    }
133    else
134    {
135      G4cerr << "ERROR - G4Tubs()::G4Tubs(): " << GetName() << G4endl
136             << "        Negative delta-Phi ! - "
137             << pDPhi << G4endl;
138      G4Exception("G4Tubs::G4Tubs()", "InvalidSetup", FatalException,
139                  "Invalid dphi.");
140    }
141  }
142 
143  // Ensure fSphi in 0-2PI or -2PI-0 range if shape crosses 0
144
145  fSPhi = pSPhi;
146
147  if ( fSPhi < 0 )
148  {
149    fSPhi = twopi - std::fmod(std::fabs(fSPhi),twopi) ;
150  }
151  else
152  {
153    fSPhi = std::fmod(fSPhi,twopi) ;
154  }
155  if (fSPhi + fDPhi > twopi )
156  {
157    fSPhi -= twopi ;
158  }
159}
160
161///////////////////////////////////////////////////////////////////////
162//
163// Fake default constructor - sets only member data and allocates memory
164//                            for usage restricted to object persistency.
165//
166G4Tubs::G4Tubs( __void__& a )
167  : G4CSGSolid(a)
168{
169}
170
171//////////////////////////////////////////////////////////////////////////
172//
173// Destructor
174
175G4Tubs::~G4Tubs()
176{
177}
178
179/////////////////////////////////////////////////////////////////////////
180//
181// Dispatch to parameterisation for replication mechanism dimension
182// computation & modification.
183
184void G4Tubs::ComputeDimensions(       G4VPVParameterisation* p,
185                                const G4int n,
186                                const G4VPhysicalVolume* pRep )
187{
188  p->ComputeDimensions(*this,n,pRep) ;
189}
190
191////////////////////////////////////////////////////////////////////////
192//
193// Calculate extent under transform and specified limit
194
195G4bool G4Tubs::CalculateExtent( const EAxis              pAxis,
196                                const G4VoxelLimits&     pVoxelLimit,
197                                const G4AffineTransform& pTransform,
198                                      G4double&          pMin, 
199                                      G4double&          pMax    ) const
200{
201
202  if ( !pTransform.IsRotated() && fDPhi == twopi && fRMin == 0 )
203  {
204    // Special case handling for unrotated solid tubes
205    // Compute x/y/z mins and maxs fro bounding box respecting limits,
206    // with early returns if outside limits. Then switch() on pAxis,
207    // and compute exact x and y limit for x/y case
208     
209    G4double xoffset, xMin, xMax ;
210    G4double yoffset, yMin, yMax ;
211    G4double zoffset, zMin, zMax ;
212
213    G4double diff1, diff2, maxDiff, newMin, newMax ;
214    G4double xoff1, xoff2, yoff1, yoff2, delta ;
215
216    xoffset = pTransform.NetTranslation().x() ;
217    xMin = xoffset - fRMax ;
218    xMax = xoffset + fRMax ;
219
220    if (pVoxelLimit.IsXLimited())
221    {
222      if ( (xMin > pVoxelLimit.GetMaxXExtent())
223        || (xMax < pVoxelLimit.GetMinXExtent()) )
224      {
225        return false;
226      }
227      else
228      {
229        if ( xMin < pVoxelLimit.GetMinXExtent() )
230        {
231          xMin = pVoxelLimit.GetMinXExtent() ;
232        }
233        if (xMax > pVoxelLimit.GetMaxXExtent() )
234        {
235          xMax = pVoxelLimit.GetMaxXExtent() ;
236        }
237      }
238    }
239    yoffset = pTransform.NetTranslation().y() ;
240    yMin    = yoffset - fRMax ;
241    yMax    = yoffset + fRMax ;
242
243    if ( pVoxelLimit.IsYLimited() )
244    {
245      if ( (yMin > pVoxelLimit.GetMaxYExtent())
246        || (yMax < pVoxelLimit.GetMinYExtent()) )
247      {
248        return false ;
249      }
250      else
251      {
252        if ( yMin < pVoxelLimit.GetMinYExtent() )
253        {
254          yMin = pVoxelLimit.GetMinYExtent() ;
255        }
256        if ( yMax > pVoxelLimit.GetMaxYExtent() )
257        {
258          yMax=pVoxelLimit.GetMaxYExtent();
259        }
260      }
261    }
262    zoffset = pTransform.NetTranslation().z() ;
263    zMin    = zoffset - fDz ;
264    zMax    = zoffset + fDz ;
265
266    if ( pVoxelLimit.IsZLimited() )
267    {
268      if ( (zMin > pVoxelLimit.GetMaxZExtent())
269        || (zMax < pVoxelLimit.GetMinZExtent()) )
270      {
271        return false ;
272      }
273      else
274      {
275        if ( zMin < pVoxelLimit.GetMinZExtent() )
276        {
277          zMin = pVoxelLimit.GetMinZExtent() ;
278        }
279        if ( zMax > pVoxelLimit.GetMaxZExtent() )
280        {
281          zMax = pVoxelLimit.GetMaxZExtent();
282        }
283      }
284    }
285    switch ( pAxis )  // Known to cut cylinder
286    {
287      case kXAxis :
288      {
289        yoff1 = yoffset - yMin ;
290        yoff2 = yMax    - yoffset ;
291
292        if ( yoff1 >= 0 && yoff2 >= 0 ) // Y limits cross max/min x => no change
293        {
294          pMin = xMin ;
295          pMax = xMax ;
296        }
297        else
298        {
299          // Y limits don't cross max/min x => compute max delta x,
300          // hence new mins/maxs
301
302          delta   = fRMax*fRMax - yoff1*yoff1;
303          diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
304          delta   = fRMax*fRMax - yoff2*yoff2;
305          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
306          maxDiff = (diff1 > diff2) ? diff1:diff2;
307          newMin  = xoffset - maxDiff;
308          newMax  = xoffset + maxDiff;
309          pMin    = (newMin < xMin) ? xMin : newMin;
310          pMax    = (newMax > xMax) ? xMax : newMax;
311        }   
312        break;
313      }
314      case kYAxis :
315      {
316        xoff1 = xoffset - xMin ;
317        xoff2 = xMax - xoffset ;
318
319        if ( xoff1 >= 0 && xoff2 >= 0 ) // X limits cross max/min y => no change
320        {
321          pMin = yMin ;
322          pMax = yMax ;
323        }
324        else
325        {
326          // X limits don't cross max/min y => compute max delta y,
327          // hence new mins/maxs
328
329          delta   = fRMax*fRMax - xoff1*xoff1;
330          diff1   = (delta>0.) ? std::sqrt(delta) : 0.;
331          delta   = fRMax*fRMax - xoff2*xoff2;
332          diff2   = (delta>0.) ? std::sqrt(delta) : 0.;
333          maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
334          newMin  = yoffset - maxDiff ;
335          newMax  = yoffset + maxDiff ;
336          pMin    = (newMin < yMin) ? yMin : newMin ;
337          pMax     =(newMax > yMax) ? yMax : newMax ;
338        }
339        break ;
340      }
341      case kZAxis:
342      {
343        pMin = zMin ;
344        pMax = zMax ;
345        break ;
346      }
347      default:
348        break;
349    }
350    pMin -= kCarTolerance ;
351    pMax += kCarTolerance ;
352    return true;   
353  }
354  else // Calculate rotated vertex coordinates
355  {
356    G4int i, noEntries, noBetweenSections4 ;
357    G4bool existsAfterClip = false ;
358    G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;
359   
360    pMin = +kInfinity ;
361    pMax = -kInfinity ;
362
363    noEntries = vertices->size() ;
364    noBetweenSections4 = noEntries - 4 ;
365    /*
366    G4cout << "vertices = " << noEntries << "\t"
367           << "v-4 = " << noBetweenSections4 << G4endl;
368    G4cout << G4endl;
369    for(i = 0 ; i < noEntries ; i++ )
370    {
371      G4cout << i << "\t" << "v.x = " << ((*vertices)[i]).x() << "\t"
372                          << "v.y = " << ((*vertices)[i]).y() << "\t"
373                          << "v.z = " << ((*vertices)[i]).z() << "\t" << G4endl;
374    }     
375    G4cout << G4endl;
376    G4cout << "ClipCrossSection" << G4endl;
377    */
378    for (i = 0 ; i < noEntries ; i += 4 )
379    {
380      // G4cout << "section = " << i << G4endl;
381      ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
382    }
383    // G4cout << "ClipBetweenSections" << G4endl;
384    for (i = 0 ; i < noBetweenSections4 ; i += 4 )
385    {
386      // G4cout << "between sections = " << i << G4endl;
387      ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax) ;
388    }
389    if (pMin != kInfinity || pMax != -kInfinity )
390    {
391      existsAfterClip = true ;
392      pMin -= kCarTolerance ; // Add 2*tolerance to avoid precision troubles
393      pMax += kCarTolerance ;
394    }
395    else
396    {
397      // Check for case where completely enveloping clipping volume
398      // If point inside then we are confident that the solid completely
399      // envelopes the clipping volume. Hence set min/max extents according
400      // to clipping volume extents along the specified axis.
401
402      G4ThreeVector clipCentre(
403             (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
404             (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
405             (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
406       
407      if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
408      {
409        existsAfterClip = true ;
410        pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
411        pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
412      }
413    }
414    delete vertices;
415    return existsAfterClip;
416  }
417}
418
419
420///////////////////////////////////////////////////////////////////////////
421//
422// Return whether point inside/outside/on surface
423
424EInside G4Tubs::Inside( const G4ThreeVector& p ) const
425{
426  G4double r2,pPhi,tolRMin,tolRMax;
427  EInside in = kOutside ;
428 
429  if (std::fabs(p.z()) <= fDz - kCarTolerance*0.5)
430  {
431    r2 = p.x()*p.x() + p.y()*p.y() ;
432
433    if (fRMin) { tolRMin = fRMin + kRadTolerance*0.5 ; }
434    else       { tolRMin = 0 ; }
435
436    tolRMax = fRMax - kRadTolerance*0.5 ;
437     
438    if (r2 >= tolRMin*tolRMin && r2 <= tolRMax*tolRMax)
439    {
440      if ( fDPhi == twopi )
441      {
442        in = kInside ;
443      }
444      else
445      {
446        // Try inner tolerant phi boundaries (=>inside)
447        // if not inside, try outer tolerant phi boundaries
448
449        pPhi = std::atan2(p.y(),p.x()) ;
450        if ((tolRMin==0)&&(p.x()==0)&&(p.y()==0))
451        {
452          in=kSurface;
453        }
454        else
455        {
456          if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; } // 0<=pPhi<2pi
457
458          if ( fSPhi >= 0 )
459          {
460            if ( (std::abs(pPhi) < kAngTolerance*0.5)
461              && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
462            { 
463              pPhi += twopi ; // 0 <= pPhi < 2pi
464            }
465            if ( (pPhi >= fSPhi + kAngTolerance*0.5)
466              && (pPhi <= fSPhi + fDPhi - kAngTolerance*0.5) )
467            {
468              in = kInside ;
469            }
470            else if ( (pPhi >= fSPhi - kAngTolerance*0.5)
471                   && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
472            {
473              in = kSurface ;
474            }
475          }
476          else  // fSPhi < 0
477          {
478            if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
479              && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
480            else if ( (pPhi <= fSPhi + twopi + kAngTolerance*0.5)
481                   && (pPhi >= fSPhi + fDPhi  - kAngTolerance*0.5) )
482            {
483              in = kSurface ;
484            }
485            else
486            {
487              in = kInside ;
488            }
489          }
490        }                   
491      }
492    }
493    else  // Try generous boundaries
494    {
495      tolRMin = fRMin - kRadTolerance*0.5 ;
496      tolRMax = fRMax + kRadTolerance*0.5 ;
497
498      if ( tolRMin < 0 )  { tolRMin = 0; }
499
500      if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
501      {
502        if ( fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis
503        {
504          in = kSurface ;
505        }
506        else // Try outer tolerant phi boundaries only
507        {
508          pPhi = std::atan2(p.y(),p.x()) ;
509
510          if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; } // 0<=pPhi<2pi
511          if ( fSPhi >= 0 )
512          {
513            if ( (std::abs(pPhi) < kAngTolerance*0.5)
514              && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
515            { 
516              pPhi += twopi ; // 0 <= pPhi < 2pi
517            }
518            if ( (pPhi >= fSPhi - kAngTolerance*0.5)
519              && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
520            {
521              in = kSurface ;
522            }
523          }
524          else  // fSPhi < 0
525          {
526            if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
527              && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
528            else
529            {
530              in = kSurface ;
531            }
532          }
533        }
534      }
535    }
536  }
537  else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5)
538  {                                          // Check within tolerant r limits
539    r2      = p.x()*p.x() + p.y()*p.y() ;
540    tolRMin = fRMin - kRadTolerance*0.5 ;
541    tolRMax = fRMax + kRadTolerance*0.5 ;
542
543    if ( tolRMin < 0 )  { tolRMin = 0; }
544
545    if ( (r2 >= tolRMin*tolRMin) && (r2 <= tolRMax*tolRMax) )
546    {
547      if (fDPhi == twopi || r2 == 0 ) // Continuous in phi or on z-axis
548      {
549        in = kSurface ;
550      }
551      else // Try outer tolerant phi boundaries
552      {
553        pPhi = std::atan2(p.y(),p.x()) ;
554
555        if ( pPhi < -kAngTolerance*0.5 )  { pPhi += twopi; }  // 0<=pPhi<2pi
556        if ( fSPhi >= 0 )
557        {
558          if ( (std::abs(pPhi) < kAngTolerance*0.5)
559            && (std::abs(fSPhi + fDPhi - twopi) < kAngTolerance*0.5) )
560          { 
561            pPhi += twopi ; // 0 <= pPhi < 2pi
562          }
563          if ( (pPhi >= fSPhi - kAngTolerance*0.5)
564            && (pPhi <= fSPhi + fDPhi + kAngTolerance*0.5) )
565          {
566            in = kSurface;
567          }
568        }
569        else  // fSPhi < 0
570        {
571          if ( (pPhi <= fSPhi + twopi - kAngTolerance*0.5)
572            && (pPhi >= fSPhi + fDPhi  + kAngTolerance*0.5) ) {;}
573          else
574          {
575            in = kSurface ;
576          }
577        }     
578      }
579    }
580  }
581  return in;
582}
583
584///////////////////////////////////////////////////////////////////////////
585//
586// Return unit normal of surface closest to p
587// - note if point on z axis, ignore phi divided sides
588// - unsafe if point close to z axis a rmin=0 - no explicit checks
589
590G4ThreeVector G4Tubs::SurfaceNormal( const G4ThreeVector& p ) const
591{ G4int noSurfaces = 0;
592  G4double rho, pPhi;
593  G4double delta   = 0.5*kCarTolerance, dAngle = 0.5*kAngTolerance;
594  G4double distZ, distRMin, distRMax;
595  G4double distSPhi = kInfinity, distEPhi = kInfinity;
596  G4ThreeVector norm, sumnorm(0.,0.,0.);
597  G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0);
598  G4ThreeVector nR, nPs, nPe;
599
600  rho = std::sqrt(p.x()*p.x() + p.y()*p.y());
601
602  distRMin = std::fabs(rho - fRMin);
603  distRMax = std::fabs(rho - fRMax);
604  distZ    = std::fabs(std::fabs(p.z()) - fDz);
605
606  if (fDPhi < twopi)   //  &&  rho ) // Protected against (0,0,z)
607  {
608    if ( rho )
609    {
610      pPhi = std::atan2(p.y(),p.x());
611   
612      if(pPhi  < fSPhi-delta)           { pPhi += twopi; }
613      else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
614
615      distSPhi = std::fabs( pPhi - fSPhi );       
616      distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 
617    }
618    else if( !fRMin )
619    {
620      distSPhi = 0.; 
621      distEPhi = 0.; 
622    }
623    nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
624    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
625  }
626  if ( rho > delta )   { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); }
627
628  if( distRMax <= delta )
629  {
630    noSurfaces ++;
631    sumnorm += nR;
632  }
633  if( fRMin && distRMin <= delta )
634  {
635    noSurfaces ++;
636    sumnorm -= nR;
637  }
638  if( fDPhi < twopi )   
639  {
640    if (distSPhi <= dAngle)  // delta)
641    {
642      noSurfaces ++;
643      sumnorm += nPs;
644    }
645    if (distEPhi <= dAngle) // delta)
646    {
647      noSurfaces ++;
648      sumnorm += nPe;
649    }
650  }
651  if (distZ <= delta) 
652  {
653    noSurfaces ++;
654    if ( p.z() >= 0.)  { sumnorm += nZ; }
655    else               { sumnorm -= nZ; }
656  }
657  if ( noSurfaces == 0 )
658  {
659#ifdef G4CSGDEBUG
660    G4Exception("G4Tube::SurfaceNormal(p)", "Notification",
661                JustWarning, "Point p is not on surface !?" );
662    G4cout.precision(20);
663    G4cout<< "G4Tubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); "
664          << G4endl << G4endl;
665#endif
666     norm = ApproxSurfaceNormal(p);
667  }
668  else if ( noSurfaces == 1 )  { norm = sumnorm; }
669  else                         { norm = sumnorm.unit(); }
670  return norm;
671}
672
673/////////////////////////////////////////////////////////////////////////////
674//
675// Algorithm for SurfaceNormal() following the original specification
676// for points not on the surface
677
678G4ThreeVector G4Tubs::ApproxSurfaceNormal( const G4ThreeVector& p ) const
679{ ENorm side ;
680  G4ThreeVector norm ;
681  G4double rho, phi ;
682  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
683
684  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
685
686  distRMin = std::fabs(rho - fRMin) ;
687  distRMax = std::fabs(rho - fRMax) ;
688  distZ    = std::fabs(std::fabs(p.z()) - fDz) ;
689
690  if (distRMin < distRMax) // First minimum
691  {
692    if ( distZ < distRMin )
693    {
694       distMin = distZ ;
695       side    = kNZ ;
696    }
697    else
698    {
699      distMin = distRMin ;
700      side    = kNRMin   ;
701    }
702  }
703  else
704  {
705    if ( distZ < distRMax )
706    {
707      distMin = distZ ;
708      side    = kNZ   ;
709    }
710    else
711    {
712      distMin = distRMax ;
713      side    = kNRMax   ;
714    }
715  }   
716  if (fDPhi < twopi  &&  rho ) // Protected against (0,0,z)
717  {
718    phi = std::atan2(p.y(),p.x()) ;
719
720    if ( phi < 0 )  { phi += twopi; }
721
722    if ( fSPhi < 0 )
723    {
724      distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ;
725    }
726    else
727    {
728      distSPhi = std::fabs(phi - fSPhi)*rho ;
729    }
730    distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
731                                     
732    if (distSPhi < distEPhi) // Find new minimum
733    {
734      if ( distSPhi < distMin )
735      {
736        side = kNSPhi ;
737      }
738    }
739    else
740    {
741      if ( distEPhi < distMin )
742      {
743        side = kNEPhi ;
744      }
745    }
746  }   
747  switch ( side )
748  {
749    case kNRMin : // Inner radius
750    {                     
751      norm = G4ThreeVector(-p.x()/rho,-p.y()/rho,0) ;
752      break ;
753    }
754    case kNRMax : // Outer radius
755    {                 
756      norm = G4ThreeVector(p.x()/rho,p.y()/rho,0) ;
757      break ;
758    }
759    case kNZ : //    + or - dz
760    {                             
761      if ( p.z() > 0 ) norm = G4ThreeVector(0,0,1)  ; 
762      else             norm = G4ThreeVector(0,0,-1) ; 
763      break ;
764    }
765    case kNSPhi:
766    {
767      norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
768      break ;
769    }
770    case kNEPhi:
771    {
772      norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
773      break;
774    }
775    default:
776    {
777      DumpInfo();
778      G4Exception("G4Tubs::ApproxSurfaceNormal()", "Notification", JustWarning,
779                  "Undefined side for valid surface normal to solid.");
780      break ;
781    }   
782  }               
783  return norm;
784}
785
786////////////////////////////////////////////////////////////////////
787//
788//
789// Calculate distance to shape from outside, along normalised vector
790// - return kInfinity if no intersection, or intersection distance <= tolerance
791//
792// - Compute the intersection with the z planes
793//        - if at valid r, phi, return
794//
795// -> If point is outer outer radius, compute intersection with rmax
796//        - if at valid phi,z return
797//
798// -> Compute intersection with inner radius, taking largest +ve root
799//        - if valid (in z,phi), save intersction
800//
801//    -> If phi segmented, compute intersections with phi half planes
802//        - return smallest of valid phi intersections and
803//          inner radius intersection
804//
805// NOTE:
806// - Precalculations for phi trigonometry are Done `just in time'
807// - `if valid' implies tolerant checking of intersection points
808
809G4double G4Tubs::DistanceToIn( const G4ThreeVector& p,
810                               const G4ThreeVector& v  ) const
811{
812  G4double snxt = kInfinity ;  // snxt = default return value
813
814  // Precalculated trig for phi intersections - used by r,z intersections to
815  //                                            check validity
816
817  G4bool seg ;        // true if segmented
818
819  G4double hDPhi, hDPhiOT, hDPhiIT, cosHDPhiOT=0., cosHDPhiIT=0. ;
820          // half dphi + outer tolerance
821
822  G4double cPhi, sinCPhi=0., cosCPhi=0. ;  // central phi
823
824  G4double tolORMin2, tolIRMax2 ;  // `generous' radii squared
825
826  G4double tolORMax2, tolIRMin2, tolODz, tolIDz ;
827
828  // Intersection point variables
829  //
830  G4double Dist, s, xi, yi, zi, rho2, inum, iden, cosPsi ; 
831
832  G4double t1, t2, t3, b, c, d ;   // Quadratic solver variables
833
834  G4double Comp ;
835  G4double cosSPhi, sinSPhi ;    // Trig for phi start intersect
836
837  G4double ePhi, cosEPhi, sinEPhi ;  // for phi end intersect
838
839  // Set phi divided flag and precalcs
840
841  if ( fDPhi < twopi )
842  {
843    seg        = true ;
844    hDPhi      = 0.5*fDPhi ;    // half delta phi
845    cPhi       = fSPhi + hDPhi ; 
846    hDPhiOT    = hDPhi + 0.5*kAngTolerance ;  // outers tol' half delta phi
847    hDPhiIT    = hDPhi - 0.5*kAngTolerance ;
848    sinCPhi    = std::sin(cPhi) ;
849    cosCPhi    = std::cos(cPhi) ;
850    cosHDPhiOT = std::cos(hDPhiOT) ;
851    cosHDPhiIT = std::cos(hDPhiIT) ;
852  }
853  else
854  {
855    seg = false  ;
856  }
857
858  // Calculate tolerant rmin and rmax
859
860  if (fRMin > kRadTolerance)
861  {
862    tolORMin2 = (fRMin - 0.5*kRadTolerance)*(fRMin - 0.5*kRadTolerance) ;
863    tolIRMin2 = (fRMin + 0.5*kRadTolerance)*(fRMin + 0.5*kRadTolerance) ;
864  }
865  else
866  {
867    tolORMin2 = 0.0 ;
868    tolIRMin2 = 0.0 ;
869  }
870  tolORMax2 = (fRMax + 0.5*kRadTolerance)*(fRMax + 0.5*kRadTolerance) ;
871  tolIRMax2 = (fRMax - 0.5*kRadTolerance)*(fRMax - 0.5*kRadTolerance) ;
872
873  // Intersection with Z surfaces
874
875  tolIDz = fDz - kCarTolerance*0.5 ;
876  tolODz = fDz + kCarTolerance*0.5 ;
877
878  if (std::fabs(p.z()) >= tolIDz)
879  {
880    if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
881    {
882      s = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ;     // Z intersect distance
883
884      if(s < 0.0)  { s = 0.0; }
885
886      xi   = p.x() + s*v.x() ;                // Intersection coords
887      yi   = p.y() + s*v.y() ;
888      rho2 = xi*xi + yi*yi ;
889
890      // Check validity of intersection
891
892      if (tolIRMin2 <= rho2 && rho2 <= tolIRMax2)
893      {
894        if (seg && rho2)
895        {
896          // Psi = angle made with central (average) phi of shape
897          //
898          inum   = xi*cosCPhi + yi*sinCPhi ;
899          iden   = std::sqrt(rho2) ;
900          cosPsi = inum/iden ;
901          if (cosPsi >= cosHDPhiIT) return s ;
902        }
903        else
904        {
905          return s ;
906        }
907      }
908    }
909    else
910    {
911      if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
912      return snxt ;  // On/outside extent, and heading away
913                     // -> cannot intersect
914    }
915  }
916
917  // -> Can not intersect z surfaces
918  //
919  // Intersection with rmax (possible return) and rmin (must also check phi)
920  //
921  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
922  //
923  // Intersects with x^2+y^2=R^2
924  //
925  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
926  //            t1                t2                t3
927
928  t1 = 1.0 - v.z()*v.z() ;
929  t2 = p.x()*v.x() + p.y()*v.y() ;
930  t3 = p.x()*p.x() + p.y()*p.y() ;
931
932  if ( t1 > 0 )        // Check not || to z axis
933  {
934    b = t2/t1 ;
935    c = t3 - fRMax*fRMax ;
936    if (t3 >= tolORMax2 && t2<0)   // This also handles the tangent case
937    {
938      // Try outer cylinder intersection
939      //          c=(t3-fRMax*fRMax)/t1;
940
941      c /= t1 ;
942      d = b*b - c ;
943
944      if (d >= 0)  // If real root
945      {
946        s = -b - std::sqrt(d) ;
947        if (s >= 0)  // If 'forwards'
948        {
949          // Check z intersection
950          //
951          zi = p.z() + s*v.z() ;
952          if (std::fabs(zi)<=tolODz)
953          {
954            // Z ok. Check phi intersection if reqd
955            //
956            if (!seg)
957            {
958              return s ;
959            }
960            else
961            {
962              xi     = p.x() + s*v.x() ;
963              yi     = p.y() + s*v.y() ;
964              cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ;
965              if (cosPsi >= cosHDPhiIT) return s ;
966            }
967          }  //  end if std::fabs(zi)
968        }    //  end if (s>=0)
969      }      //  end if (d>=0)
970    }        //  end if (r>=fRMax)
971    else 
972    {
973      // Inside outer radius :
974      // check not inside, and heading through tubs (-> 0 to in)
975
976      if (t3 > tolIRMin2 && t2 < 0 && std::fabs(p.z()) <= tolIDz)
977      {
978        // Inside both radii, delta r -ve, inside z extent
979
980        if (seg)
981        {
982          inum   = p.x()*cosCPhi + p.y()*sinCPhi ;
983          iden   = std::sqrt(t3) ;
984          cosPsi = inum/iden ;
985          if (cosPsi >= cosHDPhiIT)
986          {
987            // In the old version, the small negative tangent for the point
988            // on surface was not taken in account, and returning 0.0 ...
989            // New version: check the tangent for the point on surface and
990            // if no intersection, return kInfinity, if intersection instead
991            // return s.
992            //
993            c = t3-fRMax*fRMax; 
994            if ( c<=0.0 )
995            {
996              return 0.0;
997            }
998            else
999            {
1000              c = c/t1 ;
1001              d = b*b-c;
1002              if ( d>=0.0 )
1003              {
1004                snxt = c/(-b+std::sqrt(d)); // using safe solution
1005                                            // for quadratic equation
1006                if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
1007                return snxt ;
1008              }     
1009              else
1010              {
1011                return kInfinity;
1012              }
1013            }
1014          } 
1015        }
1016        else
1017        {   
1018          // In the old version, the small negative tangent for the point
1019          // on surface was not taken in account, and returning 0.0 ...
1020          // New version: check the tangent for the point on surface and
1021          // if no intersection, return kInfinity, if intersection instead
1022          // return s.
1023          //
1024          c = t3 - fRMax*fRMax; 
1025          if ( c<=0.0 )
1026          {
1027            return 0.0;
1028          }
1029          else
1030          {
1031            c = c/t1 ;
1032            d = b*b-c;
1033            if ( d>=0.0 )
1034            {
1035              snxt= c/(-b+std::sqrt(d)); // using safe solution
1036                                         // for quadratic equation
1037              if ( snxt<kCarTolerance*0.5 ) { snxt=0; }
1038              return snxt ;
1039            }     
1040            else
1041            {
1042              return kInfinity;
1043            }
1044          }
1045        } // end if   (seg)
1046      }   // end if   (t3>tolIRMin2)
1047    }     // end if   (Inside Outer Radius)
1048    if ( fRMin )    // Try inner cylinder intersection
1049    {
1050      c = (t3 - fRMin*fRMin)/t1 ;
1051      d = b*b - c ;
1052      if ( d >= 0.0 )  // If real root
1053      {
1054        // Always want 2nd root - we are outside and know rmax Hit was bad
1055        // - If on surface of rmin also need farthest root
1056
1057        s = -b + std::sqrt(d) ;
1058        if (s >= -0.5*kCarTolerance)  // check forwards
1059        {
1060          // Check z intersection
1061          //
1062          if(s < 0.0)  { s = 0.0; }
1063          zi = p.z() + s*v.z() ;
1064          if (std::fabs(zi) <= tolODz)
1065          {
1066            // Z ok. Check phi
1067            //
1068            if ( !seg )
1069            {
1070              return s ; 
1071            }
1072            else
1073            {
1074              xi     = p.x() + s*v.x() ;
1075              yi     = p.y() + s*v.y() ;
1076              cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ;
1077              if (cosPsi >= cosHDPhiIT)
1078              {
1079                // Good inner radius isect
1080                // - but earlier phi isect still possible
1081
1082                snxt = s ;
1083              }
1084            }
1085          }        //    end if std::fabs(zi)
1086        }          //    end if (s>=0)
1087      }            //    end if (d>=0)
1088    }              //    end if (fRMin)
1089  }
1090
1091  // Phi segment intersection
1092  //
1093  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
1094  //
1095  // o NOTE: Large duplication of code between sphi & ephi checks
1096  //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
1097  //            intersection check <=0 -> >=0
1098  //         -> use some form of loop Construct ?
1099  //
1100  if ( seg )
1101  {
1102    // First phi surface (Starting phi)
1103
1104    sinSPhi = std::sin(fSPhi) ;
1105    cosSPhi = std::cos(fSPhi) ;
1106    Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
1107                   
1108    if ( Comp < 0 )  // Component in outwards normal dirn
1109    {
1110      Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
1111
1112      if ( Dist < kCarTolerance*0.5 )
1113      {
1114        s = Dist/Comp ;
1115
1116        if (s < snxt)
1117        {
1118          if ( s < 0 )  { s = 0.0; }
1119          zi = p.z() + s*v.z() ;
1120          if ( std::fabs(zi) <= tolODz )
1121          {
1122            xi   = p.x() + s*v.x() ;
1123            yi   = p.y() + s*v.y() ;
1124            rho2 = xi*xi + yi*yi ;
1125
1126            if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1127              || ( (rho2 >  tolORMin2) && (rho2 <  tolIRMin2)
1128                && ( v.y()*cosSPhi - v.x()*sinSPhi >  0 )
1129                && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 )     )
1130              || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1131                && (v.y()*cosSPhi - v.x()*sinSPhi > 0)
1132                && (v.x()*cosSPhi + v.y()*sinSPhi < 0) )    )
1133            {
1134              // z and r intersections good
1135              // - check intersecting with correct half-plane
1136              //
1137              if ((yi*cosCPhi-xi*sinCPhi) <= 0) snxt = s ;
1138            }   
1139          }
1140        }
1141      }   
1142    }
1143     
1144    // Second phi surface (`E'nding phi)
1145
1146    ePhi    = fSPhi + fDPhi ;
1147    sinEPhi = std::sin(ePhi) ;
1148    cosEPhi = std::cos(ePhi) ;
1149    Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
1150       
1151    if (Comp < 0 )  // Component in outwards normal dirn
1152    {
1153      Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
1154
1155      if ( Dist < kCarTolerance*0.5 )
1156      {
1157        s = Dist/Comp ;
1158
1159        if (s < snxt)
1160        {
1161          if ( s < 0 )  { s = 0; }
1162          zi = p.z() + s*v.z() ;
1163          if ( std::fabs(zi) <= tolODz )
1164          {
1165            xi   = p.x() + s*v.x() ;
1166            yi   = p.y() + s*v.y() ;
1167            rho2 = xi*xi + yi*yi ;
1168            if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) )
1169                || ( (rho2 > tolORMin2)  && (rho2 < tolIRMin2)
1170                  && (v.x()*sinEPhi - v.y()*cosEPhi >  0)
1171                  && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) )
1172                || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2)
1173                  && (v.x()*sinEPhi - v.y()*cosEPhi > 0)
1174                  && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) )
1175            {
1176              // z and r intersections good
1177              // - check intersecting with correct half-plane
1178              //
1179              if ( (yi*cosCPhi-xi*sinCPhi) >= 0 )  { snxt = s; }
1180            }   
1181          }
1182        }
1183      }
1184    }         //  Comp < 0
1185  }           //  seg != 0
1186  if ( snxt<kCarTolerance*0.5 )  { snxt=0; }
1187  return snxt ;
1188}
1189 
1190//////////////////////////////////////////////////////////////////
1191//
1192// Calculate distance to shape from outside, along normalised vector
1193// - return kInfinity if no intersection, or intersection distance <= tolerance
1194//
1195// - Compute the intersection with the z planes
1196//        - if at valid r, phi, return
1197//
1198// -> If point is outer outer radius, compute intersection with rmax
1199//        - if at valid phi,z return
1200//
1201// -> Compute intersection with inner radius, taking largest +ve root
1202//        - if valid (in z,phi), save intersction
1203//
1204//    -> If phi segmented, compute intersections with phi half planes
1205//        - return smallest of valid phi intersections and
1206//          inner radius intersection
1207//
1208// NOTE:
1209// - Precalculations for phi trigonometry are Done `just in time'
1210// - `if valid' implies tolerant checking of intersection points
1211//   Calculate distance (<= actual) to closest surface of shape from outside
1212// - Calculate distance to z, radial planes
1213// - Only to phi planes if outside phi extent
1214// - Return 0 if point inside
1215
1216G4double G4Tubs::DistanceToIn( const G4ThreeVector& p ) const
1217{
1218  G4double safe=0.0, rho, safe1, safe2, safe3 ;
1219  G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
1220
1221  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1222  safe1 = fRMin - rho ;
1223  safe2 = rho - fRMax ;
1224  safe3 = std::fabs(p.z()) - fDz ;
1225
1226  if ( safe1 > safe2 ) { safe = safe1; }
1227  else                 { safe = safe2; }
1228  if ( safe3 > safe )  { safe = safe3; }
1229
1230  if (fDPhi < twopi && rho)
1231  {
1232    phiC    = fSPhi + fDPhi*0.5 ;
1233    cosPhiC = std::cos(phiC) ;
1234    sinPhiC = std::sin(phiC) ;
1235
1236    // Psi=angle from central phi to point
1237    //
1238    cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
1239
1240    if ( cosPsi < std::cos(fDPhi*0.5) )
1241    {
1242      // Point lies outside phi range
1243
1244      if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1245      {
1246        safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1247      }
1248      else
1249      {
1250        ePhi    = fSPhi + fDPhi ;
1251        safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1252      }
1253      if ( safePhi > safe )  { safe = safePhi; }
1254    }
1255  }
1256  if ( safe < 0 )  { safe = 0; }
1257  return safe ;
1258}
1259
1260//////////////////////////////////////////////////////////////////////////////
1261//
1262// Calculate distance to surface of shape from `inside', allowing for tolerance
1263// - Only Calc rmax intersection if no valid rmin intersection
1264
1265G4double G4Tubs::DistanceToOut( const G4ThreeVector& p,
1266                                const G4ThreeVector& v,
1267                                const G4bool calcNorm,
1268                                      G4bool *validNorm,
1269                                      G4ThreeVector *n    ) const
1270{ 
1271  ESide side = kNull , sider = kNull, sidephi = kNull ;
1272  G4double snxt, sr = kInfinity, sphi = kInfinity, pdist ;
1273  G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ;
1274
1275  // Vars for phi intersection:
1276
1277  G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi ;
1278  G4double cPhi, sinCPhi, cosCPhi ;
1279  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ;
1280 
1281  // Z plane intersection
1282
1283  if (v.z() > 0 )
1284  {
1285    pdist = fDz - p.z() ;
1286    if ( pdist > kCarTolerance*0.5 )
1287    {
1288      snxt = pdist/v.z() ;
1289      side = kPZ ;
1290    }
1291    else
1292    {
1293      if (calcNorm)
1294      {
1295        *n         = G4ThreeVector(0,0,1) ;
1296        *validNorm = true ;
1297      }
1298      return snxt = 0 ;
1299    }
1300  }
1301  else if ( v.z() < 0 )
1302  {
1303    pdist = fDz + p.z() ;
1304
1305    if ( pdist > kCarTolerance*0.5 )
1306    {
1307      snxt = -pdist/v.z() ;
1308      side = kMZ ;
1309    }
1310    else
1311    {
1312      if (calcNorm)
1313      {
1314        *n         = G4ThreeVector(0,0,-1) ;
1315        *validNorm = true ;
1316      }
1317      return snxt = 0.0 ;
1318    }
1319  }
1320  else
1321  {
1322    snxt = kInfinity ;    // Travel perpendicular to z axis
1323    side = kNull;
1324  }
1325
1326  // Radial Intersections
1327  //
1328  // Find intersction with cylinders at rmax/rmin
1329  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
1330  //
1331  // Intersects with x^2+y^2=R^2
1332  //
1333  // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0
1334  //
1335  //            t1                t2                    t3
1336
1337  t1   = 1.0 - v.z()*v.z() ;      // since v normalised
1338  t2   = p.x()*v.x() + p.y()*v.y() ;
1339  t3   = p.x()*p.x() + p.y()*p.y() ;
1340
1341  if ( snxt > 10*(fDz+fRMax) )  { roi2 = 2*fRMax*fRMax; }
1342  else  { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; }        // radius^2 on +-fDz
1343
1344  if ( t1 > 0 ) // Check not parallel
1345  {
1346    // Calculate sr, r exit distance
1347     
1348    if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) )
1349    {
1350      // Delta r not negative => leaving via rmax
1351
1352      deltaR = t3 - fRMax*fRMax ;
1353
1354      // NOTE: Should use rho-fRMax<-kRadTolerance*0.5
1355      // - avoid sqrt for efficiency
1356
1357      if ( deltaR < -kRadTolerance*fRMax )
1358      {
1359        b     = t2/t1 ;
1360        c     = deltaR/t1 ;
1361        d2= b*b-c;
1362        if(d2>=0.){sr    = -b + std::sqrt(d2);}
1363        else{sr=0.;};
1364        sider = kRMax ;
1365      }
1366      else
1367      {
1368        // On tolerant boundary & heading outwards (or perpendicular to)
1369        // outer radial surface -> leaving immediately
1370
1371        if ( calcNorm ) 
1372        {
1373          // if ( p.x() || p.y() )
1374          // {
1375          //  *n=G4ThreeVector(p.x(),p.y(),0);
1376          // }
1377          // else
1378          // {
1379          //  *n=v;
1380          // }
1381          *n         = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1382          *validNorm = true ;
1383        }
1384        return snxt = 0 ; // Leaving by rmax immediately
1385      }
1386    }             
1387    else if ( t2 < 0. ) // i.e.  t2 < 0; Possible rmin intersection
1388    {
1389      roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement
1390
1391      if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) )
1392      {
1393        deltaR = t3 - fRMin*fRMin ;
1394        b      = t2/t1 ;
1395        c      = deltaR/t1 ;
1396        d2     = b*b - c ;
1397
1398        if ( d2 >= 0 )   // Leaving via rmin
1399        {
1400          // NOTE: SHould use rho-rmin>kRadTolerance*0.5
1401          // - avoid sqrt for efficiency
1402
1403          if (deltaR > kRadTolerance*fRMin)
1404          {
1405            sr    = -b-std::sqrt(d2) ;
1406            sider = kRMin ;
1407          }
1408          else
1409          {
1410            if ( calcNorm ) { *validNorm = false; }  // Concave side
1411            return snxt = 0.0;
1412          }
1413        }
1414        else    // No rmin intersect -> must be rmax intersect
1415        {
1416          deltaR = t3 - fRMax*fRMax ;
1417          c     = deltaR/t1 ;
1418          d2    = b*b-c;
1419          if(d2>=0.)
1420          {
1421            sr     = -b + std::sqrt(d2) ;
1422            sider  = kRMax ;
1423          }
1424          else // Case: On the border+t2<kRadTolerance
1425               //       (v is perpendiculair to the surface)
1426          {
1427            if (calcNorm)
1428            {
1429              *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1430              *validNorm = true ;
1431            }
1432            return snxt = 0.0;
1433          }
1434        }
1435      }
1436      else if ( roi2 > fRMax*(fRMax + kRadTolerance) )
1437           // No rmin intersect -> must be rmax intersect
1438      {
1439        deltaR = t3 - fRMax*fRMax ;
1440        b      = t2/t1 ;
1441        c      = deltaR/t1;
1442        d2     = b*b-c;
1443        if(d2>=0.)
1444        {
1445          sr     = -b + std::sqrt(d2) ;
1446          sider  = kRMax ;
1447        }
1448        else // Case: On the border+t2<kRadTolerance
1449             //       (v is perpendiculair to the surface)
1450        {
1451          if (calcNorm)
1452          {
1453            *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ;
1454            *validNorm = true ;
1455          }
1456          return snxt = 0.0;
1457        }
1458      }
1459    }
1460   
1461    // Phi Intersection
1462
1463    if ( fDPhi < twopi )
1464    {
1465      sinSPhi = std::sin(fSPhi) ;
1466      cosSPhi = std::cos(fSPhi) ;
1467      ePhi    = fSPhi + fDPhi ;
1468      sinEPhi = std::sin(ePhi) ;
1469      cosEPhi = std::cos(ePhi) ;
1470      cPhi    = fSPhi + fDPhi*0.5 ;
1471      sinCPhi = std::sin(cPhi) ;
1472      cosCPhi = std::cos(cPhi) ;
1473
1474      // add angle calculation with correction
1475      // of the difference in domain of atan2 and Sphi
1476      //
1477      vphi = std::atan2(v.y(),v.x()) ;
1478
1479      if ( vphi < fSPhi - kAngTolerance*0.5  )             { vphi += twopi; }
1480      else if ( vphi > fSPhi + fDPhi + kAngTolerance*0.5 ) { vphi -= twopi; }
1481
1482
1483      if ( p.x() || p.y() )  // Check if on z axis (rho not needed later)
1484      {
1485        // pDist -ve when inside
1486
1487        pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
1488        pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
1489
1490        // Comp -ve when in direction of outwards normal
1491
1492        compS   = -sinSPhi*v.x() + cosSPhi*v.y() ;
1493        compE   =  sinEPhi*v.x() - cosEPhi*v.y() ;
1494        sidephi = kNull;
1495       
1496        if( ( (fDPhi <= pi) && ( (pDistS <= 0.5*kCarTolerance)
1497                              && (pDistE <= 0.5*kCarTolerance) ) )
1498         || ( (fDPhi >  pi) && !((pDistS >  0.5*kCarTolerance)
1499                              && (pDistE >  0.5*kCarTolerance) ) )  )
1500        {
1501          // Inside both phi *full* planes
1502         
1503          if ( compS < 0 )
1504          {
1505            sphi = pDistS/compS ;
1506           
1507            if (sphi >= -0.5*kCarTolerance)
1508            {
1509              xi = p.x() + sphi*v.x() ;
1510              yi = p.y() + sphi*v.y() ;
1511             
1512              // Check intersecting with correct half-plane
1513              // (if not -> no intersect)
1514              //
1515              if((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
1516                { sidephi = kSPhi;
1517                if (((fSPhi-0.5*kAngTolerance)<=vphi)
1518                   &&((ePhi+0.5*kAngTolerance)>=vphi))
1519                {
1520                  sphi = kInfinity;
1521                }
1522              }
1523              else if ((yi*cosCPhi-xi*sinCPhi)>=0)
1524              {
1525                sphi = kInfinity ;
1526              }
1527              else
1528              {
1529                sidephi = kSPhi ;
1530                if ( pDistS > -kCarTolerance*0.5 )
1531                {
1532                  sphi = 0.0 ; // Leave by sphi immediately
1533                }   
1534              }       
1535            }
1536            else
1537            {
1538              sphi = kInfinity ;
1539            }
1540          }
1541          else
1542          {
1543            sphi = kInfinity ;
1544          }
1545
1546          if ( compE < 0 )
1547          {
1548            sphi2 = pDistE/compE ;
1549           
1550            // Only check further if < starting phi intersection
1551            //
1552            if ( (sphi2 > -0.5*kCarTolerance) && (sphi2 < sphi) )
1553            {
1554              xi = p.x() + sphi2*v.x() ;
1555              yi = p.y() + sphi2*v.y() ;
1556             
1557              if ((std::abs(xi)<=kCarTolerance)&&(std::abs(yi)<=kCarTolerance))
1558              {
1559                // Leaving via ending phi
1560                //
1561                if(!(((fSPhi-0.5*kAngTolerance)<=vphi)
1562                    &&((ePhi+0.5*kAngTolerance)>=vphi)))
1563                {
1564                  sidephi = kEPhi ;
1565                  if ( pDistE <= -kCarTolerance*0.5 )  { sphi = sphi2 ; }
1566                  else                                 { sphi = 0.0 ; }
1567                }
1568              } 
1569              else    // Check intersecting with correct half-plane
1570
1571              if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
1572              {
1573                // Leaving via ending phi
1574                //
1575                sidephi = kEPhi ;
1576                if ( pDistE <= -kCarTolerance*0.5 ) { sphi = sphi2 ; }
1577                else                                { sphi = 0.0 ; }
1578              }
1579            }
1580          }
1581        }
1582        else
1583        {
1584          sphi = kInfinity ;
1585        }
1586      }
1587      else
1588      {
1589        // On z axis + travel not || to z axis -> if phi of vector direction
1590        // within phi of shape, Step limited by rmax, else Step =0
1591
1592        // vphi = std::atan2(v.y(),v.x()) ;//defined previosly
1593        // G4cout<<"In axis vphi="<<vphi
1594        //       <<" Sphi="<<fSPhi<<" Ephi="<<ePhi<<G4endl;
1595        // old  if ( (fSPhi < vphi) && (vphi < fSPhi + fDPhi) )
1596        // new : correction for if statement, must be '<=' 
1597               
1598        if ( ((fSPhi-0.5*kAngTolerance) <= vphi)
1599           && (vphi <=( ePhi+0.5*kAngTolerance) )) 
1600        {
1601          sphi = kInfinity ;
1602        }
1603        else
1604        {
1605          sidephi = kSPhi ; // arbitrary
1606          sphi    = 0.0 ;
1607        }
1608      }
1609      if (sphi < snxt)  // Order intersecttions
1610      {
1611        snxt = sphi ;
1612        side = sidephi ;
1613      }
1614    }
1615    if (sr < snxt)  // Order intersections
1616    {
1617      snxt = sr ;
1618      side = sider ;
1619    }
1620  }
1621  if (calcNorm)
1622  {
1623    switch(side)
1624    {
1625      case kRMax:
1626        // Note: returned vector not normalised
1627        // (divide by fRMax for unit vector)
1628        //
1629        xi = p.x() + snxt*v.x() ;
1630        yi = p.y() + snxt*v.y() ;
1631        *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ;
1632        *validNorm = true ;
1633        break ;
1634
1635      case kRMin:
1636        *validNorm = false ;  // Rmin is inconvex
1637        break ;
1638
1639      case kSPhi:
1640        if ( fDPhi <= pi )
1641        {
1642          *n         = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
1643          *validNorm = true ;
1644        }
1645        else
1646        {
1647          *validNorm = false ;
1648        }
1649        break ;
1650
1651      case kEPhi:
1652        if (fDPhi <= pi)
1653        {
1654          *n = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
1655          *validNorm = true ;
1656        }
1657        else
1658        {
1659          *validNorm = false ;
1660        }
1661        break ;
1662
1663      case kPZ:
1664        *n=G4ThreeVector(0,0,1) ;
1665        *validNorm=true ;
1666        break ;
1667
1668      case kMZ:
1669        *n         = G4ThreeVector(0,0,-1) ;
1670        *validNorm = true ;
1671        break ;
1672
1673      default:
1674        G4cout.precision(16) ;
1675        G4cout << G4endl ;
1676        DumpInfo();
1677        G4cout << "Position:"  << G4endl << G4endl ;
1678        G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1679        G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1680        G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1681        G4cout << "Direction:" << G4endl << G4endl ;
1682        G4cout << "v.x() = "   << v.x() << G4endl ;
1683        G4cout << "v.y() = "   << v.y() << G4endl ;
1684        G4cout << "v.z() = "   << v.z() << G4endl << G4endl ;
1685        G4cout << "Proposed distance :" << G4endl << G4endl ;
1686        G4cout << "snxt = "    << snxt/mm << " mm" << G4endl << G4endl ;
1687        G4Exception("G4Tubs::DistanceToOut(p,v,..)","Notification",JustWarning,
1688                    "Undefined side for valid surface normal to solid.");
1689        break ;
1690    }
1691  }
1692  if ( snxt<kCarTolerance*0.5 )  { snxt=0 ; }
1693
1694  return snxt ;
1695}
1696
1697//////////////////////////////////////////////////////////////////////////
1698//
1699// Calculate distance (<=actual) to closest surface of shape from inside
1700
1701G4double G4Tubs::DistanceToOut( const G4ThreeVector& p ) const
1702{
1703  G4double safe=0.0, rho, safeR1, safeR2, safeZ ;
1704  G4double safePhi, phiC, cosPhiC, sinPhiC, ePhi ;
1705  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
1706
1707#ifdef G4CSGDEBUG
1708  if( Inside(p) == kOutside )
1709  {
1710    G4cout.precision(16) ;
1711    G4cout << G4endl ;
1712    DumpInfo();
1713    G4cout << "Position:"  << G4endl << G4endl ;
1714    G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
1715    G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
1716    G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
1717    G4Exception("G4Tubs::DistanceToOut(p)", "Notification", JustWarning, 
1718                 "Point p is outside !?");
1719  }
1720#endif
1721
1722  if ( fRMin )
1723  {
1724    safeR1 = rho   - fRMin ;
1725    safeR2 = fRMax - rho ;
1726 
1727    if ( safeR1 < safeR2 ) { safe = safeR1 ; }
1728    else                   { safe = safeR2 ; }
1729  }
1730  else
1731  {
1732    safe = fRMax - rho ;
1733  }
1734  safeZ = fDz - std::fabs(p.z()) ;
1735
1736  if ( safeZ < safe )  { safe = safeZ ; }
1737
1738  // Check if phi divided, Calc distances closest phi plane
1739  //
1740  if ( fDPhi < twopi )
1741  {
1742    // Above/below central phi of Tubs?
1743
1744    phiC    = fSPhi + fDPhi*0.5 ;
1745    cosPhiC = std::cos(phiC) ;
1746    sinPhiC = std::sin(phiC) ;
1747
1748    if ( (p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
1749    {
1750      safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
1751    }
1752    else
1753    {
1754      ePhi    = fSPhi + fDPhi ;
1755      safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
1756    }
1757    if (safePhi < safe)  { safe = safePhi ; }
1758  }
1759  if ( safe < 0 )  { safe = 0 ; }
1760
1761  return safe ; 
1762}
1763
1764/////////////////////////////////////////////////////////////////////////
1765//
1766// Create a List containing the transformed vertices
1767// Ordering [0-3] -fDz cross section
1768//          [4-7] +fDz cross section such that [0] is below [4],
1769//                                             [1] below [5] etc.
1770// Note:
1771//  Caller has deletion resposibility
1772//  Potential improvement: For last slice, use actual ending angle
1773//                         to avoid rounding error problems.
1774
1775G4ThreeVectorList*
1776G4Tubs::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
1777{
1778  G4ThreeVectorList* vertices ;
1779  G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
1780  G4double meshAngle, meshRMax, crossAngle,
1781           cosCrossAngle, sinCrossAngle, sAngle;
1782  G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ;
1783  G4int crossSection, noCrossSections;
1784
1785  // Compute no of cross-sections necessary to mesh tube
1786  //
1787  noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;
1788
1789  if ( noCrossSections < kMinMeshSections )
1790  {
1791    noCrossSections = kMinMeshSections ;
1792  }
1793  else if (noCrossSections>kMaxMeshSections)
1794  {
1795    noCrossSections = kMaxMeshSections ;
1796  }
1797  // noCrossSections = 4 ;
1798
1799  meshAngle = fDPhi/(noCrossSections - 1) ;
1800  // meshAngle = fDPhi/(noCrossSections) ;
1801
1802  meshRMax  = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ;
1803  meshRMin = fRMin - 100*kCarTolerance ; 
1804 
1805  // If complete in phi, set start angle such that mesh will be at fRMax
1806  // on the x axis. Will give better extent calculations when not rotated.
1807
1808  if (fDPhi == pi*2.0 && fSPhi == 0 )  { sAngle = -meshAngle*0.5 ; }
1809  else                                 { sAngle =  fSPhi ; }
1810   
1811  vertices = new G4ThreeVectorList();
1812  vertices->reserve(noCrossSections*4);
1813   
1814  if ( vertices )
1815  {
1816    for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ )
1817    {
1818      // Compute coordinates of cross section at section crossSection
1819
1820      crossAngle    = sAngle + crossSection*meshAngle ;
1821      cosCrossAngle = std::cos(crossAngle) ;
1822      sinCrossAngle = std::sin(crossAngle) ;
1823
1824      rMaxX = meshRMax*cosCrossAngle ;
1825      rMaxY = meshRMax*sinCrossAngle ;
1826
1827      if(meshRMin <= 0.0)
1828      {
1829        rMinX = 0.0 ;
1830        rMinY = 0.0 ;
1831      }
1832      else
1833      {
1834        rMinX = meshRMin*cosCrossAngle ;
1835        rMinY = meshRMin*sinCrossAngle ;
1836      }
1837      vertex0 = G4ThreeVector(rMinX,rMinY,-fDz) ;
1838      vertex1 = G4ThreeVector(rMaxX,rMaxY,-fDz) ;
1839      vertex2 = G4ThreeVector(rMaxX,rMaxY,+fDz) ;
1840      vertex3 = G4ThreeVector(rMinX,rMinY,+fDz) ;
1841
1842      vertices->push_back(pTransform.TransformPoint(vertex0)) ;
1843      vertices->push_back(pTransform.TransformPoint(vertex1)) ;
1844      vertices->push_back(pTransform.TransformPoint(vertex2)) ;
1845      vertices->push_back(pTransform.TransformPoint(vertex3)) ;
1846    }
1847  }
1848  else
1849  {
1850    DumpInfo();
1851    G4Exception("G4Tubs::CreateRotatedVertices()",
1852                "FatalError", FatalException,
1853                "Error in allocation of vertices. Out of memory !");
1854  }
1855  return vertices ;
1856}
1857
1858//////////////////////////////////////////////////////////////////////////
1859//
1860// Stream object contents to an output stream
1861
1862G4GeometryType G4Tubs::GetEntityType() const
1863{
1864  return G4String("G4Tubs");
1865}
1866
1867//////////////////////////////////////////////////////////////////////////
1868//
1869// Stream object contents to an output stream
1870
1871std::ostream& G4Tubs::StreamInfo( std::ostream& os ) const
1872{
1873  os << "-----------------------------------------------------------\n"
1874     << "    *** Dump for solid - " << GetName() << " ***\n"
1875     << "    ===================================================\n"
1876     << " Solid type: G4Tubs\n"
1877     << " Parameters: \n"
1878     << "    inner radius : " << fRMin/mm << " mm \n"
1879     << "    outer radius : " << fRMax/mm << " mm \n"
1880     << "    half length Z: " << fDz/mm << " mm \n"
1881     << "    starting phi : " << fSPhi/degree << " degrees \n"
1882     << "    delta phi    : " << fDPhi/degree << " degrees \n"
1883     << "-----------------------------------------------------------\n";
1884
1885  return os;
1886}
1887
1888/////////////////////////////////////////////////////////////////////////
1889//
1890// GetPointOnSurface
1891
1892G4ThreeVector G4Tubs::GetPointOnSurface() const
1893{
1894  G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose,
1895           aOne, aTwo, aThr, aFou;
1896  G4double rRand;
1897
1898  aOne = 2.*fDz*fDPhi*fRMax;
1899  aTwo = 2.*fDz*fDPhi*fRMin;
1900  aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin);
1901  aFou = 2.*fDz*(fRMax-fRMin);
1902
1903  phi    = RandFlat::shoot(fSPhi, fSPhi+fDPhi);
1904  cosphi = std::cos(phi);
1905  sinphi = std::sin(phi);
1906
1907  rRand  = RandFlat::shoot(fRMin,fRMax);
1908 
1909  if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; }
1910 
1911  chose  = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou);
1912
1913  if( (chose >=0) && (chose < aOne) )
1914  {
1915    xRand = fRMax*cosphi;
1916    yRand = fRMax*sinphi;
1917    zRand = RandFlat::shoot(-1.*fDz,fDz);
1918    return G4ThreeVector  (xRand, yRand, zRand);
1919  }
1920  else if( (chose >= aOne) && (chose < aOne + aTwo) )
1921  {
1922    xRand = fRMin*cosphi;
1923    yRand = fRMin*sinphi;
1924    zRand = RandFlat::shoot(-1.*fDz,fDz);
1925    return G4ThreeVector  (xRand, yRand, zRand);
1926  }
1927  else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) )
1928  {
1929    xRand = rRand*cosphi;
1930    yRand = rRand*sinphi;
1931    zRand = fDz;
1932    return G4ThreeVector  (xRand, yRand, zRand);
1933  }
1934  else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) )
1935  {
1936    xRand = rRand*cosphi;
1937    yRand = rRand*sinphi;
1938    zRand = -1.*fDz;
1939    return G4ThreeVector  (xRand, yRand, zRand);
1940  }
1941  else if( (chose >= aOne + aTwo + 2.*aThr)
1942        && (chose < aOne + aTwo + 2.*aThr + aFou) )
1943  {
1944    xRand = rRand*std::cos(fSPhi);
1945    yRand = rRand*std::sin(fSPhi);
1946    zRand = RandFlat::shoot(-1.*fDz,fDz);
1947    return G4ThreeVector  (xRand, yRand, zRand);
1948  }
1949  else
1950  {
1951    xRand = rRand*std::cos(fSPhi+fDPhi);
1952    yRand = rRand*std::sin(fSPhi+fDPhi);
1953    zRand = RandFlat::shoot(-1.*fDz,fDz);
1954    return G4ThreeVector  (xRand, yRand, zRand);
1955  }
1956}
1957
1958///////////////////////////////////////////////////////////////////////////
1959//
1960// Methods for visualisation
1961
1962void G4Tubs::DescribeYourselfTo ( G4VGraphicsScene& scene ) const 
1963{
1964  scene.AddSolid (*this) ;
1965}
1966
1967G4Polyhedron* G4Tubs::CreatePolyhedron () const 
1968{
1969  return new G4PolyhedronTubs (fRMin, fRMax, fDz, fSPhi, fDPhi) ;
1970}
1971
1972G4NURBS* G4Tubs::CreateNURBS () const 
1973{
1974  G4NURBS* pNURBS ;
1975  if (fRMin != 0) 
1976  {
1977    if (fDPhi >= twopi) 
1978    {
1979      pNURBS = new G4NURBStube (fRMin,fRMax,fDz) ;
1980    }
1981    else 
1982    {
1983      pNURBS = new G4NURBStubesector (fRMin,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
1984    }
1985  }
1986  else 
1987  {
1988    if (fDPhi >= twopi) 
1989    {
1990      pNURBS = new G4NURBScylinder (fRMax,fDz) ;
1991    }
1992    else 
1993    {
1994      const G4double epsilon = 1.e-4 ; // Cylinder sector not yet available!
1995      pNURBS = new G4NURBStubesector (epsilon,fRMax,fDz,fSPhi,fSPhi+fDPhi) ;
1996    }
1997  }
1998  return pNURBS ;
1999}
Note: See TracBrowser for help on using the repository browser.