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

Last change on this file since 1198 was 1058, checked in by garnier, 15 years ago

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